FUNCTION read_filledmet_other,yr,site=site,header=header

  IF n_elements(site) EQ 0 THEN site = 'wcreek'
  dir = co2_basedir('othersites')+site+'/'
  yrstr = string(yr MOD 100,format='(i2.2)')
  fname = dir+site+yrstr+'_filledmet.xdf'
  dta = read_xdf(fname,header=header)
  return,dta
END

PRO create_resppsyn_wcreek2,yr,kind=kind,outdir=outdir
;this time use the flag in the flux file 
;and use the original flux in the flux file
;and use par/temp from resppsyn file

  site = 'wcreek'
  readfluxtext,yr,site=site,flux=f,met=m,fhead=fh,mhead=mh
  read_resppsyn,yr,altdir=site,fillednee=fno,filledheader=fnh

  flag = f[12,*]
;  lowustar = f[12,*] EQ 1
  lowustar = f[11,*] LT 0.3
;isnan fno, but notnan fn, and not lowustar

  sewind = isnan(fno[14,*]) AND isnotnan(f[8,*]) AND (~lowustar)

;  sewind = f[12,*] EQ 2
  onee1 = fno[14,*]
  onee2 = f[8,*]
  onee = onee1
  onee[where(lowustar)] = onee2[where(lowustar)]
  onee[where(sewind)] = onee2[where(sewind)]
;  onee = fno[14,*]
  temp = fno[7,*]
  par = fno[12,*]

  IF n_elements(kind) EQ 0 THEN kind = 3
  CASE kind OF
    0 : badlocs = flag EQ 20
    1 : badlocs = lowustar
    2 : badlocs = sewind
    3 : badlocs = lowustar OR sewind
    ELSE : badlocs = flag EQ 20
  ENDCASE 

  lat = 45.81
  lon = -90.08
  CalcRespPsynFits,yr,nee=onee,par=par,temp=temp,respfit=rf,psynfit=pf,fillednee=fn,badval=nan(),badlocs=badlocs,resphead=rh,psynhead=ph,filledhead=fnh,lat=lat,lon=lon,/quiet

  print,1.0368*total(fno[8,*])/48.0,1.0368*total(fno[13,*])/48.0,1.0368*total(fno[16,*])/48.0
  print,1.0368*total(fn[8,*])/48.0,1.0368*total(fn[13,*])/48.0,1.0368*total(fn[16,*])/48.0

  stop

  IF n_elements(outdir) EQ 0 THEN BEGIN 
    CASE kind OF
      0 : outdir = 'noscreen'
      1 : outdir = 'ustar'
      2 : outdir = 'SE'
      3 : outdir = 'regular'
      ELSE : outdir = 'noscreen'
    ENDCASE 
  ENDIF 

  outdr = co2_basedir('fits')+site+outdir+'/'
  create_dir,outdr
  yrstr1 = strtrim(string((fix(yr) MOD 2000),format='(i2.2)'),2)
  fname = co2_const('file_head')+yrstr1+'.filled.xdf'
  print,'Writing file: ',outdr+fname
  write_xdf,outdr+fname,fn,header=fnh

END

PRO create_resppsyn_wcreek,yr,kind=kind,outdir=outdir
;retool to read Bruce's good fills and extract out what we need
;put back in low u* or SE (not low u*)
;u* cutoff is actually 0.175

  site = 'wcreek'
  readfluxtext,yr,site=site,flux=f,met=m,fhead=fh,mhead=mh
  fm = read_filledmet_other(yr,site=site,header=fmh)
  yrstr1 = strtrim(string((fix(yr) MOD 2000),format='(i2.2)'),2)
  par = reform(fm[3,*])
  temp = reform(fm[5,*])
  stor = zapbadval(despike(ensemble_fill(despike(reform(f[6,*]),nsig=6),maxdays=60,mindays=14)))
  bs = where(isnan(stor),nbs)
  IF nbs GT 0 THEN stor[bs] = 0.0
  cflux = despike(reform(f[7,*]),nsig=6)
  nee = stor + cflux
  ustar = reform(f[11,*])
  winddir = reform(m[20,*])
  se = (winddir GE 90.0) AND (winddir LT 180.0)
  flag = reform(f[12,*])
  lat = 45.81
  lon = -90.08
  dn = daynight(yr=yr,lat=lat,lon=lon)
  night = (dn EQ 0)
  day = (dn NE 0)

  IF n_elements(kind) EQ 0 THEN kind = 0
  
;flag = 1 ustar flag = 2 winds

  CASE kind OF
    -1 : badlocs = flag EQ 5
    0 : badlocs = (ustar LT 0.3) OR (se)
    1 : badlocs = (se) OR (night AND ustar LT 0.1)
    2 : badlocs = (se) OR (night AND ustar LT 0.2)
    3 : badlocs = (se) OR (night AND ustar LT 0.3)
    4 : badlocs = (se) OR (night AND ustar LT 0.4)
    5 : badlocs = (se) OR (night AND ustar LT 0.5)
    6 : badlocs = (ustar LT 0.3)
    7 : badlocs = (se)
    8 : badlocs = (se AND night) OR (ustar LT 0.3)
    9 : badlocs = (se AND night) OR ((ustar LT 0.3) AND night)
  ENDCASE 

  CalcRespPsynFits,yr,nee=nee,par=par,temp=temp,respfit=rf,psynfit=pf,fillednee=fn,badval=nan(),badlocs=badlocs,resphead=rh,psynhead=ph,filledhead=fnh,lat=lat,lon=lon,/quiet

  IF n_elements(outdir) EQ 0 THEN outdir = strtrim(string(kind),2)
  outdr = co2_basedir('fits')+site+outdir+'/'
  create_dir,outdr
  fname = co2_const('file_head')+yrstr1+'.filled.xdf'
  print,'Writing file: ',outdr+fname
  write_xdf,outdr+fname,fn,header=fnh

  print,1.0368*total(fn[8,*])/48.0,1.0368*total(fn[13,*])/48.0,1.0368*total(fn[16,*])/48.0

