;The procedure finallagtable_wc.pro creates final lag tables for use 
;by dayflux.pro. Input is the combined output from the lag tables created
;by dayflux using the appropriate lag windows. Combining the output is done
;by the create_annual.pro.
;Use is as follows:
;                   finallagtable_wc
;
;The input files are assumed to be located in:
;/data/davis2/wwg/wcreek/output/lag/
;with filename format         wcYYYY.lag
;and the final lag files are placed in:
;/data/davis/wcreek/data/tables/lag/
;with filename format         lagtableYYYY.lag
;where     YYYY = 4-digit year
;
;Note: Input data must be complete - that is, all hours must be represented
;      for the block of dates to be processed. This should be done by
;      create_annual_wc.pro.
;Written by BWBerger 9/8/99 based heavily on the original lagcorrection.pro
;by VSinha. This code complies with file format defined by dayflux_wc.pro and
;its supporting code.

pro finallagtable_lc

; input directory
;INDIR='/data/davis2/wwg/wcreek/output/lag/'
;INDIR='/davis/s1/wang/d2wwg/lcreek/output/lag/'

INDIR='/eddy/s2/wang/d2wwg/lcreek/output/lag/'
; output directory
;OUTDIR='/data/davis/wcreek/data/tables/lag/'
;OUTDIR='/davis/s1/cheas/lcreek/data/lag/'
;OUTDIR='/davis/s2/lcreek/data/lag/'
OUTDIR='/eddy/s4/lcreek/data/lag/' ;changed on 2/18/03

badval=-999.

print,'Enter the year to be analysed:YYYY'
year = ''
read, year


;Load the input data
load,INDIR+'lc'+year +'.lag',datain

; Variable initialization
hr=datain(3,*)
jday=datain(4,*)
goodco2=datain(8,*)
goodh2o=datain(9,*)

mm=row(datain) ;total number of hours

;array to hold julian day, hour, final lags for co2 and h2o, resp.
dataout=fltarr(4,mm)

;arrays to put final lags
finalco2 = fltarr(mm)+badval
finalh2o = fltarr(mm)+badval

;********************************************************

; Missing lag times (e.g. little turbulence so no lag time could be
; computed) are replaced with a median value of the surrounding 20 days 
; of lag times.
; This loops assumes there is at least 20 days of valid lag time
; estimates.  If not, the code will not produce a useable lag time
; table.
; Also, if no good lag times exist for 20 days, no lag time will be
; filled in.
; Notes by KJ Davis, 18 Dec, 1998.
; MEDIAN_PERIOD defines in days the length of the period from which to 
; compute a median lag time
median_period = 20
; multiply by 24 hours/day to get number of lag values per median period
pts = median_period*24


for i=0,mm -1 do begin ;Loop through all hours
   ; If a valid lag time does not exist, use a median_period day median value.
   ; First patch the missing co2 lag times.
   if (goodco2(i) eq badval) then begin ;need to get median lag
      ; If PTS/2 hours of prior lag data are not available, use the first
      ; PTS hours of the record to find a median fill-in value.
      if (i-pts/2 lt 0) then begin
         g = where(goodco2(0:pts-1) ne badval,num)
         if (num gt 1) then begin
            finalco2(i) = median(goodco2(g),/EVEN)
         endif
      endif else begin
         ;If PTS/2 hours of future lag data are not available, use the last
         ; PTS hours of the record to find a median fill-in value.
         if (i+pts/2 gt mm-1) then begin
            g = where(goodco2(mm-1-pts:mm-1) ne badval,num)
            if (num gt 1) then begin
               finalco2(i) = median(goodco2(g+mm-1-pts),/EVEN)
            endif
         endif else begin
            ;If PTS/2 hours exist before and after the point I, use a
            ;time-centered MEDIAN_PERIOD day median to fill in 
            ;missing lag times.
            g = where(goodco2(i-pts/2:i+pts/2-1) ne badval,num)
            if (num gt 1) then begin
               finalco2(i) = median(goodco2(g+i-pts/2),/EVEN)
            endif
         endelse
      endelse
   endif else begin ;lag time is good
      finalco2(i)=goodco2(i)
   endelse

   ; Next patch the missing h2o lag times
   if (goodh2o(i) eq badval) then begin ;need to get median lag
      ; If PTS/2 hours of prior lag data are not available, use the first
      ; PTS hours of the record to find a median fill-in value.
      if (i-pts/2 lt 0) then begin
         g = where(goodh2o(0:pts-1) ne badval,num)
         if (num gt 1) then begin
            finalh2o(i) = median(goodh2o(g),/EVEN)
         endif
      endif else begin
         ; If PTS/2 hours of future lag data are not available, use the last
         ; PTS hours of the record to find a median fill-in value.
         if (i+pts/2 gt mm-1) then begin
            g = where(goodh2o(mm-1-pts:mm-1) ne badval,num)
            if (num gt 1) then begin
               finalh2o(i) = median(goodh2o(g+mm-1-pts),/EVEN)
            endif
         endif else begin
            ;If PTS/2 hours exist before and after the point I, use a time-
            ;centered MEDIAN_PERIOD day median to fill in missing lag times.
            g = where(goodh2o(i-pts/2:i+pts/2-1) ne badval,num)
            if (num gt 1) then begin
               finalh2o(i) = median(goodh2o(g+i-pts/2),/EVEN)
            endif
         endelse
      endelse
   endif else begin ;lag time is good
         finalh2o(i)=goodh2o(i)
   endelse

   ;print the final lags to file
   dataout(*,i)=[jday(i),hr(i),finalco2(i), finalh2o(i)]
endfor

;Do final screening for badvals in dataout
for i=0,mm-1 do begin
   if dataout(2,i) eq -999 then begin  ;co2 lag
      if i eq 0 then begin ;if first value is bad, insert first good one
         g_co2=where(dataout(2,*) ne -999.) ;there's got to be a good value!
         dataout(2,i)=dataout(2,g_co2(0))
      endif else begin ;replace badval with first previous good value
                       ;(this fills a bad block with value of left edge)
         dataout(2,i)=dataout(2,i-1)
      endelse
   endif
   if dataout(3,i) eq -999 then begin  ;h2o lag
      if i eq 0 then begin ;if first value is bad, insert first good one
         g_h2o=where(dataout(3,*) ne -999.) ;there's got to be a good value!
         dataout(3,i)=dataout(3,g_h2o(0))
      endif else begin ;replace badval with first previous good value
                       ;(this fills a bad block with value of left edge)
         dataout(3,i)=dataout(3,i-1)
      endelse
   endif
endfor

;save the final lag table
saveascii2,dataout,OUTDIR+'lagtable'+year


return
end