PRO wt_record

; This IDL program is used to generate a long-term water table record
; for Lost Creek for Faith Ann's BIOME-BGC modeling efforts.  The
; following steps are taken in the process:
;
;   - opens and plots water table record for Lost Creek and Bear River (USGS);
;   - computes and plots correlation for 2000-2005 data; and 
;   - extrapolates WT back to 1991.
;
; Also used to compare 

;Define plot colors
;   0=black
;   1=red
;   2=green
;   3=blue
;   4=cyan
;   5=magenta
;   6=yellow
;   7=white
red=[0,1,0,0,0,1,1,1]*255
green=[0,0,1,0,1,0,1,1]*255
blue=[0,0,0,1,1,1,0,1]*255
tvlct,red,green,blue

;open USGS observations
;loadsmallcommagz,'/eddy/s2/bcook/wdepth/bearcreek_obs.gz',bc_obs

;open USGS daily averaged streamflow (ft^3/s) -- 1 May 1991 through 31 Dec 2006
loadsmallcommagz,'/eddy/s2/bcook/wdepth/bearcreek_daily.gz',bc_daily
bc_jday=julday(bc_daily(1,*),bc_daily(2,*),bc_daily(0,*))-julday(12,31,bc_daily(0,*)-1)
index=where(bc_daily(3,*) GT 0)
bc_daily(3,index)=bc_daily(3,index)*0.02831685; ft^3/s --> m^3/s

;open LC#1 observations and time stamp (30 min intervals begining 1 Jan 1997)
;variables:  wdepth_ald_fill,wdepth_ald_flag,wdepth_lc_fill,wdepth_lc_flag
restore,'/eddy/s2/bcook/variables/soilt_h2o.sav'
restore,'/eddy/s2/bcook/variables/timestamp.sav'

;open LC#2 observations (30 min intervals starting June 2005)
;loadsmallgz,'/eddy/s1/cheas/micromet/wdepth/30min/wdepth2005_30min.gz',lc2
;0=year
;1=jday
;2=hrmin
;5=lc2

;open LC daily averages (starting 1 Jan 1997)
restore,'/eddy/s2/bcook/variables/daily.sav'
;0=year
;1=month
;2=jday
;46=average obs
;47=filled

;create array containing good observations that can be compared
j=0
corr=fltarr(5,4000); year,month,jday,lc,bc
FOR i=0,n_elements(daily_met(0,*))-1. DO BEGIN
   index=where(bc_daily(0,*) EQ daily_met(0,i) AND bc_jday(*) EQ daily_met(2,i) AND bc_daily(4,*) NE 3 AND daily_met(46,i) gt -80. AND bc_daily(3,*) GT 0 AND daily_met(0,i) GT 2000)
   IF index(0) NE -1 THEN BEGIN
      corr(0,j)=daily_met(0,i)
      corr(1,j)=daily_met(1,i)
      corr(2,j)=daily_met(2,i)
      corr(3,j)=daily_met(46,i)
      corr(4,j)=bc_daily(3,index)
      j=j+1
   ENDIF
ENDFOR
corr=corr(*,0:j-1)

;segmented linear fit (breakpoint=-20 cm)
ige20=where(corr(3,*) GE -20.)
fge20=linear_fit(corr(3,ige20),corr(4,ige20))
ilt20=where(corr(3,*) LT -20.)
flt20=linear_fit(corr(3,ilt20),corr(4,ilt20))

;plot data and segmented regression
set_plot,'ps'
device,file='/eddy/s2/bcook/wdepth/wtable_sflow.ps',/inches,xoffset=1,yoffset=3.5,xsize=6.5,ysize=6.5,/helvetica,/color
xequal=(flt20(0)-fge20(0))/(fge20(1)-flt20(1));x value where two equations are equal
yequal=flt20(0)+flt20(1)*xequal
;window,1,title='Water table vs. Stream flow',xsize=400,ysize=400
plot,corr(3,*),corr(4,*),psym=1,xtitle='Water table (cm above surface)',ytitle='Stream flow (m!u3!n s!u-1!n)',charthick=4,charsize=1.75,xthick=4,ythick=4
oplot,[-100,xequal],flt20(0)+flt20(1)*[-100,xequal],color=1,thick=4
oplot,[xequal,100],fge20(0)+fge20(1)*[xequal,100],color=1,thick=4
device,/close
set_plot,'x'

