;now we will find the parameters for each file

;use mcmc.pro to estimate parameters for each site
;fix LAI_MAX

;build model here based on mcmc_testmodel3
;use sipnet clim file
;photosynthesis - T,PAR,VPD,maybe precip
;respiration - T,GPPd,maybe precip
;water effects
;  precip based - running mean anomaly? 30 day? against 10 yr?
;  deltaP should have some GPP and RE effects?
;  linear for GPP
;  more complex for RE

;likelihood model should
;test for loglikely of hourly,daily,weekly,monthly,and annual

;then upscale with local met 10-yr, regional met, NHL met

FUNCTION syn_model,x,param,_extra=ex,precom=precom
;a full featured ecosystem model that runs on hourly timestep for one year
;16 parameters

;input data is climfile, unless precom set
;tair = 5
;tsoil = 6
;par = 7
;precip = 8
;vpd = 9
;precom then
;tair,tsoil,par in umol,precipm,vpd,tempda,tempds,gdd

  LUE = param[0]
  k = param[1]
  LAImin = param[2]
  LAImax = param[3]
  Tmin = param[4]
  Topt = param[5]
  VPDmax = param[6]
  VPDmin = param[7]
  alpha = param[8]
  GDDthresh = param[9]
  beta = param[10]
  TEMPthresh = param[11]
  rs = param[12]
  rv = param[13]
  b1 = param[14]
  b2 = param[15]
  b3 = param[16]
  PREmin = param[17]
  PREopt = param[18]

;let's rewrite model to be faster
  IF keyword_set(precom) THEN BEGIN
    nyears = n_elements(x[0,*,0])
    temp_all = x[*,*,0]
    temps_all = x[*,*,1]
    par_all = x[*,*,2]
    vpd_all = x[*,*,3]
    precipm = x[*,*,4] 
    tempda_all = x[*,*,5]
    tempds_all = x[*,*,6]
    gdd_all = x[*,*,7]
  ENDIF ELSE BEGIN 
    nyears = long(n_elements(x[0,*])/8760)
    temp_all = reform(x[5,*],8760,nyears)
    temps_all = reform(x[6,*],8760,nyears)
    par_all = reform(x[7,*],8760,nyears) * 277.778
    vpd_all = reform(x[9,*],8760,nyears)
    precipm = congrid(reform(mysmooth(average_arr(x[8,*],24,/tot),30,/tot,/for),365,nyears),8760,nyears,/center)
    tempd = reform(average_arr(x[5,*],24),365,nyears)
    gdd_all = congrid(total((tempd-10)>0,1,/cum),8760,nyears,/center)
    tempda_all = congrid(tempd,8760,nyears,/center)
    tempds_all = congrid(reform(average_arr(x[6,*],24),365,nyears),8760,nyears,/center)
  ENDELSE 

  output = fltarr(4,8760l*nyears)
  doy = fix(findgen(8760)/24)  
  leafshow_all = 1. / (1. + exp( (-1.0) * alpha * (gdd_all - GDDthresh)))
  fT = ((temp_all-Tmin)/(Topt-Tmin)) > 0 < 1
  fV = ((VPDmax-vpd_all)/(VPDmax-VPDmin)) > 0 < 1
  gpp1 = par_all * fT * fV
  fP = ((precipm-PREmin)/(PREopt-PREmin)) > 0 < 1
  Rh_all = fP * rs * exp ( b1 * (temps_all-15)) ;heterotrophic and maintenance

  FOR yr = 0l,nyears-1 DO BEGIN 
    leafshow = leafshow_all[*,yr]
    tempda = tempda_all[*,yr]
    tempds = tempds_all[*,yr]

    leafshow = ((leafshow - min(leafshow))/(max(leafshow)-min(leafshow))) > 0 < 1
    leafoffd = where((tempds LE tempthresh) AND (leafshow GT 0.99),nld)
    IF nld GT 0 THEN BEGIN 
      leafoff = 1. - 1. / (1. + exp( (-1.0) * beta * (doy-(leafoffd[0]/24)  )))
      leafoff = ((leafoff - min(leafoff))/(max(leafoff)-min(leafoff))) > 0 < 1
   ENDIF ELSE BEGIN
      leafoff = replicate(1.0,8760)
    ENDELSE 
 
    LAI = ((LAImax-LAImin) * leafshow * leafoff)+LAImin
    fpar = 1.0-exp(-k * LAI)
    gpp = lue * fpar * gpp1[*,yr]
    gppd = congrid(max(reform(gpp,24,365),dim=1),/center,8760)
    
    Rw = (rv * exp ( b2 * (tempda-15))) * (GPPd GT 0) ;growth for wood
    Rl = b3*shift(GPPd,-24)                         ;growth resp
    ER = Rl + Rw + Rh_all[*,yr]    
    NEE = ER - GPP

    ost = yr*8760l
    oen = ((yr+1)*8760l)-1
    
    output[0,ost:oen] = nee
    output[1,ost:oen] = er
    output[2,ost:oen] = gpp
    output[3,ost:oen] = lai
  ENDFOR 