stop
END

PRO run_wcreek_var_test

  FOR i = -1,7 DO BEGIN 
    create_Resppsyn_wcreek,2002,kind=i
  ENDFOR 

END

PRO spec_play,site,yr=yr,level=level
  IF n_elements(site) EQ 0 THEN site = 'wcreek'
  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF n_elements(level) EQ 0 THEN level = '4' ELSE level = strtrim(string(level),2)
  yrstr1 = strmid(string(yr,format='(i4.4)'),2,2)
  yrstr2 = string(yr,format='(i4.4)')

  IF site EQ 'wcreek' OR site EQ 'lcreek' THEN BEGIN 
    dir = mydocs()+'co2/data/'+site+'/specland/'+yrstr1
    days = string(fix(findgen(364)+1),format='(i3.3)')
    fnames = dir+days+'.spec'
    goodq = [0.0]
    badq = [0.0]
    goodc = goodq
    badc = badq
    
    FOR i = 0,n_elements(fnames)-1 DO BEGIN
      thefname = fnames[i]
;    print,'On fname ',thefname
      IF file_test(thefname,/read) THEN BEGIN
        dta = read_ascii(thefname,delimiter=' ')
        dta = dta.(0)
        goodc = [goodc,reform(dta[1,*])]
        badc = [badc,reform(dta[2,*])]
        goodq = [goodq,reform(dta[3,*])]
        badq = [badq,reform(dta[4,*])]
      ENDIF 
    ENDFOR 

    goodq = cdr(goodq) & badq = cdr(badq) & goodc = cdr(goodc) & badc = cdr(badc)
  ENDIF ELSE BEGIN
    IF site EQ 'sylvania' THEN BEGIN
      thefile = co2_basedir('spec')+'unfilled/'+'hl'+yrstr1+'.spec.xdf'
      dta = read_xdf(thefile,header=head)
      goodc = dta[2,*]
      badc = dta[3,*]
      goodq = dta[4,*]
      badq = dta[5,*]
    ENDIF ELSE BEGIN
      IF site EQ 'wlef' THEN BEGIN
        thefile = mydocs()+'co2/data/wlef/lef'+yrstr2+'.spec'+level
        dta = read_ascii(thefile,delimiter=' ')
        dta = dta.(0)
        goodc = dta[6,*]
        badc = dta[7,*]
        goodq = dta[8,*]
        badq = dta[9,*]
      ENDIF ELSE BEGIN
        goodq = [0.0]
        badq = [0.0]
        goodc = [0.0]
        badc = [0.0]
      ENDELSE 
    ENDELSE
  ENDELSE 

  ;stats

  validq = where(badq GE 0.0,numvalsq)
  validc = where(badc GE 0.0,numvalsc)

  gq = where(goodq GT 1.0 AND goodq LT 1.4 AND isnotnan(goodq),ngq)
  IF ngq GT 0 THEN mgq = mean(goodq[gq]) ELSE mgq = -999.0
  pergq = 100.0 * ngq/numvalsq 
  bqlo = where(badq GT 0.0 AND badq LT 1.0,nbqlo)
  IF nbqlo GT 0 THEN mbqlo = mean(badq[bqlo]) ELSE mbqlo = -999.0
  perbqlo = 100.0 * nbqlo/numvalsq
  IF site EQ 'sylvania' THEN BEGIN 
    bqhi1 = where(goodq GT 1.4,nbqhi1)
    bqhi2 = where(badq GT 1.4,nbqhi2)
    IF nbqhi1 GT 0 AND nbqhi2 GT 0 THEN mbqhi = mean([goodq[bqhi1],badq[bqhi2]]) ELSE mbqhi = -999.0
    perbqhi = 100.0 * (nbqhi1 + nbqhi2) / numvalsq
  ENDIF ELSE BEGIN 
    bqhi = where(badq GT 1.4,nbqhi)
    IF nbqhi GT 0 THEN mbqhi = mean(badq[bqhi]) ELSE mbqhi = -999.0
    perbqhi = 100.0 * nbqhi/numvalsq
  ENDELSE 
  combq1b = where(badq GT 1.4,ncombq1b)
  combq1g = where(goodq GT 1.0,ncombq1g)
  IF ncombq1b GT 0 AND ncombq1g GT 0 THEN combq1 = [goodq[combq1g],badq[combq1b]] ELSE IF ncombq1g GT 0 THEN combq1 = goodq[combq1g] ELSE IF ncombq1b GT 0 THEN combq1 = badq[combq1b] ELSE combq1 = [-999.0]
  mcombq1 = mean(combq1)
  percombq1 = 100.0 * (ncombq1b + ncombq1g) / numvalsq
  combq2b = where(badq GT 1.4 AND badq LT 2.0,ncombq2b)
  combq2g = where(goodq GT 1.0 AND goodq LT 2.0,ncombq2g)
  IF ncombq2b GT 0 AND ncombq2g GT 0 THEN combq2 = [goodq[combq2g],badq[combq2b]] ELSE IF ncombq2g GT 0 THEN combq2 = goodq[combq2g] ELSE IF ncombq2b GT 0 THEN combq2 = badq[combq2b] ELSE combq2 = [-999.0]
  mcombq2 = mean(combq2)
  percombq2 = 100.0 * (ncombq2b + ncombq2g) / numvalsq

  print,'Site '+site+' Q factor Year: '+yrstr2
  print,'    Mean Good Q: '+string(mgq,format='(f5.2)')+' Per: '+string(pergq,format='(i3)')+'%'
  print,'  Mean Bad Q Lo: '+string(mbqlo,format='(f5.2)')+' Per: '+string(perbqlo,format='(i3)')+'%'
  print,'  Mean Bad Q Hi: '+string(mbqhi,format='(f5.2)')+' Per: '+string(perbqhi,format='(i3)')+'%'
  print,' Mean Comb GT 1: '+string(mcombq1,format='(f5.2)')+' Per: '+string(percombq1,format='(i3)')+'%'
  print,'  Mean Comb 1-2: '+string(mcombq2,format='(f5.2)')+' Per: '+string(percombq2,format='(i3)')+'%'
  print,' Median:         '+string(median(combq1),format='(f5.2)')


  gc = where(goodc GT 1.0 AND goodc LT 1.3 AND isnotnan(goodc),ngc)
  IF ngc GT 0 THEN mgc = mean(goodc[gc]) ELSE mgc = -999.0
  pergc = 100.0 * ngc/numvalsc 
  bclo = where(badc GT 0.0 AND badc LT 1.0,nbclo)
  IF nbclo GT 0 THEN mbclo = mean(badc[bclo]) ELSE mbclo = -999.0
  perbclo = 100.0 * nbclo/numvalsc
  bchi = where(badc GT 1.3,nbchi)
  IF nbchi GT 0 THEN mbchi = mean(badc[bchi]) ELSE mbchi = -999.0
  perbchi = 100.0 * nbchi/numvalsc
  combc1b = where(badc GT 1.3,ncombc1b)
  combc1g = where(goodc GT 1.0,ncombc1g)
  IF ncombc1b GT 0 AND ncombc1g GT 0 THEN combc1 = [goodc[combc1g],badc[combc1b]] ELSE IF ncombc1g GT 0 THEN combc1 = goodc[combc1g] ELSE IF ncombc1b GT 0 THEN combc1 = badc[combc1b] ELSE combc1 = [-999.0]
  mcombc1 = mean(combc1)
  percombc1 = 100.0 * (ncombc1b + ncombc1g) / numvalsc
  combc2b = where(badc GT 1.3 AND badc LT 2.0,ncombc2b)
  combc2g = where(goodc GT 1.0 AND goodc LT 2.0,ncombc2g)
  IF ncombc2b GT 0 AND ncombc2g GT 0 THEN combc2 = [goodc[combc2g],badc[combc2b]] ELSE IF ncombc2g GT 0 THEN combc2 = goodc[combc2g] ELSE IF ncombc2b GT 0 THEN combc2 = badc[combc2b] ELSE combc2 = [-999.0]
  mcombc2 = mean(combc2)
  percombc2 = 100.0 * (ncombc2b + ncombc2g) / numvalsc

  print,'Site '+site+' C factor Year: '+yrstr2
  print,'    Mean Good C: '+string(mgc,format='(f5.2)')+' Per: '+string(pergc,format='(i3)')+'%'
  print,'  Mean Bad C Lo: '+string(mbclo,format='(f5.2)')+' Per: '+string(perbclo,format='(i3)')+'%'
  print,'  Mean Bad C Hi: '+string(mbchi,format='(f5.2)')+' Per: '+string(perbchi,format='(i3)')+'%'
  print,' Mean Comb GT 1: '+string(mcombc1,format='(f5.2)')+' Per: '+string(percombc1,format='(i3)')+'%'
  print,'  Mean Comb 1-2: '+string(mcombc2,format='(f5.2)')+' Per: '+string(percombc2,format='(i3)')+'%'
  print,' Median:         '+string(median(combc1),format='(f5.2)')