;extrapolate water table measurements and output data
wt=fltarr(6,n_elements(bc_daily(0,*)))-999.
wt(0:4,*)=bc_daily
FOR i=0,n_elements(bc_daily(0,*))-1 DO BEGIN
   IF wt(3,i) LT yequal THEN wt(3,i)=(wt(3,i)-flt20(0))/flt20(1) $
   ELSE wt(3,i)=(wt(3,i)-fge20(0))/fge20(1)
   wt(4,i)=julday(wt(1,i),wt(2,i),wt(0,i))-julday(12,31,wt(0,i)-1)
   index=where(daily_met(0,*) EQ wt(0,i) AND daily_met(2,*) EQ wt(4,i) AND daily_met(47,*) Gt -80)
   IF index(0) NE -1 THEN wt(5,i)=daily_met(47,index)
ENDFOR
index=where(bc_daily(3,*) EQ -999.)
IF index(0) NE -1 THEN wt(3,index)=-999.

;plot long-term stream flow
xyear=findgen(n_elements(bc_daily(3,*)))/365.25; cumulative years since start of stream flow record
startdate=1991.+((julday(bc_daily(1,0),bc_daily(2,0),bc_daily(0,0))-julday(12,31,bc_daily(0,0)-1.))/365.25)
x=xyear+startdate
index=where(bc_daily(3,*) NE -999. AND xyear GE 10)
sf_trend=linfit(xyear(index),bc_daily(3,index),yfit=sf_line)
window,2,title='Stream flow trend',xsize=800,ysize=400
plot,x(index),bc_daily(3,index);,/ylog,yrange=[3,1000]
oplot,x(index),sf_line,linestyle=1,color=5

;plot extrapolated long-term water table with trendline
set_plot,'ps'
device,file='/eddy/s2/bcook/wdepth/wt_longterm.ps',/inches,xoffset=0,yoffset=5,xsize=8,ysize=5,/helvetica,/color
index=where(wt(3,*) GT -100)
wt_trend=linfit(xyear(index),wt(3,index),yfit=wt_line)
plot,x(index),wt(3,index),thick=2,charsize=1.5,xthick=2,ythick=2,charthick=2,xtitle='Year',ytitle='Water table (cm above surface)',yrange=[-100,80]
oplot,x(index),smooth(wt(3,index),3*365.,/edge_truncate),color=3,thick=2
oplot,x(index),wt_line,color=1,thick=2
xyouts,.2,.2,'Water depth = '+strmid(strtrim(string(wt_trend(0)),2),0,5)+' + ('+strmid(strtrim(string(wt_trend(1)),2),0,5)+' * year)',/normal,charsize=1,color=1,charthick=2
legend,['Extrapolated','3-year moving average','Linear trendline'],color=[0,3,1],linestyle=[0,0,0],/right,box=0,charthick=1,thick=2
device,/close
set_plot,'x'

;plot long-term water table (observations and extrapolation)
window,3,title='Extrapolated water table',xsize=800,ysize=400
plot,wt(3,*),yrange=[-100,80],xtitle='Days since 1 June 1991',ytitle='Water table (cm above surface)',min_val=-998.
oplot,wt(5,*),min_val=-998,color=1
legend,['Extrapolated','Gap-Filled Observations'],color=[0,1],linestyle=[0,0],/right,box=0

;Output composite of extrapolated and gapfilled data
index=where(wt(5,*) GT -999.)
wt(3,index)=wt(5,index)
openw,lun1,filepath('wt_record_lc',root_dir='/eddy/s2/bcook/wdepth/'),/get_lun
printf,lun1,float(wt(0:3,*))
close,lun1 
free_lun,lun1,/force

stop
return
end