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

;using model biome vlass

;

; required command line input parameters:

;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl

;

; histogram normalized by rain and compute correleration

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

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"

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

procedure set_line(lines:string,nline:integer,newlines:string) 

begin

; add line to ascci/html file

  nnewlines = dimsizes(newlines)

  if(nline+nnewlines-1.ge.dimsizes(lines))

    print("set_line: bad index, not setting anything.") 

    return

  end if 

  lines(nline:nline+nnewlines-1) = newlines

;  print ("lines = " + lines(nline:nline+nnewlines-1))

  nline = nline + nnewlines

  return 

end

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

; Main code.

begin

  plot_type     = "ps"

  plot_type_new = "png"

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

; model name and grid       

  model_grid = "T42"

  model_name  = "cn"

  model_name1 = "i01.06cn"

  model_name2 = "i01.10cn"

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

; get biome data: model

  biome_name_mod = "Model PFT Class"

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

  film = "class_pft_"+model_grid+".nc"

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

  classmod = fm->CLASS_PFT

  delete (fm)

; model data has 17 land-type classes

  nclass_mod = 17

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

; get model data: landfrac and area

  dirm = "/fis/cgd/cseg/people/jeff/surface_data/" 

  film = "lnd_T42.nc"

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

  landmask = fm->landmask

  landfrac = fm->landfrac

  area     = fm->area

  delete (fm)

; change area from km**2 to m**2

  area = area * 1.e6             

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

; read data: time series, model

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

 film = model_name2 + "_Fire_C_1979-2004_monthly.nc"

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

 data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:)

 delete (fm)

; Units for these variables are:

; g C/m^2/s

; change unit to g C/m^2/month

  nsec_per_month = 60*60*24*30

  data_mod = data_mod * nsec_per_month 

  data_mod@unit = "gC/m2/month"

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

; read data: time series, observed

 dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"

 film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"

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

 data_ob = fm->FIRE_C(0:7,:,:,:)

 delete (fm)

 ob_name = "GFEDv2"

; Units for these variables are:

; g C/m^2/month

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

; take into account landfrac

 data_mod = data_mod * conform(data_mod, landfrac, (/2,3/))

 data_ob  = data_ob  * conform(data_ob,  landfrac, (/2,3/))

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

; get time-mean

  x          = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:))

  data_mod_m = dim_avg_Wrap(       x(lat|:,lon|:,month|:))

  delete (x)

  x          = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:))

  data_ob_m  = dim_avg_Wrap(       x(lat|:,lon|:,month|:))

  delete (x)

; printVarSummary(y)

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

; compute correlation coef

  landmask_1d = ndtooned(landmask)

  data_mod_1d = ndtooned(data_mod_m)

  data_ob_1d  = ndtooned(data_ob_m)

  good = ind(landmask_1d .gt. 0.)

; print (dimsizes(good))

  cc = esccr(data_mod_1d(good),data_ob_1d(good),0)

; print (cc)

  delete (landmask_1d)

  delete (data_mod_1d)

  delete (data_ob_id)

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

; compute M_global

  score_max = 1.

  Mscore = cc * cc * score_max

  M_global = sprintf("%.2f", Mscore)

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

; global res

  resg                      = True             ; Use plot options

  resg@cnFillOn             = True             ; Turn on color fill

  resg@gsnSpreadColors      = True             ; use full colormap

  resg@cnLinesOn            = False            ; Turn off contourn lines

  resg@mpFillOn             = False            ; Turn off map fill

  resg@cnLevelSelectionMode = "ManualLevels"   ; Manual contour invtervals

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

; global contour: model vs ob

  plot_name = "global_model_vs_ob"

  wks = gsn_open_wks (plot_type,plot_name)   

  gsn_define_colormap(wks,"gui_default")     

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

  resg@gsnFrame             = False          ; Do not draw plot 

  resg@gsnDraw              = False          ; Do not advance frame

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

; plot correlation coef

  gRes               = True

  gRes@txFontHeightF = 0.02

  gRes@txAngleF      = 90

  correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"

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

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

; plot ob

  data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue)

  title     = ob_name

  resg@tiMainString  = title

  resg@cnMinLevelValF       = 1.             

  resg@cnMaxLevelValF       = 10.             

  resg@cnLevelSpacingF      = 1.

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

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

; plot model

  data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue)

  title     = "Model "+ model_name

  resg@tiMainString  = title

  resg@cnMinLevelValF       = 1.             

  resg@cnMaxLevelValF       = 10.             

  resg@cnLevelSpacingF      = 1.

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

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

; plot model-ob

     resg@cnMinLevelValF  = -8.           

     resg@cnMaxLevelValF  =  2.            

     resg@cnLevelSpacingF =  1.

  zz = data_ob_m

  zz = data_mod_m - data_ob_m

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

  resg@tiMainString    = title

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

; plot panel

  pres                            = True        ; panel plot mods desired

  pres@gsnMaximize                = True        ; fill the page

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

exit

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

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

  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  clear (wks)

  delete (plot)

  delete (data_ob_m)

  delete (data_mod_m)

  delete (zz)

  resg@gsnFrame             = True          ; Do advance frame 

  resg@gsnDraw              = True          ; Do draw plot

end