;finalflux_wc.pro takes concatenated output from dayflux and calculates final
;fluxes for Willow Creek. This involves correcting the fluxes for spectral 
;degradation and converting the kinematic units of the fluxes into more 
;common units.
;
;Use is as follows:
;                         finalflux_wc
;
;The input files and where they come from are as follows, where YYYY is the 
;4-digit year:
;            lefYYYY.flx    |Created by create_annual.pro from .flx and .stats
;            lefYYYY.stats  |files output from dayflux.pro.
;    location: Usually in same location as output of dayflux_wc 
;                        (see compileflux_wc_yourname)
;              Operator must set this location at the load statements below.
;------------------------------------------------------------------------------
;      spectableYYYY.L       |Created by finalspectable.pro from .spec
;                            |files output from dayflux.pro.
;    location: /data/davis/cheas/wcreek/data/spec/
;------------------------------------------------------------------------------
; 
;Note: The output of this code must be modified to place files where operator
;      wants to put it.
;
;The final flux files have the following column format:
;yyyy mm dd hh jday fjday cflux(umol/(m^2*s)) cflux_corr 
;qflux(W/m^2) qflux_corr tflux(W/m^2) tflux_corr u*(m/s) mflux(Pa)
;
;A row exists for each hour. The number of days and hours is based on the
;common block of time present in the lefYYYY.flx, lefYYYY.stats and 
;spectableYYYY files.
;
;Written by BWBerger 9/21/99 for wlef. Modified to be used for willow creek
; ~12/99 by WWang.
;Modified by BWB 2/23/00 to make some input files come from
;          /data/davis/cheas/wcreek/data/
;instead of /data/davis2/wcreek/data/tables/
; wwg, Sep 12 2000, change directies

pro finalflux_lc

badval=-999.

print,'Enter 4-digit year (yyyy) to process:'
yyyy=''
read,yyyy

;get output filename
;outfile='/data/davis2/wwg/wcreek/output/flx/wc'+yyyy+'.finflx'
outfile='/eddy/s2/wang/d2wwg/lcreek/output/flx/lc'+yyyy+'.finflx'

;load the data files
;load,'/data/davis2/wwg/wcreek/output/flx/wc'+yyyy+'.flx',flx
;load,'/data/davis/cheas/wcreek/data/spec/spectable'+yyyy,spec
;load,'/data/davis2/wwg/wcreek/output/stats/wc'+yyyy+'.stats',stats
; changed to use in PSU,

load,'/eddy/s2/wang/d2wwg/lcreek/output/flx/lc'+yyyy+'.flx',flx
;change s2 to s4 on 02/21/2003
load,'/eddy/s4/lcreek/data/spec/spectable'+yyyy,spec
;load,'/davis/s2/lcreek/data/spec/spectable'+yyyy,spec
load,'/eddy/s2/wang/d2wwg/lcreek/output/stats/lc'+yyyy+'.stats',stats

;get min and max dates in input files
minflx=jdate2(flx(0,0),flx(1,0),flx(2,0))
maxflx=jdate2(flx(0,row(flx)-1),flx(1,row(flx)-1),flx(2,row(flx)-1))
minspec=spec(0,0)
maxspec=spec(0,row(spec)-1)
minstats=jdate2(stats(0,0),stats(1,0),stats(2,0))
maxstats=jdate2(stats(0,row(stats)-1),stats(1,row(stats)-1),stats(2,row(stats)-1))

mindate=max([minflx,minspec,minstats]) ;common start date
maxdate=min([maxflx,maxspec,maxstats]) ;common end date
ndays=maxdate-mindate+1

;Get hourly pressure and mixing ratio data for Willow Creek.
;airpress and qmixratio are in mb and g h2o/kg dry air, respectively.
startyy=strmid(yyyy,2,2); string of 2-digit year
days=strtrim(string(ndays),2); string of number of days to be processed
hrmet_lc_psu,startyy,mindate,days,airpress,qmixratio

;trim data arrays as necessary - based on  mindate and maxdate
flx=flx(*,(mindate-minflx)*24*2L:(mindate-minflx)*24*2L+ndays*24*2L-1)
spec=spec(*,(mindate-minspec)*24*2L:(mindate-minspec)*24*2L+ndays*24*2L-1)
stats=stats(*,(mindate-minstats)*24*2L:(mindate-minstats)*24*2L+ndays*24*2L-1)

n=ndays*24*2L ;the number of hours in period common to all data sets

