PRO syn_iavplot
;20081120 is the original one
;20090906 is the second one

  dr1 = mydocs()+'synthesis/filled/'
  IF file_test(dr1+'plotdata.sav',/read) THEN BEGIN
    restore,dr1+'plotdata.sav'
  ENDIF ELSE BEGIN 
;    syn_readiav,dt='20081120',/readonly,obs=f,model=a
    syn_readiav,dt='20100125',/readonly,obs=f,model=a
;    syn_readiav,dt='20090906',/readonly,model=s
    syn_readiav,dt='20100127',/readonly,model=s,/spatial
    save,f,a,s,filename=dr1+'plotdata.sav'
  ENDELSE 

  stop

;filling frac from NACP
  fillfrac = [0.506,0.386,0.480,0.321,0.276]

;error bars for an_nee! - look at Barr NACP files?
  r_err = make_array(5,10,/float,value=nan())
  u_err = r_err
  anee_err = r_err
  u_err[0,3:8] = [1.73,11.14,19.54,12.43,nan(),3.52]
  u_err[1,2:9] = [10.97,8.33,12.95,6.99,11.62,11.63,6.27,3.36]
  u_err[2,[5,9]] = [3.24,6.87]
  u_err[3,4:9] = [8.37,3.60,5.05,1.51,1.55,2.13]
  u_err[4,0:4] = [1.58,4.04,4.96,2.91,6.04]
  FOR i = 0,4 DO u_err[i,where(~finite(u_err[i,*]))] = mean(u_err[i,*],/nan)
  r_err[0,3:8] = [12.91,49.90,25.57,nan(),nan(),21.1877]
  r_err[1,2:9] = [23.08,17.57,18.76,18.63,20.01,19.17,21.02,18.94]
  r_err[2,[5,9]] = [18.17,16.16]
  r_err[3,4:9] = [11.18,7.40,8.63,10.35,8.04,8.478]
  r_err[4,0:4] = [17.12,27.30,25.28,16.10,18.93] 
  FOR i = 0,4 DO r_err[i,where(~finite(r_err[i,*]))] = mean(r_err[i,*],/nan)
  anee_err = sqrt(r_err^2 + u_err^2)

;CUP calculations

  cup = make_array(5,10,6,/float,value=nan())
  cupdoy = findgen(365)+1
  FOR i = 0,4 DO BEGIN
    FOR y = 0,9 DO BEGIN 
      cupnee_mod = smooth(average_arr(a.iln[i,(y*730):((y+1)*730)-1],2,/nan),8,/nan,/edge)
      cupnee_obs = smooth(average_arr(f.nee[i,(y*730):((y+1)*730)-1],2,/nan),8,/nan,/edge)
      cupnee_mloc = where(cupnee_MOD LT 0 AND cupdoy GT 100 AND cupdoy LT 300,nvalm)
      cupnee_oloc = where(cupnee_obs LT 0 AND cupdoy GT 100 AND cupdoy LT 300,nvalo)
      IF nvalm GT 2 THEN BEGIN 
        cup[i,y,0] = cupnee_mloc[0]+1
        cup[i,y,2] = cupnee_mloc[nvalm-1]+1
        cup[i,y,4] = cup[i,y,2]-cup[i,y,0]
      ENDIF 
      IF nvalo GT 2 THEN BEGIN 
        cup[i,y,1] = cupnee_oloc[0]+1
        cup[i,y,3] = cupnee_oloc[nvalo-1]+1
        cup[i,y,5] = cup[i,y,3]-cup[i,y,1]
      ENDIF 
    ENDFOR 
  ENDFOR 

;plot one per site
;show each year, cup_start and cup_end, and model means, sort by
;latitude (umbs, wcr, pfa, los, syv)
;also show modis evi? NO
;but temperature

;WC, UM, SY, LO, PF
;r2 = 0.71, 0.74, 0.14, 0.56, 0.02  with MAT and CUP_l
;r2 = 0.87, 0.76, 0.64, 0.50, 0.01 with MAT and CUP_Start
;poor for cup_end
;Mean annual PAR not a good predictor
;MAP when Temp > 0, r2 = 0.45,0.57,0.16,0.03,0.01,   MAP and CUP_len
;parpp and tempp correl at 0.67 (-0.82)
;role of remote sensing discussion of fpar, EVI, and temp

;UMBS 145   274
;WC 138    287
;PFA 121  280
;LOS 137  267
;SYL 130 224

;trend in fall
;spatial var: fall more variables (25) than spring (9) in mean CUP date
;spring temporal mean var is 14, fall temporal mean var is 13, but larger at SYL 

  loadct,0
  snames = ['US-WCr','US-UMB','US-Syv','US-Los','US-PFa'] 
  lstyle = [0,2,0,2,1]
  cols = [0,0,180,180,100]
  thcks = [2,2,3,3,4]
  doy = findgen(365)+1
  yrs = findgen(10)+1997

;MAKE FIGURE

;Fig 8
  
