;creates ED input files to run model (also useful for biomebgc)

;HCN Minocqua Dam 475516-02
;Spooner 478027-01




;forestry
;primary and secondary e-folding factors

;plantation fraction

;land use
;transition years

;mortality
;frost mortality rate, winter decay, fraction leaf C retained at leaf
;drop
;height of treefal for disturbance
;allometry: leaf, structural, height
;seed, sapwood, leaf, root, stored-leaf, strucutral decay rates
;VM0
;SLA
;FIA density independent mortality rate
;need soil carbon pools, strucutral C (for testing)



;----

;time to do real work:

;1. read in daily tmin, tmax, precip from minocqua dam
;process one line at a time, copy right number of days into happy
;array
;skip 1905, fill missing with ensemble

PRO ed_readminoc

  basedir = './data/'
  flmax_nm = basedir+'minocqua_daily_tmaxf.csv'
  flmin_nm = basedir+'minocqua_daily_tminf.csv'
  flpre_nm = basedir+'minocqua_daily_precip.csv'
  
  maxt = read_ascii(flmax_nm,delimiter=',')
  mint = read_ascii(flmin_nm,delimiter=',')
  precip = read_ascii(flpre_nm,delimiter=',')
  
  maxt = float(maxt.(0))
  mint = float(mint.(0))
  precip = float(precip.(0))

  yrs = numgen(min(maxt[0,*]),max(maxt[0,*]),1.0)
  nyrs = long(n_elements(yrs))
  outarr = make_array(5,nyrs*365l,/float,value=-999.)
  outhdr = ['year','yearday','Tmax','Tmin','Prcp']
  units = ['','','(deg C)','(deg C)','(cm)']
  outarr[0,*] = expand_arr(yrs,365l)
  outarr[1,*] = replicate_arr(numgen(1.0,365.0),nyrs)

  workarr = make_array(31,12,nyrs,3,/float,value=-999)

  nrows = [n_elements(maxt[0,*]),n_elements(mint[0,*]),n_elements(precip[0,*])]

;clean up array (has missing months, etc...)
  FOR i = 0,max(nrows)-1 DO BEGIN
    IF i LT nrows[0] THEN BEGIN
      theyr = long(maxt[0,i])-yrs[0]
      themo = long(maxt[1,i])-1
      workarr[*,themo,theyr,0] = maxt[3:33,i]
    ENDIF 

    IF i LT nrows[1] THEN BEGIN
      theyr = long(mint[0,i])-yrs[0]
      themo = long(mint[1,i])-1
      workarr[*,themo,theyr,1] = mint[3:33,i]
    ENDIF 
    
    IF i LT nrows[2] THEN BEGIN
      theyr = long(precip[0,i])-yrs[0]
      themo = long(precip[1,i])-1
      workarr[*,themo,theyr,2] = precip[3:33,i]
    ENDIF 

  ENDFOR

;change badval to nan
  bad = where(workarr EQ -999,nbad)
  IF nbad GT 0 THEN workarr[bad] = nan()

;convert units
  workarr[*,*,*,0] = (workarr[*,*,*,0]-32.0)/1.8
  workarr[*,*,*,1] = (workarr[*,*,*,1]-32.0)/1.8
  workarr[*,*,*,2] = (workarr[*,*,*,2]/100.0)*2.54

;fill badvals

;first interpolate by month
  IF nbad GT 0 THEN BEGIN 
    FOR i = 0,nyrs-1 DO BEGIN
      FOR j = 0,11 DO BEGIN
        FOR k = 0,2 DO BEGIN 
          thearr = reform(workarr[*,j,i,k])
          bd = where(isnan(thearr),nbd)
          IF nbd GT 0 AND nbd LE 15 THEN workarr[bd,j,i,k] = (zapbadval(thearr))[bd]
        ENDFOR 
      ENDFOR 
    ENDFOR

;then interpolate by year if still missing
    bad2 = where(isnan(workarr),nbad2)
    IF nbad2 GT 0 THEN BEGIN
      FOR i = 0,11 DO BEGIN 
        FOR j = 0,30 DO BEGIN 
          FOR k = 0,2 DO BEGIN 
            thearr = reform(workarr[j,i,*,k])
            bd = where(isnan(thearr),nbd)
            IF nbd GT 0 THEN workarr[j,i,bd,k] = (zapbadval(thearr))[bd]
        ENDFOR 
      ENDFOR 
      ENDFOR  
    ENDIF
  ENDIF 

;now flatten, and remove extra monthdays

  ndays = [31,28,31,30,31,30,31,31,30,31,30,31]
  ndcum = [0,total(ndays,/cum)]

  locs1 = replicate(0,12)
  locs2 = ndays-1

  FOR i = 0l,nyrs-1l DO BEGIN
    FOR j = 0l,11l DO BEGIN 
      row1 = (i*365l)+ndcum[j]
      row2 = row1+locs2[j]
      outarr[2:4,row1:row2] = transpose(reform(workarr[locs1[j]:locs2[j],j,i,*]))
    ENDFOR 
  ENDFOR 

;clean up weirdness
  mxmnflip = where(outarr[3,*] GT outarr[2,*],nflip)

  IF nflip GT 0 THEN BEGIN 
    temp = outarr[3,mxmnflip]
    outarr[3,mxmnflip] = outarr[2,mxmnflip]
    outarr[2,mxmnflip] = temp
  ENDIF 

;save output

  write_xdf,basedir+'minocqua.xdf',outarr,header=outhdr
  openw,outfl,basedir+'minocqua.mtcin',/get_lun
  printf,outfl,outhdr
  printf,outfl,units
  printf,outfl,outarr,format='(i4.4,1x,i3,1x,f6.2,1x,f6.2,1x,f5.2)'
  free_lun,outfl