;initialize array to hold final fluxes
flux=fltarr(14,n)

;**************************************************************
;Loop through each hour to:
;Convert kinematic units to common units
;Apply spectral correction

for i=0L,n-1 do begin

   ;get necessary quantities
   yyyy=flx(0,i)
   mm=flx(1,i)
   dd=flx(2,i)
   hh=flx(3,i)
   jday=flx(4,i)
   fjday=flx(5,i)
   cflx=flx(6,i)
   qflx=flx(7,i)
   if (cflx ne badval) and (spec(2,i) ne badval) then begin
      cflx_corr=flx(6,i)*spec(2,i) ;spectrally corrected co2 flux
   endif else begin
      cflx_corr=badval
   endelse
   if (qflx ne badval) and (spec(3,i) ne badval) then begin
      qflx_corr=flx(7,i)*spec(3,i) ;spectrally corrected h2o flux
   endif else begin
      qflx_corr=badval
   endelse
   tvflx=flx(8,i)
   mflx=flx(9,i)
   tvmean=stats(14,i)
   qmean=stats(18,i)
   airp=airpress(i)
   qmix=qmixratio(i)

   ;calculate u* from mflx
   if (mflx ne badval) then begin
      ustar=sqrt(abs(mflx))
   endif else begin
      ustar=badval
   endelse

   ;Compute temperature flux rather than virtual temperature flux.

   ;First convert mean virtual temperature (VIRTEMP-K) for the period 
   ;of the flux calculation to mean temperature (AIRTEMP-K).
   if (tvmean ne badval) then begin
      virtemp=tvmean+273.15 ;convert virtual temp. to K
   end else begin
      virtemp=badval
   endelse
   if qmix eq badval then begin ;See if hrmet2 mixing ratio is bad.
      if qmean ne badval then begin
         qmix=qmean
         print,'H2O mixing ratio bad, jday:',jday,' hr:',hh
         print,'Replaced with Licor value.'
         print,'qmixat=',qmix
      endif else begin ;If bad, replace by licor value.
         print,'H2O mixing ratio bad, jday:',jday,' hr:',hh
         qmix=badval
      endelse
   endif
   virtotemp,qmix,virtemp,airtemp,badval ;output is airtemp
   ;Note: Unless otherwise specified, badvals are returned by routines
   ;      if badvals are input.

   ; Next compute temperature flux (TFLX-K*m/s) from virtual temperature
   ; flux TVFLX - K*m/s. (First to get non-spectrally-corrected fluxes,then
   ; to get spectrally-corrected fluxes.)
   virtotempflx,qmix,virtemp,qflx,tvflx,tflx,badval ;output is tflx
   virtotempflx,qmix,virtemp,qflx_corr,tvflx,tflx_corr,badval

   ; Compute dry air density (DRHO - kg/m^3) averaged over the time of
   ; the flux calculation
   ddensity,airp,qmix,airtemp,drho,badval ;output is drho

   ; Next compute water vapor density (QRHO - kg/m^3) averged over the
   ; time of the flux calculation.
   qdensity,airp,qmix,airtemp,qrho,badval  ;output is qrho

   ; Convert flux units (First for non-spectrally-corrected fluxes,then
   ; for spectrally-corrected fluxes.)
   convert,drho,qrho,airtemp,cflx,qflx,tflx,mflx,$ ;output is cflux,qflux
           cflux,qflux,tflux,mflux,badval          ;tflux,mflux
   convert,drho,qrho,airtemp,cflx_corr,qflx_corr,tflx_corr,mflx,$
           cflux_corr,qflux_corr,tflux_corr,mflux_corr,badval
           ;mflux_corr should equal mflux

   ;load-up final flux array
   flux(*,i)=[yyyy,mm,dd,hh,jday,fjday,cflux,cflux_corr,qflux,qflux_corr,tflux,tflux_corr,ustar,mflux]

endfor

;write final flux array to file
openw,fout,outfile,/get_lun
;printf,fout,'yyyy mm dd hh jday fjday cflux(umol/(m^2*s)) cflux_corr qflux(W/m^2) qflux_corr tflux(W/m^2) tflux_corr u*(m/s) mflux(Pa)'
for i=0L,n-1 do begin
   printf,fout,flux(*,i),format='(I4,1x,I2,1x,I2,1x,f4.1,1x,I3,1x,e15.7,1x,8(e15.7,1x))'
endfor
close,fout
free_lun,fout


return
end