;  !p.multi = [0,2,2]

  plothist,combq1,bin=0.1,xrange=[0.95,5],peak=1.0,yrange=[0,1.05],xticks=8,xtickv=numgen(1.0,5.0,0.5),title='Q factor hist '+site+' '+yrstr1 
  oplot,[1.4,1.4],[-10,10]
  plothist,combc1,bin=0.01,xrange=[0.99,1.3],peak=1.0,yrange=[0,1.05],xticks=6,xtickv=numgen(1.0,1.3,0.05),title='C factor hist '+site+' '+yrstr1
;  !p.multi = 0

  stop

END

PRO energystorage_other,yr

;calculate energy storage for LH, H and G for wcreek

  readfluxtext,yr,flux=flux,met=met,fhead=fh,mhead=mh,/wcreek

;replace RH with Q soon (34-38)
  press = zapbadval(reform(met[50,*]))*10.0
  met[34,*] = ensemble_fill(mixing_ratio(met[25,*],press,met[34,*]/100.0))
  met[35,*] = ensemble_fill(mixing_ratio(met[26,*],press,met[35,*]/100.0))
  met[36,*] = ensemble_fill(mixing_ratio(met[27,*],press,met[36,*]/100.0))
  met[37,*] = ensemble_fill(mixing_ratio(met[28,*],press,met[37,*]/100.0))
  met[38,*] = ensemble_fill(mixing_ratio(met[29,*],press,met[38,*]/100.0))

  qheights = [2,25,40,60,80,97]*0.3048
  qarrlocs = [34,35,36,37,38,33]
  theights1 = [0.25,0.5,0.75,1.0,2.0] 
  theights2 = [25,40,60,80,97]*0.3048
  theights = [theights1,theights2]
  tarrlocs = [21,22,23,24,25,26,27,28,29,30]

  sheight = .075
  sarrloc = 40

  lh = reform(flux[9,*])
  sh = reform(flux[10,*])
  g = reform(met[68,*])
  nelem = n_elements(lh)

  airt = ensemble_fill(reform(met[30,*]))+273.15
  qmix = zapbadval(reform(met[33,*]))
  soilm = zapbadval(reform(met[45,*]))/100.0
  soilt = reform(met[40,numgen(0,nelem,4.0)])
  soilt = interpol(soilt,numgen(0,nelem,4.0),numgen(0,nelem,1.0))

  tstr = make_array(nelem,/float,value=nan())
  qstr = tstr
  gstr = tstr



  FOR i = 0,17518 DO BEGIN
    loc = i
    loc2 = i + 1
    IF isnotnan(lh[loc]) THEN BEGIN 
      deltaq = zapbadval(reform(met[[qarrlocs],loc2]-met[[qarrlocs],loc]))/1800.0
      intq = int_tabulated(qheights,deltaq)
    ENDIF ELSE intq = nan()
    IF isnotnan(sh[loc]) THEN BEGIN 
      deltat = zapbadval(reform(met[[tarrlocs],loc2]-met[[tarrlocs],loc]))/1800.0
      intt = int_tabulated(theights,deltat)
    ENDIF ELSE intt = nan()
    IF isnotnan(g[loc]) THEN BEGIN 
      deltas = (soilt[loc2]-soilt[loc])/1800.0
      ints = deltas * sheight
    ENDIF ELSE ints = nan()
    drho = ddensity(press[loc],qmix[loc],airt[loc])
    qrho = qdensity(press[loc],qmix[loc],airt[loc])
    convert_flux,drho,qrho,airt[loc],qflx=intq,tflx=intt,gflx=ints,lwc=soilm[loc],qflux=qstor,tflux=tstor,gflux=gstor
    tstr[loc] = tstor
    gstr[loc] = gstor
    qstr[loc] = qstor