END

PRO ed_fixminoc
  dta = read_xdf('data/temp/minocqua_old.xdf',header=hdr)
  dta[2,*]*=0.98
  dta[2,*]-=0.70
  dta[3,*]*=0.95
  dta[3,*]+=3.3
  
  mxmnflip = where(dta[3,*] GT dta[2,*],nflip)

  IF nflip GT 0 THEN BEGIN 
    temp = dta[3,mxmnflip]
    dta[3,mxmnflip] = dta[2,mxmnflip]
    dta[2,mxmnflip] = temp
  ENDIF 

  write_xdf,'data/temp/minocqua.xdf',dta,header=hdr
  openw,outfl,'data/temp/minocqua.mtcin',/get_lun
  printf,outfl,hdr
  units = ['','','(deg C)','(deg C)','(cm)']
  printf,outfl,units
  printf,outfl,dta,format='(i4.4,1x,i3,1x,f6.2,1x,f6.2,1x,f5.2)'
  free_lun,outfl
END

;4b. create harmonized daily data set 1905-2004
PRO ed_createdailymet
  t1 = read_ascii('data/temp/ParkFalls4804.mtc43',data_start=4)
  t2 = read_Ascii('data/temp/minocqua.mtc43',data_start=4)
  
  t1=t1.(0)
  t2=t2.(0)

  t2 = t2[*,where(t2[0,*] LT 1948)]
  t1 = t1[*,where(t1[0,*] GE 1948)]
  t1[5,where(t1[0,*] GE 2003)] *= 100.

  output = [[t2],[t1]]

  write_xdf,'data/dailymet.xdf',output,header=['year','yday','Tmax','Tmin','Tday','prcp','VPD','srad','daylen']

END

;5. read hourly met data for wlef - from a variety of places
;a. wlef1997 - wlef2003
;b. co2/WLEF1995 1996
;c. fill wlef1995 and 1996 (ensembling)

PRO ed_ncdcfix
  ncdc = (read_ascii('data/temp/ncdctemps.csv',delimiter=',')).(0)
;1995-2004 in UTC

  d = where(ncdc EQ -999,nd)
  IF nd GT 0 THEN ncdc[d] = nan()

  ncdc=double(ncdc)

  ncdc[6:7,*]-=32d
  ncdc[6:7,*]/=1.8d

  Es = 6.112d * exp( (17.67d * ncdc[6,*]) / (243.5d + ncdc[6,*]))
  Ea = 6.112d * exp( (17.67d * ncdc[7,*]) / (243.5d + ncdc[7,*]))

  tmsin = ncdc[0,*]+(ncdc[3,*]/365d)+(ncdc[5,*]/(365d*24d))
  vpd = (Es-Ea)*100.

  outarr = make_array(5,87600,/double,value=nan())
  outhdr = ['Year','DOY','hhmm','T','VPD']
  outarr[0,*] = expand_arr(numgen(1995.,2004.),8760)
  outarr[1,*] = replicate_arr(expand_arr(numgen(1.,365.),24.),10)
  outarr[2,*] = replicate_arr(numgen(0.,23.),3650)
  
  locin = round((tmsin-1995.)*8760)-6
  
  goodt = where(isnotnan(ncdc[6,*]) AND locin GE 0,ngt)
  IF ngt GT 0 THEN outarr[3,locin[goodt]] = ncdc[6,goodt]
  goodvpd = where(isnotnan(vpd) AND locin GE 0,ngv)
  IF ngv GT 0 THEN outarr[4,locin[goodvpd]] = vpd[goodvpd]

  write_xdf,'data/temp/ncdc.xdf',outarr,header=outhdr
END