;filled square = model, open triangle = obs
;gray line = model mean

  !p.multi = [0,2,5,0,1]
  plot,[0,0],[0,0],/nodata,xrange=[1996.5,2006.5],xticks=9,xtickv=yrs,yrange=[90,340],ytitle='Day',xtitle='Year'
  oplot,yrs,cup[1,*,0],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[1,*,2],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[1,*,1],psym=symbol('TRIANGLE',size=3)
  oplot,yrs,cup[1,*,3],psym=symbol('TRIANGLE',size=3)
  xyouts,1997,310,snames[1]

  plot,[0,0],[0,0],/nodata,xrange=[1996.5,2006.5],xticks=9,xtickv=yrs,yrange=[90,340],ytitle='Day',xtitle='Year'
  oplot,yrs,cup[0,*,0],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[0,*,2],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[0,*,1],psym=symbol('TRIANGLE',size=3)
  oplot,yrs,cup[0,*,3],psym=symbol('TRIANGLE',size=3)
  xyouts,1997,310,snames[0]

  plot,[0,0],[0,0],/nodata,xrange=[1996.5,2006.5],xticks=9,xtickv=yrs,yrange=[90,340],ytitle='Day',xtitle='Year'
  oplot,yrs,cup[4,*,0],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[4,*,2],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[4,*,1],psym=symbol('TRIANGLE',size=3)
  oplot,yrs,cup[4,*,3],psym=symbol('TRIANGLE',size=3)
  xyouts,1997,310,snames[4]

  plot,[0,0],[0,0],/nodata,xrange=[1996.5,2006.5],xticks=9,xtickv=yrs,yrange=[90,340],ytitle='Day',xtitle='Year'
  oplot,yrs,cup[3,*,0],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[3,*,2],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[3,*,1],psym=symbol('TRIANGLE',size=3)
  oplot,yrs,cup[3,*,3],psym=symbol('TRIANGLE',size=3)
  xyouts,1997,310,snames[3]

  plot,[0,0],[0,0],/nodata,xrange=[1996.5,2006.5],xticks=9,xtickv=yrs,yrange=[90,340],ytitle='Day',xtitle='Year'
  oplot,yrs,cup[2,*,0],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[2,*,2],psym=-symbol('SQUARE',/fill,size=2),color=150
  oplot,yrs,cup[2,*,1],psym=symbol('TRIANGLE',size=3)
  oplot,yrs,cup[2,*,3],psym=symbol('TRIANGLE',size=3)
  xyouts,1997,310,snames[2]


stop
;Fig. 1 - map done

;Fig. New
;Mean cum nee for each site and each model

;f.nee[5,7300]
;a.hln[5,7300]
;a.iln[5,7300]
;s.iln[5,7300]

  mean_nee_obs = total(reform(f.nee,5,730,10),3,/nan)/total(reform(finite(f.nee),5,730,10),3)
  aahln = a.hln
  aailn = a.iln
  ssiln = s.iln
  bfnee = where(~finite(f.nee))
  aahln[bfnee] = nan()
  aailn[bfnee] = nan()
  ssiln[bfnee] = nan()
  mean_nee_ah = total(reform(aahln,5,730,10),3,/nan)/total(reform(finite(f.nee),5,730,10),3)
  mean_nee_ai = total(reform(aailn,5,730,10),3,/nan)/total(reform(finite(f.nee),5,730,10),3)
  mean_nee_s = total(reform(ssiln,5,730,10),3,/nan)/total(reform(finite(f.nee),5,730,10),3)
  doy = findgen(365)+1

  !p.multi = [0,2,5,0,1]
  plot,doy,total(average_arr(mean_nee_obs[0,*],2,/nan),/cum,/nan),yrange=[-300,100],xticks=12,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D',' ']
  oplot,doy,total(average_arr(mean_nee_ah[0,*],2,/nan),/cum,/nan),linestyle=1,thick=2
  oplot,doy,total(average_arr(mean_nee_ai[0,*],2,/nan),/cum,/nan),linestyle=2,thick=2
  oplot,doy,total(average_arr(mean_nee_s[0,*],2,/nan),/cum,/nan),color=200,thick=3

  oplot,[10,30],[-100,-100]
  oplot,[10,30],[-150,-150],linestyle=1,thick=2
  oplot,[10,30],[-200,-200],linestyle=2,thick=2
  oplot,[10,30],[-250,-250],color=200,thick=3

  plot,doy,total(average_arr(mean_nee_obs[1,*],2,/nan),/cum,/nan),yrange=[-250,100],xticks=12,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D',' ']
  oplot,doy,total(average_arr(mean_nee_ah[1,*],2,/nan),/cum,/nan),linestyle=1,thick=2
  oplot,doy,total(average_arr(mean_nee_ai[1,*],2,/nan),/cum,/nan),linestyle=2,thick=2
  oplot,doy,total(average_arr(mean_nee_s[1,*],2,/nan),/cum,/nan),color=200,thick=3

  plot,doy,total(average_arr(mean_nee_obs[2,*],2,/nan),/cum,/nan),yrange=[-100,50],xticks=12,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D',' ']
  oplot,doy,total(average_arr(mean_nee_ah[2,*],2,/nan),/cum,/nan),linestyle=1,thick=2
  oplot,doy,total(average_arr(mean_nee_ai[2,*],2,/nan),/cum,/nan),linestyle=2,thick=2
  oplot,doy,total(average_arr(mean_nee_s[2,*],2,/nan),/cum,/nan),color=200,thick=3

  plot,doy,total(average_arr(mean_nee_obs[3,*],2,/nan),/cum,/nan),yrange=[-100,50],xticks=12,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D',' ']
  oplot,doy,total(average_arr(mean_nee_ah[3,*],2,/nan),/cum,/nan),linestyle=1,thick=2
  oplot,doy,total(average_arr(mean_nee_ai[3,*],2,/nan),/cum,/nan),linestyle=2,thick=2
  oplot,doy,total(average_arr(mean_nee_s[3,*],2,/nan),/cum,/nan),color=200,thick=3

  plot,doy,total(average_arr(mean_nee_obs[4,*],2,/nan),/cum,/nan),yrange=[-50,100],xticks=12,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D',' ']
  oplot,doy,total(average_arr(mean_nee_ah[4,*],2,/nan),/cum,/nan),linestyle=1,thick=2
  oplot,doy,total(average_arr(mean_nee_ai[4,*],2,/nan),/cum,/nan),linestyle=2,thick=2
  oplot,doy,total(average_arr(mean_nee_s[4,*],2,/nan),/cum,/nan),color=200,thick=3

  !p.multi = 0