;  alpha - leaf expansion coeff
;  LUE - light use efficiency
;  k   - light attenuation
;  LAImax - max LAI
;  Tmin, Topt
;  alpha - leaf expansion coeff
;  GDDtresh - GDD to show leaves
; Rs, Rv, b1, b2, b3, pmin, popt
;vpdmax = 0, vpdmin = 1 (backward ramp)

; maybe a precip factor?  linear 0-1 precip: pmin to popt - running
; mean - 1 month? for Rh and GPP

;  stop
  return,output
END

FUNCTION syn_model_convclim,x

  nyears = long(n_elements(x[0,*])/8760)
  temp_all = reform(x[5,*],8760,nyears)
  temps_all = reform(x[6,*],8760,nyears)
  par_all = reform(x[7,*],8760,nyears) * 277.778
  vpd_all = reform(x[9,*],8760,nyears)
  precipm = congrid(reform(mysmooth(average_arr(x[8,*],24,/tot),30,/tot,/for),365,nyears),8760,nyears,/center)
  tempd = reform(average_arr(x[5,*],24),365,nyears)
  gdd_all = congrid(total((tempd-10)>0,1,/cum),8760,nyears,/center)
  tempda_all = congrid(tempd,8760,nyears,/center)
  tempds_all = congrid(reform(average_arr(x[6,*],24),365,nyears),8760,nyears,/center)

  newx = fltarr(8760,nyears,8)
  newx[*,*,0] = temp_all
  newx[*,*,1] = temps_all
  newx[*,*,2] = par_all
  newx[*,*,3] = vpd_all
  newx[*,*,4] = precipm
  newx[*,*,5] = tempda_all
  newx[*,*,6] = tempds_all
  newx[*,*,7] = gdd_all

  return,newx

END

FUNCTION syn_model_likelihood,x,y,param,model,valid,dat=dat,_EXTRA=ex,modout=modout,honly=honly
;default is to call model with x,param, compare output[0,*] to y[0,*]
;set dat to choose different columns to compare

 
;  IF n_elements(dat) LE 1 THEN dat = [0,0]
;  modout = call_FUNCTION(model,x,param,_EXTRA=ex)
  modout = reform((syn_model(x,param,/precom))[0,*])

  gv = where(valid AND finite(modout),nvalid)
  IF nvalid GT 1 THEN BEGIN 
    sq = total( (modout[gv] - y[gv])^2  )
    sigma = sqrt(sq/nvalid) > 1e-16
 ;   loglike = (nvalid*alog(sigma))  +  (sq/(2*(sigma^2)))
    loglike = nvalid * alog(4.1327314 * sigma)

    IF ~keyword_set(honly) THEN BEGIN 
    
      modout_t = modout
      y_t = y
      bv = where(~finite(y_t) or ~finite(modout_t),nbv)
      IF nbv gt 0 then begin 
        modout_t[bv] = !values.f_nan
        y_t[bv] = !values.f_nan
      ENDIF        
      modout_d = average_arr(modout_t,24,/nan)
      modout_m = average_arr(modout_t,730,/nan)
      modout_y = average_arr(modout_t,8760,/nan)
      y_d = average_arr(y_t,24,/nan)
      y_m = average_arr(y_t,730,/nan)
      y_y = average_arr(y_t,8760,/nan)
      gvd = where(finite(modout_d) AND finite(y_d),nv_d)
      gvm = where(finite(modout_m) AND finite(y_m),nv_m)
      gvy = where(finite(modout_y) AND finite(y_y),nv_y)

      IF nv_d GT 1 THEN BEGIN 
        sq_d = total((modout_d[gvd]-y_d[gvd])^2)*24
        sigma_d = sqrt(sq_d/(24*nv_d)) > 1e-16
 ;       loglike += ((24*nv_d)*alog(sigma_d)) + (sq_d/(2*(sigma_d^2)))
        loglike += 24 * nv_d * alog(4.1327314 * sigma_d)
      ENDIF 
      
      IF nv_m GT 1 THEN BEGIN 
        sq_m = total((modout_m[gvm]-y_m[gvm])^2)*730
        sigma_m = sqrt(sq_m/(730*nv_m)) > 1e-16
 ;       loglike += ((730*nv_m)*alog(sigma_m)) + (sq_m/(2*(sigma_m^2)))
        loglike += 730 * nv_m * alog(4.1327314 * sigma_m)
      ENDIF 
      
      IF nv_y GT 1 THEN BEGIN 
        sq_y = total((modout_y[gvy]-y_y[gvy])^2)*8760
        sigma_y = sqrt(sq_y/(8760*nv_y)) > 1e-16
 ;       loglike += ((8760*nv_y)*alog(sigma_y)) + (sq_y/(2*(sigma_y^2)))
        loglike += 8760 * nv_y * alog(4.1327314 * sigma_y)
      ENDIF 
    ENDIF 