PRO ed_readlefmet
  w95 = (read_ascii('data/temp/WLEF1995_met.txt',data_start=1)).(0)
  bv = where(w95 EQ -999,nbv)
  IF nbv GT 0 THEN w95[bv] = nan()
  temp95 = reform(w95[21,*])
  temp95 = merge_array(temp95,fit_array(reform(w95[20,*]),temp95))
  rh95 = reform(w95[34,*])
  rh95 = merge_array(rh95,fit_array(reform(w95[33,*]),rh95))
  es95 = esat(temp95)*100.
  vpd95 = es95-((rh95/100.)*es95)
  par95 = merge_array(reform(w95[36,*]),reform(w95[37,*]))
  dn = daynight(yr=1995,interval=24.,lat=45.946,lon=-90.272)
  nt = where(dn EQ 0)
  par95[nt] = 0.0
  par95 = zapbadval(par95)

  w96 = (read_ascii('data/temp/WLEF1996_met.txt',data_start=1)).(0)
  bv = where(w96 EQ -999,nbv)
  IF nbv GT 0 THEN w96[bv] = nan()
  temp96 = reform(w96[21,*])
  temp96 = merge_array(temp96,fit_array(reform(w96[20,*]),temp96))
  rh96 = reform(w96[34,*])
  rh96 = merge_array(rh96,fit_array(reform(w96[33,*]),rh96))
  es96 = esat(temp96)*100.
  vpd96 = es96-((rh96/100.)*es96)
  par96 = merge_array(reform(w96[36,*]),reform(w96[37,*]))
  dn = daynight(yr=1996,interval=24.,lat=45.946,lon=-90.272)
  nt = where(dn EQ 0)
  par96[nt] = 0.0
  par96 = zapbadval(par96)

  w97 = (read_ascii('data/temp/wlef1997')).(0)
  bv = where(w97 EQ -999,nbv)
  IF nbv GT 0 THEN w97[bv] = nan()
  temp97 = reform(w97[3,*])
  btmp = where(w97[4,*] GE 7,nbt)
  IF nbt GT 0 THEN temp97[btmp] = nan()
  temp97 = average_arr(temp97,2,/nan)
  vpd97 = reform(w97[5,*])
  bvpd = where(w97[6,*] GE 7,nbv)
  IF nbv GT 0 THEN vpd97[bvpd] = nan()
  vpd97 = average_arr(vpd97,2,/nan)*1000.

  w98 = (read_ascii('data/temp/wlef1998')).(0)
  bv = where(w98 EQ -999,nbv)
  IF nbv GT 0 THEN w98[bv] = nan()
  temp98 = reform(w98[3,*])
  btmp = where(w98[4,*] GE 7,nbt)
  IF nbt GT 0 THEN temp98[btmp] = nan()
  temp98 = average_arr(temp98,2,/nan)
  vpd98 = reform(w98[5,*])
  bvpd = where(w98[6,*] GE 7,nbv)
  IF nbv GT 0 THEN vpd98[bvpd] = nan()
  vpd98 = average_arr(vpd98,2,/nan)*1000.

  w99 = (read_ascii('data/temp/wlef1999')).(0)
  bv = where(w99 EQ -999,nbv)
  IF nbv GT 0 THEN w99[bv] = nan()
  temp99 = reform(w99[3,*])
  btmp = where(w99[4,*] GE 7,nbt)
  IF nbt GT 0 THEN temp99[btmp] = nan()
  temp99 = average_arr(temp99,2,/nan)
  vpd99 = reform(w99[5,*])
  bvpd = where(w99[6,*] GE 7,nbv)
  IF nbv GT 0 THEN vpd99[bvpd] = nan()
  vpd99 = average_arr(vpd99,2,/nan)*1000.

  w00 = (read_ascii('data/temp/wlef2000')).(0)
  bv = where(w00 EQ -999,nbv)
  IF nbv GT 0 THEN w00[bv] = nan()
  temp00 = reform(w00[3,*])
  btmp = where(w00[4,*] GE 7,nbt)
  IF nbt GT 0 THEN temp00[btmp] = nan()
  temp00 = average_arr(temp00,2,/nan)
  vpd00 = reform(w00[5,*])
  bvpd = where(w00[6,*] GE 7,nbv)
  IF nbv GT 0 THEN vpd00[bvpd] = nan()
  vpd00 = average_arr(vpd00,2,/nan)*1000.

  w01 = (read_ascii('data/temp/wlef2001')).(0)
  bv = where(w01 EQ -999,nbv)
  IF nbv GT 0 THEN w01[bv] = nan()
  temp01 = reform(w01[3,*])
  btmp = where(w01[4,*] GE 7,nbt)
  IF nbt GT 0 THEN temp01[btmp] = nan()
  temp01 = average_arr(temp01,2,/nan)
  vpd01 = reform(w01[5,*])
  bvpd = where(w01[6,*] GE 7,nbv)
  IF nbv GT 0 THEN vpd01[bvpd] = nan()
  vpd01 = average_arr(vpd01,2,/nan)*1000.

  w02 = (read_ascii('data/temp/wlef2002')).(0)
  bv = where(w02 EQ -999,nbv)
  IF nbv GT 0 THEN w02[bv] = nan()
  temp02 = reform(w02[3,*])
  btmp = where(w02[4,*] GE 7,nbt)
  IF nbt GT 0 THEN temp02[btmp] = nan()
  temp02 = average_arr(temp02,2,/nan)
  vpd02 = reform(w02[5,*])
  bvpd = where(w02[6,*] GE 7,nbv)
  IF nbv GT 0 THEN vpd02[bvpd] = nan()
  vpd02 = average_arr(vpd02,2,/nan)*1000.

  w03 = (read_ascii('data/temp/wlef2003')).(0)
  bv = where(w03 EQ -999,nbv)
  IF nbv GT 0 THEN w03[bv] = nan()
  temp03 = reform(w03[3,*])
  btmp = where(w03[4,*] GE 7,nbt)
  IF nbt GT 0 THEN temp03[btmp] = nan()
  temp03 = average_arr(temp03,2,/nan)
  vpd03 = reform(w03[5,*])
  bvpd = where(w03[6,*] GE 7,nbv)
  IF nbv GT 0 THEN vpd03[bvpd] = nan()
  vpd03 = average_arr(vpd03,2,/nan)*1000.

  w04 = (read_ascii('data/temp/wlef2004')).(0)
  bv = where(w04 EQ -999,nbv)
  IF nbv GT 0 THEN w04[bv] = nan()
  temp04 = reform(w04[3,*])
  btmp = where(w04[4,*] GE 7,nbt)
  IF nbt GT 0 THEN temp04[btmp] = nan()
  temp04 = average_arr(temp04,2,/nan)
  vpd04 = reform(w04[5,*])
  bvpd = where(w04[6,*] GE 7,nbv)
  IF nbv GT 0 THEN vpd04[bvpd] = nan()
  vpd04 = average_arr(vpd04,2,/nan)*1000.
  par04 = reform(w04[10,*])
  par04 = average_arr(par04,2,/nan)

