; ****************************************************

 ; combine scatter, histogram, global and zonal plots

 ; compute all correlation coef and M score

 ; ****************************************************

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

 ;************************************************

 begin

  plot_type     = "ps"

  plot_type_new = "png"

 ;************************************************

 ; read in model data

 ;************************************************

 ;film = "b30.061n_1995-2004_ANN_climo_lnd.nc"

 ;model_name = "b30.061n"

 ;model_grid = "T31"

 ;film = "newcn05_ncep_1i_ANN_climo_lnd.nc"

 ;model_name = "newcn"

 ;model_grid = "1.9"

 ;--------------------------------------------

  film = "i01.10cn_1948-2004_ANN_climo.nc"

  model_name = "10cn"

  model_grid = "T42"

  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

  fm   = addfile (dirm+film,"r")

 ;unit: gC/m2  

  data1  = fm->LIVESTEMC

  data2  = fm->DEADSTEMC

  data3  = fm->LEAFC

  datamod0 = data1(0,:,:)

  datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)

 ;-------------------------------------------

 ;film = "i01.10casa_1948-2004_ANN_climo.nc"

 ;model_name = "10casa"

 ;model_grid = "T42"

 ;dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

 ;fm   = addfile (dirm+film,"r")

 ;unit: gC/m2

 ;factor_WOODC = 0.7  

 ;data1  = fm->WOODC

 ;data2  = fm->LEAFC

 ;datamod0 = data1(0,:,:)

 ;datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)

 ;----------------------------------------------------

  html_name = "table.html." + model_name

  html_new  = html_name +".new"

  system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \

          "mv -f "+html_new+" "+html_name)

  xm       = fm->lon  

  ym       = fm->lat

 ;************************************************

 ; read in ob data

 ;************************************************

  ob_name = "LC15_Amazon_Biomass"

  diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"

  filo = "amazon_biomass_"+model_grid+".nc"

  fo   = addfile (diro+filo,"r")

  dataob   = fo->BIOMASS

  xo       = fo->lon  

  yo       = fo->lat

 ;************************************************

 ; Units for these variables are:

 ; dataob   : MgC/ha

 ; datamod0 : gC/m2

 ; We want to convert these to KgC/m2

 ; ha = 100m*100m = 10,000 m2

 ; MgC/ha*1000/10,000 = KgC/m2

   factor_aboveground = 0.5

   factor_unit_ob     = 0.1

   factor_unit_mod    = 0.001         

   dataob   = dataob * factor_aboveground * factor_unit_ob

   datamod0 = datamod0 * factor_unit_mod 

   dataob@units      = "KgC/m2"

   datamod0@units    = "KgC/m2"

   dataob@long_name      = "Amazon Biomass"

   datamod0@long_name    = "Amazon Biomass"

 ;********************************************************

 ; get subset of datamod0 that match dataob

   nlon = dimsizes(xo)

   nlat = dimsizes(yo)

   ind_lonL = ind(xm .eq. xo(0))

   ind_lonR = ind(xm .eq. xo(nlon-1))

   ind_latS = ind(ym .eq. yo(0))

   ind_latN = ind(ym .eq. yo(nlat-1))

   datamod  = dataob

   datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)

 ; print (datamod)

 ;---------------------------------------------------------------------- 

 ; global contour 

 ;res

   resg                     = True             ; Use plot options

   resg@cnFillOn            = True             ; Turn on color fill

   resg@gsnSpreadColors     = True             ; use full colormap

 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

 ; resg@lbLabelAutoStride   = True

   resg@cnLinesOn           = False            ; Turn off contourn lines

   resg@mpFillOn            = False            ; Turn off map fill

   resg@gsnAddCyclic        = False

   resg@gsnSpreadColors      = True            ; use full colormap

   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   resg@cnMinLevelValF       = 0.              ; Min level

   resg@cnMaxLevelValF       = 30.             ; Max level

   resg@cnLevelSpacingF      = 2.              ; interval

   resg@mpMinLatF            = -21.1      ; range to zoom in on

   resg@mpMaxLatF            =  13.8

   resg@mpMinLonF            =  277.28

   resg@mpMaxLonF            =  326.43

 ;------------------------------------------------------------------------

 ;global contour ob

   plot_name = "global_ob"

   title     = ob_name+" "+ model_grid

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   plot = gsn_csm_contour_map_ce(wks,dataob,resg)   

   frame(wks)

   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

          "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

          "rm "+plot_name+"-*."+plot_type_new+";"+ \

          "rm "+plot_name+"."+plot_type)

   clear (wks)

 ;------------------------------------------------------------------------

 ;global contour model

   plot_name = "global_model"

   title     = "Model "+ model_name 

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   plot = gsn_csm_contour_map_ce(wks,datamod,resg)   

   frame(wks)

   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

          "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \

          "rm "+plot_name+"-*."+plot_type_new+";"+ \

          "rm "+plot_name+"."+plot_type)

   clear (wks)

 ;------------------------------------------------------------------------

 ;global contour model vs ob

   plot_name = "global_model_vs_ob"

   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   delete (plot)

   plot=new(3,graphic)                        ; create graphic array

   resg@gsnFrame             = False          ; Do not draw plot 

   resg@gsnDraw              = False          ; Do not advance frame

 ;(d) compute correlation coef and M score

   uu = ndtooned(datamod)

   vv = ndtooned(dataob)

   good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))

   ug = uu(good)

   vg = vv(good)

   ccrG = esccr(ug,vg,0)

 ; print (ccrG)

   MG   = (ccrG*ccrG)* 5.0

 ; new eq

   bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))

   MG   = (1. - (bias/dimsizes(ug)))*5.

   M_biomass = sprintf("%.2f", MG)

   system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \

          "mv -f "+html_new+" "+html_name)

   print (M_biomass)

 ; plot correlation coef

   gRes  = True

   gRes@txFontHeightF = 0.02

   gRes@txAngleF = 90

   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"

   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

 ;--------------------------------------------------------------------

 ;(a) ob

   title     = ob_name+" "+ model_grid

   resg@tiMainString  = title

   plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)       

 ;(b) model

   title     = "Model "+ model_name

   resg@tiMainString  = title

   plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg) 

 ;(c) model-ob

   zz = datamod

   zz = datamod - dataob

   title = "Model_"+model_name+" - Observed"

   resg@cnMinLevelValF  = -10.          ; Min level

   resg@cnMaxLevelValF  =  10.          ; Max level

   resg@cnLevelSpacingF =  2.           ; interval

   resg@tiMainString    = title

   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 

   pres                            = True        ; panel plot mods desired

 ; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

                                                 ; indiv. plots in panel

   pres@gsnMaximize                = True        ; fill the page

   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot

   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

          "rm "+plot_name+"."+plot_type)

   clear (wks)

   delete (plot)

 ;-------------------------------------------------------------------

   temp_name = "biomass." + model_name

   system("mkdir -p " + temp_name+";"+ \

          "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \

          "mv *.png " + temp_name +";"+ \ 

          "tar cf "+ temp_name +".tar " + temp_name)

 ;------------------------------------------------------------------- 

 end