;    IF isnan(tstor) THEN tstor = 0.0
;    IF isnan(qstor) THEN qstor = 0.0
;    IF isnan(gstor) THEN gstor = 0.0
;    lh[loc] = lh[loc] + qstor
;    sh[loc] = sh[loc] + tstor
;    g[loc] = g[loc] + (gstor * (-1.0))
  ENDFOR
  
  lh = lh + ensemble_fill(qstr)
  sh = sh + ensemble_fill(tstr)
  g = g - ensemble_fill(gstr)
  flux[9,*] = lh
  flux[10,*] = sh
  met[68,*] = g
  energyplot,flux=flux,met=met,site='wcreek'
  stop

END

PRO split_chen_data,site,yr

;Read 2002 chen flux files
;fill into year long file with firs five columns same as Sylvania
;rest of columns fill with data
;make sane header
;update readfluxtext to work with these files (mhw,yhw,mrp,yrp,pd)

  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr LT 2000 THEN yr = yr + 2000
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)
  
  ndays = days_in_year(yr)
  nrows = ndays * 48.0

  sitename = co2_basedir('othersites')+'chen/'+yrstr+'/'+site+'_'+yrstr+'.csv'
  outdir = co2_basedir('othersites')+'chen/'
  fluxfile = outdir+site+yrstr+'_flux.xdf'

  dta = read_ascii(sitename,delimiter=',',data_start=1,header=hdr)
  dta = dta.(0)
  nfcols = n_elements(dta[*,0])+4
  get_fluxmetheader,fhead=newhdr,site=site

  fluxarr = make_array(nfcols,nrows,/float,value=nan())
  fluxarr[0,*] = yr
  fluxarr[4,*] = (indgen(nrows) / 48) + 1
  fluxarr[5,*] = (indgen(nrows) MOD 48) / 48.0
  fluxarr[3,*] = fluxarr[5,*] * 24.0
  FOR mo=1,12 DO BEGIN
    doys = doy_mo(mo,yr=yr)
    ndoys = n_elements(doys)
    fluxarr[1,((doys[0]-1)*48):(((doys[ndoys-1]-1)*48)+47)] = mo
    fluxarr[2,((doys[0]-1)*48):(((doys[ndoys-1]-1)*48)+47)] = (indgen(ndoys*48) / 48) + 1
  ENDFOR

  FOR i = 0,n_elements(dta[0,*])-1 DO BEGIN
    IF i MOD 1000 EQ 0 THEN print,'on line ',i
    loc = round(((dta[0,i]-1) + dta[1,i])*48)
    fluxarr[6:45,loc] = dta[2:41,i]
  ENDFOR 
  bval = where(fluxarr EQ -999.,nbv)
  IF nbv GT 0 THEN fluxarr[bval] = nan()
  fluxarr[6:10,*] = fluxarr[6:10,*] * 22.7272

  write_xdf,fluxfile,float(fluxarr),header=newhdr
END

PRO split_all_chen_Data,yr
  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr EQ 2002 THEN sites = ['mhw','yhw','mrp','yrp','pd'] ELSE sites = ['mhw','ihw','mrp','irp','pb']
  FOR i = 0,n_elements(sites)-1 DO BEGIN
    print,'On site ',sites[i]
    split_chen_data,sites[i],yr
  ENDFOR 
END

;need for 2003 to split met_fix data before fill_met

PRO split_chen_metfix,site
  yr = 2003
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)
  readfluxtext,yr,site=site,flux=f,fhead=fh

  fnew = f

  dta = read_ascii(co2_basedir('othersites')+'chen/'+yrstr+'/'+site+'_03_met_fix.csv',delimiter=',',data_start=1,header=hdr)
  dta = dta.(0)

  locs = reform(round((dta[0,*] + (dta[1,*]))*48.0))
  glocs = where(locs GE 0 AND locs LE 17519 AND isnotnan(locs),ngl)
;stop
  IF ngl GT 0 THEN BEGIN 