;have no badvals
  par97 = reform(((read_ascii('data/solar/prefpar1997.txt')).(0))[3,*])
  par98 = reform(((read_ascii('data/solar/prefpar1998.txt')).(0))[3,*])
  par99 = reform(((read_ascii('data/solar/prefpar1999.txt')).(0))[3,*])
  par00 = reform(((read_ascii('data/solar/prefpar2000.txt')).(0))[3,*])
  par01 = reform(((read_ascii('data/solar/prefpar2001.txt')).(0))[3,*])
  par02 = reform(((read_ascii('data/solar/prefpar2002.txt')).(0))[3,*])
  par03 = reform(((read_ascii('data/solar/prefpar2003.txt')).(0))[3,*])

  par96 = par96[0:8759]
  temp96 = temp96[0:8759]
  vpd96 = vpd96[0:8759]
  par00 = par00[0:8759]
  temp00 = temp00[0:8759]
  vpd00 = vpd00[0:8759]
  par04 = par04[0:8759]
  temp04 = temp04[0:8759]
  vpd04 = vpd04[0:8759]

;create hourly arrays of t, par, vpd

  temp = [temp95,temp96,temp97,temp98,temp99,temp00,temp01,temp02,temp03,temp04]
  vpd = [vpd95,vpd96,vpd97,vpd98,vpd99,vpd00,vpd01,vpd02,vpd03,vpd04]
  par = [par95,par96,par97,par98,par99,par00,par01,par02,par03,par04]
  vl0 = vpd LT 0
  m = where(vl0,nm)
  IF nm GT 0 THEN vpd[m] = nan()

  ncdc = read_xdf('data/temp/ncdc.xdf',header=nchd)
  tempnc = ncdc[3,*]
  vpdnc = ncdc[4,*]

  vpd = merge_array(vpd,fit_array(vpdnc,vpd))
  m2 = where(vl0 AND isnan(vpd),nm2)
  IF nm2 GT 0 THEN vpd[m2] = 0.0

  temp = merge_array(temp,fit_array(tempnc,temp))

  temp = zapbadval(ensemble_fill(temp,ptsperday=24.))
  vpd = zapbadval(ensemble_fill(vpd,ptsperday=24.))
  par = zapbadval(ensemble_fill(par,ptsperday=24.))

  outarr = make_array(7,87600,/float,value=nan())
  outhdr = ['Year','DOY','Hr','T','VPD','PAR','CO2']
  outarr[0,*] = expand_arr(numgen(1995.,2004.),8760)
  outarr[1,*] = replicate_arr(expand_arr(numgen(1.,365.),24.),10)
  outarr[2,*] = replicate_arr(numgen(0.,23.),3650)
  outarr[3,*] = temp
  outarr[4,*] = vpd
  outarr[5,*] = par

  write_xdf,'data/hourlymet1995.xdf',outarr,header=outhdr

END

;6. create hourly simulator for 1905-
;for each day in dailymet
;find set of wlef hours that looks like day
;find days with similar tmax-tmin (within 5) and tmean (within 5) vpd
;and srad (within 20)
;pick a random one - save

PRO ed_hourlymet
  daily = read_xdf('data/dailymet.xdf',header=dhdr)
  deltat = (daily[2,*]-daily[3,*])
  temp = daily[4,*]
  par = daily[7,*] * 2.07 * daily[8,*] / 86400.
  vpd = daily[6,*] * daily[8,*] / 86400.
  
  h95 = read_xdf('data/hourlymet1995.xdf',header=hhdr)
  lef_temp = average_arr(h95[3,*],24)
  lef_par = average_arr(h95[5,*],24)
  lef_vpd = average_arr(h95[4,*],24)
  lef_maxtemp = average_arr(h95[3,*],24,/max)
  lef_mintemp = average_arr(h95[3,*],24,/min)
  lef_deltat = lef_maxtemp-lef_mintemp

  outarr = make_array(6,876000l,/float,value=nan())
  outhdr = ['Year','DOY','Hr','T','VPD','PAR']
  outarr[0,*] = expand_arr(numgen(1905.,2004.),8760)
  outarr[1,*] = replicate_arr(expand_arr(numgen(1.,365.),24.),100)
  outarr[2,*] = replicate_arr(numgen(0.,23.),36500l)

  outarr[*,788400l:*] = h95[0:5,*]

  seed = systime(/seconds)
  
  yr = 1905
  FOR i = 0l,32849l DO BEGIN
    gvals = where(abs(lef_deltat-deltat[i]) LE 2.5 AND abs(lef_temp-temp[i]) LE 2.5 AND $
                  abs(lef_par-par[i]) LE 20. AND abs(lef_vpd-vpd[i]) LE 150.,ngvals)
    IF i MOD 365 EQ 0 THEN print,'On year: '+string(yr++)
    IF ngvals GT 0 THEN BEGIN
      IF ngvals EQ 1 THEN theval = long(gvals[0]) ELSE BEGIN 
        IF ngvals GE 10 THEN gvals = gvals[sort(sqrt((lef_deltat[gvals]-deltat[i])^2+(lef_temp[gvals]-temp[i])^2+(lef_par[gvals]-par[i])^2+(lef_vpd[gvals]-vpd[i])^2))]
        ngvals<=10
        theval = long(gvals[fix(randomu(seed)*ngvals) < (ngvals-1) > 0 ])
      ENDELSE 
 ;correct T
;      spancorr = deltat[i]/lef_deltat[theval]
;      tcorr = temp[i]-(lef_temp[theval] * spancorr)
;      parcorr = par[i]/lef_par[theval]
;      vpdcorr = (vpd[i]-lef_vpd[theval])
      spancorr = 1.0 & tcorr = 0.0 & vpdcorr = 0.0 & parcorr = 1.0
      outarr[3,(i*24l):((i*24l)+23l)] = (h95[3,(theval*24l):((theval*24l)+23l)]*spancorr)+tcorr
      outarr[4,(i*24l):((i*24l)+23l)] = (h95[4,(theval*24l):((theval*24l)+23l)]+vpdcorr)>0
      outarr[5,(i*24l):((i*24l)+23l)] = (h95[5,(theval*24l):((theval*24l)+23l)]*parcorr)>0
    ENDIF ELSE BEGIN