stop


;Fig. 2  IAV all sites - in std dev + err
  anee_s = fltarr(5)
  FOR i = 0,4 DO anee_s[i] = stddev(f.an_nee[i,*],/nan)
  an_nee2 = f.an_nee
  FOR i = 0,4 DO an_nee2[i,*] /= anee_s[i]
  an_nee_err = anee_err
  FOR i = 0,4 DO an_nee_err[i,*] /= anee_s[i]

  !p.multi = [0,2,2]
  plot,[0,0],[0,0],/nodata,yrange=[-3,3],xrange=[1996.5,2006.5],xtitle='Year',ytitle='NEE anomaly (sigma)'
  oplot,[0,10000],[0,0]
  FOR i = 0,4 DO oplot,findgen(10)+1997,an_nee2[i,*],psym=0,thick=thcks[i],linestyle=lstyle[i],color=cols[i]
  FOR i = 0,4 DO FOR y = 0,9 DO oplot,[y+1997,y+1997],[an_nee2[i,y]-an_nee_err[i,y],an_nee2[i,y]+an_nee_err[i,y]]
  FOR i = 0,4 DO oplot,[1997,1998],[2.75-(i*0.4),2.75-(i*0.4)],thick=thcks[i],linestyle=lstyle[i],color=cols[i]
  FOR i = 0,4 DO xyouts,1998.1,2.75-(i*0.4),snames[i]

stop

;Fig. 3 - a model IL correlation - honly and IAV

  !p.multi = [0,2,2]
;Fig 3a 
  plot,[0,0],[0,0],/nodata,yrange=[-200,300],xrange=[-200,300],title=correl(f.an_nee,a.an_hla)^2
  oplot,[-1000,1000],[-1000,1000]
  FOR i = 0,4 DO FOR y = 0,9 DO oplot,[f.an_nee[i,y]-anee_err[i,y],f.an_nee[i,y]+anee_err[i,y]],[a.an_hla[i,y],a.an_hla[i,y]],color=0
  FOR i = 0,4 DO oplot,f.an_nee[i,*],a.an_hla[i,*],psym=symbol(i,/fill,size=(i+2)/1.5),color=50.*i
  FOR i = 0,4 DO FOR y = 0,9 DO oplot,[f.an_nee[i,y],f.an_nee[i,y]],[a.an_hla_min[i,y],a.an_hla_max[i,y]],color=0
  FOR i = 0,4 DO plots,250,105-(i*60),psym=symbol(i,/fill,size=(i+2)/1.5),color=50.*i
  FOR i = 0,4 DO xyouts,260,105-(i*60),snames[i]

;Fig 3b
  plot,[0,0],[0,0],/nodata,yrange=[-200,300],xrange=[-200,300],title=correl(f.an_nee,a.an_ila)^2
  oplot,[-1000,1000],[-1000,1000]
  FOR i = 0,4 DO FOR y = 0,9 DO oplot,[f.an_nee[i,y]-anee_err[i,y],f.an_nee[i,y]+anee_err[i,y]],[a.an_ila[i,y],a.an_ila[i,y]],color=0
  FOR i = 0,4 DO oplot,f.an_nee[i,*],a.an_ila[i,*],psym=symbol(i,/fill,size=(i+2)/1.5),color=50.*i
  FOR i = 0,4 DO FOR y = 0,9 DO oplot,[f.an_nee[i,y],f.an_nee[i,y]],[a.an_ila_min[i,y],a.an_ila_max[i,y]],color=0
  FOR i = 0,4 DO plots,250,105-(i*60),psym=symbol(i,/fill,size=(i+2)/1.5),color=50.*i
  FOR i = 0,4 DO xyouts,260,105-(i*60),snames[i]



  stop

