forrest@0: ; ****************************************************
forrest@0: ; combine scatter, histogram, global and zonal plots
forrest@0: ; compute all correlation coef and M score
forrest@0: ; ****************************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ;************************************************
forrest@0: begin
forrest@0: 
forrest@0:  plot_type     = "ps"
forrest@0:  plot_type_new = "png"
forrest@0: ;************************************************
forrest@0: ; read in model data
forrest@0: ;************************************************
forrest@0: ;film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
forrest@0: ;model_name = "b30.061n"
forrest@0: ;model_grid = "T31"
forrest@0: 
forrest@0: ;film = "newcn05_ncep_1i_ANN_climo_lnd.nc"
forrest@0: ;model_name = "newcn"
forrest@0: ;model_grid = "1.9"
forrest@0: ;--------------------------------------------
forrest@0: ;film = "i01.10cn_1948-2004_ANN_climo.nc"
forrest@0: ;model_name = "10cn"
forrest@0: ;model_grid = "T42"
forrest@0: 
forrest@0: ;dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0: ;fm   = addfile (dirm+film,"r")
forrest@0: 
forrest@0: ;unit: gC/m2  
forrest@0: ;data1  = fm->LIVESTEMC
forrest@0: ;data2  = fm->DEADSTEMC
forrest@0: ;data3  = fm->LEAFC
forrest@0: ;datamod0 = data1(0,:,:)
forrest@0: ;datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
forrest@0: ;-------------------------------------------
forrest@0:  film = "i01.10casa_1948-2004_ANN_climo.nc"
forrest@0:  model_name = "10casa"
forrest@0:  model_grid = "T42"
forrest@0: 
forrest@0:  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:  fm   = addfile (dirm+film,"r")
forrest@0: 
forrest@0: ;unit: gC/m2
forrest@0:  factor_WOODC = 0.7  
forrest@0:  data1  = fm->WOODC
forrest@0:  data2  = fm->LEAFC
forrest@0:  datamod0 = data1(0,:,:)
forrest@0:  datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
forrest@0: ;----------------------------------------------------
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:  xm       = fm->lon  
forrest@0:  ym       = fm->lat
forrest@0: ;************************************************
forrest@0: ; read in ob data
forrest@0: ;************************************************
forrest@0:  ob_name = "LC15_Amazon_Biomass"
forrest@0: 
forrest@0:  diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
forrest@0:  filo = "amazon_biomass_"+model_grid+".nc"
forrest@0:  fo   = addfile (diro+filo,"r")
forrest@0:  
forrest@0:  dataob   = fo->BIOMASS
forrest@0:  xo       = fo->lon  
forrest@0:  yo       = fo->lat
forrest@0: ;************************************************
forrest@0: ; Units for these variables are:
forrest@0: ; dataob   : MgC/ha
forrest@0: ; datamod0 : gC/m2
forrest@0: ; We want to convert these to KgC/m2
forrest@0: ; ha = 100m*100m = 10,000 m2
forrest@0: ; MgC/ha*1000/10,000 = KgC/m2
forrest@0: 
forrest@0:   factor_aboveground = 0.5
forrest@0:   factor_unit_ob     = 0.1
forrest@0:   factor_unit_mod    = 0.001         
forrest@0: 
forrest@0:   dataob   = dataob * factor_aboveground * factor_unit_ob
forrest@0:   datamod0 = datamod0 * factor_unit_mod 
forrest@0: 
forrest@0:   dataob@units      = "KgC/m2"
forrest@0:   datamod0@units    = "KgC/m2"
forrest@0: 
forrest@0:   dataob@long_name      = "Amazon Biomass"
forrest@0:   datamod0@long_name    = "Amazon Biomass"
forrest@0: ;********************************************************
forrest@0: ; get subset of datamod0 that match dataob
forrest@0:   
forrest@0:   nlon = dimsizes(xo)
forrest@0:   nlat = dimsizes(yo)
forrest@0: 
forrest@0:   ind_lonL = ind(xm .eq. xo(0))
forrest@0:   ind_lonR = ind(xm .eq. xo(nlon-1))
forrest@0:   ind_latS = ind(ym .eq. yo(0))
forrest@0:   ind_latN = ind(ym .eq. yo(nlat-1))
forrest@0: 
forrest@0:   datamod  = dataob
forrest@0:   datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
forrest@0: 
forrest@0: ; print (datamod)
forrest@0:     
forrest@0: ;---------------------------------------------------------------------- 
forrest@0: ; global contour 
forrest@0: 
forrest@0: ;res
forrest@0:   resg                     = True             ; Use plot options
forrest@0:   resg@cnFillOn            = True             ; Turn on color fill
forrest@0:   resg@gsnSpreadColors     = True             ; use full colormap
forrest@0: ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color
forrest@0: ; resg@lbLabelAutoStride   = True
forrest@0:   resg@cnLinesOn           = False            ; Turn off contourn lines
forrest@0:   resg@mpFillOn            = False            ; Turn off map fill
forrest@0:   resg@gsnAddCyclic        = False
forrest@0: 
forrest@0:   resg@gsnSpreadColors      = True            ; use full colormap
forrest@0:   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
forrest@0:   resg@cnMinLevelValF       = 0.              ; Min level
forrest@0:   resg@cnMaxLevelValF       = 30.             ; Max level
forrest@0:   resg@cnLevelSpacingF      = 2.              ; interval
forrest@0: 
forrest@0:   resg@mpMinLatF            = -21.1      ; range to zoom in on
forrest@0:   resg@mpMaxLatF            =  13.8
forrest@0:   resg@mpMinLonF            =  277.28
forrest@0:   resg@mpMaxLonF            =  326.43
forrest@0: ;------------------------------------------------------------------------
forrest@0: ;global contour ob
forrest@0:   
forrest@0:   plot_name = "global_ob"
forrest@0:   title     = ob_name+" "+ model_grid
forrest@0:   resg@tiMainString  = title
forrest@0: 
forrest@0:   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0:   gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0: 
forrest@0:   plot = gsn_csm_contour_map_ce(wks,dataob,resg)   
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:   system("rm "+plot_name+"-1."+plot_type_new)
forrest@0:   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0: 
forrest@0:   clear (wks)
forrest@0: ;------------------------------------------------------------------------
forrest@0: ;global contour model
forrest@0: 
forrest@0:   plot_name = "global_model"
forrest@0:   title     = "Model "+ model_name 
forrest@0:   resg@tiMainString  = title
forrest@0: 
forrest@0:   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0:   gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0: 
forrest@0:   plot = gsn_csm_contour_map_ce(wks,datamod,resg)   
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:   system("rm "+plot_name+"-1."+plot_type_new)
forrest@0:   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0: 
forrest@0:   clear (wks)
forrest@0: ;------------------------------------------------------------------------
forrest@0: ;global contour model vs ob
forrest@0: 
forrest@0:   plot_name = "global_model_vs_ob"
forrest@0: 
forrest@0:   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0:   gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0: 
forrest@0:   delete (plot)
forrest@0:   plot=new(3,graphic)                        ; create graphic array
forrest@0: 
forrest@0:   resg@gsnFrame             = False          ; Do not draw plot 
forrest@0:   resg@gsnDraw              = False          ; Do not advance frame
forrest@0: 
forrest@0: ;(d) compute correlation coef and M score
forrest@0: 
forrest@0:   uu = ndtooned(datamod)
forrest@0:   vv = ndtooned(dataob)
forrest@0:  
forrest@0:   good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))
forrest@0: 
forrest@0:   ug = uu(good)
forrest@0:   vg = vv(good)
forrest@0: 
forrest@0:   ccrG = esccr(ug,vg,0)
forrest@0: ; print (ccrG)
forrest@0: 
forrest@0:   MG   = (ccrG*ccrG)* 5.0
forrest@0: ; new eq
forrest@0:   bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
forrest@0:   MG   = (1. - (bias/dimsizes(ug)))*5.
forrest@0: 
forrest@0:   M_biomass = sprintf("%.2f", MG)
forrest@0:   system("sed s#M_biomass#"+M_biomass+"# table.html > table.html.new")
forrest@0:   system("mv -f table.html.new table.html")
forrest@0:   print (M_biomass)
forrest@0: 
forrest@0: ; plot correlation coef
forrest@0: 
forrest@0:   gRes  = True
forrest@0:   gRes@txFontHeightF = 0.02
forrest@0:   gRes@txAngleF = 90
forrest@0: 
forrest@0:   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
forrest@0: 
forrest@0:   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0: ;--------------------------------------------------------------------
forrest@0:   
forrest@0: ;(a) ob
forrest@0: 
forrest@0:   title     = ob_name+" "+ model_grid
forrest@0:   resg@tiMainString  = title
forrest@0: 
forrest@0:   plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)       
forrest@0: 
forrest@0: ;(b) model
forrest@0: 
forrest@0:   title     = "Model "+ model_name
forrest@0:   resg@tiMainString  = title
forrest@0: 
forrest@0:   plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg) 
forrest@0: 
forrest@0: ;(c) model-ob
forrest@0: 
forrest@0:   zz = datamod
forrest@0:   zz = datamod - dataob
forrest@0:   title = "Model_"+model_name+" - Observed"
forrest@0: 
forrest@0:   resg@cnMinLevelValF  = -10.          ; Min level
forrest@0:   resg@cnMaxLevelValF  =  10.          ; Max level
forrest@0:   resg@cnLevelSpacingF =  2.           ; interval
forrest@0:   resg@tiMainString    = title
forrest@0: 
forrest@0:   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 
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:   gsn_panel(wks,plot,(/3,1/),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: ; system("rm "+plot_name+"-1."+plot_type_new)
forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
forrest@0: 
forrest@0:   frame (wks)
forrest@0:   clear (wks)
forrest@0:   delete (plot)
forrest@0: ;------------------------------------------------------------------------
forrest@0:   temp_name = "biomass." + 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: end