; doy time H P  RN PAR Precip winddir windspd
; 0   1    2 3  4   5   6       7      8
;         14 36 37 26  38     39      40 
    fnew[14,locs[glocs]] = dta[2,glocs]
    fnew[40,locs[glocs]] = dta[3,glocs]
    fnew[41,locs[glocs]] = dta[4,glocs]
    fnew[30,locs[glocs]] = dta[5,glocs]
    fnew[42,locs[glocs]] = dta[6,glocs]
    fnew[43,locs[glocs]] = dta[7,glocs]
    fnew[44,locs[glocs]] = dta[8,glocs]
  ENDIF 
  c = where(fnew EQ -999,nc)
  IF nc GT 0 THEN fnew[c] = nan()
  write_xdf,co2_basedir('othersites')+'chen/'+site+yrstr+'_flux.xdf',fnew,header=fh
END

PRO make_fill_met_chen,yr

  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr EQ 2002 THEN sites = ['mhw','yhw','mrp','yrp','pd'] ELSE sites = ['mhw','ihw','mrp','irp','pb']
  IF yr LT 2000 THEN yr = yr + 2000
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)
  temp = make_array(17520,/float,value=nan())
  par = temp

  FOR i = 0,n_elements(Sites)-1 DO BEGIN 
    readfluxtext,yr,site=sites[i],flux=f
    temp = merge_array(temp,reform(f[22,*]))
    par = merge_array(par,reform(f[30,*]))
stop
  ENDFOR 

  k = where(par LT 0.0,nk)
  IF nk GT 0 THEN par[k] = 0.0

  temp[5712:*] = ensemble_fill(temp[5712:*])
  par[5712:*] = ensemble_fill(par[5712:*])

  outp = [transpose(temp),transpose(par)]
  write_xdf,co2_basedir('othersites')+'chen/fill_met_'+yrstr+'.xdf',outp,header=['temp','par']
END

PRO fill_chen_data,site,yr

;read site xdf file
;do ustar cutoff 
;mhw 0.3 yhw 0.15 mrp 0.15 yrp 0.2 pd 0.15
;use flag field in col 45 1 = good
;find lat/lon for each site
;data is in cst as far as i can tell
;par is 30 - not filled - fill - after 5/20 only?
;nee is 9
;ts10 10 cm soil temp is 22
;cutoff - fill only after 5/20 
;ustar is 33
;winddir is 45 - no screening here yet

  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr LT 2000 THEN yr = yr + 2000
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)
  IF n_elements(site) EQ 0 THEN site = 'mhw'

  readfluxtext,yr,site=site,flux=f,fhead=fh
  fillmet = read_xdf(co2_basedir('othersites')+'chen/fill_met_'+yrstr+'.xdf',header=fmh)

  par = merge_array(f[30,*],fillmet[1,*])
  temp = merge_array(f[22,*],fillmet[0,*])
  k = where(par LT 0.0,nk)
  IF nk GT 0 THEN par[k] = 0.0

  nee = f[9,*]
  ustar = f[33,*]
  flag = f[45,*]

  CASE site OF 
    'mhw' : BEGIN
      ustarcut = 0.3
      lat = 46.63546
      lon = -91.10068
      age = 65
    END 
    'yhw' : BEGIN 
      ustarcut = 0.2
      lat = 46.63546
      lon = -91.10058
      age = 14
    END 
    'ihw' : BEGIN
      ustarcut = 0.2
      lat = 46.63546
      lon = -91.10068
      age = 0
    END 
    'mrp' : BEGIN 
      ustarcut = 0.3
      lat = 46.71889
      lon = -91.17928
      age = 63
    END 
    'yrp' : BEGIN 
      ustarcut = 0.2
      lat = 46.71889
      lon = -91.17928
      age = 8
    END 
    'irp' : BEGIN
      ustarcut = 0.2
      lat = 46.69770
      lon = -91.122080
      age = 26
    END 
    'pd' : BEGIN
      ustarcut = 0.15
      lat = 46.71889
      lon = -91.17928
      age = 0
    END 
    'pb' : BEGIN
      ustarcut = 0.2
      lat = 46.63546
      lon = -91.10068
      age = 0
    END 
  ENDCASE 
  
  badvals = where(isnan(nee) OR (ustar LE ustarcut) OR (flag EQ 0) OR (nee GT 20.0) OR (nee LT -40.0),nbv)
  IF nbv GT 0 THEN nee[badvals] = nan()
stop

  IF site EQ 'ihw' OR site EQ 'irp' OR site EQ 'pb' THEN co = [120,300] ELSE co = [120,365]
  IF site EQ 'pb' THEN minpoints = 100
  calcresppsynfits,yr,nee=nee,par=par,temp=temp,latitude=lat,longitude=lon,cutoff=co,fillednee=fn,filledhead=fnh,badval=nan(),minpoints=minpoints

  outdr = co2_basedir('fits')+site+'/'
  create_dir,outdr
  outfl = co2_const('file_head')+yrstr+'.filled.xdf'
  write_xdf,outdr+outfl,fn,header=fnh
END

PRO fill_all_chen_Data,yr
  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr EQ 2002 THEN sites = ['mhw','yhw','mrp','yrp','pd'] ELSE sites = ['mhw','ihw','mrp','irp','pb']
  FOR i = 0,n_elements(sites)-1 DO BEGIN
    print,'On site ',sites[i]
    fill_chen_data,sites[i],yr
  ENDFOR 
END

PRO split_yjp,yr