;FIG 4. - WC lon loff
wc_loff = [nan(),nan(),268,270,269,273,277,278,269,268]
wc_lon = [nan(),nan(),147,149,147,115,159,163,158,152]
wc_loff2 = [nan(),nan(),293,295,291,266,299,301,301,288]
wc_lon2 = [nan(),nan(),115,120,119,107,131,131,130,116]
wc_loff3 = (wc_loff+wc_loff2)/2.0
wc_lon3 = (wc_lon+wc_lon2)/2.0
wc_loff -= mean(wc_loff,/nan)
wc_lon -= mean(wc_lon,/nan)
wc_loff2 -= mean(wc_loff2,/nan)
wc_lon2 -= mean(wc_lon2,/nan)
wc_loff3 -= mean(wc_loff3,/nan)
wc_lon3 -= mean(wc_lon3,/nan)

s_loff = s.leafoff[0,*]-mean(s.leafoff[0,2:*])
s_lon = s.leafon[0,*]-mean(s.leafon[0,2:*])
;s_lon = (total(s.leafon,1)/5.0)-mean(s.leafon[0,2:*])
;s_loff = (total(s.leafoff,1)/5.0)-mean(s.leafoff[0,2:*])
a_loff = a.leafoff[0,*]-mean(a.leafoff[0,2:*])
a_lon = a.leafon[0,*]-mean(a.leafon[0,2:*])

!p.multi = [0,2,2]
plot,wc_loff2,a_loff,yrange=[-20,20],xrange=[-20,20],psym=5,xtitle='Observed (days)',ytitle='Modeled (days)'
oplot,wc_lon2,a_lon,psym=6
oplot,[-20,20],[-20,20]
plots,-18,18,psym=6
xyouts,-17,18,'Leaf On (r2='+string(correl(wc_lon2,a_lon)^2,format='(f0.2)')+')'
plots,-18,14,psym=5
xyouts,-17,14,'Leaf Off (r2='+string(correl(wc_loff2,a_loff)^2,format='(f0.2)')+')'
!p.multi = 0

stop

;Fig. 5 - s model IL correlation
  !p.multi = [0,2,2]
  plot,[0,0],[0,0],/nodata,yrange=[-200,300],xrange=[-200,300],title=correl(f.an_nee,s.an_ila)^2
  oplot,[-1000,1000],[-1000,1000]
  FOR i = 0,4 DO FOR y = 0,9 DO oplot,[f.an_nee[i,y]-anee_err[i,y],f.an_nee[i,y]+anee_err[i,y]],[s.an_ila[i,y],s.an_ila[i,y]],color=0
  FOR i = 0,4 DO oplot,f.an_nee[i,*],s.an_ila[i,*],psym=symbol(i,/fill,size=(i+2)/1.5),color=50.*i
  FOR i = 0,4 DO FOR y = 0,9 DO oplot,[f.an_nee[i,y],f.an_nee[i,y]],[s.an_ila_min[i,y],s.an_ila_max[i,y]],color=0
  FOR i = 0,4 DO plots,250,105-(i*60),psym=symbol(i,/fill,size=(i+2)/1.5),color=50.*i
  FOR i = 0,4 DO xyouts,260,105-(i*60),snames[i]

  stop

;Fig. 5 NEE sensitivity of asyn
;  corrleaf = fltarr(3,5)
;  FOR i = 0,4 DO corrleaf[*,i] = [correl(s.an_leafon[i,*],s.an_ila[i,*]),correl(s.an_leafoff[i,*],s.an_ila[i,*]),correl(s.an_gsl[i,*],s.an_ila[i,*])]
;  slope_lon = make_array(1,5,/float,value=nan())
;  slope_loff = slope_lon
;  slope_gsl = slope_lon
;  sigma_lon = slope_lon
;  sigma_loff = slope_lon
;  sigma_gsl = slope_lon
;  FOR i = 0,4 DO BEGIN
;    ft = myfit(s.an_leafon[i,*],s.an_ila[i,*],pval=p,sigma=sigma_l,measure_err=s.ila_s[i,*]*2)
;    IF p LT 1.0 THEN slope_lon[0,i] = ft[1]
;    sigma_lon[0,i] = sigma_l[1]
;;    ft = myfit(s.an_leafon[i,*],s.an_ilg[i,*],pval=p)
;;    IF p LT 1.0 THEN slope_lon[1,i] = ft[1]
;;    ft = myfit(s.an_leafon[i,*],s.an_ilar[i,*],pval=p)
;;    IF p LT 1.0 THEN slope_lon[2,i] = ft[1]
;    ft = myfit(s.an_leafoff[i,*],s.an_ila[i,*],pval=p,sigma=sigma_l,measure_err=s.ila_s[i,*]*2)
;    IF p LT 1.0 THEN slope_loff[0,i] = ft[1]
;   sigma_loff[0,i] = sigma_l[1]
;;    ft = myfit(s.an_leafoff[i,*],s.an_ilg[i,*],pval=p)
;;    IF p LT 1.0 THEN slope_loff[1,i] = ft[1]
;;    ft = myfit(s.an_leafoff[i,*],s.an_ilar[i,*],pval=p)
;;    IF p LT 1.0 THEN slope_loff[2,i] = ft[1]
;    ft = myfit(s.an_gsl[i,*],s.an_ila[i,*],pval=p,sigma=sigma_l,measure_err=s.ila_s[i,*]*2)
;    IF p LT 1.0 THEN slope_gsl[0,i] = ft[1]
;    sigma_gsl[0,i] = sigma_l[1]
;;    ft = myfit(s.an_gsl[i,*],s.an_ilg[i,*],pval=p)
;;    IF p LT 1.0 THEN slope_gsl[1,i] = ft[1]
;;    ft = myfit(s.an_gsl[i,*],s.an_ilar[i,*],pval=p)
;;    IF p LT 1.0 THEN slope_gsl[2,i] = ft[1]
;  ENDFOR 
;  sigma_lon*=2
;  sigma_loff*=2
;  sigma_gsl*=2

  !p.multi = [0,3,3]