;true loglike += nvalid*alog(sqrt(2*!pi)) ; loglike*=-1
  ENDIF ELSE BEGIN
    loglike = 1e31
  ENDELSE 
;stop
  return,loglike

END

FUNCTION syn_testmodel2,x,param,_extra=ex
;a bit more complex NEE model that includes dormant and growing season difference
;10 parameters - leafon,leafoff,4 dormant, 4 growing season
  leafon = ((fix(param[0])-1)*24) > 1
  leafoff = ((fix(param[1])-1)*24) > leafon+1 < 8759
  nyrs = n_elements(x[0,*,0])
  nee = fltarr(1,8760,nyrs)
;winter
  FOR i = 0,nyrs-1 DO BEGIN 
    nee[0,0:(leafon-1),i] = (param[2] * exp(param[3]*(x[0:(leafon-1),i,1]-15.0))) -(( param[4]*x[0:(leafon-1),i,2]) / (param[5]+x[0:(leafon-1),i,2]))
;summer
    nee[0,leafon:(leafoff-1),i] = (param[6] * exp(param[7]*(x[leafon:(leafoff-1),i,1]-15.0))) -(( param[8]*x[leafon:(leafoff-1),i,2]) / (param[9]+x[leafon:(leafoff-1),i,2]))
;fall  
    nee[0,leafoff:8759,i] = (param[2] * exp(param[3]*(x[leafoff:8759,i,1]-15.0))) -(( param[4]*x[leafoff:8759,i,2]) / (param[5]+x[leafoff:8759,i,2]))
  ENDFOR 

  return,reform(nee,1,8760l*nyrs)
END


PRO syn_estimate,fillclim=fillclim,filldat=filldat,site=site,fillvalid=fillvalid,honly=honly,altmod=altmod,fname=fname
  IF n_elements(fillclim) EQ 0 OR n_elements(filldat) EQ 0 OR n_elements(fillvalid) EQ 0 THEN BEGIN
    outdir_f = mydocs()+'synthesis/filled/'
    restore,outdir_f+'filled.sav'
  ENDIF 

 name = ['yhw', 'ihw', 'mhw', 'wcr', 'umb', 'syl', 'pb1', 'rpc', 'yrp', $
         'irp', 'mrp', 'yjp', 'mjp', 'lcr', 'wfl', 'sfk', 'lef']
 syr =   [ 2002, 2003, 2002, 2000, 1999, 2002, 2002, 2005, 2002, 2003, 2002, 2004, 2004, 2001, 2005, 2005, 1997]
 eyr =   [ 2002, 2003, 2004, 2006, 2003, 2006, 2002, 2005, 2002, 2003, 2005, 2005, 2004, 2006, 2006, 2006, 2005]
 spt = (syr-1997l)*8760l
 ept = ((eyr-1997l)*8760l)+8759l