;Read 2002 chen flux files
;fill into year long file with firs five columns same as Sylvania
;rest of columns fill with data
;make sane header
;update readfluxtext to work with these files (mhw,yhw,mrp,yrp,pd)

  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr LT 2000 THEN yr = yr + 2000
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)
  
  ndays = days_in_year(yr)
  nrows = ndays * 48.0

  sitename = co2_basedir('othersites')+'yjp/yjp_'+yrstr+'.csv'
  outdir = co2_basedir('othersites')+'yjp/'
  fluxfile = outdir+'yjp'+yrstr+'_flux.xdf'

  dta = read_ascii(sitename,delimiter=',',data_start=1,header=hdr)
  dta = dta.(0)
  nfcols = n_elements(dta[*,0])+6-3
  get_fluxmetheader,fhead=newhdr,site='yjp'

  fluxarr = make_array(nfcols,nrows,/float,value=nan())
  fluxarr[0,*] = yr
  fluxarr[4,*] = (indgen(nrows) / 48) + 1
  fluxarr[5,*] = (indgen(nrows) MOD 48) / 48.0
  fluxarr[3,*] = fluxarr[5,*] * 24.0
  FOR mo=1,12 DO BEGIN
    doys = doy_mo(mo,yr=yr)
    ndoys = n_elements(doys)
    fluxarr[1,((doys[0]-1)*48):(((doys[ndoys-1]-1)*48)+47)] = mo
    fluxarr[2,((doys[0]-1)*48):(((doys[ndoys-1]-1)*48)+47)] = (indgen(ndoys*48) / 48) + 1
  ENDFOR

  FOR i = 0,n_elements(dta[0,*])-1 DO BEGIN
    IF i MOD 1000 EQ 0 THEN print,'on line ',i
    loc = round(((dta[1,i]-1)+hhmm_To_fjday(dta[2,i])-(3.0/24.0))*48)
    IF loc GE 0 AND loc LE 17519 THEN fluxarr[6:*,loc] = dta[3:*,i]
  ENDFOR 
  bval = where(fluxarr EQ -999.,nbv)
  IF nbv GT 0 THEN fluxarr[bval] = nan()
  fluxarr[6,*] = fluxarr[6,*] * 22.7272
  fluxarr[9:10,*] = fluxarr[9:10,*] * 22.7272
  write_xdf,fluxfile,float(fluxarr),header=newhdr
END

PRO split_umbs,yr
;convert eastern to central time (sub one)
;double up hours - every data point gets interp on it and 30 minute
;after

  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr LT 2000 THEN yr = yr + 2000
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)
  
  ndays = days_in_year(yr)
  nrows = ndays * 48.0

  sitename = co2_basedir('othersites')+'umbs/umbs_'+yrstr+'.csv'
  outdir = co2_basedir('othersites')+'umbs/'
  fluxfile = outdir+'umb'+yrstr+'_met.xdf'

  dta = read_ascii(sitename,delimiter=',',data_start=1,header=hdr)
  dta = dta.(0)
  nfcols = n_elements(dta[*,0])+6-5
  get_fluxmetheader,mhead=newhdr,site='umb'

  fluxarr = make_array(nfcols,nrows,/float,value=nan())
  fluxarr[0,*] = yr
  fluxarr[4,*] = (indgen(nrows) / 48) + 1
  fluxarr[5,*] = (indgen(nrows) MOD 48) / 48.0
  fluxarr[3,*] = fluxarr[5,*] * 24.0
  FOR mo=1,12 DO BEGIN
    doys = doy_mo(mo,yr=yr)
    ndoys = n_elements(doys)
    fluxarr[1,((doys[0]-1)*48):(((doys[ndoys-1]-1)*48)+47)] = mo
    fluxarr[2,((doys[0]-1)*48):(((doys[ndoys-1]-1)*48)+47)] = (indgen(ndoys*48) / 48) + 1
  ENDFOR

  FOR i = 0,n_elements(dta[0,*])-1 DO BEGIN
    IF i MOD 1000 EQ 0 THEN print,'on line ',i
    loc = round(((dta[1,i]-1)+((dta[2,i]-1.0)/24.0))*48.)
    IF loc GE 0 AND loc LE 17518 THEN BEGIN 
      fluxarr[6:*,loc] = dta[5:*,i]
      fluxarr[6:*,loc+1] = fluxarr[6:*,loc]
    ENDIF 
  ENDFOR 
  bval = where(fluxarr EQ -999.,nbv)
  IF nbv GT 0 THEN fluxarr[bval] = nan()
  write_xdf,fluxfile,float(fluxarr),header=newhdr
END

PRO fill_umbs,yr
;use 2 cm soil temp and PAR - ensemble fill them
;already have gap filled dataset - compare the two
;u* 0.35 cutoff
;test wind - ok

;nee = 22
;soil temp = 13
;PAR = 14
;ustar = 8

;northern hardwood, conifer, age = 90, 20 m canopy

  IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr LT 2000 THEN yr = yr + 2000
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)

  readfluxtext,yr,met=f,mhead=fh,site='umb'
  nee = f[22,*]
  temp = ensemble_fill(f[13,*])
  par = f[14,*]
  k1 = where(par LT -10,nk1)
  IF nk1 GT 0 THEN par[k1] = nan()
  par = ensemble_fill(par)
  k = where(fix(par) LE 0,nk)
  IF nk GT 0 THEN par[k] = 0.0
  par = par * 4.15 ;unit conversion W/m2 to umol/m2/s (approx)
  ustar = f[8,*]
  badvals = where(isnan(nee) OR (nee GT 20.) OR (nee LT -40.0) OR (ustar LT 0.35),nbv)
  IF nbv GT 0 THEN nee[badvals] = nan()

  lat = 45.55984
  lon = -84.71382

  calcresppsynfits,yr,nee=nee,par=par,temp=temp,latitude=lat,longitude=lon,fillednee=fn,filledhead=fnh,badval=nan()
  
  site = 'umbs'
  outdr = co2_basedir('fits')+site+'/'
  create_dir,outdr
  outfl = co2_const('file_head')+yrstr+'.filled.xdf'
  write_xdf,outdr+outfl,fn,header=fnh
stop
;maybe replace with their model in future?

END

PRO fill_yjp,yr
;need to fill PAR with Sylvania - compare
;Soil temp - need to model - maybe use Sylvania - compare
;need u* and wind curofffs