;;   plot,[0,0],[0,0],/nodata,yrange=[-15,10],xrange=[-0.05,0.95],ytitle='Slope (gC m-2 day-2)',xtitle='r2'
;;   FOR i = 0,4 DO plots,corrleaf[0,i]^2,slope_lon[0,i],psym=symbol(i,/fill,size=(i+2)/2.0),color=50.*i
;;   FOR i = 0,4 DO plots,[corrleaf[0,i]^2,corrleaf[0,i]^2],[slope_lon[0,i]+sigma_lon[0,i],slope_lon[0,i]-sigma_lon[0,i]]
;;   oplot,[0.3025,0.3025],[-100,100],linestyle=2
;;   oplot,[-1,1],[0,0]
;;   plot,[0,0],[0,0],/nodata,yrange=[-15,10],xrange=[-0.05,0.95],ytitle='Slope (gC m-2 day-2)',xtitle='r2'
;;   FOR i = 0,4 DO plots,corrleaf[1,i]^2,slope_loff[0,i],psym=symbol(i,/fill,size=(i+2)/2.0),color=50.*i
;;   FOR i = 0,4 DO plots,[corrleaf[1,i]^2,corrleaf[1,i]^2],[slope_loff[0,i]+sigma_loff[0,i],slope_loff[0,i]-sigma_loff[0,i]]
;;   oplot,[0.3025,0.3025],[-100,100],linestyle=2
;;   oplot,[-1,1],[0,0]
;;   plot,[0,0],[0,0],/nodata,yrange=[-15,15],xrange=[-0.05,0.95],ytitle='Slope (gC m-2 day-2)',xtitle='r2'
;;   FOR i = 0,4 DO plots,corrleaf[2,i]^2,slope_gsl[0,i],psym=symbol(i,/fill,size=(i+2)/2.0),color=50.*i
;;   FOR i = 0,4 DO plots,[corrleaf[2,i]^2,corrleaf[2,i]^2],[slope_gsl[0,i]+sigma_gsl[0,i],slope_gsl[0,i]-sigma_gsl[0,i]]
;;   oplot,[0.3025,0.3025],[-100,100],linestyle=2
;;   oplot,[-1,1],[0,0]

  corrleaf = fltarr(3,5)
  FOR i = 0,4 DO corrleaf[*,i] = [correl(a.an_leafon[i,*],a.an_ila[i,*]),correl(a.an_leafoff[i,*],a.an_ila[i,*]),correl(a.an_gsl[i,*],a.an_ila[i,*])]
  slope_lon = make_array(1,5,/float,value=nan())
  slope_loff = slope_lon
  slope_gsl = slope_lon
  sigma_lon = slope_lon
  sigma_loff = slope_lon
  sigma_gsl = slope_lon
  FOR i = 0,4 DO BEGIN
    ft = myfit(a.an_leafon[i,*],a.an_ila[i,*],pval=p,sigma=sigma_l,measure_err=a.ila_s[i,*]*2)
    IF p LT 1.0 THEN slope_lon[0,i] = ft[1]
    sigma_lon[0,i] = sigma_l[1]
;    ft = myfit(a.an_leafon[i,*],a.an_ilg[i,*],pval=p)
;    IF p LT 1.0 THEN slope_lon[1,i] = ft[1]
;    ft = myfit(a.an_leafon[i,*],a.an_ilar[i,*],pval=p)
;    IF p LT 1.0 THEN slope_lon[2,i] = ft[1]
    ft = myfit(a.an_leafoff[i,*],a.an_ila[i,*],pval=p,sigma=sigma_l,measure_err=a.ila_s[i,*]*2)
    IF p LT 1.0 THEN slope_loff[0,i] = ft[1]
   sigma_loff[0,i] = sigma_l[1]
