;wlef_process_ch4

;; process_dr_ch4
;;     unzip the folder  (unzip -d DIR file) to a temp folder
;;     for each file (find across subirs all .dat files, process file_ch4
;;        file_search(dir_sepc,recur_pattern)
;;     delete files

;; process_file_ch4
;;     modify all fast to have a new fast file - in it store:  CH4, CO2, Temp, Press
;;     open User file, extract timestamp, (year, doy, hr, mn, sec)
;;     option1: read "user" - interpolate to 10 Hz, option2: user_sync
;;     call standard fast processing    


;convert user->user_sync
;open user file, interpolate to 10 hz, save

;new picarro sensor columns:
                                ;YMD HMS FracDays FracHrs Epoch
                                ;Alarm C UX UY UZ CH4_dry CH4 CO2_dry
                                ;CO2 P T DT ET GPS_L GPS_Lo GPF_Ft
                                ;GPS_RLat GPS_RLon GPS_Tm H2O_s MPV
                                ;Outlet Warmbox

;old picarro sensor colmns
;YMD HMS FracDay FracHrs Epoch Alarm species solenoid DT CP Mode CO2
;CH4 H2O





PRO wlef_processfast_ch4_old,dr,newch4=newch4,nounzip=nounzip,fixold=fixold
;go through directory

  ddr = '/air/incoming/wlef/ch4/'+dr+'/User_Sync/'
  print,'WLEF process fast CH4 for directory '+ddr
  curdoy = -1
  curyr = -1
  wlef_fl = -1

  IF ~keyword_set(nounzip) THEN BEGIN 
    zfile = file_search(ddr+'*.zip',/test_read,count=zc)
    IF zc GE 1 THEN BEGIN 
      IF keyword_set(fixold) AND zc GT 1 THEN BEGIN
        FOR i = 1,zc-1 DO spawn,'unzip -d /air/incoming/wlef/ch4/work/tmp '+zfile[i]
      ENDIF ELSE BEGIN 
        FOR i = 0,zc-1 DO spawn,'unzip -d /air/incoming/wlef/ch4/work/tmp '+zfile[i]
      ENDELSE 
    ENDIF 
    fnames = file_search('/air/incoming/wlef/ch4/work/tmp/','*_User_Sync.dat',count=c,/fully_qualify)
  ENDIF ELSE BEGIN
    fnames = file_search(ddr,'*_User_Sync.dat',count=c,/fully_qualify)
  ENDELSE 

  IF c GT 0 THEN BEGIN
    FOR i = 0,c-1 DO BEGIN
      fname = fnames[i]
      print,'Processing '+fname
      wlef_processfile,fname,12,17,3,curdoy,curyr,wlef_fl,/ch4,newch4=newch4,idch4='Flux'
      IF ~keyword_set(nounzip) THEN spawn,'rm -f '+fname
    ENDFOR 
    IF ~keyword_set(nounzip) THEN spawn,'rm -rf /air/incoming/wlef/ch/work/tmp/*'
  ENDIF 

  wlef_closefast,wlef_fl
END

PRO wlef_processfast_ch4,yr,mo,dy,loaner=loaner,force=force,nosync=nosync,old=old

  if ~keyword_set(nosync) then prefix='User_Sync' else prefix = 'User'
  IF keyword_set(old) THEN ddr_prefix='/air/incoming/' ELSE ddr_prefix='/air/incoming/'

  yrstr = string(yr,format='(i4.4)')
  mostr = string(mo,format='(i2.2)')
  dystr = string(dy,format='(i2.2)')
  IF keyword_set(loaner) THEN BEGIN
    ddr = '/air/incoming/WLEFPicarro/Loaner/DataLog_'+prefix+'/'+yrstr+'/'+mostr+'/'+dystr+'/'
    procfile = '/air/incoming/WLEFPicarro/Loaner/DataLog_'+prefix+'/processed_'+yrstr+'_'+mostr+'_'+dystr
  ENDIF ELSE BEGIN
    ddr = ddr_prefix+'/WLEFPicarro/DataLog_'+prefix+'/'+yrstr+'/'+mostr+'/'+dystr+'/'
    procfile = ddr_prefix+'WLEFPicarro/DataLog_'+prefix+'/processed_'+yrstr+'_'+mostr+'_'+dystr
  ENDELSE 
  IF ~keyword_set(force) AND file_test(procfile,/read) THEN BEGIN
    openr,fl,procfile,/get_lun
    s = ''
    readf,fl,s
    free_lun,fl
    print,'Skipping day ',ddr,' which was already split on ',s
  ENDIF ELSE BEGIN 
    curdoy = -1
    curyr = -1
    wlef_fl = -1
    fnames = file_search(ddr,'*_'+prefix+'.dat',count=c,/fully_qualify)
    print,'WLEF process fast CH4 for directory '+ddr+' with ',c,' files'
    IF c GT 0 THEN BEGIN
      spawn,'scp -i /home/adesai/ameriflux_keys/us-pfa '+ddr+'*_'+prefix+'.dat'+' fluxnet@dtn01.nersc.gov:'
      FOR i = 0,c-1 DO BEGIN
        fname = fnames[i]
        id = strmid((reverse(strsplit(fnames[0],'/',/extract)))[0],5,4)
        print,'  Processing '+fname
        wlef_processfile,fname,12,17,3,curdoy,curyr,wlef_fl,/ch4,idch4=id
      ENDFOR 
    ENDIF 
    wlef_closefast,wlef_fl
    openw,fl,procfile,/get_lun
    printf,fl,systime()
    free_lun,fl
  ENDELSE
END
 
PRO wlef_ch4_usersync,dr

;for 11Picarro139

  IF n_elements(dr) EQ 0 THEN dr = '11Picarro139'

;convert user dir to user sync dir

  ddr = '/air/incoming/wlef/ch4/'+dr+'/User/'
  ddr2 =  '/air/incoming/wlef/ch4/'+dr+'/User_Sync/'
  create_dir,ddr2
;  zfile = file_search(ddr+'*.zip',/test_read,count=zc)
;  IF zc GE 1 THEN BEGIN 
;    zfile = zfile[0]
;    spawn,'unzip -d /air/incoming/wlef/ch4/work/tmp '+zfile
  fnames = file_search(ddr,'Flux-*_User.dat',count=c,/fully_qualify)
  secPerDay = Long(24L*60L*60L)
  day1 = julday(1,1,1970)

  IF c GT 0 THEN BEGIN
    FOR i = 0,c-1 DO BEGIN
      fname = fnames[i]
      print,'Converting '+fname
      openr,fl,fname,/get_lun
      s = ''
      readf,fl,s
      fln = file_lines(fname)-1
      indata = make_Array(18,fln,/double,value=nan())
      ln = 0
      WHILE ~eof(fl) DO BEGIN
        s = ''
        readf,fl,s
        indata[*,ln] = [strsplit(strmid(s,0,20),'-',/extract),strsplit(strmid(s,20),' :',/extract)]
        ln++
      ENDWHILE 
      free_lun,fl
      etime = reform(indata[8,*])
      etimei = numgen(long64(etime[0]*10l)/10.0d,round(etime[n_elements(etime)-1]*10l,/l64)/10.0d,0.1d)
      
      outdata = make_array(18,n_elements(etimei),/double,value=nan())
      jday_foo = double(day1) + etimei/double(secPerDay) - double(12./24.)
      caldat,jday_foo,mm,dd,yyyy,hh,min,ss
      outdata[0:5,*] = [transpose(yyyy),transpose(mm),transpose(dd),transpose(hh),transpose(min),transpose(ss)]
      FOR j = 6,17 DO outdata[j,*] = interpol(indata[j,*],etime,etimei)
      
      newfname = strsplit(fname,'/',/extract)
      newfname = newfname[n_elements(newfname)-1]
      newfname = strmid(newfname,0,strlen(newfname)-4)+'_Sync.dat'
      print,'Writing ',ddr2+newfname
      openw,fl2,ddr2+newfname,/get_lun
      printf,fl2,'DATE                      TIME                      FRAC_DAYS_SINCE_JAN1      FRAC_HRS_SINCE_JAN1       EPOCH_TIME                ALARM_STATUS              species                   solenoid_valves           das_temp                  cavity_pressure           mode_id                   co2_conc_sync             ch4_conc_sync             h2o_conc_sync ',format='(a0)'
      FOR k = 0,n_elements(etimei)-1 DO printf,fl,outdata[*,k],format='(i4,"-",i2.2,"-",i2.2," ",i2.2,":",i2.2,":",f-6.2," ",12f20.6)'
      free_lun,fl2

    ENDFOR
  ENDIF
END 

PRO wlef_ch4_makeslow,yr,startd=startd,endd=endd
;; make_wlef_slow_ch4 - read all available fast files for a year, stitch into CO2 and CH4 time series - hourly average   lefch4_YYYY.txt lefch4YYYY.xdf 
;;    CH4_122_p CO2_122_p, then later CH4 from the three levels (8 columns)
;;   Year DOY HR 

;we need: hourly (yearly)

  outdr = '/air/incoming/wlef/ch4/profile/'

  IF n_elements(yr) EQ 0 THEN yr = 2011
  ystr = string(yr,format='(i4.4)')

  diy = long(days_in_year(yr))
  IF n_elements(startd) EQ 0 THEN IF yr EQ 2010 THEN startd = 270l ELSE startd = 1l
  IF n_elements(endd) EQ 0 THEN endd = diy

  slowfile = outdr+'lefch4_'+ystr+'.xdf'
  IF file_test(slowfile,/read) THEN output = read_xdf(slowfile,header=header) ELSE BEGIN 
    header=['year','doy','hr','co2_122_p','ch4_122_p','ch4_30','ch4_122','ch4_396','co2_30','co2_122','co2_396','co2_30_n','co2_122_n','co2_396_n']
    output = make_array(n_elements(header),8784l,/float,value=nan())
    output[0,*] = yr
    output[1,*] = (lindgen(8784)/24)+1
    output[2,*] = lindgen(8784) MOD 24
  ENDELSE  

  FOR d = startd,endd DO BEGIN
    print,'Year ',yr,' DOY ',d
    dta = wlef_readfast(yr,d,/ch4,header=fh)
    sloc = (d-1)*24l
    eloc = sloc+23l
;['Hr','Min','Sec','T','P','Mode','CO2','CH4','H2O']
    output[3,sloc:eloc] = average_arr(dta[6,*],36000.,/nan)
    output[4,sloc:eloc] = average_arr(dta[7,*],36000.,/nan)

    prof = wlef_readlosgatos(d,yr,header=ph)
    profavg = average_cols(prof,60.,/nan)
    output[5:13,sloc:eloc]=profavg
  ENDFOR
 
  bco2 = wherE(finite(output[3,*]) AND (output[3,*] LT 300 OR output[3,*] GT 600),nbco2)
  IF nbco2 GT 0 THEN output[3,bco2]=nan()
  FOR j = 8,13 DO BEGIN 
    bco2 = wherE(finite(output[j,*]) AND (output[j,*] LT 300 OR output[j,*] GT 600),nbco2)
    IF nbco2 GT 0 THEN output[j,bco2]=nan()
  ENDFOR 
  FOR j = 4,7 DO BEGIN 
    bch4 = where(finite(output[j,*]) AND (output[j,*] LT 1.0 OR output[j,*] GT 3.0),nbch4)
    IF nbch4 GT 0 THEN output[j,bch4]=nan()
  ENDFOR
 
  write_xdf,slowfile,output,header=header
  outf = outdr+'lefch4_'+ystr+'.txt'
  openw,fl,outf,/get_lun
  printf,fl,'YEAR DOY HR CO2_122_p CH4_122_p   CH4_30    CH4_122   CH4_396   CO2_30    CO2_122   CO2_396   CO2_30_N  CO2_122_N CO2_396_N'
  bval = where(~finite(output),nbv)
  IF nbv GT 0 THEN output[bval] = -999.0
  printf,fl,output,format='(i4," ",i3," ",i2.2," ",f9.3," ",f9.3," ",f9.3," ",f9.3," ",f9.3," ",f9.3," ",f9.3," ",f9.3," ",f9.3," ",f9.3," ",f9.3)'
  free_lun,fl
;stop

END

PRO wlef_processflux_ch4,yr,doy,flux=flux,fheader=fheader,diag=diag,dheader=dheader,debug=debug,final=final,finrot=finrot,error=error,eheader=eheader,f_err=f_err
;MAIN PROGRAM

;TEST CODE
  IF n_elements(yr) EQ 0 THEN yr = 2011
  IF n_elements(doy) EQ 0 THEN doy = 160

  print,'WLEF CH4 fluxes for ',yr,' ',doy

;Output headers
  IF n_elements(fheader) EQ 0 THEN BEGIN
    fheader = ['Year','MO','DD','HH','DOY','fDOY',$ ;0-5
               'NEE_CO2_122','NEE_CH4_122','NEE_Q_122']         ;6-8
  ENDIF

  IF n_elements(dheader) EQ 0 THEN BEGIN
    dheader = ['Year','MO','DD','HH','DOY','fDOY',$ ;0-5
               'FLAG1','FLAG2','FLAG3','FLAG4','FLAG5','FLAG6',$ ;6-11 for future use
               'Cstor_122','CH4stor_122','Qstor_122',$ ;12-14
               'Qmix_122','Tv_122','Tair_122',$ ;15-17
               'Press_122','Drho_122','Qrho_122',$ ;18-20
               'Theta_122','Phi_122',$ ;21-22
               'Clag_122','CH4lag_122','Qlag_122',$ ;23-25
               'Cflx_122','CH4flx_122','Qflx_122',$ ;26-28
               'Cspec_122','CH4spec_122','Qspec_122'] ;29-31
  ENDIF

  IF n_elements(eheader) EQ 0 THEN BEGIN
    eheader = ['Year','MO','DD','HH','DOY','fDOY',$ ;0-5
              'Cflux_err','CH4flux_err','LE_err'] ;6-8
  ENDIF 

;Setup output arrays
  flux = make_array(n_elements(fheader),24,value=nan())
  IF ~keyword_set(final) THEN diag =  make_Array(n_elements(dheader),24,value=nan())
  yrstr = string(yr,format='(i4.4)')
  yrdr = yrstr+'/'
  IF keyword_set(final) THEN yrstr = 'final_'+yrstr
  IF n_elements(diag_c) EQ 0 THEN diag_c = read_xdf('/air/incoming/wlef/flux/'+yrdr+'diag_'+yrstr+'.xdf',header=dcheader)
  diag_c = diag_c[*,where(diag_c[4,*] EQ doy)]
  flux[0,*] = yr
  flux[4,*] = doy
  mmdd = jd_to_dy(doy,yr=yr)
  flux[1,*] = float(strmid(mmdd,4,2))
  flux[2,*] = float(strmid(mmdd,6,2))
  flux[3,*] = findgen(24)
  flux[5,*] = doy+(findgen(24)/24.)
  diag[0:5,*] = flux[0:5,*]
  f_err = flux
  ngood = 22000l

  print,'  ',yr,' ',doy,' Reading data' 

;read fast file - 10 Hz data
  fast = wlef_readfast(yr,doy,header=fhead)
  bfast = where(fast LE -999.0,nbf)
  IF nbf gt 0 THEN fast[bfast] = !values.f_nan
  bfast2 = where(fast EQ -99.99,nbf2)
  IF nbf2 gt 0 THEN fast[bfast2] = !values.f_nan

;read Picarro file - 10 Hz
  fast_p =  wlef_readfast(yr,doy,header=fhead_p,/ch4)
  bfast_p = where(fast_p LE -999.0,nbf_p)
  IF nbf_p gt 0 THEN fast_p[bfast_p] = !values.f_nan
  
;read One minute Los Gatos
  lgch4 = wlef_readlosgatos(doy,yr,header=lghead)
  blg = where(lgch4 LE -999,nblg)
  IF nblg GT 0 THEN lgch4[blg] = !values.f_nan

;FAST HEADER Picarro: ['Hr','Min','Sec','T','P','Mode','CO2','CH4','H2O']
;    header = ['Hr','Min','Sec',$ ;0,1,2
;              'LiC_6','LiQ_6','LiT_6','LiC_4','LiQ_4','LiT_4','LiC_2','LiQ_2','LiT_2',$ ;3-11 (trailer)
;              'ST_6','ST_4','ST_2','SP_6','SP_4','SP_2',$ ;12-17 (trailer)
;              'U_6i','V_6i','W_6i','Ts_6i','LiC_6i','LiQ_6i','LiT_6i','ST_6i','SP_6i',$ ;18-26 (top)
;              'U_4i','V_4i','W_4i','Ts_4i','LiC_4i','LiQ_4i','LiT_4i','ST_4i','SP_4i',$ ;27-35 (mid)
;              'U_2i','V_2i','W_2i','Ts_2i','LiC_2i','LiQ_2i','LiT_2i','ST_2i','SP_2i'] ;36-44 (bot)


;Extract variables
  c_122 = fast_p[6,*] ;co2 in ppm
  ch4_122 = fast_p[7,*]*1000. ;convert to ppb
  q_122 = fast_p[8,*]*10.*18./28.97 ;g/kg

;despike u,v,w,Ts - 3 levels
;  print,'  ',yr,' ',doy,' Extracting variables and despiking'
  
  IF (yr EQ 2010 AND doy GE 288) OR (yr EQ 2011 AND doy LT 200) OR (yr EQ 2019 AND doy GT 150) OR (yr GT 2019) THEN BEGIN
    badsonic_122 = where(abs(fast[29,*]) GE 10.0,nbs_122)
  ENDIF ELSE BEGIN 
    badsonic_122 = where((abs(fast[27,*]) GE 35.0) OR (abs(fast[28,*]) GE 35.0) OR (abs(fast[29,*]) GE 10.0) OR (abs(fast[30,*]) GE 45.0),nbs_122)
  ENDELSE 
  IF nbs_122 GT 0 THEN fast[27:30,badsonic_122]=nan()

  u_122 = reform(despike(fast[27,*]))
  v_122 = reform(despike(fast[28,*]))
  w_122 = reform(despike(fast[29,*]))
  t_122 = reform(despike(fast[30,*]))

  c_122 = reform(c_122)
  ch4_122 = reform(ch4_122)
  q_122 = reform(q_122)
  IF total(abs(q_122),/nan) EQ 0 THEN q_122[*] = nan()

;hourly mixing ratio, virtual temperature, air temp, pressure, air density,vapor density
;just read these in from DIAG!

;READ DIAG FILE, extract variables
;qmix_122,tv_122,tair_122,press_122,drho_122,qrho_122
  qmix_122 = diag_c[22,*]
  tv_122 = diag_c[25,*]
  tair_122 = diag_c[28,*]
  press_122 = diag_c[31,*]
  drho_122 = diag_c[34,*]
  qrho_122 = diag_c[37,*]
  diag[15:20,*] = [qmix_122,tv_122,tair_122,press_122,drho_122,qrho_122]

;Storage
  cstor = diag_c[13,*]
  qstor = diag_c[16,*]
  ch4_fl = lgch4[0:2,*]*1000.
  ch4_122f = average_arr(ch4_122,600.,/nan)
  FOR i = 0,2 DO ch4_fl[i,*] = merge_array(ch4_fl[i,*],ch4_122f)
  ch4stor = (wlef_storage(ch4_fl))[1,*]

  co2_fl =   lgch4[3:5,*]
  co2_122f = average_arr(c_122,600.,/nan)
  FOR i = 0,2 DO co2_fl[i,*] = merge_array(co2_fl[i,*],co2_122f)
  co2stor = (wlef_storage(co2_fl))[1,*]
  IF yr EQ 2011 AND doy GE 158 AND doy LE 180 THEN cstor = merge_array(cstor,co2stor) ELSE cstor = merge_Array(co2stor,cstor)

  q_fl =   make_array(3,1440,/float,value=nan())
  q_122f = average_arr(q_122,600.,/nan)
  FOR i = 0,2 DO q_fl[i,*] = q_122f
  qstor2 = (wlef_storage(q_fl))[1,*]
  qstor = merge_Array(qstor,qstor2)

  diag[12:14,*] = [cstor,ch4stor,qstor]

;For each hour
  allbad = replicate(!values.f_nan,36000l)
;  FOR i = 18l,18l DO BEGIN 
  FOR i = 0l,23l DO  BEGIN 
    hstart = i*36000l
    hend = ((i+1l)*36000l)-1l

    print,'  ',yr,' ',doy,' Processing hour ',i
;    print,'  ',yr,' ',doy,' Rotating U,V,W for hour ',i
;Rotate u,v,w (3 levels)
;  calc_rot,u,v,w,ur=ur,vr=vr,wr=wr,finrot=finrot,rot_tp=rot_tp,rot_ntp=rot_ntp,therot=therot,rotflag=rotflag
;fix for bad u in Oct 14, 2010 122m sonic
    IF (yr EQ 2010 AND doy GE 288) OR (yr GE 2011 AND doy LT 200) THEN BEGIN 
      tmp = where(finite(w_122[hstart:hend]),nuvw122) 
    ENDIF ELSE BEGIN  
      tmp = where(finite(u_122[hstart:hend]) AND finite(v_122[hstart:hend]) AND finite(w_122[hstart:hend]),nuvw122)
    ENDELSE 

;FINAL CODE HERE (finrot subset)
    rot_122 = [nan(),nan()]
    ur_122 = make_array(18000,/float,value=nan())

    IF keyword_set(final) AND n_elements(finrot) NE 0 THEN finrot122 = finrot[0:5]

    IF nuvw122 GT ngood THEN BEGIN
      IF (yr EQ 2010 AND doy GE 288) OR (yr EQ 2011 AND doy LT 200) THEN BEGIN 
        wr_122 = w_122[hstart:hend] ;future - use long term rotation
        ur_122 = u_122[hstart:hend]
        vr_122 = v_122[hstart:hend]
      ENDIF ELSE BEGIN 
        calc_rot,u_122[hstart:hend],v_122[hstart:hend],w_122[hstart:hend],ur=ur_122,vr=vr_122,wr=wr_122,therot=rot_122,/quiet,finrot=finrot122
      ENDELSE 
    ENDIF 
    diag[21:22,i] = [rot_122]

;    print,'  ',yr,' ',doy,' Detrending for hour ',i    
;detrend data and despike hourly data
    IF nuvw122 GT ngood THEN BEGIN
      dur_122 = despike(detrend(ur_122),/nodetrend)
      dvr_122 = despike(detrend(vr_122),/nodetrend)
      dwr_122 = despike(detrend(wr_122),/nodetrend)
    ENDIF ELSE BEGIN
     dur_122 = allbad & dvr_122 = allbad & dwr_122 = allbad
   ENDELSE

   dc_122 = despike(detrend(c_122[hstart:hend]),/nodetrend)
   dq_122 = despike(detrend(q_122[hstart:hend]),/nodetrend)
   dch4_122 = despike(detrend(ch4_122[hstart:hend]),/nodetrend)
   dt_122 = despike(detrend(t_122[hstart:hend]),/nodetrend)

   tmp = where(finite(dc_122),ndc_122)
   tmp = where(finite(dch4_122),ndch4_122)
   tmp = where(finite(dq_122),ndq_122)
   tmp = where(finite(dt_122),ndt_122)

;FINAL CODE HERE (don't calc lags, take lags)
;    print,'  ',yr,' ',doy,' Computing lags hour ',i
;Compute lags for c and q (5 levels)
    lag_122 = [nan(),nan()]
    lag_122_ch4 = [nan(),nan()]

    IF keyword_set(final) THEN BEGIN
      lag_122 = diag[[23,25],i]
      lag_122_ch4 = diag[24:25,i]
      wlef_applylag,lag_122,dwr_122,dc_122,dq_122,dt_122,level=122,ldwr_c=ldwr_c_122,ldt_c=ldt_c_122,ldc=ldc_122,ldwr_q=ldwr_q_122,ldt_q=ldt_q_122,ldq=ldq_122
      wlef_applylag,lag_122_ch4,dwr_122,dch4_122,dq_122,dt_122,level=122,ldwr_c=ldwr_ch4_122,ldt_c=ldt_ch4_122,ldc=ldch4_122
    ENDIF ELSE BEGIN 
      lag_122 = wlef_lag(dwr_122,dc_122,dq_122,dt_122,level=122,ldwr_c=ldwr_c_122,ldt_c=ldt_c_122,ldc=ldc_122,ldwr_q=ldwr_q_122,ldt_q=ldt_q_122,ldq=ldq_122)
      lag_122_ch4 = wlef_lag(dwr_122,dch4_122,dq_122,dt_122,level=122,ldwr_c=ldwr_ch4_122,ldt_c=ldt_ch4_122,ldc=ldch4_122)
      diag[23:25,i] = [lag_122[0],lag_122_ch4[0],lag_122[1]]
    ENDELSE
    
;    print,'  ',yr,' ',doy,' Computing fluxes for hour ',i
    uflx_30 = nan()
    qflx_122 = uflx_30 
    cflx_122 = uflx_30 & ch4flx_122 = uflx_30
    tvflx_122 = uflx_30

;Compute H and momentum flux (3 levels) - no lags
    dwrs_122 = shift(reform(dwr_122),1)
    IF ndt_122 GT ngood THEN tvflx_122 = covar(dt_122,dwrs_122)
    diag[6,i] = tvflx_122

;compute CO2, CH4, H2O flux
    IF ndq_122 GT ngood THEN qflx_122 = covar(ldq_122,ldwr_q_122)
    IF ndc_122 GT ngood THEN cflx_122 = covar(ldc_122,ldwr_c_122)
    IF ndch4_122 GT ngood THEN ch4flx_122 = covar(ldch4_122,ldwr_ch4_122)

    diag[26:28,i] = [cflx_122,ch4flx_122,qflx_122]


    IF keyword_set(error) THEN BEGIN
;      print,'Computing flux error'
       IF nuvw122 GT ngood THEN u_bar = mean(abs(ur_122),/nan) ELSE u_bar = 10.0
      IF ~finite(u_bar) THEN u_bar = 10.0
      cvar = ldc_122*ldwr_c_122
      fcvar = where(finite(cvar),ncvar)
      IF ncvar GT ngood THEN cflx_err = fluxerr(cvar[fcvar],freq=10,z=122.0,U_bar=U_bar,avg_t=3600.0,nfilters=10) ELSE cflx_err = nan()
      c4var = ldch4_122*ldwr_ch4_122
      fc4var = where(finite(c4var),nc4var)
      IF nc4var GT ngood THEN ch4flx_err = fluxerr(c4var[fc4var],freq=10,z=122.0,U_bar=U_bar,avg_t=3600.0,nfilters=10) ELSE ch4flx_err = nan()
      qvar = ldq_122*ldwr_q_122
      fqvar = where(finite(qvar),nqvar)
      IF nqvar GT ngood THEN qflx_err = fluxerr(qvar[fqvar],freq=10,z=122.0,U_bar=U_bar,avg_t=1800.0,nfilters=10) ELSE qflx_err = nan()
      cflux_err = nan() & ch4flux_err = nan() & qflux_err = nan()
      convert_flux,drho_122[i],qrho_122[i],tair_122[i],cflx=cflx_err,qflx=qflx_Err,cflux=cflux_err,qflux=qflux_err
      convert_flux,drho_122[i],qrho_122[i],tair_122[i],cflx=ch4flx_err,cflux=ch4flux_err
      f_err[6:8,i] = [cflux_err,ch4flux_err,qflux_err]
 ;     stop
    ENDIF 


;FINAL CODE HERE - USE SPECS
;    print,'  ',yr,' ',doy,' Spectral correction for hour ',i
;spectral corrections - need to fix constants to make this work
    IF isallnan(ldt_c_122) THEN ldt_c_122=ldc_122
    IF isallnan(ldt_q_122) THEN ldt_q_122=ldq_122
    IF isallnan(ldt_ch4_122) THEN ldt_ch4_122=ldch4_122

    IF keyword_set(final) THEN BEGIN
      spec122 = diag[[29,31],i]
      spec122_ch4 = diag[30:31,i]
    ENDIF ELSE BEGIN 
      spec122 = wlef_spec(ldq_122,ldc_122,ldt_q_122,ldt_c_122,t_122,ldwr_q_122,ldwr_c_122,lag_122[0]*(-1.0),level=122.,tmp=tv_122[i]-273.15)
      spec122_ch4 = wlef_spec(ldq_122,ldch4_122,ldt_q_122,ldt_ch4_122,t_122,ldwr_q_122,ldwr_ch4_122,lag_122_ch4[0]*(-1.0),level=122.,tmp=tv_122[i]-273.15)
      diag[29:31,i] = [spec122[0],spec122_ch4[0],spec122[1]]
    ENDELSE 

    IF finite(spec122[0]) THEN cflx_122*=spec122[0]
    IF finite(spec122[1]) THEN qflx_122*=spec122[1]
    IF finite(spec122_ch4[0]) THEN ch4flx_122*=spec122_ch4[0]

;Compute NEE
;    print,'  ',yr,' ',doy,' Outputting for hour ',i
    qnee_122 = qflx_122 +  qstor[i]
    nee_122 = cflx_122 +  cstor[i]
    nee_ch4_122 = ch4flx_122 +  ch4stor[i]

;convert units to umol/m2/s and W/m2
    convert_flux,drho_122[i],qrho_122[i],tair_122[i],cflx=nee_122,qflx=qnee_122,cflux=cflux_122,qflux=qflux_122
    convert_flux,drho_122[i],qrho_122[i],tair_122[i],cflx=nee_ch4_122,cflux=ch4flux_122

    flux[6:8,i] = [cflux_122,ch4flux_122,qflux_122]

;    IF keyword_set(debug) THEN stop

  ENDFOR 

  if keyword_set(debug) then stop

END

PRO wlef_processfluxyear_ch4,yr,error=error,startday=startday,endday=endday,email=email
  
  yrstr = string(yr,format='(i4.4)')
  yrdr = yrstr+'/'
  IF keyword_set(final) THEN yrstr='final_'+yrstr

  IF n_elements(yr) EQ 0 THEN yr=2010
  diy = long(days_in_year(yr))
  
  IF n_elements(startday) EQ 0 THEN startday = 1l
  IF n_elements(endday) EQ 0 THEN endday = diy

  IF file_test('/air/incoming/wlef/flux/'+yrdr+'fluxch4_'+yrstr+'.xdf',/read) THEN f = read_xdf('/air/incoming/wlef/flux/'+yrdr+'fluxch4_'+yrstr+'.xdf',header=fheader)
  IF file_test('/air/incoming/wlef/flux/'+yrdr+'diagch4_'+yrstr+'.xdf',/read) THEN d = read_xdf('/air/incoming/wlef/flux/'+yrdr+'diagch4_'+yrstr+'.xdf',header=dheader)
  IF file_test('/air/incoming//wlef/flux/'+yrdr+'fluxch4err_'+yrstr+'.xdf',/read) AND keyword_set(error) THEN e = read_xdf('/air/incoming//wlef/flux/'+yrdr+'fluxch4err_'+yrstr+'.xdf',header=eheader)

  FOR doy = startday,endday DO BEGIN 
    wlef_processflux_ch4,yr,doy,flux=flux,fheader=fheader,diag=diag,dheader=dheader,error=error,eheader=eheader,f_err=f_err
    IF n_elements(f) EQ 0 THEN f = make_array(n_elements(fheader),diy*24l,/float,value=nan())
    IF n_elements(d) EQ 0 THEN d = make_array(n_elements(dheader),diy*24l,/float,value=nan())
    IF keyword_set(error) AND n_elements(e) EQ 0 THEN e = make_array(n_elements(eheader),diy*24l,/float,value=nan())
 
    lstart = (doy-1l)*24l
    lend = (doy*24l)-1l
    f[*,lstart:lend]=flux
    d[*,lstart:lend]=diag
    IF keyword_set(error) THEN e[*,lstart:lend]=f_err
    IF doy MOD 30 EQ 0 THEN BEGIN 
      write_xdf,'/air/incoming/wlef/flux/'+yrdr+'fluxch4_'+yrstr+'.xdf',f,header=fheader
      write_xdf,'/air/incoming/wlef/flux/'+yrdr+'diagch4_'+yrstr+'.xdf',d,header=dheader
      IF keyword_set(error) THEN write_xdf,'/air/incoming/wlef/flux/'+yrdr+'fluxch4err_'+yrstr+'.xdf',e,header=eheader
    ENDIF 

  ENDFOR
 
  write_xdf,'/air/incoming/wlef/flux/'+yrdr+'fluxch4_'+yrstr+'.xdf',f,header=fheader
  write_xdf,'/air/incoming/wlef/flux/'+yrdr+'diagch4_'+yrstr+'.xdf',d,header=dheader
  IF keyword_set(error) THEN write_xdf,'/air/incoming/wlef/flux/'+yrdr+'fluxch4err_'+yrstr+'.xdf',e,header=eheader

  IF keyword_set(email) THEN BEGIN
;    openw,fl,'/air/incoming/wlef/flux/lastday_ch4.txt',/get_lun
    openw,fl,'/air/incoming/wlef/lastday.txt',/get_lun,/append
;    printf,fl,'CC: Jonathan E. Thom <jthom@ssec.wisc.edu>'
;    printf,fl,'Subject: WLEF CH4 Diagnostic '+string(yr,format='(i4.4)')+' DOY '+string(Endday,format='(i3.3)')
;    printf,fl,'MIME-Version: 1.0'
;    printf,fl,'Content-Type: text/html'
;    printf,fl,'Content-Disposition: inline'
;    printf,fl,'<html>'
;    printf,fl,'<body>'
;    printf,fl,'<pre style="font-family: consolas,monospace;font-size:8pt">'
    printf,fl,''
    printf,fl,'WLEF FLUX CH4'
    totvar = float(lend-lstart)+1.0
    skiphead = ['NEE_Q_122']
    FOR i = 0,n_elements(fheader)-1 DO IF ~in(skiphead,fheader[i]) THEN IF total(~finite(f[i,lstart:lend])) GT (totvar/1.8) THEN printf,fl,fheader[i],' ',100*total(~finite(f[i,lstart:lend]))/totvar,'% missing'
;    printf,fl,fheader,format='(a4," ",a2," ",a2," ",a2," ",a3," ",a7," ",2(a10," "),a10)'
;    FOR i=lstart,lend DO printf,fl,f[*,i],format='(i4," ",i2," ",i2," ",i2," ",i3," ",f7.3," ",2(f10.1," "),f10.1)'
    printf,fl,''
    printf,fl,'WLEF DIAG CH4'
    skiphead = ['FLAG2','FLAG3','FLAG4','FLAG5','FLAG6','Qlag_122','Qflx_122']
    FOR i = 0,n_elements(dheader)-1 DO IF ~in(skiphead,dheader[i]) THEN IF total(~finite(d[i,lstart:lend])) GT (totvar/1.8) THEN printf,fl,dheader[i],' ',100*total(~finite(d[i,lstart:lend]))/totvar,'% missing'
;    printf,fl,dheader,format='(a4," ",a2," ",a2," ",a2," ",a3," ",a7," ",25(a12," "),a12)'
;    FOR i=lstart,lend DO printf,fl,d[*,i],format='(i4," ",i2," ",i2," ",i2," ",i3," ",f7.3," ",25(f12.2," "),f12.2)'
;    printf,fl,'</pre>'
;    printf,fl,'</body>'
;    printf,fl,'</html>'
    free_lun,fl
;    spawn,'/usr/sbin/sendmail professorankurdesai@gmail.com < /air/incoming/wlef/flux/lastday_ch4.txt'
  ENDIF

END

PRO wlef_finflx_ch4,yr,error=error
  yrstr = string(yr,format='(i4.4)')
  yrdr = yrstr+'/'
  yrstrf = 'final_'+yrstr
  IF n_elements(yr) EQ 0 THEN yr=2011
  diy = long(days_in_year(yr))

  f = read_xdf('/air/incoming/wlef/flux/'+yrdr+'fluxch4_'+yrstr+'.xdf',header=fh)
  d = read_xdf('/air/incoming/wlef/flux/'+yrdr+'diagch4_'+yrstr+'.xdf',header=dh)
  IF keyword_set(error) THEN BEGIN
    ename = '/air/incoming/wlef/flux/'+yrdr+'fluxch4err_'+yrstr+'.xdf'
    IF file_test(ename,/read) THEN BEGIN
      e = read_xdf('/air/incoming/wlef/flux/'+yrdr+'fluxch4err_'+yrstr+'.xdf',header=eh)
    ENDIF ELSE BEGIN
      eh = ['Year','MO','DD','HH','DOY','fDOY','Cflux_err','CH4flux_err','LE_err']  
      e = make_array(n_elements(eh),diy*24l,/float,value=nan())
    ENDELSE 
  ENDIF 

  ff = f
  df = d
  IF keyword_set(error) THEN ef = e

;1. final rotation
  finrot122 = wlef_finalrot(d[21,*],d[22,*])

;2. Screen and fill lag times - despike, then median per day

  FOR j = 23,25 DO BEGIN
    daily = average_arr(d[j,*],24,/nan,/med)
    daily = merge_Array(daily,amedian(daily,7))
    bd = where(~finite(daily),nbd)
    IF nbd GT 0 THEN daily[bd] = median(daily)
    df[j,*] = expand_arr(daily,24)
  ENDFOR 

;3. Screen and fill spectral corrections - despike
  FOR j = 29,31 DO BEGIN
    bd = where(d[j,*] EQ 1 OR ~finite(d[j,*]),nbd)
    IF nbd GT 0 THEN df[j,bd] = nan()
    daily = average_arr(df[j,*],24,/nan,/med)
    daily = merge_Array(daily,amedian(daily,7))
    bd = where(~finite(daily),nbd)
    IF nbd GT 0 THEN daily[bd] = median(daily)
    df[j,*] = expand_arr(daily,24)>1
  ENDFOR 

  FOR doy = 1l,diy DO BEGIN
    lstart = (doy-1l)*24l
    lend = (doy*24l)-1l
    diag = df[*,lstart:lend]
    wlef_processflux_ch4,yr,doy,flux=flux,diag=diag,finrot=finrot,/final,error=error,f_err=f_err
    ff[*,lstart:lend]=flux
    df[*,lstart:lend]=diag
    IF keyword_set(error) THEN ef[*,lstart:lend]=f_err
  ENDFOR 

;write data
  write_xdf,'/air/incoming/wlef/flux/'+yrdr+'fluxch4_'+yrstrf+'.xdf',ff,header=fh
  write_xdf,'/air/incoming/wlef/flux/'+yrdr+'diagch4_'+yrstrf+'.xdf',df,header=dh
  IF keyword_set(error) THEN write_xdf,'/air/incoming/wlef/flux/'+yrdr+'fluxch4err_'+yrstrf+'.xdf',ef,header=eh
END


PRO wlef_ch4_roxana
;read hourly fill NEE and CH4 flux, convert to useful things
  f = read_xdf('/air/incoming/wlef/flux/2011/prefnee_2011.xdf',header=hf)

  nee = f[10,*]
  gpp = f[11,*]
  re = f[12,*]
  temp = f[14,*]
  par = f[15,*]
  ch4 = f[20,*]

  par[where(~finite(par))]=0.0

  cfill = nee
  warm = where(temp GT 0,complement=cold)
  cfill[warm] = cfill[warm]*.594416 + 2.95902
  cfill[cold] = cfill[cold]*0.0972571 + 2.52475

  ch4f = merge_array(ch4,cfill)

  outf1 = make_array(12,8760,/float,value=nan())
  outf1[0:5,*] = f[0:5,*]
  outf1[6,*] = temp
  outf1[7,*] = par
  outf1[8,*] = nee
  outf1[9,*] = gpp
  outf1[10,*] = re
  outf1[11,*] = ch4f

  outf2 = outf1
  outf2[8:11,*] = outf2[8:11,*]*12.0*1e-6*3600.0

  outd = make_array(8,365,/float,value=nan())
  outd[0,*] = f[0,0]
  outd[1,*] = lindgen(365)+1
  outd[2,*] = average_arr(outf2[6,*],24,/nan)
  outd[3,*] = average_arr(outf2[7,*],24,/nan)
  outd[4,*] = average_arr(outf2[8,*],24,/tot,/nan)
  outd[5,*] = average_arr(outf2[9,*],24,/tot,/nan)
  outd[6,*] = average_arr(outf2[10,*],24,/tot,/nan)
  outd[7,*] = average_arr(outf2[11,*],24,/tot,/nan)

  openw,fl,'~/outf1.txt',/get_lun
  FOR i = 0,8759 DO printf,fl,outf1[*,i],format='(i4," ",i2," ",i2," ",i2," ",i3," ",f7.3," ",6(f10.1," "))'
  free_lun,fl
  openw,fl,'~/outf2.txt',/get_lun
  FOR i = 0,8759 DO printf,fl,outf2[*,i],format='(i4," ",i2," ",i2," ",i2," ",i3," ",f7.3," ",6(f10.1," "))'
  free_lun,fl
  openw,fl,'~/outd.txt',/get_lun
  FOR i = 0,364 DO printf,fl,outd[*,i],format='(i4," ",i3," ",6(f10.1," "))'
  free_lun,fl

  stop

;gap fill CH4 (T-fit?)


;output csv in:
;umol/m2/s nmol/m2/s

;g/m2/hour mg/m2/hour

;daily file in g/m2/day  mg/m2/day

;Timestamp + T PAR NEE_CO2(umol/m2/s) NEE_CH4(nmol/m2/s) NEE



END