;age 13, canopy is 2.3 m - so young
  
  lat = 46.6465
  lon = -88.5194

 IF n_elements(yr) EQ 0 THEN yr = 2002
  IF yr LT 2000 THEN yr = yr + 2000
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)

  readfluxtext,yr,flux=f,fhead=fh,site='yjp'
  fm = read_fill_met(yr,header=fmh)

  nee = f[10,*]
  ustar = f[7,*]
  par = merge_array(f[51,*]*2.0,fm[10,*])
  k = where(fix(par) LE 0,nk)
  IF nk GT 0 THEN par[k] = 0.0
  temp = fm[6,*]
  temp2 = merge_array(f[55,*],fm[8,*])

  badvals = where(isnan(nee) OR (nee GT 20.) OR (nee LT -40.0) OR (ustar LT 0.3),nbv)
  IF nbv GT 0 THEN nee[badvals] = nan()

  calcresppsynfits,2002,nee=nee,par=par,temp=temp,latitude=lat,longitude=lon,fillednee=fn,filledhead=fnh,badval=nan(),cutoff=[90,290]
  
  site = 'yjp'
  outdr = co2_basedir('fits')+site+'/'
  create_dir,outdr
  outfl = co2_const('file_head')+yrstr+'.filled.xdf'
  write_xdf,outdr+outfl,fn,header=fnh

END

;make WLEF 5 yr avg

PRO make_wlef_avg
;read 5 year XDF and bin_avg by fjday , binsize is 30 minutes 0 to 3652330
;check with Dan for something more robust - future

  dta = read_xdf(co2_basedir('fits')+'wlef/neefilledall.xdf',header=hdr)
  read_resppsyn,filledhead=fnh,2002,fillednee=fn
  x = where(dta eq -999,nx)
  IF nx GT 0 THEN dta[x] = nan()

  fn[4:16,*] = nan()

  FOR i = 0,17519 DO BEGIN
    locs = where(dta[2,*] EQ fn[2,i] AND dta[1,*] EQ fn[1,i],nlocs)
    IF nlocs GT 0 THEN FOR j = 4,16 DO fn[j,i] = mean(dta[j,locs],/nan)
  ENDFOR 

  for j = 4,16 do fn[j,*] = zapbadval(fn[j,*])

  write_xdf,co2_basedir('fits')+'wlef/wlefavg.xdf',fn,header=fnh
END

PRO rotate_chen,site,yr
;replace Cflux, Qflux and u* data in flux file with rot file
;create new arrays and split into them then replace columns

;xxx_rot.csv = DOY, Hour, u*, LEwpl, FCwpl
;remember to convert units of FCwpl to umol/m2/s (* 22.7273)
;replace columns 7 (flux wpl), 11 (LE)
;update column 9 (NEE) (bad storage = 0)
;column 33 is ustar
;refill... - check u* vals again - plot
;keep backups - DONE

;will be different for 2003!

  IF n_elements(yr) EQ 0 THEN yr = 2002
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)
  readfluxtext,yr,site=site,flux=f,fhead=fh

  fnew = f
  fnew[7,*] = nan()
  fnew[9,*] = nan()
  fnew[11,*] = nan()
  fnew[33,*] = nan()
  stor = fnew[8,*]
  bstor = where(isnan(stor),nbs)
  IF nbs GT 0 THEN stor[bstor] = 0.0

  dta = read_ascii(co2_basedir('othersites')+'chen/'+yrstr+'/'+site+'_rot.csv',delimiter=',',data_start=1,header=hdr)
  dta = dta.(0)

  locs = reform(round((dta[0,*] + (dta[1,*]/24.0))*48.0))
  glocs = where(locs GE 0 AND locs LE 17519 AND isnotnan(locs),ngl)
  IF ngl GT 0 THEN BEGIN 
    fnew[7,locs[glocs]] = dta[4,glocs]
    fnew[11,locs[glocs]] = dta[3,glocs]
    fnew[33,locs[glocs]] = dta[2,glocs]
  ENDIF 
  c = where(fnew[7,*] LE (-999) OR fnew[7,*] GE 999,nc)
  IF nc GT 0 THEN fnew[7,c] = nan()
  c = where(fnew[11,*] LE (-999) OR fnew[11,*] GE 999,nc)
  IF nc GT 0 THEN fnew[11,c] = nan()
  c = where(fnew[33,*] LE (-999) OR fnew[33,*] GE 999,nc)
  IF nc GT 0 THEN fnew[33,c] = nan()
  fnew[7,*] = fnew[7,*] * 22.7273
  fnew[9,*] = fnew[7,*] + stor

  sum2 = where(f[1,*] GE 6 AND f[1,*] LE 8 AND daynight() EQ 0 AND fnew[9,*] GE 0)
  plot,l,bin_avg(fnew[33,sum2],fnew[9,sum2],binsize=0.05,locations=l),psym=1,yrange=[0,10],xrange=[0,1.6]

  write_xdf,co2_basedir('othersites')+'chen/'+site+yrstr+'_flux.xdf',fnew,header=fh
;stop
END