;    ft = myfit(a.an_leafoff[i,*],a.an_ilg[i,*],pval=p)
;    IF p LT 1.0 THEN slope_loff[1,i] = ft[1]
;    ft = myfit(a.an_leafoff[i,*],a.an_ilar[i,*],pval=p)
;    IF p LT 1.0 THEN slope_loff[2,i] = ft[1]
    ft = myfit(a.an_gsl[i,*],a.an_ila[i,*],pval=p,sigma=sigma_l,measure_err=a.ila_s[i,*]*2)
    IF p LT 1.0 THEN slope_gsl[0,i] = ft[1]
    sigma_gsl[0,i] = sigma_l[1]
;    ft = myfit(a.an_gsl[i,*],a.an_ilg[i,*],pval=p)
;    IF p LT 1.0 THEN slope_gsl[1,i] = ft[1]
;    ft = myfit(a.an_gsl[i,*],a.an_ilar[i,*],pval=p)
;    IF p LT 1.0 THEN slope_gsl[2,i] = ft[1]
  ENDFOR 
  sigma_lon*=2
  sigma_loff*=2
  sigma_gsl*=2

  plot,[0,0],[0,0],/nodata,yrange=[-15,10],xrange=[-0.05,0.95],ytitle='Slope (gC m-2 day-2)',xtitle='r2'
  FOR i = 0,4 DO plots,corrleaf[0,i]^2,slope_lon[0,i],psym=symbol(i,/fill,size=(i+2)/2.0),color=50.*i
  FOR i = 0,4 DO plots,[corrleaf[0,i]^2,corrleaf[0,i]^2],[slope_lon[0,i]+sigma_lon[0,i],slope_lon[0,i]-sigma_lon[0,i]]
  oplot,[0.3025,0.3025],[-100,100],linestyle=2
  oplot,[-1,1],[0,0]
  plot,[0,0],[0,0],/nodata,yrange=[-15,10],xrange=[-0.05,0.95],ytitle='Slope (gC m-2 day-2)',xtitle='r2'
  FOR i = 0,4 DO plots,corrleaf[1,i]^2,slope_loff[0,i],psym=symbol(i,/fill,size=(i+2)/2.0),color=50.*i
  FOR i = 0,4 DO plots,[corrleaf[1,i]^2,corrleaf[1,i]^2],[slope_loff[0,i]+sigma_loff[0,i],slope_loff[0,i]-sigma_loff[0,i]]
  oplot,[0.3025,0.3025],[-100,100],linestyle=2
  oplot,[-1,1],[0,0]
  plot,[0,0],[0,0],/nodata,yrange=[-15,15],xrange=[-0.05,0.95],ytitle='Slope (gC m-2 day-2)',xtitle='r2'
  FOR i = 0,4 DO plots,corrleaf[2,i]^2,slope_gsl[0,i],psym=symbol(i,/fill,size=(i+2)/2.0),color=50.*i
  FOR i = 0,4 DO plots,[corrleaf[2,i]^2,corrleaf[2,i]^2],[slope_gsl[0,i]+sigma_gsl[0,i],slope_gsl[0,i]-sigma_gsl[0,i]]
  oplot,[0.3025,0.3025],[-100,100],linestyle=2
  oplot,[-1,1],[0,0]

; oplot,[0.399,0.399],[-100,100],linestyle=1
;  oplot,[0.513,0.513],[-100,100],linestyle=3

stop

;another possible fig - LAI time series for WC vs model (restore leaf data)


;Table 1 - sites - DONE - need mean flux and uncertainty
print,'Table 1'
print,'Site NEE NEE_err IAV'
FOR i = 0,4 DO print,snames[i],round(mean(f.anee[i,*],/nan)),round(mean(anee_err[i,*],/nan)),round(stddev(f.anee[i,*],/nan)),format='(a6," ",i0," ",i0," ",i0)'
print

;Table 2 -  param names and priors - IN DOC - need fix
print,'Table 2'
print,'laimax: 5.3 3.7 4.1 4.9 3.7'
print,'laimin: 0.0 0.0 0.5 0.0 0.5'
print
pname = ['LUE','k','LAImin','LAImax','Tmin','Topt','VPDmax','VPDmin','alpha','GDDthresh','beta','TEMPthresh','rs','rv','b1','b2','b3']
baseparam = [0.25,0.5,0.0,5.0,0.0,15.0,3000.0,100.0,0.05,200.0,0.1,4.0,2.0,2.0,0.025,0.025,0.05]
changeable = [1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1]
maxparam = [1.0,0.75,10.0,10.0,10.0,40.0,20000.0,2000.0,0.5,400.0,0.5,20.0,5.0,5.0,0.5,0.25,0.25]
minparam = [0.001,0.001,0.0,0.0,-15.0,5.0,0.0,0.0,0.05,10.0,0.05,0.0,0.1,0.1,0.0,0.0,0.0]
print,['Name','Default','Min','Max'],format='(4a12)'
FOR i = 0,n_elements(pname)-1 DO print,pname[i],baseparam[i],minparam[i],maxparam[i],format='(a12,3f12.2)'