;LUE,k,LAImin,LAImax,Tmin,Topt,VPDmax,VPDmin,alpha,GDDthresh
;beta,THEMPthresh,rs,rv,b1,b2,b3,PREmin,PREopt


  pname = ['LUE','k','LAImin','LAImax','Tmin','Topt','VPDmax','VPDmin','alpha','GDDthresh','beta','TEMPthresh','rs','rv','b1','b2','b3','PREmin','PREopt']
  baseparam = [0.025,0.5,0.0,5.0,0.0,15.0,5000.0,1000.0,0.05,200.0,0.1,4.0,2.0,2.0,0.025,0.025,0.05,0.0,20.0]
  changeable = [1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
  maxparam = [1.0,0.75,10.0,10.0,10.0,50.0,10000.0,2500.0,0.5,500.0,0.5,5.0,20.0,20.0,0.3,0.3,0.3,100.0,200.0]
  minparam = [0.001,0.001,0.0,0.0,-15.0,5.0,0.0,0.0,0.001,0.0,0.001,-10.0,0.1,0.1,0.0,0.0,0.0,0.0,0.0]

  laimax = [1.19,3.0,3.86,5.3,3.7,4.06,0.2,0.2,0.52,0.8,2.5, 0.93,2.5,4.9,6.3, 6.3, 3.7]
  laimin = [0.0, 0.0,0.0, 0.0,0.0,0.53,0.1,0.0,0.27,0.4,1.25,0.93,2.5,0.0,6.3, 0.0, 0.5]

;let's use an alternative model
  IF keyword_set(altmod) THEN BEGIN
    model = 'syn_testmodel2'
    like = 'mcmc_likelihood'
    ;x0 = temp, x1 = par
    pname = ['L_on','L_off','ra','rb','pc','pd','ra_s','rb_s','pc_s','pd_s']
    baseparam = [150,300,1.0,0.05,0.5,10000.0,3.0,0.1,20.0,1000.0]
    changeable = [1,1,1,1,1,1,1,1,1,1]
    maxparam = [180,310,10.0,0.5,100.0,100000.0,20.0,0.5,100.0,100000.0]
    minparam = [90,181,0.0,0.0,0.0,10.0,0.0,0.0,0.0,10.0]
  ENDIF 

  IF n_elements(site) EQ 0 THEN site = findgen(17)
  fluxes = fltarr(4,8760l*10l,n_elements(site))
  bestparam = make_array(n_elements(baseparam),n_elements(site),/float,value=nan())
  allparam = make_array(n_elements(baseparam),10000,n_elements(site),/float,value=nan())
  allll = make_array(10000,n_elements(site),/float,value=nan())

  FOR i = 0,n_elements(site)-1 DO BEGIN
    print,'On site ',name[site[i]]
    paramval = baseparam
    IF ~keyword_Set(altmod) THEN BEGIN 
      paramval[2] = laimin[site[i]]
      paramval[3] = laimax[site[i]]
    ENDIF 

    param = {name:pname, value:paramval, max:maxparam, min:minparam, changeable:changeable}

;one speed up - subset x for period of record
    x = syn_model_convclim(fillclim[*,spt[site[i]]:ept[site[i]],site[i]])
    y = reform(filldat[1,spt[site[i]]:ept[site[i]],site[i]] / 0.0432)
    v = reform(fillvalid[1,spt[site[i]]:ept[site[i]],site[i]])
    IF keyword_set(altmod) THEN BEGIN 
      mcmc,x,y,param,/fast,model=model,likelihood=like,outputll=outputll,outputvalue=outputvalue,outputy=outputy,validdata=v,honly=honly       
   ENDIF ELSE BEGIN 
     mcmc,x,y,param,/medium,/precom,model='syn_model',likelihood='syn_model_likelihood',outputll=outputll,outputvalue=outputvalue,outputy=outputy,validdata=v,honly=honly
;     mcmc,x,y,param,/fast,/precom,model='syn_model',likelihood='syn_model_likelihood',outputll=outputll,outputvalue=outputvalue,outputy=outputy,validdata=v,honly=honly

    ENDELSE 
    bestparam[*,i] = outputvalue[*,0]
    npar = (n_elements(outputvalue[0,*])-1)<9999
    allparam[*,0:npar,i] = outputvalue[*,0:npar]
    allll[0:npar,i] = outputll[0:npar]
    IF keyword_Set(altmod) THEN BEGIN
      fluxes[0,*,i] = syn_testmodel2(syn_model_convclim(fillclim[*,*,site[i]]),bestparam[*,i])
    ENDIF ELSE BEGIN 
      fluxes[*,*,i] = syn_model(fillclim[*,*,site[i]],bestparam[*,i])    
    ENDELSE 
;    stop
;saves as you go in case of error
    
    IF n_elements(fname) EQ 0 THEN BEGIN 
      fname = 'synmodel_param_'
      IF keyword_set(honly) THEN fname+='honly_' ELSE fname+='iav_'
      fname+=strjoin(string(site,format='(i2.2)'))
      fname+='.sav'
    ENDIF 

    save,filename=mydocs()+'synthesis/filled/'+fname,bestparam,allparam,allll,fluxes
  ENDFOR
 
;tests - fixed GDD info
;regional vs local meteo
;MCMC likelihood


END

PRO syn_upscale_nhl,prob=prob

  names = ['yhw','ihw','mhw','wcr','umb','syl','pb1','rpc','yrp','irp','mrp','yjp','mjp','lcr','wfl','sfk','lef']

;mhw is problem (2), also mrp 10, wfl, sfk 14,15

;read output parameters
  restore,mydocs()+'synthesis/filled/synmodel_param_honly.sav'
  hbestparam = bestparam
  hallparam = allparam
  restore,mydocs()+'synthesis/filled/synmodel_param.sav'
;  bestparam[*,2] = hbestparam[*,2]
  bestparam[*,10] = hbestparam[*,10]
  bestparam[*,14] = hbestparam[*,14]
  bestparam[*,15] = hbestparam[*,15]
  allparam[*,*,10] = hallparam[*,*,10]
  allparam[*,*,14] = hallparam[*,*,14]
  allparam[*,*,15] = hallparam[*,*,15]
  nruns = n_elements(allll[*,0])

;read met file
  nhlclim = read_xdf(mydocs()+'synthesis/lter/nhlclim.xdf')

;apply each param to met file
;scale NEE by % cover for each
;con,hard,mix,shrub,grass,crop,forwet,herbwet

  sites = findgen(17)
  cover = make_array(n_elements(sites),8)
  
  cover[*,0] = [0,0,0,0,0,0,0,0.01,0.09,0.17,0.68,0.02,0.04,0,0,0,0]
  cover[*,1] = [0.22,0.31,0.47,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
  cover[*,2] = [0.06,0.04,0,0.4,0.4,0.1,0,0,0,0,0,0,0,0,0,0,0]
  cover[*,3] = [0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0]
  cover[*,4] = [0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0.5]
  cover[*,5] = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
  cover[*,6] = [0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0]
  cover[*,7] = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,0]
;make sure all add to 1
  FOR i = 0,7 DO cover[*,i]/=total(cover[*,i])
  
  outnee = make_array(8,n_elements(nhlclim[0,*]),/float,value=0)
  outgpp = outnee
  outre = outnee
  
  sitenee = make_array(n_elements(sites),n_elements(nhlclim[0,*]),/float,value=0)
  sitegpp = sitenee
  sitere = sitenee

  IF keyword_set(prob) THEN BEGIN 
    psitenee = make_array(n_elements(sites),n_elements(nhlclim[0,*])/730.,nruns,/float,value=0)
    psitegpp = psitenee
    psitere = psitenee
    poutnee = make_array(8,n_elements(nhlclim[0,*])/730.,nruns,/float,value=0)
    poutgpp = poutnee
    poutre = poutnee
  ENDIF 
  
;for each site, run syn model save NEE
  
  FOR i = 0,n_elements(sites)-1 DO BEGIN
    print,'On site ',names[i]
    output_s = syn_model(nhlclim,bestparam[*,i])
    sitenee[i,*] = output_s[0,*]
    sitere[i,*] = output_s[1,*]
    sitegpp[i,*] = output_s[2,*]
    IF keyword_set(prob) THEN BEGIN 
      FOR j = 0,nruns-1 DO BEGIN
        IF j MOD 1000 EQ 0 THEN print,' Prob run ',j+1
        output_p = syn_model(nhlclim,allparam[*,j,i])
        psitenee[i,*,j] = average_arr(output_p[0,*],730,/tot)/24.*1.0368
        psitere[i,*,j] = average_arr(output_p[1,*],730,/tot)/24.*1.0368
        psitegpp[i,*,j] = average_arr(output_p[2,*],730,/tot)/24.*1.0368
      ENDFOR 
    ENDIF
  ENDFOR 

  IF keyword_set(prob) THEN BEGIN
    psitenee/=10000.
    psitere/=10000.
    psitegpp/=10000.
  ENDIF 
  
 ; stop

  FOR i = 0,7 DO BEGIN
    print,'On cover ',i
    IF i NE 5 THEN BEGIN 
      FOR j = 0,n_elements(sites)-1 DO BEGIN
        IF cover[j,i] NE 0 THEN BEGIN 
          outnee[i,*]+=sitenee[j,*]*cover[j,i]
          outgpp[i,*]+=sitegpp[j,*]*cover[j,i]
          outre[i,*]+=sitere[j,*]*cover[j,i]
          IF keyword_set(prob) THEN BEGIN 
            FOR k = 0,nruns-1 DO BEGIN
              IF k MOD 1000 EQ 0 THEN print,' Prob run ',j+1
              poutnee[i,*,j] += poutnee[j,*,k] * cover[j,i]
              poutgpp[i,*,j] += poutgpp[j,*,k] * cover[j,i]
              poutre[i,*,j] += poutre[j,*,k] * cover[j,i]
            ENDFOR 
          ENDIF
        ENDIF 
      ENDFOR 
    ENDIF ELSE BEGIN 
      q10 = 2.0
      r10 = 4.5
      qy = 0.025
      amax = 25.0
      soilt = nhlclim[6,*]
      wpar = nhlclim[7,*]*1e6/3600
      er_ag= r10 * q10^((soilt-10)/10)
      gpp_ag = (qy * wpar * amax) / ((qy * wpar) + amax)
      outnee[i,*] = er_ag-gpp_ag
      outgpp[i,*] = gpp_ag
      outre[i,*] = er_ag
      IF keyword_set(prob) THEN BEGIN
        FOR k = 0,nruns-1 DO BEGIN
          poutnee[i,*,k] = outnee[i,*]
          poutgpp[i,*,k] = outgpp[i,*]
          poutre[i,*,k] = outre[i,*]
        ENDFOR 
      ENDIF 
    ENDELSE 
  ENDFOR
   
;  stop

;total cover
;water 13.2, dev 4.8
;cover: 7.8,24.5,19.8,0.8,1.6,1.1,23.4,4.2

;con,hard,mix,shrub,grass,crop,forwet,herbwet

;plot, mean NEE for each covertype

;scaled NEE - upland
  totcov = [7.8,24.5,19.8,0.8,1.6,1.1,23.4,4.2]
  totcov/=total(totcov)
  tnee = outnee[0,*]

  tnee[*] = 0.0
  tgpp = tnee
  tre = tnee

  FOR i = 0,7 DO tnee+=outnee[i,*]*totcov[i]
  FOR i = 0,7 DO tgpp+=outgpp[i,*]*totcov[i]
  FOR i = 0,7 DO tre+=outre[i,*]*totcov[i]

  stop

;output monthly and annual



  save,filename=mydocs()+'synthesis/filled/nhlscaled.sav',outnee,outre,outgpp,tnee,tgpp,tre,cover,totcov,sitenee,sitere,sitegpp

  stop

;plots
;mean flux for each type
;mean net flux
  types = ['con','hard','mix','shrub','grass','crop','forwet','herbwet','NET']
  mnee = fltarr(n_elements(types))
  FOR i = 0,n_elements(types)-1 DO IF i LT 8 THEN mnee[i]=mean(outnee[i,*]) ELSE mnee[i]=mean(tnee)
  mnee*=1.0368*365
  



END


PRO syn_plots



  restore,mydocs()+'synthesis/filled/filled.sav'
  restore,mydocs()+'synthesis/filled/synmodel_param_honly.sav'
  hfluxes=fluxes
  hbestparam = bestparam
  hallll = allll
  hallparam = allparam
  restore,mydocs()+'synthesis/filled/synmodel_param.sav'

;  stop

;observed IAV of all sites

;1:1 of all data? - diurnal and interannual

  nee_o = reform(filldat[1,*,*]/0.0432)
  nee_h = reform(hfluxes[0,*,*])
  nee_i = reform(fluxes[0,*,*])

  aneeo = average_cols(transpose(nee_o),8760,/nan)*1.0368*365
  aneeh = average_cols(transpose(nee_h),8760)*1.0368*365
  aneei = average_cols(transpose(nee_i),8760)*1.0368*365

;change sfk and wfk to new data

  mneeo = average_cols(transpose(nee_o),730,/nan)*1.0368*730./24.
  mneeh = average_cols(transpose(nee_h),730)*1.0368*730./24.
  mneei = average_cols(transpose(nee_i),730)*1.0368*730./24.

  aneeo2 = aneeo
  aneeh2 = aneeh
  aneei2 = aneei
  FOR i = 0,16 DO aneeo2[i,*]-=mean(aneeo[i,*],/nan)
  FOR i = 0,16 DO aneeh2[i,*]-=mean(aneeh[i,*],/nan)
  FOR i = 0,16 DO aneei2[i,*]-=mean(aneei[i,*],/nan)

  
;  iav = [2,10,11,3,4,5,13,14,15,16] ;mhw mrp yjp wcr umb syl lcr wfl
;  sfk lef

  iav = [10,11,3,4,5,13,14,15,16] 

  yr = findgen(10)+1997
stop
  sym = symbol(0,/fill)
  colors = [0,fsc_color('red',100),fsc_color('blue',101),fsc_Color('green',102),fsc_color('brown',103),fsc_color('cyan',104),fsc_color('pink',105),fsc_color('orange',106),fsc_Color('purple',107),fsc_Color('magenta',108)]
  plot,[0,0],[0,0],/nodata,xrange=[1997,2007],yrange=[-300,300],xtitle='Year',ytitle='NEE-mean(NEE)',xticks=10
  FOR i = 0,n_elements(iav)-1 DO oplot,yr,aneeo2[iav[i],*],color=colors[i],thick=3,psym=-sym
  FOR i = 0,n_elements(iav)-1 DO oplot,[1998,1998.5],[-50-i*25,-50-i*25],color=colors[i],thick=3
  FOR i = 0,n_elements(iav)-1 DO xyouts,1998.75,-50-i*25,names[iav[i]]
  oplot,[1900,2100],[0,0]

stop
  plot,[-600,600],[-600,600],xrange=[-600,600],yrange=[-600,600],xtitle='Observed',ytitle='Modeled w/ Hourly',title='Regular cost function'
  FOR i = 0,n_elements(iav)-1 DO oplot,aneeo2[iav[i],*],aneeh2[iav[i],*],psym=sym,color=colors[i]
  FOR i = 0,n_elements(iav)-1 DO oplot,[200,250],[-80-i*35,-80-i*35],color=colors[i],thick=3
  FOR i = 0,n_elements(iav)-1 DO xyouts,275,-80-i*35,names[iav[i]]


stop
  plot,[-600,600],[-600,600],xrange=[-600,600],yrange=[-600,600],xtitle='Observed',ytitle='Modeled w/ Interannual',title='IAV cost Function'
  FOR i = 0,n_elements(iav)-1 DO oplot,aneeo2[iav[i],*],aneei2[iav[i],*],psym=sym,color=colors[i]
 FOR i = 0,n_elements(iav)-1 DO oplot,[200,250],[-80-i*35,-80-i*35],color=colors[i],thick=3
  FOR i = 0,n_elements(iav)-1 DO xyouts,275,-80-i*35,names[iav[i]]

stop
;IAV 3,4,5, 13, 14, 15, 16
;per site

;make a best param set for nhl

  ch_h = fltarr(17)
  ch_i = ch_h
  ca_h = ch_h
  ca_i = ch_h
  cm_h = ch_h
  cm_i = ch_h

  FOR i = 0,16 DO ch_h[i] = correl(filldat[1,*,i],hfluxes[0,*,i])^2
  FOR i = 0,16 DO ch_i[i] = correl(filldat[1,*,i],fluxes[0,*,i])^2
  FOR i = 0,16 DO ca_h[i] = correl(aneeo[i,*],aneeh[i,*])^2
  FOR i = 0,16 DO ca_i[i] = correl(aneeo[i,*],aneei[i,*])^2
  FOR i = 0,16 DO cm_h[i] = correl(mneeo[i,*],mneeh[i,*])^2
  FOR i = 0,16 DO cm_i[i] = correl(mneeo[i,*],mneei[i,*])^2

  !p.multi = [0,1,3]
  plot,ch_i,psym=sym,symsize=1.2,xticks=16,xtickn=names,yrange=[0,1],title='Hourly',ytitle='r2'
  oplot,ch_h,psym=sym,color=200
  plot,cm_i,psym=sym,symsize=1.2,xticks=16,xtickn=names,yrange=[0,1],title='Monthly',ytitle='r2'
  oplot,cm_h,psym=sym,color=200
  plot,ca_i,psym=sym,symsize=1.2,xticks=16,xtickn=names,yrange=[0,1],title='Annual',ytitle='r2'
  oplot,ca_h,psym=sym,color=200


  !p.multi = 0


stop

;next test - remove GDD

;UMBS good case study , mean NEE, ER, GPP the same
;most params same, 
;IAV case better able to find:
;optimum temp for psyn
;alpha and GDDthresh big diff (sigmoidal function), LUE diff
;precip doesn't explain much
  tm = fillclim[1,*,0]+fillclim[2,*,0]/365.+fillclim[3,*,0]/8760.
  spt = (syr-1997)*8760l
  ept = (eyr-1997)*8760l+8759l
  jfluxes = fluxes
  jhfluxes = hfluxes
  jhfluxes[0,5l*8760l:6l*8760l,16]=nan()
  jfluxes[0,5l*8760l:6l*8760l,16]=nan()

  jsites = [3,4,5,13,16]
  FOR kk = 0,n_elements(jsites)-1 DO BEGIN 

    jsite = jsites[kk]
    
    !p.multi = [0,1,2]

    cum_o = total(filldat[1,spt[jsite]:ept[jsite],jsite],/cum,/nan)
    cum_h = total(jhfluxes[0,spt[jsite]:ept[jsite],jsite],/cum,/nan)/24.*1.0368
    cum_i = total(jfluxes[0,spt[jsite]:ept[jsite],jsite],/cum,/nan)/24.*1.0368

    rng = [min([cum_o,cum_h,cum_i],/nan),max([cum_o,cum_h,cum_i],/nan)]

    plot,tm[spt[jsite]:ept[jsite]],cum_o,title=names[jsite]+' red=Hourly blue=IAV',yrange=rng
    oplot,tm[spt[jsite]:ept[jsite]],cum_h,color=colors[1]
    oplot,tm[spt[jsite]:ept[jsite]],cum_i,color=colors[2]
    
;  plot,1997+findgen(120)/12.,average_arr(filldat[1,spt[jsite]:ept[jsite],jsite],730,/tot),title='WLEF red=Hourly blue=IAV'
;  oplot,1997+findgen(120)/12.,average_arr(fluxes[0,spt[jsite]:ept[jsite],jsite],730,/tot)/24.*1.0368,color=colors[1]
;  oplot,1997+findgen(120)/12.,average_arr(hfluxes[0,spt[jsite]:ept[jsite],jsite],730,/tot)/24.*1.0368,color=colors[2]
    
    lefann = average_arr(filldat[1,spt[jsite]:ept[jsite],jsite],8760,/nan,/tot)
    ann_h = average_arr(hfluxes[0,spt[jsite]:ept[jsite],jsite],/nan,8760,/tot)/24.*1.0368  
    ann_i = average_arr(fluxes[0,spt[jsite]:ept[jsite],jsite],/nan,8760,/tot)/24.*1.0368
    IF jsite EQ 16 THEN lefann[5] = nan()
    rng2 = [min([lefann,ann_h,ann_i],/nan),max([lefann,ann_h,ann_i],/nan)]

    plot,findgen(eyr[jsite]-syr[jsite]+1)+syr[jsite],lefann,title=names[jsite]+' red=Hourly blue=IAV',psym=sym,xticks=eyr[jsite]-syr[jsite],yrange=rng2,symsize=1.5 
    oplot,findgen(eyr[jsite]-syr[jsite]+1)+syr[jsite],ann_h,color=colors[1],psym=sym
    oplot,findgen(eyr[jsite]-syr[jsite]+1)+syr[jsite],ann_i,color=colors[2],psym=sym
    !p.multi = 0
;    stop

  ENDFOR 

;params
  pname = ['LUE','k','LAImin','LAImax','Tmin','Topt','VPDmax','VPDmin','alpha','GDDthresh','beta','TEMPthresh','rs','rv','b1','b2','b3','PREmin','PREopt']
  !p.multi = [0,2,2]
  FOR jj = 0,18 DO BEGIN
    bparam_h = hbestparam[jj,*]
    bparam_h_max = max(hallparam[jj,*,*],dim=2,/nan)
    bparam_h_min = min(hallparam[jj,*,*],dim=2,/nan)

    bparam_i = bestparam[jj,*]
    bparam_i_max = max(allparam[jj,*,*],dim=2,/nan)
    bparam_i_min = min(allparam[jj,*,*],dim=2,/nan)

    rng3 = [min([bparam_h_min,bparam_i_min],/nan),max([bparam_h_max,bparam_i_max],/nan)]
    plot,[0,0],[0,0],/nodata,xrange=[0,n_elements(names)-1],yrange=rng3,xtickn=names,xticks=n_elements(names)-1,title=pname[jj]
    oplot,bparam_h,color=colors[1],psym=sym
    FOR ii = 0,n_elements(names)-1 DO oplot,[ii,ii],[bparam_h_min[ii],bparam_h_max[ii]],color=colors[1]
    oplot,bparam_i,color=colors[2],psym=sym
    FOR ii = 0,n_elements(names)-1 DO oplot,[ii,ii],[bparam_i_min[ii],bparam_i_max[ii]],color=colors[2]

  ;  stop

  ENDFOR 
  !p.multi = 0

  print,'Ready for WLEF'
  stop
 

  tmy = findgen(8760)/730.
;  plot,[0,0],yrange=[-400,400],xrange=[0,12]
  !p.multi = [0,2,2]
  FOR j = 0l,8l DO BEGIN
    jsite=16
    IF j NE 5 THEN BEGIN
      plot,[0,0],yrange=[-200,200],xrange=[0,12],title='Year '+string(j+1997)
      cum_o = total(filldat[1,(j*8760l):(((j+1l)*8760l)-1l),jsite],/cum,/nan)
      cum_h = total(jhfluxes[0,(j*8760l):(((j+1l)*8760l)-1l),jsite],/cum,/nan)/24.*1.0368
      cum_i = total(jfluxes[0,(j*8760l):(((j+1l)*8760l)-1l),jsite],/cum,/nan)/24.*1.0368      
      oplot,tmy,cum_o,thick=2
      oplot,tmy,cum_i,color=colors[2],thick=2
      oplot,tmy,cum_h,color=colors[1],thick=2
    ENDIF 
  ENDFOR
  !p.multi = 0
 
;WLEF case
;much stronger GDD effect and Topt psyn
;very small diff in GPP/NEE/ER

END