PRO compare_chen_AND_all,yr,mos,totals=totals
;add UMBS -orange and YJP - purple
;also 5 years WLEF - brown
;11 sites
;split graphs 3,4,5 / 6,7,8 / 9,10,11 - 3 pages

  IF n_elements(Yr) EQ 0 THEN yr = 2002
  yrstr = strmid(string(yr,format='(i4.4)'),2,2)

  IF yr EQ 2002 THEN BEGIN 
    sites=['sylnew_std','wcreek','mhw','mrp','umbs','yhw','yrp','yjp','pba','lcreek','wlef']
    sy = ['SY','WC','MH','MR','UM','YH','YR','YJ','PB','LC','WL']
  ENDIF ELSE BEGIN
    sites=['sylnew_std','wcreek','mhw','mrp','null','ihw','irp','yjp','pbb','lcreek','wlef']
    sy = ['SY','WC','MH','MR','  ','IH','IR','YJ','PB','LC','WL','  ']
  ENDELSE 

  colors = [fsc_color('black',100),fsc_color('green',101),fsc_color('blue',102),fsc_color('red',103),fsc_color('orange',104),fsc_color('cyan',105),fsc_color('pink',106),fsc_color('yellow',107),fsc_color('purple',108),fsc_color('gray',109),fsc_color('brown',110)]

  allf = make_array(n_elements(sites),17,17520)
  FOR i = 0,n_elements(sites)-1 DO BEGIN
    IF sites[i] EQ 'null' THEN allf[i,*,*] = 0.0 ELSE BEGIN 
      read_resppsyn,yr,altdir=sites[i],fillednee=fn,filledhead=fnh
      allf[i,*,*] = fn
    ENDELSE 
  ENDFOR 

;  mos = [5,6,7,8,9,10]
  IF n_elements(mos) EQ 0 THEN BEGIN
    mos = [6,7,8]
    motext = '6-8'
  ENDIF ELSE motext = strtrim(string(min(mos)),2)+'-'+strtrim(string(max(mos)),2)

  snee = make_array(n_elements(sites),n_elements(mos),/float,value=nan())
  sresp = snee
  sgep = snee
  FOR i = 0,n_elements(mos)-1 DO BEGIN
    indoys = doy_mo(mos[i])
    FOR j = 0,n_elements(sites)-1 DO BEGIN
      locs = where(allf[j,1,*] GE min(indoys) AND allf[j,1,*] LE max(indoys),nlocs)
      IF nlocs GT 0 THEN BEGIN
        snee[j,i] = mean(allf[j,16,locs],/nan)*48.0*1.0368
        sresp[j,i] = mean(allf[j,8,locs],/nan)*48.0*1.0368
        sgep[j,i] = mean(allf[j,8,locs]-allf[j,16,locs],/nan)*48.0*1.0368
      ENDIF 
    ENDFOR 
  ENDFOR

  k = where(isnan(snee),nk)
  IF nk GT 0 THEN snee[k] = 0
  k = where(isnan(sresp),nk)
  IF nk GT 0 THEN sresp[k] = 0
  k = where(isnan(sgep),nk)
  IF nk GT 0 THEN sgep[k] = 0

  totnee = total(snee,2)
  toter = total(sresp,2)
  totgep = total(sgep,2)

  IF keyword_set(totals) THEN BEGIN 
    plot,[0,0],[0,0],/nodata,xrange=[0,n_elements(sites)-1],yrange=[-200,1200],title='20'+yrstr+' -NEE '+motext,xticks=1,xtickn=[' ',' '],yticks=7,ytitle='gC m!U-2!N'
    oplot,[0,20],[0,0]
    bar_plot,totnee*(-1),/over,colors=colors,/outline,barnames=sy
    plot,[0,0],[0,0],/nodata,xrange=[0,n_elements(sites)-1],yrange=[0,1600],title='20'+yrstr+' ER '+motext,xticks=1,xtickn=[' ',' '],yticks=8,ytitle='gC m!U-2!N'
    bar_plot,toter,/over,colors=colors,/outline,barnames=sy
    plot,[0,0],[0,0],/nodata,xrange=[0,n_elements(sites)-1],yrange=[0,1600],title='20'+yrstr+' GEP '+motext,xticks=1,xtickn=[' ',' '],yticks=8,ytitle='gC m!U-2!N'
    bar_plot,totgep,/over,colors=colors,/outline,barnames=sy
    stop

  ENDIF ELSE BEGIN 
    !p.multi = [0,1,3]
    
    plot,[0,0],[0,0],/nodata,xrange=[mos[0],mos[n_elements(mos)-1]+1],yrange=[-400,100],xticklen=0.5,xminor=1,title='20'+yrstr+' NEE '+motext,xticks=n_elements(mos)
    oplot,[0,12],[0,0]
    bwid = 1.0/(n_elements(sites)+1)
    bspc = 1.0 - bwid
    FOR i = 0,n_elements(sites)-1 DO bar_plot,snee[i,*],/over,colors=replicate(colors[i],n_elements(mos)),barwidth=bwid,baroffset=1.0+(bwid*i*(n_elements(sites)+1)),barspace=bspc
    
    plot,[0,0],[0,0],/nodata,xrange=[mos[0],mos[n_elements(mos)-1]+1],yrange=[0,600],xticklen=0.5,xminor=1,title='20'+yrstr+' ER '+motext,xticks=n_elements(mos)
    oplot,[0,12],[0,0]
    bwid = 1.0/(n_elements(sites)+1)
    bspc = 1.0 - bwid
    FOR i = 0,n_elements(sites)-1 DO bar_plot,sresp[i,*],/over,colors=replicate(colors[i],n_elements(mos)),barwidth=bwid,baroffset=1.0+(bwid*i*(n_elements(sites)+1)),barspace=bspc
    
    plot,[0,0],[0,0],/nodata,xrange=[mos[0],mos[n_elements(mos)-1]+1],yrange=[0,600],xticklen=0.5,xminor=1,title='20'+yrstr+' GEP '+motext,xticks=n_elements(mos)
    oplot,[0,12],[0,0]
    bwid = 1.0/(n_elements(sites)+1)
    bspc = 1.0 - bwid
    FOR i = 0,n_elements(sites)-1 DO bar_plot,sgep[i,*],/over,colors=replicate(colors[i],n_elements(mos)),barwidth=bwid,baroffset=1.0+(bwid*i*(n_elements(sites)+1)),barspace=bspc
    
    !p.multi = 0

  ENDELSE 

stop
END