;Table 3 - r2 and MAE for each site, daily and annual, for 
;NEED P-VALUES
nyrs = [7,5,5,6,8]

hcorr = [correl(f.nee[0,*],a.hln[0,*])^2,correl(f.nee[1,*],a.hln[1,*])^2,correl(f.nee[2,*],a.hln[2,*])^2,correl(f.nee[3,*],a.hln[3,*])^2,correl(f.nee[4,*],a.hln[4,*])^2]
acorr = [correl(f.nee[0,*],a.iln[0,*])^2,correl(f.nee[1,*],a.iln[1,*])^2,correl(f.nee[2,*],a.iln[2,*])^2,correl(f.nee[3,*],a.iln[3,*])^2,correl(f.nee[4,*],a.iln[4,*])^2]
scorr = [correl(f.nee[0,*],s.iln[0,*])^2,correl(f.nee[1,*],s.iln[1,*])^2,correl(f.nee[2,*],s.iln[2,*])^2,correl(f.nee[3,*],s.iln[3,*])^2,correl(f.nee[4,*],s.iln[4,*])^2]

hcorr_y = [correl(f.an_nee[0,*],a.an_hla[0,*])^2,correl(f.an_nee[1,*],a.an_hla[1,*])^2,correl(f.an_nee[2,*],a.an_hla[2,*])^2,correl(f.an_nee[3,*],a.an_hla[3,*])^2,correl(f.an_nee[4,*],a.an_hla[4,*])^2]
acorr_y = [correl(f.an_nee[0,*],a.an_ila[0,*])^2,correl(f.an_nee[1,*],a.an_ila[1,*])^2,correl(f.an_nee[2,*],a.an_ila[2,*])^2,correl(f.an_nee[3,*],a.an_ila[3,*])^2,correl(f.an_nee[4,*],a.an_ila[4,*])^2]
scorr_y = [correl(f.an_nee[0,*],s.an_ila[0,*])^2,correl(f.an_nee[1,*],s.an_ila[1,*])^2,correl(f.an_nee[2,*],s.an_ila[2,*])^2,correl(f.an_nee[3,*],s.an_ila[3,*])^2,correl(f.an_nee[4,*],s.an_ila[4,*])^2]

hcorr_p = fltarr(5)
FOR i = 0,4 DO hcorr_p[i] = corr_ttest(corr=abs(hcorr_y[i]),n=nyrs[i])
acorr_p = fltarr(5)
FOR i = 0,4 DO acorr_p[i] = corr_ttest(corr=abs(acorr_y[i]),n=nyrs[i])
scorr_p = fltarr(5)
FOR i = 0,4 DO scorr_p[i] = corr_ttest(corr=abs(scorr_y[i]),n=nyrs[i])

print,'Table 3'
print,'Daily r2'
print,'Hsync ',hcorr,format='(a10,5f10.2)'
print,'Async ',acorr,format='(a10,5f10.2)'
print,'Sync ',scorr,format='(a10,5f10.2)'
print,'Daily MAE'
print,'Hsync ',mean(abs(f.nee[0,*]-a.hln[0,*]),/nan),$
      mean(abs(f.nee[1,*]-a.hln[1,*]),/nan),$
      mean(abs(f.nee[2,*]-a.hln[2,*]),/nan),$
      mean(abs(f.nee[3,*]-a.hln[3,*]),/nan),$
      mean(abs(f.nee[4,*]-a.hln[4,*]),/nan),format='(a10,5f10.2)'
print,'Async ',mean(abs(f.nee[0,*]-a.iln[0,*]),/nan),$
      mean(abs(f.nee[1,*]-a.iln[1,*]),/nan),$
      mean(abs(f.nee[2,*]-a.iln[2,*]),/nan),$
      mean(abs(f.nee[3,*]-a.iln[3,*]),/nan),$
      mean(abs(f.nee[4,*]-a.iln[4,*]),/nan),format='(a10,5f10.2)'
print,'Sync ',mean(abs(f.nee[0,*]-s.iln[0,*]),/nan),$
      mean(abs(f.nee[1,*]-s.iln[1,*]),/nan),$
      mean(abs(f.nee[2,*]-s.iln[2,*]),/nan),$
      mean(abs(f.nee[3,*]-s.iln[3,*]),/nan),$
      mean(abs(f.nee[4,*]-s.iln[4,*]),/nan),format='(a10,5f10.2)'
