forrest@0: ; ***********************************************
forrest@0: ; using gsn_table for all
forrest@0: ; output: line plot for each site (4 fields)
forrest@0: ;         table for M_score
forrest@0: ; ***********************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ;load "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"
forrest@0: ;************************************************
forrest@0: begin
forrest@0: 
forrest@0:   plot_type     = "ps"
forrest@0:   plot_type_new = "png"
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; read model data
forrest@0: ;************************************************
forrest@0:   model_name = "10cn"
forrest@0:   film = "i01.10cn_1948-2004_MONS_climo.nc"
forrest@0: 
forrest@0: ; model_name = "10casa"
forrest@0: ; film = "i01.10casa_1948-2004_MONS_climo.nc"
forrest@0: 
forrest@0:  system("sed s#model_name#"+model_name+"# table.html > table.html.new")
forrest@0:  system("mv -f table.html.new table.html")
forrest@0: ;--------------------------------------------
forrest@0:   dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:   fm    = addfile(dirm+film,"r")
forrest@0: 
forrest@0:   xm    = fm->lon
forrest@0:   ym    = fm->lat
forrest@0:   nlat = dimsizes(ym)
forrest@0:   nlon = dimsizes(xm)
forrest@0: 
forrest@0:   nmon      = 12
forrest@0:   nfield    = 4
forrest@0: 
forrest@0:   data_mod0 = new ((/nfield,nmon,nlat,nlon/),float)
forrest@0: 
forrest@0: ; change to unit of observed (u mol/m2/s)
forrest@0: ; Model_units [=] gC/m2/s
forrest@0: ; 12. = molecular weight of C
forrest@0: ; u mol = 1e-6 mol 
forrest@0:   data = fm->NEE
forrest@0: 
forrest@0:   factor = 1e6 /12.
forrest@0: 
forrest@0:   data_mod0(0,:,:,:) = data(:,:,:) * factor
forrest@0:   delete (data)
forrest@0: 
forrest@0:   data = fm->NETRAD
forrest@0:   data_mod0(1,:,:,:) = data(:,:,:) 
forrest@0:   delete (data)
forrest@0: 
forrest@0:   data = fm->LATENT
forrest@0:   data_mod0(2,:,:,:) = data(:,:,:) 
forrest@0:   delete (data)
forrest@0: 
forrest@0:   data = fm->SENSIBLE
forrest@0:   data_mod0(3,:,:,:) = data(:,:,:) 
forrest@0:   delete (data)
forrest@0: 
forrest@0: ; printVarSummary (data_mod0)
forrest@0: ;************************************************
forrest@0: ; read in data: observed
forrest@0: ;************************************************
forrest@0:  station = (/"BOREAS_NSA_OBS" \
forrest@0:             ,"CastelPorziano" \
forrest@0:             ,"Hyytiala" \
forrest@0:             ,"Kaamanen" \
forrest@0:             ,"LBA_Tapajos_KM67" \
forrest@0:             ,"Lethbridge" \
forrest@0:             ,"Tharandt" \
forrest@0:             ,"Vielsalm" \
forrest@0:             /)
forrest@0: 
forrest@0:  year_ob = (/"1994-2004" \
forrest@0:             ,"1997-2003" \
forrest@0:             ,"1996-2003" \
forrest@0:             ,"2000-2003" \
forrest@0:             ,"2002-2005" \
forrest@0:             ,"1999-2004" \
forrest@0:             ,"1996-2003" \
forrest@0:             ,"1998-2003" \
forrest@0:             /)
forrest@0: 
forrest@0:  field   = (/"CO2 Flux" \
forrest@0:             ,"Net Radiation" \
forrest@0:             ,"Latent Heat" \
forrest@0:             ,"Sensible Heat" \
forrest@0:             /)
forrest@0: 
forrest@0:  nstation  = dimsizes(station)
forrest@0:  nmon      = 12
forrest@0:  nfield    = dimsizes(field)
forrest@0: 
forrest@0:  data_ob   = new ((/nstation, nfield, nmon/),float)
forrest@0:  lat_ob    = new ((/nstation/),float)
forrest@0:  lon_ob    = new ((/nstation/),float)
forrest@0: 
forrest@0:  diri_root  = "/fis/cgd/cseg/people/jeff/clamp_data/fluxnet/"
forrest@0: 
forrest@0:  do n = 0,nstation-1
forrest@0:     diri = diri_root + station(n)+"/"
forrest@0:     fili = station(n)+"_"+year_ob(n)+"_monthly.nc"
forrest@0:     g     = addfile (diri+fili,"r")
forrest@0:  
forrest@0:     lon_ob(n) = g->lon 
forrest@0:     lat_ob(n) = g->lat
forrest@0: 
forrest@0:     data      = g->CO2_FLUX
forrest@0:     data_ob(n,0,:) = dim_avg(data(month|:,year|:))
forrest@0:     delete (data)
forrest@0: 
forrest@0:     data      = g->RAD_FLUX
forrest@0:     data_ob(n,1,:) = dim_avg(data(month|:,year|:))
forrest@0:     delete (data)
forrest@0: 
forrest@0:     data      = g->LH_FLUX
forrest@0:     data_ob(n,2,:) = dim_avg(data(month|:,year|:))
forrest@0:     delete (data)
forrest@0: 
forrest@0:     data      = g->SH_FLUX
forrest@0:     data_ob(n,3,:) = dim_avg(data(month|:,year|:))
forrest@0:     delete (data)
forrest@0: 
forrest@0:     delete (g)
forrest@0:  end do
forrest@0:  
forrest@0: ;print (lat_ob)
forrest@0: ;print (lon_ob)
forrest@0: ;printVarSummary (data_ob)
forrest@0: 
forrest@0: ;************************************************************
forrest@0: ; interpolate model data into observed station
forrest@0: ; note: model is 0-360E, 90S-90N
forrest@0: ;************************************************************
forrest@0: 
forrest@0: ; to be able to handle observation at (-89.98,-24.80)
forrest@0: ; print (ym(0))
forrest@0:   ym(0) = -90.  
forrest@0: 
forrest@0:   yy = linint2_points_Wrap(xm,ym,data_mod0,True,lon_ob,lat_ob,0)
forrest@0: 
forrest@0:   delete (data_mod0)
forrest@0:   yy!0 = "field"
forrest@0:   data_mod = yy(pts|:,field|:,time|:)
forrest@0: ; printVarSummary (data_mod)
forrest@0: 
forrest@0: ;************************************************************
forrest@0: ; compute correlation coef and M score
forrest@0: ;************************************************************
forrest@0: 
forrest@0:  score_max = 5.
forrest@0: 
forrest@0:  ccr     = new ((/nstation, nfield/),float)
forrest@0:  M_score = new ((/nstation, nfield/),float) 
forrest@0: 
forrest@0:  do n=0,nstation-1
forrest@0:  do m=0,nfield-1   
forrest@0:     ccr(n,m) = esccr(data_ob(n,m,:),data_mod(n,m,:),0)
forrest@0:     bias = sum(abs(data_mod(n,m,:)-data_ob(n,m,:))/(abs(data_mod(n,m,:))+abs(data_ob(n,m,:))))
forrest@0:     M_score(n,m) = (1. -(bias/nmon)) * score_max
forrest@0:  end do
forrest@0:  end do
forrest@0: 
forrest@0: ;*******************************************************************
forrest@0: ; for station line plot
forrest@0: ;*******************************************************************
forrest@0: 
forrest@0: ; for x-axis in xyplot
forrest@0:   mon = ispan(1,12,1)
forrest@0:   mon@long_name = "month"
forrest@0: 
forrest@0:   res                   = True               ; plot mods desired
forrest@0:   res@xyLineThicknesses = (/2.0,2.0/)        ; make 2nd lines thicker
forrest@0:   res@xyLineColors      = (/"blue","red"/)   ; line color (ob,model)
forrest@0: ;-------------------------------------------------------------------------
forrest@0: ; Add a boxed legend using the more simple method
forrest@0: 
forrest@0:   res@pmLegendDisplayMode    = "Always"
forrest@0: ; res@pmLegendWidthF         = 0.1
forrest@0:   res@pmLegendWidthF         = 0.08
forrest@0:   res@pmLegendHeightF        = 0.06
forrest@0: ; res@pmLegendOrthogonalPosF = -1.17
forrest@0: ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
forrest@0:   res@pmLegendOrthogonalPosF = -0.30  ;(downward)
forrest@0: 
forrest@0: ; res@pmLegendParallelPosF   =  0.18
forrest@0:   res@pmLegendParallelPosF   =  0.23  ;(rightward)
forrest@0: 
forrest@0: ; res@lgPerimOn             = False
forrest@0:   res@lgLabelFontHeightF     = 0.015
forrest@0:   res@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0: ;-------------------------------------------------------------------
forrest@0: ; for panel plot
forrest@0:   res@gsnFrame     = False                   ; Do not draw plot 
forrest@0:   res@gsnDraw      = False                   ; Do not advance frame
forrest@0: 
forrest@0:   pres                            = True     ; panel plot mods desired
forrest@0:   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
forrest@0:                                              ; indiv. plots in panel
forrest@0:   pres@gsnMaximize                = True     ; fill the page
forrest@0: ;-------------------------------------------------------------------
forrest@0: 
forrest@0:   plot_data   = new((/2,12/),float)
forrest@0:   plot_data!0 = "case"
forrest@0:   plot_data!1 = "month"
forrest@0: 
forrest@0:   do n = 0,nstation-1
forrest@0: ;----------------------------
forrest@0: ; for observed
forrest@0: 
forrest@0:     plot_name = station(n)+"_ob"    
forrest@0:     title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
forrest@0:     res@tiMainString = title
forrest@0: 
forrest@0:     wks = gsn_open_wks (plot_type,plot_name)
forrest@0:     plot=new(4,graphic)                        ; create graphic array   
forrest@0:                            
forrest@0:     plot_data(0,:) = (/data_ob (n,0,:)/)
forrest@0:     plot_data@long_name = field(0)   
forrest@0:     plot(0)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 1
forrest@0: 
forrest@0:     plot_data(0,:) = (/data_ob (n,1,:)/)
forrest@0:     plot_data@long_name = field(1)
forrest@0:     plot(1)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 2
forrest@0: 
forrest@0:     plot_data(0,:) = (/data_ob (n,2,:)/)
forrest@0:     plot_data@long_name = field(2)   
forrest@0:     plot(2)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 3
forrest@0: 
forrest@0:     plot_data(0,:) = (/data_ob (n,3,:)/)
forrest@0:     plot_data@long_name = field(3)
forrest@0:     plot(3)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 4
forrest@0: 
forrest@0:     gsn_panel(wks,plot,(/2,2/),pres)                 ; create panel plot
forrest@0: 
forrest@0:     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:     system("rm "+plot_name+"."+plot_type)
forrest@0:     clear (wks)  
forrest@0:     delete (plot)
forrest@0: ;----------------------------
forrest@0: ; for model_vs_ob
forrest@0: 
forrest@0:     plot_name = station(n)+"_model_vs_ob"
forrest@0:     title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"    
forrest@0:     res@tiMainString = title
forrest@0: 
forrest@0:     wks = gsn_open_wks (plot_type,plot_name)
forrest@0:     plot=new(4,graphic)                        ; create graphic array   
forrest@0:                            
forrest@0:     plot_data(0,:) = (/data_ob (n,0,:)/)
forrest@0:     plot_data(1,:) = (/data_mod(n,0,:)/)
forrest@0:     plot_data@long_name = field(0)   
forrest@0:     plot(0)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 1
forrest@0: 
forrest@0:     plot_data(0,:) = (/data_ob (n,1,:)/)
forrest@0:     plot_data(1,:) = (/data_mod(n,1,:)/)
forrest@0:     plot_data@long_name = field(1)
forrest@0:     plot(1)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 2
forrest@0: 
forrest@0:     plot_data(0,:) = (/data_ob (n,2,:)/)
forrest@0:     plot_data(1,:) = (/data_mod(n,2,:)/)
forrest@0:     plot_data@long_name = field(2)   
forrest@0:     plot(2)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 3
forrest@0: 
forrest@0:     plot_data(0,:) = (/data_ob (n,3,:)/)
forrest@0:     plot_data(1,:) = (/data_mod(n,3,:)/)
forrest@0:     plot_data@long_name = field(3)
forrest@0:     plot(3)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 4
forrest@0: 
forrest@0:     gsn_panel(wks,plot,(/2,2/),pres)                 ; create panel plot
forrest@0: 
forrest@0:     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:     system("rm "+plot_name+"."+plot_type)
forrest@0:     clear (wks)  
forrest@0:     delete (plot)
forrest@0:  end do
forrest@0: 
forrest@0: ;*******************************************************************
forrest@0: ; for table of site score
forrest@0: ;*******************************************************************
forrest@0: 
forrest@0:   table_length = 0.8 
forrest@0: 
forrest@0: ; table header name
forrest@0:   table_header_name = "Site" 
forrest@0: 
forrest@0: ; column (not including header column)
forrest@0:   col_header = (/"Latitude","Longitude","observed" \
forrest@0:                 ,"CO2 Flux","Net Radiation","Latent Heat" \
forrest@0:                 ,"Sensible Heat","Average" \
forrest@0:                 /)
forrest@0:   ncol = dimsizes(col_header) 
forrest@0: 
forrest@0: ; row (not including header row) 
forrest@0:   nrow       = nstation + 1
forrest@0:   row_header = new ((/nrow/),string )
forrest@0:   row_header(0:nstation-1) = station(:)
forrest@0:   row_header(nrow-1)       = "All Sites" 
forrest@0: 
forrest@0: ; Table header
forrest@0:   ncr1  = (/1,1/)               ; 1 row, 1 column
forrest@0:   x1    = (/0.005,0.15/)        ; Start and end X
forrest@0:   y1    = (/0.800,0.895/)       ; Start and end Y
forrest@0:   text1 = table_header_name
forrest@0:   res1               = True
forrest@0:   res1@txFontHeightF = 0.03
forrest@0:   res1@gsFillColor   = "CornFlowerBlue"
forrest@0: 
forrest@0: ; Column header (equally space in x2)
forrest@0:   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns
forrest@0:   x2    = (/x1(1),0.995/)       ; start from end of x1
forrest@0:   y2    = y1                    ; same as y1
forrest@0:   text2 = col_header
forrest@0:   res2               = True
forrest@0:   res2@txFontHeightF = 0.010
forrest@0:   res2@gsFillColor   = "Gray"
forrest@0: 
forrest@0: ; Row header (equally space in y2)
forrest@0:   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns
forrest@0:   x3    = x1                         ; same as x1
forrest@0:   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1
forrest@0:   text3 = row_header
forrest@0:   res3               = True
forrest@0:   res3@txFontHeightF = 0.010
forrest@0:   res3@gsFillColor   = "Gray"
forrest@0: 
forrest@0: ; Main table body
forrest@0:   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns
forrest@0:   x4    = x2                      ; Start and end x
forrest@0:   y4    = y3                      ; Start and end Y
forrest@0:   text4 = new((/nrow,ncol/),string)
forrest@0: 
forrest@0:   color_fill4           = new((/nrow,ncol/),string)
forrest@0:   color_fill4           = "white"
forrest@0:   color_fill4(:,ncol-1) = "grey"
forrest@0:   color_fill4(nrow-1,:) = "green"
forrest@0: 
forrest@0:   res4               = True       ; Set up resource list
forrest@0: ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.
forrest@0:   res4@txFontHeightF = 0.015
forrest@0:   res4@gsFillColor   = color_fill4
forrest@0: 
forrest@0:   delete (color_fill4)
forrest@0: ;-------------------------------------------------------------------
forrest@0: ; for table value
forrest@0: 
forrest@0:   do n = 0,nrow-2
forrest@0:      text4(n,0) = sprintf("%5.2f", lat_ob(n))  
forrest@0:      text4(n,1) = sprintf("%5.2f", lon_ob(n))  
forrest@0:      text4(n,2) = year_ob(n)  
forrest@0:      text4(n,3) = sprintf("%5.2f", M_score(n,0))      ; CO2 Flux
forrest@0:      text4(n,4) = sprintf("%5.2f", M_score(n,1))      ; Net Radiation
forrest@0:      text4(n,5) = sprintf("%5.2f", M_score(n,2))      ; Latent Heat
forrest@0:      text4(n,6) = sprintf("%5.2f", M_score(n,3))      ; Sensible Heat
forrest@0:      text4(n,7) = sprintf("%5.2f", avg(M_score(n,:))) ; avg all fields    
forrest@0:   end do
forrest@0: 
forrest@0: ;    for the last row
forrest@0: 
forrest@0:      M_co2 = avg(M_score(:,0))
forrest@0:      M_rad = avg(M_score(:,1))
forrest@0:      M_lh  = avg(M_score(:,2))
forrest@0:      M_sh  = avg(M_score(:,3))
forrest@0:      M_all = M_co2+ M_rad +M_lh + M_sh
forrest@0:      
forrest@0:      text4(nrow-1,0) = "-"
forrest@0:      text4(nrow-1,1) = "-"
forrest@0:      text4(nrow-1,2) = "-"
forrest@0:      text4(nrow-1,3) = sprintf("%5.2f", M_co2) ; avg all sites
forrest@0:      text4(nrow-1,4) = sprintf("%5.2f", M_rad) ; avg all sites
forrest@0:      text4(nrow-1,5) = sprintf("%5.2f", M_lh ) ; avg all sites
forrest@0:      text4(nrow-1,6) = sprintf("%5.2f", M_sh ) ; avg all sites
forrest@0:      text4(nrow-1,7) = sprintf("%5.2f", M_all) ; avg all sites
forrest@0: 
forrest@0:      print (M_co2)
forrest@0:      print (M_rad)
forrest@0:      print (M_lh )
forrest@0:      print (M_sh)
forrest@0:      print (M_all)
forrest@0: 
forrest@0:   M_energy_co2 = sprintf("%.2f", M_co2)
forrest@0:   system("sed s#M_energy_co2#"+M_energy_co2+"# table.html > table.html.new")
forrest@0:   system("mv -f table.html.new table.html")
forrest@0:   print (M_energy_co2)
forrest@0:  
forrest@0:   M_energy_rad = sprintf("%.2f", M_rad)
forrest@0:   system("sed s#M_energy_rad#"+M_energy_rad+"# table.html > table.html.new")
forrest@0:   system("mv -f table.html.new table.html")
forrest@0:   print (M_energy_rad)
forrest@0:  
forrest@0:   M_energy_lh = sprintf("%.2f", M_lh)
forrest@0:   system("sed s#M_energy_lh#"+M_energy_lh+"# table.html > table.html.new")
forrest@0:   system("mv -f table.html.new table.html")
forrest@0:   print (M_energy_lh)
forrest@0:  
forrest@0:   M_energy_sh = sprintf("%.2f", M_sh)
forrest@0:   system("sed s#M_energy_sh#"+M_energy_sh+"# table.html > table.html.new")
forrest@0:   system("mv -f table.html.new table.html")
forrest@0:   print (M_energy_sh) 
forrest@0: ;---------------------------------------------------------------------------
forrest@0:  
forrest@0:   plot_name = "table_site_score"
forrest@0:  
forrest@0:   wks = gsn_open_wks (plot_type,plot_name)
forrest@0: ;------------------------------------------
forrest@0: ; for table title
forrest@0: 
forrest@0:   gRes               = True
forrest@0:   gRes@txFontHeightF = 0.02
forrest@0: ; gRes@txAngleF      = 90
forrest@0: 
forrest@0:   title_text = "Model  " + model_name +  "  M_Score"
forrest@0: 
forrest@0:   gsn_text_ndc(wks,title_text,0.50,0.95,gRes)
forrest@0: ;------------------------------------------   
forrest@0: 
forrest@0:   gsn_table(wks,ncr1,x1,y1,text1,res1)
forrest@0:   gsn_table(wks,ncr2,x2,y2,text2,res2)
forrest@0:   gsn_table(wks,ncr3,x3,y3,text3,res3)
forrest@0:   gsn_table(wks,ncr4,x4,y4,text4,res4) 
forrest@0: 
forrest@0:   frame(wks)
forrest@0: 
forrest@0:   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:   system("rm "+plot_name+"."+plot_type)
forrest@0: ;-------------------------------------------------------------------
forrest@0:   temp_name = "energy." + model_name
forrest@0:   system("mkdir -p " + temp_name)
forrest@0:   system("cp table.html " + temp_name)
forrest@0:   system("mv *.png " + temp_name)
forrest@0:   system("tar cf "+ temp_name +".tar " + temp_name)
forrest@0: ;------------------------------------------------------------------- 
forrest@0: 
forrest@0: end