;use 10yr diurnal avg for that day, correct for offsets
      baseloc = (i MOD 365l)*24l
      locs = numgen(baseloc,baseloc+86400l,8760l)
      FOR k = 0l,23l DO outarr[3:5,(i*24l)+k] = total(h95[3:5,locs+k],/nan,2)/10.
      newdeltat = max(outarr[3,(i*24l):((i*24l)+23l)])-min(outarr[3,(i*24l):((i*24l)+23l)])
      outarr[3,(i*24l):((i*24l)+23l)] *= deltat[i]/newdeltat
      outarr[3,(i*24l):((i*24l)+23l)] += temp[i]-(mean(outarr[3,(i*24l):((i*24l)+23l)]))
      outarr[4,(i*24l):((i*24l)+23l)] += vpd[i]-(mean(outarr[4,(i*24l):((i*24l)+23l)]))
      outarr[4,(i*24l):((i*24l)+23l)] >= 0
      outarr[5,(i*24l):((i*24l)+23l)] *= par[i]/mean(outarr[5,(i*24l):((i*24l)+23l)])
      outarr[5,(i*24l):((i*24l)+23l)] >= 0
    ENDELSE
  ENDFOR
  write_xdf,'data/hourlymet.xdf',outarr,header=outhdr
END


;if none, use +/-15 day diurnal avg T to get pattern
;get T and PAR pattern
;compute RH from vpd

;need: T, PAR, rh

;7 do stuff for co2 from 1750-2004
;read maunaloa.co2 - monthly from 1950s onward
;read trace_gas.txt - annual to past
;make trace_gas monthlies like maunaloa for pre 195?
;make maunaloa like hourly wlef (use wlef to make monthly diurnal
;pattern)
;create hourly co2 1750-

PRO sinefit,x,a,f,pder
  ;a sin (bx + c) + d
  inx = sin((a[1]*x)+a[2])
  f = (a[0]*inx)+a[3]  
  pder = [[inx],[a[0]*x*cos((a[1]*x)+a[2])],[a[0]*cos((a[1]*x)+a[2])],[replicate(1,n_elements(x))]]
END


PRO ed_makemonthco2
  ann = (read_Ascii('data/co2/trace_gas.txt',data_start=1)).(0)
  mon = (read_ascii('data/co2/maunaloa.co2c')).(0)
  j = where(mon LT 250.,nj)
  IF nj GT 0 THEN mon[j] = nan()
  monsub = zapbadval(mon[1:12,*])
  monann = mon[14,*]
  monyr = mon[0,*]

  monrng = monsub - ts_Smooth(monsub[*],24)

  rngf = findgen(45)
  FOR i = 0,44 DO rngf[i] = stddev(monrng[*,i])
  trnd = linfit(monann,rngf)

  monrng2 = monrng
  FOR i = 0,44 DO monrng2[*,i]/=rngf[i]

  mons = numgen(1.,12.)+0.5
  monrng3 = total(monrng2,2)/45.

  gdan = where(ann[0,*] GE 1650 AND ann[0,*] LT 1959,ngdan)
  ayr = ann[0,gdan]
  aco2 = ann[1,gdan]
  mon = numgen(ayr[0],ayr[n_elements(ayr)-1]+(11./12.),(1.0/12.0))
  mco2 = interpol(aco2,ayr+0.5,mon)

  ampl = replicate_arr(monrng3,n_elements(ayr))
  muly = expand_arr(trnd[0]+trnd[1]*aco2,12)
  newamp = ampl*muly
  
  mco2+=newamp

  mco2 = [mco2,monsub[*]]

;now for hourly
;use wlef

  yrs = ['1995','1996','1997','1998','1999','2000','2001','2002','2003','2004']
  fname = 'data/co2/WLEF'+yrs+'_met.txt'
  hrco2 = make_array(8760l,10,/float,value=nan())
  FOR i = 0,n_elements(fname)-1 DO BEGIN
    print,'Reading '+fname[i]
    mdta = ((read_Ascii(fname[i])).(0))[6:8,0:8759]
    k = where(mdta LT 250. OR mdta GT 450.,nk)
    IF nk GT 0 THEN mdta[k] = nan()
    hrco2[*,i] = merge_array(mdta[1,*],mdta[0,*],mdta[2,*])
  ENDFOR

  hrco2 = despike(hrco2)
  lefmo = zapbadval(average_arr(hrco2[*],730,/nan))
  hrco3 = hrco2-expand_arr(lefmo,730)
  ens = total(hrco3,/nan,2)/10.
  fillco2 = expand_arr(lefmo,730)+replicate_arr(ens,10)

  hrcof = hrco2
  bhc = where(isnan(hrcof),nbhc)
  IF nbhc GT 0 THEN hrcof[bhc] = fillco2[bhc]

;create composites for daily cycle (% of mean monthly value)
  lefmo2 = average_arr(hrcof[*],730,/nan)
  endev = (hrcof-expand_arr(lefmo,730))/expand_arr(lefmo,730)

  endevens = total(endev,/nan,2)/10.
  hrens = make_array(24,12,/float,value=nan())
  mos = long([31,28,31,30,31,30,31,31,30,31,30,31])
  cummo = [0l,total(mos,/cum)]*24.

  FOR i = 0,11 DO hrens[*,i] = smooth(total(reform(endevens[cummo[i]:(cummo[i+1]-1)],24,mos[i]),2)/mos[i],3,/edge)

  hrensyr = reform(congrid(hrens,24l,365l,/center,/interp),8760l)

  mco3 = mco2[2:*]
  lval = mco3[n_elements(mco3)-1]
  mco3 = [mco3,lval+0.88,lval+0.88+0.92]
  msmoo = ts_Smooth(mco3[*],24)
  mco3 -= msmoo