print,'Annual r2'
print,'Hsync ',hcorr_y,format='(a10,5f10.2)'
print,'Async ',acorr_y,format='(a10,5f10.2)'
print,'Sync ',scorr_y,format='(a10,5f10.2)'
print,'Annual p'
print,'Hsync ',hcorr_p,format='(a10,5f10.2)'
print,'Async ',acorr_p,format='(a10,5f10.2)'
print,'Sync ',scorr_p,format='(a10,5f10.2)'
print,'Annual MAE'
print,'Hsync ',mean(abs(f.an_nee[0,*]-a.an_hla[0,*]),/nan),$
      mean(abs(f.an_nee[1,*]-a.an_hla[1,*]),/nan),$
      mean(abs(f.an_nee[2,*]-a.an_hla[2,*]),/nan),$
      mean(abs(f.an_nee[3,*]-a.an_hla[3,*]),/nan),$
      mean(abs(f.an_nee[4,*]-a.an_hla[4,*]),/nan),format='(a10,5f10.0)'
print,'Async ',mean(abs(f.an_nee[0,*]-a.an_ila[0,*]),/nan),$
      mean(abs(f.an_nee[1,*]-a.an_ila[1,*]),/nan),$
      mean(abs(f.an_nee[2,*]-a.an_ila[2,*]),/nan),$
      mean(abs(f.an_nee[3,*]-a.an_ila[3,*]),/nan),$
      mean(abs(f.an_nee[4,*]-a.an_ila[4,*]),/nan),format='(a10,5f10.0)'
print,'Sync ',mean(abs(f.an_nee[0,*]-s.an_ila[0,*]),/nan),$
      mean(abs(f.an_nee[1,*]-s.an_ila[1,*]),/nan),$
      mean(abs(f.an_nee[2,*]-s.an_ila[2,*]),/nan),$
      mean(abs(f.an_nee[3,*]-s.an_ila[3,*]),/nan),$
      mean(abs(f.an_nee[4,*]-s.an_ila[4,*]),/nan),format='(a10,5f10.0)'
print

  pname = ['LUE','k','LAImin','LAImax','Tmin','Topt','VPDmax','VPDmin','alpha','GDDthresh','beta','TEMPthresh','rs','rv','b1','b2','b3']
  baseparam = [0.25,0.5,0.0,5.0,0.0,15.0,3000.0,100.0,0.05,200.0,0.1,4.0,2.0,2.0,0.025,0.025,0.05]
  changeable = [1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1]
  maxparam = [1.0,0.75,10.0,10.0,10.0,40.0,20000.0,2000.0,0.5,400.0,0.5,20.0,5.0,5.0,0.5,0.25,0.25]
  minparam = [0.001,0.001,0.0,0.0,-15.0,5.0,0.0,0.0,0.05,10.0,0.05,0.0,0.1,0.1,0.0,0.0,0.0]
  pform = ['f0.3','f0.3','f0.3','f0.3','f0.1','f0.1','i0','i0','f0.3','f0.1','f0.3','f0.1','f0.2','f0.2','f0.4','f0.4','f0.4']

;Table 5 - a model HL - params all sites (best and range)
;order  alpga,GDD,beta,TEMPt, , LUE Tmin Topt VPDmax VPDmin, ,rs,rv,b1,b2,b3
  porder = [8,9,10,11,-1,0,4,5,6,7,-1,12,13,14,15,16]

  print,'Table 4'
  print,['Pname ',snames],format='(6a12)'
  FOR i = 0,n_elements(porder)-1 DO BEGIN
    IF porder[i] NE -1 THEN BEGIN
      tp = porder[i]
      out = fltarr(15)
      FOR j = 0,4 DO out[(j*3):((j*3)+2)] =[a.hlp[tp,j],min(a.hlap[tp,*,j],/nan),max(a.hlap[tp,*,j],/nan)]
      print,pname[tp],out,format='(a12,5(" ",'+pform[tp]+',"_(",'+pform[tp]+',"-",'+pform[tp]+',")"))'
    ENDIF ELSE print
  ENDFOR 
  print

;Table 4 - a model IL
  print,'Table 5'
  print,['Pname ',snames],format='(6a12)'
  FOR i = 0,n_elements(porder)-1 DO BEGIN
    IF porder[i] NE -1 THEN BEGIN
      tp = porder[i]
      out = fltarr(15)
      FOR j = 0,4 DO out[(j*3):((j*3)+2)] =[a.ilp[tp,j],min(a.ilap[tp,*,j],/nan),max(a.ilap[tp,*,j],/nan)]
      print,pname[tp],out,format='(a12,5(" ",'+pform[tp]+',"_(",'+pform[tp]+',"-",'+pform[tp]+',")"))'
    ENDIF ELSE print
  ENDFOR 
  print

;Table 6 - s model - params all sites (best and range)
  print,'Table 6'
  print,['Pname ',snames],format='(6a12)'
  FOR i = 0,n_elements(porder)-1 DO BEGIN
    IF porder[i] NE -1 THEN BEGIN
      tp = porder[i]
      out = fltarr(15)
      FOR j = 0,4 DO out[(j*3):((j*3)+2)] =[s.ilp[tp,j],min(s.ilap[tp,*,j],/nan),max(s.ilap[tp,*,j],/nan)]
      print,pname[tp],out,format='(a12,5(" ",'+pform[tp]+',"_(",'+pform[tp]+',"-",'+pform[tp]+',")"))'
    ENDIF ELSE print
  ENDFOR 
  print

  stop
END