forrest@0: ;********************************************************
forrest@0: ;using model biome vlass
forrest@0: ;
forrest@0: ; required command line input parameters:
forrest@0: ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
forrest@0: ;
forrest@0: ; histogram normalized by rain and compute correleration
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: procedure set_line(lines:string,nline:integer,newlines:string) 
forrest@0: begin
forrest@0: ; add line to ascci/html file
forrest@0:     
forrest@0:   nnewlines = dimsizes(newlines)
forrest@0:   if(nline+nnewlines-1.ge.dimsizes(lines))
forrest@0:     print("set_line: bad index, not setting anything.") 
forrest@0:     return
forrest@0:   end if 
forrest@0:   lines(nline:nline+nnewlines-1) = newlines
forrest@0: ;  print ("lines = " + lines(nline:nline+nnewlines-1))
forrest@0:   nline = nline + nnewlines
forrest@0:   return 
forrest@0: end
forrest@0: ;**************************************************************
forrest@0: ; Main code.
forrest@0: begin
forrest@0:  
forrest@0:   plot_type     = "ps"
forrest@0:   plot_type_new = "png"
forrest@0: 
forrest@0: ;---------------------------------------------------
forrest@0: ; model name and grid       
forrest@0: 
forrest@0:   model_grid = "T42"
forrest@0: 
forrest@0:   model_name  = "cn"
forrest@0:   model_name1 = "i01.06cn"
forrest@0:   model_name2 = "i01.10cn"
forrest@0: 
forrest@0: ;---------------------------------------------------
forrest@0: ; get biome data: model
forrest@0: 
forrest@0:   biome_name_mod = "Model PFT Class"
forrest@0: 
forrest@0:   dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:   film = "class_pft_"+model_grid+".nc"
forrest@0:   fm   = addfile(dirm+film,"r")
forrest@0:  
forrest@0:   classmod = fm->CLASS_PFT
forrest@0: 
forrest@0:   delete (fm)
forrest@0: 
forrest@0: ; model data has 17 land-type classes
forrest@0: 
forrest@0:   nclass_mod = 17
forrest@0: 
forrest@0: ;--------------------------------------------------
forrest@0: ; get model data: landfrac and area
forrest@0: 
forrest@0:   dirm = "/fis/cgd/cseg/people/jeff/surface_data/" 
forrest@0:   film = "lnd_T42.nc"
forrest@0:   fm   = addfile (dirm+film,"r")
forrest@0:   
forrest@0:   landmask = fm->landmask
forrest@0:   landfrac = fm->landfrac
forrest@0:   area     = fm->area
forrest@0: 
forrest@0:   delete (fm)
forrest@0: 
forrest@0: ; change area from km**2 to m**2
forrest@0:   area = area * 1.e6             
forrest@0: 
forrest@0: ;----------------------------------------------------
forrest@0: ; read data: time series, model
forrest@0: 
forrest@0:  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:  film = model_name2 + "_Fire_C_1979-2004_monthly.nc"
forrest@0:  fm   = addfile (dirm+film,"r")
forrest@0: 
forrest@0:  data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:)
forrest@0: 
forrest@0:  delete (fm)
forrest@0: 
forrest@0: ; Units for these variables are:
forrest@0: ; g C/m^2/s
forrest@0: 
forrest@0: ; change unit to g C/m^2/month
forrest@0: 
forrest@0:   nsec_per_month = 60*60*24*30
forrest@0:  
forrest@0:   data_mod = data_mod * nsec_per_month 
forrest@0: 
forrest@0:   data_mod@unit = "gC/m2/month"
forrest@0: ;----------------------------------------------------
forrest@0: ; read data: time series, observed
forrest@0: 
forrest@0:  dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
forrest@0:  film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"
forrest@0:  fm   = addfile (dirm+film,"r")
forrest@0: 
forrest@0:  data_ob = fm->FIRE_C(0:7,:,:,:)
forrest@0: 
forrest@0:  delete (fm)
forrest@0: 
forrest@0:  ob_name = "GFEDv2"
forrest@0: 
forrest@0: ; Units for these variables are:
forrest@0: ; g C/m^2/month
forrest@0: 
forrest@0: ;---------------------------------------------------
forrest@0: ; take into account landfrac
forrest@0: 
forrest@0:  data_mod = data_mod * conform(data_mod, landfrac, (/2,3/))
forrest@0:  data_ob  = data_ob  * conform(data_ob,  landfrac, (/2,3/))
forrest@0: 
forrest@0: ;---------------------------------------------------
forrest@0: ; get time-mean
forrest@0:   
forrest@0:   x          = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:))
forrest@0:   data_mod_m = dim_avg_Wrap(       x(lat|:,lon|:,month|:))
forrest@0:   delete (x)
forrest@0: 
forrest@0:   x          = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:))
forrest@0:   data_ob_m  = dim_avg_Wrap(       x(lat|:,lon|:,month|:))
forrest@0:   delete (x)
forrest@0: 
forrest@0: ; printVarSummary(y)
forrest@0: 
forrest@0: ;----------------------------------------------------
forrest@0: ; compute correlation coef
forrest@0: 
forrest@0:   landmask_1d = ndtooned(landmask)
forrest@0:   data_mod_1d = ndtooned(data_mod_m)
forrest@0:   data_ob_1d  = ndtooned(data_ob_m)
forrest@0: 
forrest@0:   good = ind(landmask_1d .gt. 0.)
forrest@0: ; print (dimsizes(good))
forrest@0: 
forrest@0:   cc = esccr(data_mod_1d(good),data_ob_1d(good),0)
forrest@0: ; print (cc)
forrest@0: 
forrest@0:   delete (landmask_1d)
forrest@0:   delete (data_mod_1d)
forrest@0:   delete (data_ob_id)
forrest@0: 
forrest@0: ;----------------------------------------------------
forrest@0: ; compute M_global
forrest@0: 
forrest@0:   score_max = 1.
forrest@0: 
forrest@0:   Mscore = cc * cc * score_max
forrest@0: 
forrest@0:   M_global = sprintf("%.2f", Mscore)
forrest@0:  
forrest@0: ;----------------------------------------------------
forrest@0: ; global res
forrest@0: 
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@cnLinesOn            = False            ; Turn off contourn lines
forrest@0:   resg@mpFillOn             = False            ; Turn off map fill
forrest@0:   resg@cnLevelSelectionMode = "ManualLevels"   ; Manual contour invtervals
forrest@0:       
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)   
forrest@0:   gsn_define_colormap(wks,"gui_default")     
forrest@0: 
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: ;----------------------
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("%.2f", cc)+")"
forrest@0: 
forrest@0:   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0: 
forrest@0: ;-----------------------  
forrest@0: ; plot ob
forrest@0: 
forrest@0:   data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue)
forrest@0: 
forrest@0:   title     = ob_name
forrest@0:   resg@tiMainString  = title
forrest@0: 
forrest@0:   resg@cnMinLevelValF       = 1.             
forrest@0:   resg@cnMaxLevelValF       = 10.             
forrest@0:   resg@cnLevelSpacingF      = 1.
forrest@0: 
forrest@0:   plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg)       
forrest@0: 
forrest@0: ;-----------------------
forrest@0: ; plot model
forrest@0: 
forrest@0:   data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue)
forrest@0: 
forrest@0:   title     = "Model "+ model_name
forrest@0:   resg@tiMainString  = title
forrest@0: 
forrest@0:   resg@cnMinLevelValF       = 1.             
forrest@0:   resg@cnMaxLevelValF       = 10.             
forrest@0:   resg@cnLevelSpacingF      = 1.
forrest@0: 
forrest@0:   plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg) 
forrest@0: 
forrest@0: ;-----------------------
forrest@0: ; plot model-ob
forrest@0: 
forrest@0:      resg@cnMinLevelValF  = -8.           
forrest@0:      resg@cnMaxLevelValF  =  2.            
forrest@0:      resg@cnLevelSpacingF =  1.
forrest@0: 
forrest@0:   zz = data_ob_m
forrest@0:   zz = data_mod_m - data_ob_m
forrest@0:   title = "Model_"+model_name+" - Observed"
forrest@0:   resg@tiMainString    = title
forrest@0: 
forrest@0:   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 
forrest@0: 
forrest@0: ; plot panel
forrest@0: 
forrest@0:   pres                            = True        ; panel plot mods desired
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: exit
forrest@0: 
forrest@0: ; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: ;        "rm "+plot_name+"."+plot_type)
forrest@0:   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0: 
forrest@0:   clear (wks)
forrest@0:   delete (plot)
forrest@0: 
forrest@0:   delete (data_ob_m)
forrest@0:   delete (data_mod_m)
forrest@0:   delete (zz)
forrest@0: 
forrest@0:   resg@gsnFrame             = True          ; Do advance frame 
forrest@0:   resg@gsnDraw              = True          ; Do draw plot
forrest@0: 
forrest@0: end
forrest@0: