;The procedure finalspectable_wc.pro creates final spectral correction tables ;for use by finalflux_wc.pro. Input is the combined .spec tables created ;by dayflux_wc.pro where the combined file is made by create_annual_wc.pro. ;Use is as follows: ; finalspectable_wc ; ;The input file location must be specified by altering the load command below. ;Input filename format: lefYYYY.spec ;Usually these files will be located where the operator has placed the spec ;file output of dayflux_wc. ; ;The final lag files are placed in: ;/data/davis/cheas/wcreek/data/spec/ ;with filename format spectableYYYY ;where YYYY = 4-digit year ; ;Columns of spectable contain the following: ; jday hr c_cfactor q_cfactor ;Where ; jday = julian day ; hr = hour ; c_cfactor = co2 correction factor ; q_cfactor = h2o correction factor ; ;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.pro. ; ;Operating principles of this code - ;Spectral correction factors are checked to see if they are good. If they ;aren't, they are filled using the following methods: ;a)co2 ;Good corrections are found by looking for co2 flags that are equal to zero. ;An array is created that is complete for all hours, each hour holding a good ;correction or a badval. Each element of this array is then checked. Good ;values are retained and badvals are replaced with the average value for that ;particular hour over a selected number of days (mean_period) centered about ;that hour. (If the hour is less than half of mean_period to the start or end ;of the series, the first or last mean_period days are used to calculate the ;average.) If there are no good values in the mean_period, as a last resort, ;the diurnal average is calculated over the mean_period and is spline fit; the ;value from the spline fit is then used for the missing hour. ;b)h2o ;An array for all hours and containing good h2o spectral corrections or badvals ;is created for h2o similarly to co2. Badvals are then replaced with the ;average value for all hours over the selected number of days (mean_period). ;Like co2, if the hour to be filled is less than half of mean_period from the ;start or end of the series, the first or last mean_period days are used to ;calculate the average. ; ;Note: The reason that co2 uses an average for a particular hour where h2o ;does not is because the co2 corrections show a marked diurnal pattern and ;h2o does not. ; ;Written by BWBerger 9/8/99 for wlef. This code complies with file format ;defined by dayflux_wc.pro and its supporting code. ;Modified by WWang ~12/99 to work for Willow Creek. ;Modified by BWB 2/23/00 to make output go to ; /data/davis/wcreek/cheas/data/spec/ ;instead of /data/davis2/wcreek/data/tables/spec/ ;modified by WWG, Sep 12 2000, change directies, ; INDIR='/davis/s1/wang/d2wwg/wcreek/output/spec/' ; OUTDIR='/davis/s1/cheas/wcreek/data/spec/' ; ; finalspectable_lc_cherrey uses kcherrey specific path ; pro finalspectable_lc_cherrey ; input directory ;INDIR='/data/davis2/wwg/wcreek/output/spec/' ; INDIR='/davis/s1/wang/d2wwg/lcreek/output/spec/' INDIR='/abl/s0/users/kcherrey/lcreek/output/spec/' ; output directory ;OUTDIR='/data/davis/cheas/wcreek/data/spec/' ; OUTDIR='/davis/s1/cheas/lcreek/data/spec/' ; OUTDIR='/davis/s2/lcreek/data/spec/' ;change s2 into s4 on 2/21/2003 OUTDIR='/eddy/s4/lcreek/data/spec/' badval=-999. print,'Enter the year to be analysed:YYYY' ;used for filename year = '' read, year ;Load the input data load,INDIR+'lc'+year +'.spec',datain ;*********************** ; Variable initialization n=row(datain) ;total number of half-hours hr=datain(3,*) jday=datain(4,*) ccfac=datain(6,*) c_flag=datain(7,*) qcfac=datain(8,*) q_flag=datain(9,*) good_ccfac=fltarr(n)-999. ;array to hold good ccfactors good_qcfac=fltarr(n)-999. ;array to hold good qcfactors ;arrays to put final lags into final_ccfac = fltarr(n)+badval final_qcfac = fltarr(n)+badval ;The final table. ;Array to hold julian day, hour, final spec corrections for co2 and h2o, resp. dataout=fltarr(4,n) ;************************************* ;Build good arrays from the input data. ;flags are equal to zero when the correction is good. rc=where(c_flag eq 0,cct) if cct ne 0 then begin good_ccfac(rc)=ccfac(rc) ;leaves badvals where not good endif rq=where(q_flag eq 0,qct) if qct ne 0 then begin good_qcfac(rq)=qcfac(rq) ;leaves badvals where not good endif ;******************************************************** ;Missing co2 spectral corrections are replaced with a mean value ;for the specific hour of the surrounding 'mean_period' days of corrections. ;This means that missing values are replaced with a diurnal average value. mean_period = 45L ;number of days to do average over (an averaging window) ; multiply by 48 half-hours/day to get number max no. of values per mean period ;pts = mean_period*24 mean_period = 60L ; inital calculation due to not much data; general set to 45L pts = mean_period*48;for wcreek, 48 half-an-hour per day for i=0L,n-1 do begin ;Loop through all hours ; If a valid ccfac does not exist, use a mean_period day mean value ; for that particular hour. if (good_ccfac(i) eq badval) then begin ;need to get mean correction ;********************************************* ;If PTS/2 hours of prior correction data are not available, use the first ;PTS hours of the record to find a mean fill-in value. if (i-pts/2 lt 0) then begin ; print,'beginning' g = where((good_ccfac(0:pts-1) ne badval) and (hr(0:pts-1) eq hr(i)),num) ;get indexes for particular hour within window. if (num ge 1) then begin final_ccfac(i) = meanbadval(good_ccfac(g)) endif else begin ;need to get spline fit to fill missing hour ;because no good data exists for that hour in window ; print,'spline - beginning' ; print,'hr(i)=',hr(i) ccfac_diavg=fltarr(24*2) for j=0,47 do begin ;get average for all hours r_ccfac=where((c_flag(0:pts-1) eq 0.) and (hr(0:pts-1) eq j),ct) if ct ne 0 then begin ccfac_diavg(j)=meanbadval(ccfac(r_ccfac)) endif else begin ccfac_diavg(j)=-999. endelse endfor tm=findgen(48) ;a time vector for a day rcg=where(ccfac_diavg ne -999.) ;where diurnal average is good rcb=where(ccfac_diavg eq -999.) ;where diurnal average is bad ccfac_sp=spline(tm(rcg),ccfac_diavg(rcg),tm,8) ;fit di_avg final_ccfac(i)=ccfac_sp(hr(i)) ;extract missing hour from fit ;make plots of fits for troubleshooting purposes ;also print out how many hours are missing in the diurnal average. ; plot,tm(rcg),ccfac_diavg(rcg),xrange=[0,23],xstyle=1,$ ; title='ccfactor DIURNAL AVG., LEVEL='+level,$ ; xtitle='Hr (UTC)',ytitle='ccfactor',/ynozero ; oplot,tm,ccfac_sp,linestyle=2,color=7800000 ; refline,1,0,rcb,1 ; print,'c_cfactor diurnal avg not calculated at hours:' ; print,rcb ; print,'Press any key to continue' ; ent='' ; read,ent endelse endif else begin ;***************************************************** ;If PTS/2 hours of future correction data are not available, use the ;last PTS hours of the record to find a mean fill-in value. if (i+pts/2 gt n-1) then begin ; print,'last' g = where((good_ccfac(n-1-pts:n-1) ne badval) and (hr(n-1-pts:n-1) eq hr(i)),num) if (num ge 1) then begin final_ccfac(i) = meanbadval(good_ccfac(g+n-1-pts)) endif else begin ;need to get spline fit to fill missing hour ;because no good data exists for that hour in window ; print,'spline - last' ; print,'hr(i)=',hr(i) ccfac_diavg=fltarr(24*2) for j=0,47 do begin ;get average for all hours r_ccfac=where((c_flag(n-1-pts:n-1) eq 0.) and (hr(n-1-pts:n-1) eq j),ct) if ct ne 0 then begin ccfac_diavg(j)=meanbadval(ccfac(r_ccfac+n-1-pts)) endif else begin ccfac_diavg(j)=-999. endelse endfor tm=findgen(24*2) ;a time vector for a day rcg=where(ccfac_diavg ne -999.) ;where diurnal average is good rcb=where(ccfac_diavg eq -999.) ;where diurnal average is bad ccfac_sp=spline(tm(rcg),ccfac_diavg(rcg),tm,8) ;fit di_avg final_ccfac(i)=ccfac_sp(hr(i)) ;extract missing hour from fit ;make plots of fits for troubleshooting purposes ;also print out how many hours are missing in diurnal average. ; plot,tm(rcg),ccfac_diavg(rcg),xrange=[0,23],xstyle=1,$ ; title='ccfactor DIURNAL AVG., LEVEL='+level,$ ; xtitle='Hr (UTC)',ytitle='ccfactor',/ynozero ; oplot,tm,ccfac_sp,linestyle=2,color=7800000 ; refline,1,0,rcb,1 ; print,'ccfactor diurnal avg not calculated at hours:' ; print,rcb ; print,'Press any key to continue' ; ent='' ; read,ent endelse endif else begin ;If PTS/2 hours exist before and after the point I, use a ;time-centered MEAN_PERIOD day mean to fill in ;missing qcfacs. ; print,'middle' g = where((good_ccfac(i-pts/2:i+pts/2-1) ne badval) and (hr(i-pts/2:i+pts/2-1) eq hr(i)),num) if (num ge 1) then begin final_ccfac(i) = meanbadval(good_ccfac(g+i-pts/2)) endif else begin ;need to get spline fit to fill missing hour ;because no good data exists for that hour in window ; print,'spline - middle' ; print,'hr(i)=',hr(i) ccfac_diavg=fltarr(24*2) for j=0,47 do begin ;get average for all hours r_ccfac=where((c_flag(i-pts/2:i+pts/2-1) eq 0.) and (hr(i-pts/2:i+pts/2-1) eq j),ct) if ct ne 0 then begin ccfac_diavg(j)=meanbadval(ccfac(r_ccfac+i-pts/2)) endif else begin ccfac_diavg(j)=-999. endelse endfor tm=findgen(24*2) ;a time vector for a day rcg=where(ccfac_diavg ne -999.) ;where diurnal average is good rcb=where(ccfac_diavg eq -999.) ;where diurnal average is bad ccfac_sp=spline(tm(rcg),ccfac_diavg(rcg),tm,8) ;fit di_avg final_ccfac(i)=ccfac_sp(hr(i)) ;extract missing hour from fit ;make plots of fits for troubleshooting purposes ;also print out how many hours are missing in diurnal average. ; plot,tm(rcg),ccfac_diavg(rcg),xrange=[0,23],xstyle=1,$ ; title='ccfactor DIURNAL AVG., LEVEL='+level,$ ; xtitle='Hr (UTC)',ytitle='ccfactor',/ynozero ; oplot,tm,ccfac_sp,linestyle=2,color=7800000 ; refline,1,0,rcb,1 ; print,'ccfactor diurnal avg not calculated at hours:' ; print,rcb ; print,'Press any key to continue' ; ent='' ; read,ent endelse endelse endelse endif else begin ;qcfac is good ; print,'good' final_ccfac(i)=good_ccfac(i) endelse ;print,'i=',i endfor ;******************************************************** ;Missing h2o spectral corrections are replaced with a mean value ;of the surrounding 20 days of corrections. ;If there are no valid corrections in the 20 day interval, a badval is placed ;in that hour. for i=0L,n -1 do begin ;Loop through all hours ; If a valid qcfac does not exist, use a mean_period day mean value. if (good_qcfac(i) eq badval) then begin ;need to get mean correction ;If PTS/2 hours of prior correction data are not available, use the first ;PTS hours of the record to find a mean fill-in value. if (i-pts/2 lt 0) then begin g = where(good_qcfac(0:pts-1) ne badval,num) if (num gt 1) then begin final_qcfac(i) = meanbadval(good_qcfac(g)) endif endif else begin ;If PTS/2 hours of future correction data are not available, use the ;last PTS hours of the record to find a mean fill-in value. if (i+pts/2 gt n-1) then begin g = where(good_qcfac(n-1-pts:n-1) ne badval,num) if (num gt 1) then begin final_qcfac(i) = meanbadval(good_qcfac(g+n-1-pts)) endif endif else begin ;If PTS/2 hours exist before and after the point I, use a ;time-centered MEAN_PERIOD day mean to fill in ;missing qcfacs. g = where(good_qcfac(i-pts/2:i+pts/2-1) ne badval,num) if (num gt 1) then begin final_qcfac(i) = meanbadval(good_qcfac(g+i-pts/2)) endif endelse endelse endif else begin ;qcfac is good final_qcfac(i)=good_qcfac(i) endelse endfor ;build final correction table dataout=[jday,hr,transpose(final_ccfac),transpose(final_qcfac)] ;save final correction table saveascii2,dataout,OUTDIR+'spectable'+year return end