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

; 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,:,:)

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

 system("sed s#model_name#"+model_name+"# table.html > table.html.new")

 system("mv -f table.html.new table.html")

 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)

; system("rm "+plot_name+"."+plot_type)

  system("rm "+plot_name+"-1."+plot_type_new)

  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  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)

; system("rm "+plot_name+"."+plot_type)

  system("rm "+plot_name+"-1."+plot_type_new)

  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  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+"# table.html > table.html.new")

  system("mv -f table.html.new table.html")

  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)

; system("rm "+plot_name+"."+plot_type)

; system("rm "+plot_name+"-1."+plot_type_new)

; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  frame (wks)

  clear (wks)

  delete (plot)

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

  temp_name = "biomass." + model_name

  system("mkdir -p " + temp_name)

  system("cp table.html " + temp_name)

  system("mv *.png " + temp_name)

  system("tar cf "+ temp_name +".tar " + temp_name)

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

end