; The program dayflux.pro does the following operations using ; WCREEK data (some are by request only): ; 1) Calculates lag times for flux calculations ; 2) Calculates rotation angles for flux calculations ; 3) Calculates fluxes ; 4) Calculates spectral corrections ; 4) Plots hourly and/or daily data ; Use is as follows: ; dayflux ; ; Input: Necessary information is prompted during use. ; Final rotation, lag and lag window files may be used if they exist ; ; Output: See file 'compileflux' for output file locations. ; See openclose_files_wc.pro and save_wcreek.pro for ; descriptions of output data. ; ;------------------------------------------------------------------------------ ; Written by BWBerger and Dana Carrington 8/99 ; Based on dayflux.pro and combineddayflux.pro ; Modified by BWB on 2/17/00 to add printed message files. ; Modified by BWB on 3/2/00 to make yyddd date format more consistent. ; This format was being used in output file names up to now, but it was ; being done by accident! Now yyddd is the sanctioned format and the ; input date that is requested from the operator is also yyddd. ;Modified by W WAng, Sep 20, 2000 for Lost Creek pro dayflux_lc badval=-999. ;**************************************************************** ;Determine how operator wants analysis to run print,'ENTER STARTING DAY FOR DATA PROCESSING IN yyddd FORMAT:' yyddd='' read,yyddd print,'DO YOU WANT TO RUN A WHOLE DAY: (y/n)?' wholeday='' read,wholeday if (wholeday ne 'y') then begin print,'ENTER HALF HOUR (hhmm format) TO BE EXAMINED (e.g., 0830):' halfhourstring='' read,halfhourstring days='1' endif else begin print,'ENTER NUMBER OF DAYS TO PROCESS (Do not overlap Jan. 1):' days='' read,days endelse print,'IS INPUT DATA COMPRESSED: (y/n)?' print,'(Yes will make code read compressed data files.)' comp='' read, comp print,'DO YOU WANT TO CALCULATE FLUXES: (y/n)?' print,'(Answering yes overwrites .flx and .sflx files for dates chosen.)' calc_flux='' read, calc_flux print,'DO YOU WANT TO USE LONG TERM ROTATION CORRECTIONS: (y/n)?' print,'(Note: .rot files are overwritten for dates chosen regardless of choice.)' usetable_rot = '' read, usetable_rot print, 'DO YOU WANT TO USE FINAL LAG TIME VALUE TABLES: (y/n)?' print,'(Fixed lag times are used if you answer no.)' print,'(Note: .lag files are not affected by this choice.)' ;If lag tables are not used then fixed values are used. See gettablelags.pro usetable_lag='' read, usetable_lag print, 'DO YOU WANT TO PERFORM SPECTRAL CORRECTIONS: (y/n)?' print,'(Answering yes overwrites .spec files for dates chosen.)' calc_speccorr='' read, calc_speccorr print,'DO YOU WANT TO CALCULATE LAG TIMES: (y/n)?' print,'(Answering yes overwrites .lag files for dates chosen.)' calc_lag = '' read, calc_lag ;Lags are determined from cross-correlations and put into the lag files. These ;files are used to create the final lag tables. The determination of the lags ;is improved if a window that limits possible choices is used. Since the lag ;is where the cross-correlation peaks, a window allows only peaks in a ;reasonable time range to be selected. This helps keep more good lags than ;just screening lags calculated without a window. The windows are needed most ;during low-turbulence periods when peaks are less distinguished. ;***Window widths are either selected by the operator ;***or set to a wide default value. ;lagwindow is array as follows: ;ww & wt lag max, wc lag max, wq lag max ;ww & wt lag min, wc lag min, wq lag min lagwindow=[[5,10,10],[-5,-10,-10]] if calc_lag eq 'y' then begin print,'LAG WINDOWS ARE AS FOLLOWS:' print,'ww & wt lag max, wc lag max, wq lag max (typically zero or negative)' print,'ww & wt lag min, wc lag min, wq lag min (typically negative)' print,lagwindow print,'Change windows by editing dayflux_lc.pro line: ~97' print,'' endif print, 'DO YOU WANT TO USE LICOR CALIBRATION TABLES: (y/n)?' usetable_cal='' read,usetable_cal print,'DO YOU WANT TO GENERATE HALF HOURLY PLOTS: (y/n)?' print,'(Answering yes overwrites existing hourly plots for dates chosen.)' plot_halfhourly='' read,plot_halfhourly if wholeday eq 'y' then begin print,'DO YOU WANT TO GENERATE DAILY PLOTS: (y/n)?' print,'(Answering yes overwrites existing daily plots for dates chosen.)' plot_daily='' read,plot_daily endif else begin plot_daily='n' endelse print,'Do you want spec plot(y/n)' specplot='' read,specplot ;*************************************************************************** ;Create and open a file to place printed messages from a run of this code. ;This file is closed at the bottom of this code. ;Note that logical unit number 8 should not be used anywhere ;including openclose_files_wc.pro. ;BWBerger 2/17/00 ;filename for messages msgs = !MSGS+'msgs'+yyddd+'.'+datestring(0) openw,8,msgs printf,8,' *LOST CREEK* FLUX CODE MESSAGES' printf,8,'RUNNING DAYFLUX_lC.PRO FOR START DATE: '+yyddd printf,8,'------------------------------------------------------------------' printf,8,'Run Date and Time: '+systime(0) printf,8,'------------------------------------------------------------------' printf,8,'Filename of this file:' printf,8,msgs printf,8,'------------------------------------------------------------------' printf,8,'RUN SELECTIONS:' printf,8,'(Note: A blank response to a question =n by default)' printf,8,'------' printf,8,'ENTER STARTING DAY FOR DATA PROCESSING IN yyddd FORMAT:' printf,8,yyddd printf,8,'DO YOU WANT TO RUN A WHOLE DAY: (y/n)?' printf,8,wholeday if (wholeday ne 'y') then begin printf,8,'ENTER HALF HOUR (hhmm format) TO BE EXAMINED (e.g., 0830):' printf,8,halfhourstring endif else begin printf,8,'ENTER NUMBER OF DAYS TO PROCESS (Do not overlap Jan. 1):' printf,8,days endelse printf,8,'IS INPUT DATA COMPRESSED: (y/n)?' printf,8,'(Yes will make code read compressed data files.)' printf,8,comp printf,8,'DO YOU WANT TO CALCULATE FLUXES: (y/n)?' printf,8,'(Answering yes overwrites .flx and .sflx files for dates chosen.)' printf,8,calc_flux printf,8,'DO YOU WANT TO USE LONG TERM ROTATION CORRECTIONS: (y/n)?' printf,8,'(Note: .rot files are overwritten for dates chosen regardless of choice.)' printf,8,usetable_rot printf,8, 'DO YOU WANT TO USE FINAL LAG TIME VALUE TABLES: (y/n)?' printf,8,'(Fixed lag times are used if you answer no.)' printf,8,'(Note: .lag files are not affected by this choice.)' printf,8,usetable_lag printf,8, 'DO YOU WANT TO PERFORM SPECTRAL CORRECTIONS: (y/n)?' printf,8,'(Answering yes overwrites .spec files for dates chosen.)' printf,8,calc_speccorr printf,8,'DO YOU WANT TO CALCULATE LAG TIMES: (y/n)?' printf,8,'(Answering yes overwrites .lag files for dates chosen.)' printf,8,calc_lag if calc_lag eq 'y' then begin printf,8,'LAG WINDOWS ARE AS FOLLOWS:' printf,8,'ww & wt lag max, wc lag max, wq lag max (typically zero or negative)' printf,8,'ww & wt lag min, wc lag min, wq lag min (typically negative)' printf,8,lagwindow printf,8,'Change windows by editing dayflux_wc.pro line: ~97' printf,8,'' endif printf,8, 'DO YOU WANT TO USE LICOR CALIBRATION TABLES: (y/n)?' printf,8,usetable_cal printf,8,'DO YOU WANT TO GENERATE HALF HOURLY PLOTS: (y/n)?' printf,8,'(Answering yes overwrites existing hourly plots for dates chosen.)' printf,8,plot_halfhourly if wholeday eq 'y' then begin printf,8,'DO YOU WANT TO GENERATE DAILY PLOTS: (y/n)?' printf,8,'(Answering yes overwrites existing daily plots for dates chosen.)' printf,8,plot_daily endif printf,8,'*****************************************************************' ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;Get date and time elements - yyddd and days are initial string inputs ;make sure input strings are trimmed startyyddd=strtrim(string(yyddd),2) ;string of start date yyddd format days=strtrim(string(days),2) ;string of number of days to process ;get strings of date elements startyy = strmid(startyyddd,0,2) ;string of 2-digit year if (fix(startyy) le 99) and (fix(startyy) ge 50) then begin startyyyy='19'+startyy ;string of 4-digit year endif else begin startyyyy='20'+startyy ;string of 4-digit year endelse startjday=fix(strmid(startyyddd,2,3)) ;integer of starting jdate ;get month and day from jday startyyyymmdd=jdatetodate(startjday,fix(startyyyy)) ;integer startyyyymmdd startmm = strmid(strtrim(string(startyyyymmdd),2),4,2) ;string of 2-digit month startdd = strmid(strtrim(string(startyyyymmdd),2),6,2) ;string of 2-digit day ;get integers of date elements based on start date ;month and day get incremented as necessary if number of days ndays > 1 day = fix(startdd)-1 ;gets incremented in loop below to correct day ;day=0 is ok for getnextdaydate.pro below. month = fix(startmm) year = fix(startyyyy) ;get half hour being processed based on whether operator wants whole day ;or a specific half hour halfhours=['0000','0030','0100','0130','0200','0230','0300','0330',$ '0400','0430','0500','0530','0600','0630','0700','0730',$ '0800','0830','0900','0930','1000','1030','1100','1130',$ '1200','1230','1300','1330','1400','1430','1500','1530',$ '1600','1630','1700','1730','1800','1830','1900','1930',$ '2000','2030','2100','2130','2200','2230','2300','2330'] if wholeday ne 'y' then begin vindex=where(halfhours eq strtrim(halfhourstring,2)) vindex=vindex(0) ;vindex must be a scalar, not an array or ;hrflux_wc will crash ;vindex keeps track of where the selected half hour is located in ;airp, mixr and table lags tbllag_c, tbllag_q. halfhours=[strtrim(halfhourstring,2)] endif ;&&&&&&&&&&&&&&&&&&&&&&&&&& if calc_flux eq 'y' then begin ; Get half hourly pressure and mixing ratio data from qmix ; files. Airp and Mixr are in mb and g h2o/kg dry air, respectively. ; hrmet_wc,startyy,startjday,days,airp,mixr hrmet_lc_psu,startyy,startjday,days,airp,mixr endif ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ;Get lags for the days requested for processing -from table or use fixed ;values depending on usetable_lag (y/n). ;************************ gettablelags_wc,usetable_lag,startyyyy,startjday,days,tbllag_c,tbllag_q ;******************************************************** ;############################################################# ;************************************************************* ;NOW LOOP THROUGH EACH DAY REQUESTED. ;Day and month are incremented as necessary. ;************************************************************* numofdays=0 ;Looping continues while numofdays<=days repeat begin filenotfound: ;a label: where to start again in loop if a file isn't found numofdays=numofdays+1 ;Keeps track of the number of days processed. ;For terminating daily loop after all days are processed ;and for selecting correct element in lag array for hrflux. ;get next days date - year,month,day are input; month and day are ;incremented as necessary and returned in same variables: month and day. ;On first entry into loop day should be set to startday-1 then ;getnextdaydate will increment it to the correct value for the start day. ;Caution: year is not incremented do not overlap Dec 31 and Jan 1. getnextdaydate,year,month,day ;works when day=0 (next day is day=1 but the ;month doesn't change) ; compute julian day julian,year,month,day,jday ;&&&&&&&&&&&&&&&&&&&&&&&&&&& ;Use date being processed in this loop to get date elements to ;build up string YYDDD which is used for filenames, plot titles, etc. ;It is also concatenated later with HHMM to create YYDDDHHMM which is ;also useful. ;get 3-digit jday string if jday lt 10 then begin sjday='00'+strtrim(string(jday),2) endif else begin if jday lt 100 then begin sjday='0'+strtrim(string(jday),2) endif else begin sjday=strtrim(string(jday),2) endelse endelse yyddd=startyy+sjday ;&&&&&&&&&&&&&&&&&&&&&&&&&& ;Create output filenames and open them open=1 openclose_files_wc,open,yyddd,calc_flux,calc_lag,calc_speccorr ;&&&&&&&&&&&&&&&&&&&&&&&&& ;***Get calibration matrix CALMAT for the day ;***(For calculating h2o and h2o mixing ratios). if (usetable_cal eq 'y') then begin getcalmat_lc, year,month,day,calmat endif else begin ;use crude approximations ; calmat=[[.011],$ ;mq ; [-1.74],$ ;bq ; [.20],$ ;mc ; [-49.88]] ;bc ; default value in Sep. calmat=[[.009986],$ ;mq [-12],$ ;bq [.22],$ ;mc [-200]] ;bc endelse print,'calmat=' print,calmat printf,8,'calmat=' printf,8,calmat ;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ;########################################################################### ;*************************************************************************** ;Loop through each half hour - the following major loop is used to load each ; hour, calculate fluxes and statistics, save the ; results, and create postscript files of hour-long ; variables and cross-correlations if requested. ;*************************************************************************** for i=0,n_elements(halfhours)-1 do begin yydddhhmm=yyddd+halfhours(i) ; a string of date & half hour ;in yydddhhmm format ; Get hhmm to extract the proper mean air pressure and mixing ratio data ; (if needed) and for other uses such as getting correct lag out of ; lag array. hhmm=fix(halfhours(i)) ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ;Get data matricies ;********************** print,'-----> reading: ',yydddhhmm printf,8,'-----> reading: ',yydddhhmm ;Read and repair data ; Each variable is returned as a column vector 18000 pts long ; cv and qv are in mV, tl in K, pl in hPa=mb getdata_lc,yydddhhmm,comp,u,v,w,t,cv,qv,tl,pl ;**** ;Clean up data ;screen the sonic data screenu_lc,u,v,w,t ;screen the licor data screenl_lc,cv,qv,tl,pl ;Despike the sonic data despuvwt2,u,v,w,t,6.,badval ;Despike the licor despike,cv,6.,badval despike,qv,6.,badval despike,pl,6.,badval despike,tl,6.,badval ;*** ;Get c and q from licor values using calibrations qlicor_wc,calmat,badval,qv,tl,pl,q clicor_wc,calmat,badval,q,cv,tl,pl,c ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ;Get rotations ;************************* ;Get rotated velocities - table values or calculated values are used ;for rotation based on usetable_rot (usetable_rot = 'y' uses tabular). ;Also get the rotation angle array, rot, to be printed to file later. getrotated_lc,year,jday,usetable_rot,u,v,w,ur,vr,wr,rot,usedrot ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ;Detrend the time series ;********************** detren, ur, dur, mur, trendur, badval detren, vr, dvr, mvr, trendvr, badval detren, wr, dwr, mwr, trendwr, badval detren, t, dt, mt, trendt, badval detren, c, dc, mc, trendc, badval detren, q, dq, mq, trendq, badval ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ;Calculate lag times for lag files (not used for flux calculations) ;output is lag,lagcov,time - the lag time for the peak covariance and ; the lagcov and time arrays ;************************ lagcalc_wc,calc_lag,lagwindow,dur,dvr,dwr,dt,dc,dq,lag,lagcov,lagcovtime ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ;Get fluxes, statistics and determine spectral corrections. ;Fluxes and spectral corrections naturally go well together ;because both require use of lagged time series. ;*********************** ;********* if wholeday eq 'y' then begin ind=i+(numofdays-1)*48L ;index is location in airp,mixr,lagc,lagq ;vectors for specific half hour given by loop index i. ;These vectors have lengths equal to 48x(number of days to be ;processed) ie 48 half hours/day. endif else begin ind=vindex endelse if calc_flux eq 'y' then begin ; Get air pressure and h2o mixing ratio for this half hour. ; These quantities are used to convert flux units airpress=airp(ind) mixrat=mixr(ind) ;specific values for half hour being processed endif lag_c=tbllag_c(ind) lag_q=tbllag_q(ind) hrflux_lc,year,jday,calc_flux,calc_speccorr,airpress,mixrat,$ ur,vr,wr,t,c,q,pl,dur,dvr,dwr,dt,dc,dq,lag_c,lag_q,$ flx,sflx,stats,speccorr ; draw spec plots IF specplot EQ 'y' THEN BEGIN specplot_lc,dq,dc,dt,dwr,yydddhhmm ENDIF ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ;Save necessary quantities to files - ;eg fluxes, rotations, lags, statistics save_wc,calc_lag,calc_flux,calc_speccorr,yydddhhmm,calmat,$ rot,usedrot,lag,lag_c,lag_q,flx,sflx,stats,speccorr ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ if plot_halfhourly eq 'y' then begin ;Plot variables plothourly_lc,yydddhhmm,calc_lag,ur,vr,wr,t,c,q,tl,pl,$ lag,lagcov,lagcovtime endif ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ if plot_daily eq 'y' then begin ;Add hour to daily variables for printing day-long time series later. ;The daily variables are reduced to a data point every 40 points in ;the 5 Hz data. This routine concatinates and uses hr=0 to know if ;input time series is first and therefore isn't concatinated with ;a non-existent previous hour. getdailydata_wc,hhmm,u,v,wr,t,vt,vrh,c,q,tl,pl,$ day_u,day_v,day_wr,day_t,day_c,day_q,day_tl,day_pl endif endfor ;******************************* ;End of hourly loop ;******************************* ;###################################################################### ;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ;Save calmat file for i=0,row(calmat)-1 do begin printf,7,calmat(*,i),format='(5(e15.7,1x))' endfor ;$$$$$$$$$$$$$$$$$$$$ if plot_daily eq 'y' then begin ;Plot day-long variables plotvariables_lc,yyddd,day_u,day_v,day_wr,day_t,day_c,day_q, $ day_tl, day_pl endif ;$$$$$$$$$$$$$$$$$$$$ ;Close all files open=0 openclose_files_wc,open,yyddd,calc_flux,calc_lag,calc_speccorr endrep until(numofdays ge (days)) printf,8,'***END OF DAY (or single half hour)********************************' close,8 ;close code run message file return end