;meed to redo fit

;  mco4 = (mco3*1.60283)-0.0443622
  mco4 = (mco3*1.55063)-0.0857359
  mco4 += msmoo + 4.0

  finhr = congrid(mco4,3101040l,/center)
  addhr = finhr * replicate_arr(hrensyr,n_elements(mco4)/12l)
  finhr += addhr
  addhr = 0b

  plot,finhr[3022200l:*],hrcof[0:78839],psym=1
  oplot,[0,1000],[0,1000],color=180,thick=4

  finhr[3022200l:*] = hrcof[0:78839]
  finhr = [finhr,hrcof[78840:*]]

  write_xdf,'data/hourlyco2.xdf',transpose(finhr),header=['CO2']
;include future? not now

END

;8 create output met for 1750-
;randomly select periods of 10-25 yrs from 1905-1975 and replicate
;to save sapce, only include pointers to 1905 met

PRO ed_createoldmet
;randomly select chunks of 10 to 25 years

;two random sets:
;1. size of chunks - max needed 26 - range 10-25
;2. start date of year
;from 1905-1975
;will be some overlap
  bdir = '~/apps/ed/createinputs/'

  tm1 = round((randomu(systime(/seconds),26)*15)+10)
  tm2 = round((randomu(systime(/seconds)/200.*tm1[tm1[0]],26)*70)+1905)
  nel = where(total(tm1,/cum) LE 280,nne1)  
  yrs=[0]
  FOR i = 0,nne1-1 DO yrs=[yrs,numgen(tm2[i],tm2[i]+tm1[i]-1)]
  yrs = cdr(yrs)
  yrs = transpose(yrs[0:254])
  plot,yrs
  write_xdf,bdir+'data/oldmet.xdf',yrs,header=['Years']
END

;9. run mech files on it! insanity!

PRO ed_createmech,noread=noread,calcmech=calcmech
;read look up table (600 MB ram)
  bdir = '~/apps/ed/createinputs/'
  print,'Reading lut'
  IF NOT keyword_set(noread) THEN arrays=farq_readarrays()
;read yrs array for recreating the past
;  yrsarr = long(read_xdf(bdir+'data/oldmet.xdf'))
;read reconstructed met array (23 mb)
  print,'Reading hourly met'
  hrmet = read_xdf(bdir+'data/hourlymet.xdf',header=mhd)
  temp = reform(hrmet[3,*])
  vpd = reform(hrmet[4,*])
  par = reform(hrmet[5,*])
;  hrmet = 0b
  es = esat(temp)*100.0
  rh = ((1.0-(vpd/es))*100.0)>10.0<100.0
;  es = 0b
;  vpd = 0b

;read reconstructed CO2 (11 mb)
  print,'Reading co2'
  ppm = reform(read_xdf(bdir+'data/hourlyco2.xdf'))

;remember to convert vpd to rh
;  yrs = numgen(1650l,2004l)
  yrs = numgen(1800l,1905l)

  FOR i = 0l,n_elements(yrs)-1l DO BEGIN
    theyear = yrs[i]
    print,'Making mech file for year ',theyear
    IF theyear LT 1905l THEN metyear = (theyear+100)>1905<2004 ELSE metyear = theyear ;yrarrs[i]
    mloc1 = (metyear-1905l)*365l*24l
    mloc2 = ((metyear-1905l+1l)*365l*24l)-1l
    cloc1 = (theyear-1650l)*365l*24l
    cloc2 = ((theyear-1650l+1l)*365l*24l)-1l 
;    tm = systime(/sec)
;    farq_runyear,par[mloc1:mloc2],ppm[cloc1:cloc2],temp[mloc1:mloc2],rh[mloc1:mloc2],arrays=arrays,year=theyear,/nowrite,vmrange=[25.]
;    print,'LUT version: ',tm-systime(/sec)
    tm = systime(/sec)
    farq_runyear,par[mloc1:mloc2],ppm[cloc1:cloc2],temp[mloc1:mloc2],rh[mloc1:mloc2],arrays=arrays,year=theyear,calcmech=calcmech
    print,'Time required: ',tm-systime(/sec)    
  ENDFOR 

  ptr_free,arrays
END

;next steps:
;create monthly outputs
;monthly mean temp, DGD, CHD, min_temp, precip, soil_temp, 

PRO ed_fixdaily
;par: 0.792291 78.8151
;temp: 2.56440 1.00759
;tmax: 0.656 0.98262
;tmin: 1.349 0.96713
;vpd: 0.871512 72.775

  daily = read_xdf('data/dailymetold.xdf',header=dhdr)
  old = daily

  daily[2,*] = (daily[2,*]-0.656)/0.983
  daily[3,*] = (daily[3,*]-1.349)/0.967
  daily[4,*] = (daily[4,*]-2.564)/1.008
  jj = where(daily[3,*] GE daily[4,*],njj)
  IF njj GT 0 THEN daily[3,jj] = (2.0*daily[4,jj])-daily[2,jj]
  
  b = where(daily[0,*] EQ 2002)
  daily[5,b]*=100.

  par = daily[7,*] * daily[8,*] * 2.07 / 86400.
  par = ((par-78.815)/0.792) > 10.0
  daily[7,*] = (par / 2.07) * 86400. / daily[8,*]
  
  vpd = daily[6,*] * daily[8,*] / 86400.
  vpd = (vpd-72.775)/0.872
  daily[6,*] = vpd * 86400. / daily[8,*]
  
  stop
  write_xdf,'data/dailymet.xdf',daily,header=dhdr

;rerun hourlymet
;rerun mech :(

END

PRO ed_stemp
  yrs = ['1998','1999','2000','2001','2002','2003','2004']
  fname = 'data/temp/wc'+yrs+'_filled'
  stemp = make_array(8760l,7,/float,value=nan())
  FOR i = 0,n_elements(fname)-1 DO BEGIN
    print,'Reading '+fname[i]
    mdta = ((read_Ascii(fname[i])).(0))[5,*]
    stemp[*,i] = (average_arr(reform(mdta),2,/nan))[0:8759]
  ENDFOR 
  stop
  write_xdf,'data/temp/soilt.xdf',stemp
END

PRO ed_writefl,dr,vname,dta,yrstart=yrstart,nyrs=nyrs,xdr=xdr,format=format
  IF n_elements(yrstart) EQ 0 THEN yrstart = 1650
  IF n_elements(nyrs) EQ 0 THEN nyrs = 355
  IF n_elements(xdr) EQ 0 THEN xdr = 'outputs/xdf/'
  IF n_elements(format) EQ 0 THEN format='f-0.3'
  yrs = numgen(yrstart,yrstart+nyrs-1)
  fnames = dr+vname+'.'+string(yrs,format='(i4.4)')+'.dat'
  favg = dr+vname+'.avg.dat'
  xname = xdr+vname+'.xdf'
  odta = reform(dta,12,nyrs)
  write_xdf,xname,odta,header=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
  lathead = [45.5,-90.5]
  latfor = 'f-0.1,1x,f-0.1,1x'
  FOR i = 0,nyrs-1 DO BEGIN
    openw,fl,fnames[i],/get_lun
    printf,fl,[lathead,reform(odta[*,i])],format='('+latfor+','+'11('+format+',1x),'+format+')'
    free_lun,fl
  ENDFOR 
  avgdta = total(odta,2)/nyrs
  openw,fl,favg,/get_lun
  printf,fl,[lathead,reform(avgdta)],format='('+latfor+','+'11('+format+',1x),'+format+')'
  free_lun,fl
END


PRO ed_monthlies
;format is typically - one file per year

;file name:
;val.year.dat

  fhead = 'outputs/'
  xhead = fhead+'xdf/'

;lat lon 12 values space delimited
;lat is xx.x lon is -xx.x (left just)

;first save xdf file, then output text files

;1650-2004

  daily = read_xdf('data/dailymet.xdf',header=dhdr)
  hourly = read_xdf('data/hourlymet.xdf',header=hhdr)
  oldyrs = read_xdf('data/oldmet.xdf')
  
  pickyrs = long([reform(oldyrs),numgen(1905.,2004.)])
  pyrs = pickyrs-1905l

  modays = long([31,28,31,30,31,30,31,31,30,31,30,31])
  mocum = [0l,total(modays,/cum)]
  hrcum = mocum*24l

  nayrs = 100l
  nyrs = n_elements(pickyrs)

  allmoday1 = [0l,total(replicate_arr(modays,nayrs),/cum)]
  allmoday = lonarr(nyrs*12l)
  FOR i = 0l,11l DO allmoday[i:*:12l] = allmoday1[((pickyrs-1905l)*12l)+i]
  allmohr = allmoday*24l

  allmos = (expand_arr(pyrs,12)*12l)+replicate_arr(numgen(0l,11l),nyrs)

;mean temp  temp/temp.xxxx.dat
  dtemp = reform(daily[4,*],365,nayrs)
  mtemp = fltarr(12,nyrs)
  FOR i = 0,11 DO mtemp[i,*] = total(dtemp[mocum[i]:(mocum[i+1]-1),pyrs],1)/modays[i]
  ed_writefl,fhead+'temp/','temp',mtemp,format='f-6.2'

;min temp min_temp/min_temp.xxxx.dat
;(threshold for mortality < 5)  
;is it mean min temp or min min temp or min mean temp?
  dmtemp = reform(daily[3,*],365,nayrs)
  mintemp = fltarr(12,nyrs)
  FOR i = 0,11 DO mintemp[i,*] = total(dmtemp[mocum[i]:(mocum[i+1]-1),pyrs],1)/modays[i]
  ed_writefl,fhead+'min_temp/','min_temp',mintemp,format='f-6.2'

;ED:
;45.5 -90.5 -18.5 -16.9 -9.2 -1.4 4.7 10.025 13 11.575 6.65 1 -5.875 -14.425

;DGD temp/temp.dgd.xxxx.dat
;use hourly values
  thresh = 5.0
  dgd = reform((average_arr((hourly[3,*]-thresh)>0,730)*30.4167)[allmos],12,nyrs)
  ed_writefl,fhead+'temp/','temp.dgd',dgd,format='f-6.2'

;might be neat to see if you we can relate DGD to start of psyn?
;from ED database:
;45.5 -90.5      0.052   0.560   13.516  80.144  238.858 384.649
;487.451 474.411 322.601 153.737 25.456 1.516
;41.5 -91.5 0 0 0 109.446 338.9881 451.707 552.73 562.5291 416.601 132.4289 0 0

;CHD  temp/temp.chd.xxxx.dat
;CHD (# days with min < 0 C) - from daily met = sum(tmin lt 0)
  chd = fltarr(12,nyrs)
  cthresh = 0.0  ;should be 5 according to ED!
  FOR i = 0,11 DO chd[i,*] = total(dmtemp[mocum[i]:(mocum[i+1]-1),pyrs] LT cthresh,1)
  ed_writefl,fhead+'temp/','temp.chd',chd,format='i-2'

;45.5 -90.5 31 28 31 30 0 0 0 0 0 31 30 31
;45.5 -90.5	30.982	27.964	27.000	13.327	2.145	0.000	0.000	0.000	0.291	7.164	23.509	30.455

;precip precip/precip.xxxx.dat
;monthly total precip (compare a few products here) in mm
;use wisc n.central precip - better regional avg
  mpre = (read_ascii('data/precip/monthlyprecip.csv',data_start=1,delimiter=',')).(0)
  mpreearl = mpre[1:12,0:9]
  precip = reform((reform(mpre[1:12,10:*],1200))[allmos],12,nyrs)
  precip[*,245:254] = mpreearl ;copy in 1895-1904 observed data
  ed_writefl,fhead+'precip/','precip',precip,format='f-6.2'

;soiltemp soil_temp/soil_temp.xxxx.dat
;soil temp (use correction system - compare to actual LEF?)
  st1 = read_xdf('data/temp/soilt1.xdf')
  st2 = read_xdf('data/temp/soilt.xdf')
  st = [st1[*],st2[*]]
  stempo = ts_Smooth(zapbadval(average_arr(st[*],730,/nan)),3)
  mts = ts_smooth(mtemp[*],12)
  stg = mtemp[*]-mts
  stemp = (((((ts_smooth(stg,5)*0.82) > (-8.0))+mts)*1.1)+1.5) > (-1.0)
  stemp[4140:*] = stempo
  stemp = reform(stemp,12,nyrs)
  ed_writefl,fhead+'soil_temp/','soil_temp',stemp,format='f-6.2'

;phenology - play with a value, maybe tune with bloom data
;later - create explicit phenology database?

;NOTES
;monthly dgd base 0C - also have this from historical record - on avg
;can compute from daily
;median spring thaw: 5/16 5/02-5/31
;median fall freeze: 9/23 9/09-10/07
;monthly cdd (# days < 0 C) - also WI stat climatlo - on avg
;can compute from daily
;n days, max < 0 C
; 27.4, 19.5, 10.1, 1.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 11.6, 24.8
;n days, min < 0 C
;31.0, 27.7, 27.8, 18.6, 4.1, 0.1, 0.0, 0.0, 2.2, 13.0, 26.4, 30.8
;monthly total precip - have from WI
;monthly soil temp (some f of temp?) in K
; could do this: Jan-Mar +5.5, Apr-May +1.5, Jun-Aug -1.5, Sep +0.25, Oct-Nov
; +3.0, Dec +5.5

;tweak phenology to leaf out around day 130 = may 10
;phenology notes:
;Botta GCB 2000 - , equation is a + b*e(c*CHD)
;a = -68, b = 638, c =-0.01
;maybe move a to 20 (play around)
;get leafout dates
;= GDD required for leaf out, base 5 in botta, start jan 1
;CHD is base 0 or 5 depending on biome mean T, start nov 1
;leaf off is a function of photoperiod and critical soil temp or
;superciritical soil temp


;in BGC: onset crit: exp 4.795 + 0.129 mean annual T
;prephenology.c and pheneology.c
;leaf off is where soilt (11-day dampded running avg t avg) is less
;than critical soil t = mean fall t (244-304) and daylength less than
;39300 s.


END


;make a separate wetland run?

;mortality
;regular rates
;from FIA get volume rates m3 / m3
;filter on growth rate?
;also SORTIE rates for density-dep mort

;also: imposed mortality?

;vm0

;sla - have check units

;allometry - landis, std papers

;phenology (maybe use remote sensing)

;soil respiration pool decay rates
;= flux/pool (eg. curtis)


PRO ed_makelinks
;spin up run
;symlink 1400-1649 as f of 1650-1680 (30 yr climate)
;directories need to do that:
;data/ed/run1/inputs

;min_temp/min_temp.yyyy.dat
;soil_temp/soil_temp.yyyy.dat
;temp/temp.yyyy.dat
;temp/temp.chd.yyyy.dat
;temp/temp.dgd.yyyy.dat
;mech/c3/lat45.5lon-90.5.yxxxx.c3.mech
;mech/c4/lat45.5lon-90.5.yxxxx.c4.mech
;update /mech/c3/mech.years
;precip/precip.yyyy.dat

  yrin = string(numgen(1650,1680,1),format='(i4.4)')
  yrout = string(numgen(1400,1649,1),format='(i4.4)')

  bdir = '~/data/ed/run1/inputs/'
  prefs = ['min_temp/min_temp.','soil_temp/soil_temp.','temp/temp.','temp/temp.chd.','temp/temp.dgd.','mech/c3/lat45.5lon-90.5.y','mech/c4/lat45.5lon-90.5.y','precip/precip.']
  sufs = ['.dat','.dat','.dat','.dat','.dat','.c3.mech','.c4.mech','.dat']

  prefs = ['mech/c4/xdf/lat45.5lon-90.5.y']
  sufs = ['.c4.mech.xdf']

  FOR i = 0,n_elements(prefs)-1 DO BEGIN
    flout = bdir+prefs[i]+yrout+sufs[i]
    flin = bdir+(replicate_arr(prefs[i]+yrin+sufs[i],20))[0:(n_elements(flout)-1)]
    FOR j = 0,n_elements(flout)-1 DO BEGIN 
      file_delete,flout[j],/allow_nonexistent
      file_link,flin[j],flout[j],/verbose
    ENDFOR 
  ENDFOR 
END





;do a run!