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: ; edit table.html of current model for movel1_vs_model2
forrest@0: 
forrest@0:  if (isvar("compare")) then
forrest@0:     html_name2 = compare+"/table.html"  
forrest@0:     html_new2  = html_name2 +".new"
forrest@0:  end if
forrest@0: 
forrest@0: ;------------------------------------------------------
forrest@0: ; edit table.html for current model
forrest@0: 
forrest@0:  html_name = model_name+"/table.html"  
forrest@0:  html_new  = html_name +".new"
forrest@0: 
forrest@0: ;------------------------------------------------------ 
forrest@0: ; get biome data: model
forrest@0: 
forrest@0:   biome_name_mod = "Model PFT Class"
forrest@0: 
forrest@0:   film_c   = "class_pft_"+ model_grid +".nc"
forrest@0:   fm_c     = addfile (dirs+film_c,"r") 
forrest@0:   classmod = fm_c->CLASS_PFT
forrest@0: 
forrest@0:   delete (fm_c)
forrest@0: 
forrest@0: ; model data has 17 land-type classes
forrest@0:   nclass_mod = 17
forrest@0: 
forrest@0: ;--------------------------------
forrest@0: ; get model data: landmask, landfrac and area
forrest@0:  
forrest@0:   film_l   = "lnd_"+ model_grid +".nc"
forrest@0:   fm_l     = addfile (dirs+film_l,"r")  
forrest@0:   landmask = fm_l->landmask
forrest@0:   landfrac = fm_l->landfrac
forrest@0:   area     = fm_l->area
forrest@0: 
forrest@0:   delete (fm_l)
forrest@0: 
forrest@0: ; change area from km**2 to m**2
forrest@0:   area = area * 1.e6
forrest@0:              
forrest@0: ; take into account landfrac
forrest@0:   area     = area * landfrac
forrest@0: 
forrest@0: ;--------------------------------
forrest@0: ; read data: time series, model
forrest@0: 
forrest@0:  fm   = addfile (dirm+film7,"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 gC/m2/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@units = "gC/m2/month"
forrest@0: 
forrest@0: ;----------------------------------------------------
forrest@0: ; read data: time series, observed
forrest@0: 
forrest@0:  dir_f   = diro + "fire/"
forrest@0:  fil_f   = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"
forrest@0:  fm      = addfile (dir_f+fil_f,"r")
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: gC/m2/month
forrest@0: 
forrest@0:  data_ob@units = "gC/m2/month"
forrest@0: 
forrest@0: ;-------------------------------------------------------------
forrest@0: ; html table1 data
forrest@0: 
forrest@0: ; column (not including header column)
forrest@0: 
forrest@0:   col_head  = (/"Observed Fire_Flux (PgC/yr)" \
forrest@0:                ,"Model Fire_Flux (PgC/yr)" \
forrest@0:                ,"Correlation Coefficient" \
forrest@0:                ,"Ratio model/observed" \
forrest@0:                ,"M_score" \
forrest@0:                ,"Timeseries plot" \
forrest@0:                /)
forrest@0: 
forrest@0:   ncol = dimsizes(col_head)
forrest@0: 
forrest@0: ; row (not including header row)                   
forrest@0: 
forrest@0: ; using model biome class:  
forrest@0:   row_head  = (/"Not Vegetated" \
forrest@0:                ,"Needleleaf Evergreen Temperate Tree" \
forrest@0:                ,"Needleleaf Evergreen Boreal Tree" \
forrest@0: ;              ,"Needleleaf Deciduous Boreal Tree" \
forrest@0:                ,"Broadleaf Evergreen Tropical Tree" \
forrest@0:                ,"Broadleaf Evergreen Temperate Tree" \
forrest@0:                ,"Broadleaf Deciduous Tropical Tree" \
forrest@0:                ,"Broadleaf Deciduous Temperate Tree" \
forrest@0: ;              ,"Broadleaf Deciduous Boreal Tree" \
forrest@0: ;              ,"Broadleaf Evergreen Shrub" \
forrest@0:                ,"Broadleaf Deciduous Temperate Shrub" \
forrest@0:                ,"Broadleaf Deciduous Boreal Shrub" \
forrest@0:                ,"C3 Arctic Grass" \
forrest@0:                ,"C3 Non-Arctic Grass" \
forrest@0:                ,"C4 Grass" \
forrest@0:                ,"Corn" \
forrest@0: ;              ,"Wheat" \                      
forrest@0:                ,"All Biomes" \                
forrest@0:                /)  
forrest@0:   nrow = dimsizes(row_head)                  
forrest@0: 
forrest@0: ; arrays to be passed to table. 
forrest@0:   text = new ((/nrow, ncol/),string ) 
forrest@0: 
forrest@0: ;*****************************************************************
forrest@0: ; (A) get time-mean
forrest@0: ;*****************************************************************
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: ;----------------------------------------------------
forrest@0: ; compute correlation coef: space
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:   area_1d     = ndtooned(area)
forrest@0:   landfrac_1d = ndtooned(landfrac)
forrest@0: 
forrest@0:   good = ind(landmask_1d .gt. 0.)
forrest@0: 
forrest@0:   global_mod = sum(data_mod_1d(good)*area_1d(good)) * 1.e-15 * 12.
forrest@0:   global_ob  = sum(data_ob_1d(good) *area_1d(good)) * 1.e-15 * 12.
forrest@0: ; print (global_mod)
forrest@0: ; print (global_ob)  
forrest@0: 
forrest@0:   global_area= sum(area_1d)
forrest@0:   global_land= sum(area_1d(good))
forrest@0: ; print (global_area)
forrest@0: ; print (global_land)
forrest@0: 
forrest@0:   cc_space = esccr(data_mod_1d(good)*landfrac_1d(good),data_ob_1d(good)*landfrac_1d(good),0)
forrest@0: 
forrest@0:   delete (landmask_1d)
forrest@0:   delete (landfrac_1d)
forrest@0: ; delete (area_1d)
forrest@0:   delete (data_mod_1d)
forrest@0:   delete (data_ob_1d)
forrest@0:   delete (good)
forrest@0: 
forrest@0: ;----------------------------------------------------
forrest@0: ; compute M_global
forrest@0: 
forrest@0:   score_max = 1.
forrest@0: 
forrest@0:   Mscore1 = cc_space * cc_space * score_max
forrest@0: 
forrest@0:   M_global = sprintf("%.2f", Mscore1)
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_space)+")"
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: ; change from gC/m2/month to gC/m2/yr
forrest@0:   month_to_year = 12. 
forrest@0:  
forrest@0:   data_ob_m@units  = "gC/m2/yr"
forrest@0:   data_mod_m@units = "gC/m2/yr"
forrest@0: 
forrest@0:   data_ob_m = data_ob_m * month_to_year
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       = 10.             
forrest@0:   resg@cnMaxLevelValF       = 100.             
forrest@0:   resg@cnLevelSpacingF      = 10.
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 = data_mod_m * month_to_year
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       = 10.             
forrest@0:   resg@cnMaxLevelValF       = 100.             
forrest@0:   resg@cnLevelSpacingF      = 10.
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  = -80.           
forrest@0:   resg@cnMaxLevelValF  =  20.            
forrest@0:   resg@cnLevelSpacingF =  10.
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:   delete (wks)
forrest@0:   delete (plot)
forrest@0: 
forrest@0:   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0:         "rm "+plot_name+"."+plot_type)
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: ;*******************************************************************
forrest@0: ; (B) Time series : per biome
forrest@0: ;*******************************************************************
forrest@0: 
forrest@0:  data_n = 2
forrest@0: 
forrest@0:  dsizes = dimsizes(data_mod)
forrest@0:  nyear  = dsizes(0)
forrest@0:  nmonth = dsizes(1)
forrest@0:  ntime  = nyear * nmonth
forrest@0: 
forrest@0:  year_start = 1997
forrest@0:  year_end   = 2004
forrest@0:                 
forrest@0: ;-------------------------------------------
forrest@0: ; Calculate "nice" bins for binning the data
forrest@0: 
forrest@0: ; using model biome class
forrest@0:   nclass = nclass_mod
forrest@0: 
forrest@0:   range  = fspan(0,nclass,nclass+1)
forrest@0: 
forrest@0: ; Use this range information to grab all the values in a
forrest@0: ; particular range, and then take an average.
forrest@0: 
forrest@0:   nx = dimsizes(range) - 1
forrest@0: 
forrest@0: ;-------------------------------------------
forrest@0: ; put data into bins
forrest@0: 
forrest@0: ; using observed biome class
forrest@0: ; base  = ndtooned(classob)
forrest@0: ; using model biome class
forrest@0:   base  = ndtooned(classmod)
forrest@0: 
forrest@0: ; output
forrest@0: 
forrest@0:   area_bin = new((/nx/),float)
forrest@0:   yvalues  = new((/ntime,data_n,nx/),float)
forrest@0: 
forrest@0: ; Loop through each range, using base.
forrest@0: 
forrest@0:   do i=0,nx-1
forrest@0: 
forrest@0:      if (i.ne.(nx-1)) then
forrest@0:         idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
forrest@0:      else
forrest@0:         idx = ind(base.ge.range(i))
forrest@0:      end if
forrest@0: ;---------------------
forrest@0: ;    for area  
forrest@0: 
forrest@0:      if (.not.any(ismissing(idx))) then 
forrest@0:         area_bin(i) = sum(area_1d(idx))
forrest@0:      else
forrest@0:         area_bin(i) = area_bin@_FillValue
forrest@0:      end if
forrest@0: 
forrest@0: ;#############################################################
forrest@0: ; using model biome class:
forrest@0: ;     set the following 4 classes to _FillValue:
forrest@0: ;     (3)Needleleaf Deciduous Boreal Tree,
forrest@0: ;     (8)Broadleaf Deciduous Boreal Tree,
forrest@0: ;     (9)Broadleaf Evergreen Shrub,
forrest@0: ;     (16)Wheat
forrest@0: 
forrest@0:      if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
forrest@0:         area_bin(i) = area_bin@_FillValue
forrest@0:      end if
forrest@0: ;#############################################################  
forrest@0: 
forrest@0: ;---------------------
forrest@0: ; for data_mod and data_ob
forrest@0: 
forrest@0:   do n = 0,data_n-1
forrest@0: 
forrest@0:      t = -1
forrest@0:      do m = 0,nyear-1
forrest@0:      do k = 0,nmonth-1
forrest@0:     
forrest@0:         t = t + 1 
forrest@0: 
forrest@0:         if (n.eq.0) then
forrest@0:            data = ndtooned(data_ob(m,k,:,:))
forrest@0:         end if
forrest@0: 
forrest@0:         if (n.eq.1) then
forrest@0:            data = ndtooned(data_mod(m,k,:,:))
forrest@0:         end if
forrest@0: 
forrest@0: ;       Calculate average
forrest@0:  
forrest@0:         if (.not.any(ismissing(idx))) then 
forrest@0:            yvalues(t,n,i) = sum(data(idx)*area_1d(idx))
forrest@0:         else
forrest@0:            yvalues(t,n,i) = yvalues@_FillValue
forrest@0:         end if
forrest@0: 
forrest@0: ;#############################################################
forrest@0: ; using model biome class:
forrest@0: ;     set the following 4 classes to _FillValue:
forrest@0: ;     (3)Needleleaf Deciduous Boreal Tree,
forrest@0: ;     (8)Broadleaf Deciduous Boreal Tree,
forrest@0: ;     (9)Broadleaf Evergreen Shrub,
forrest@0: ;     (16)Wheat
forrest@0: 
forrest@0:         if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
forrest@0:            yvalues(t,n,i) = yvalues@_FillValue
forrest@0:         end if
forrest@0: ;#############################################################  
forrest@0: 
forrest@0:      end do
forrest@0:      end do
forrest@0: 
forrest@0:      delete(data)
forrest@0:   end do 
forrest@0: 
forrest@0:     delete(idx)
forrest@0:   end do
forrest@0: 
forrest@0:   delete (base)
forrest@0:   delete (data_mod)
forrest@0:   delete (data_ob)
forrest@0: 
forrest@0:   global_bin = sum(area_bin)
forrest@0: ; print (global_bin)
forrest@0: 
forrest@0: ;----------------------------------------------------------------
forrest@0: ; get area_good
forrest@0: 
forrest@0:   good = ind(.not.ismissing(area_bin))
forrest@0: 
forrest@0:   area_g = area_bin(good)  
forrest@0: 
forrest@0:   n_biome = dimsizes(good)
forrest@0: 
forrest@0:   global_good = sum(area_g)
forrest@0: ; print (global_good)
forrest@0: 
forrest@0: ;----------------------------------------------------------------
forrest@0: ; data for tseries plot
forrest@0: 
forrest@0:   yvalues_g = new((/ntime,data_n,n_biome/),float)
forrest@0: 
forrest@0:   yvalues_g@units = "TgC/month"
forrest@0: 
forrest@0: ; change unit to Tg C/month
forrest@0: ; change unit from g to Tg (Tera gram)
forrest@0:   factor_unit = 1.e-12
forrest@0: 
forrest@0:   yvalues_g = yvalues(:,:,good) * factor_unit
forrest@0: 
forrest@0:   delete (good)
forrest@0: 
forrest@0: ;-------------------------------------------------------------------
forrest@0: ; general settings for line plot
forrest@0: 
forrest@0:   res                   = True               
forrest@0:   res@xyDashPatterns    = (/0,0/)          ; make lines solid
forrest@0:   res@xyLineThicknesses = (/2.0,2.0/)      ; make lines thicker
forrest@0:   res@xyLineColors      = (/"blue","red"/) ; line color
forrest@0: 
forrest@0:   res@trXMinF   = year_start
forrest@0:   res@trXMaxF   = year_end + 1
forrest@0: 
forrest@0:   res@vpHeightF = 0.4                 ; change aspect ratio of plot
forrest@0: ; res@vpWidthF  = 0.8
forrest@0:   res@vpWidthF  = 0.75   
forrest@0: 
forrest@0:   res@tiMainFontHeightF = 0.025       ; size of title 
forrest@0: 
forrest@0:   res@tmXBFormat  = "f"               ; not to add trailing zeros
forrest@0: 
forrest@0: ; res@gsnMaximize = True
forrest@0: 
forrest@0: ;----------------------------------------------
forrest@0: ; Add a boxed legend using the simple method
forrest@0: 
forrest@0:   res@pmLegendDisplayMode    = "Always"
forrest@0: ; res@pmLegendWidthF         = 0.1
forrest@0:   res@pmLegendWidthF         = 0.08
forrest@0:   res@pmLegendHeightF        = 0.06
forrest@0:   res@pmLegendOrthogonalPosF = -1.17
forrest@0: ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
forrest@0: ; res@pmLegendOrthogonalPosF = -0.30  ;(downward)
forrest@0: 
forrest@0: ; res@pmLegendParallelPosF   =  0.18
forrest@0:   res@pmLegendParallelPosF   =  0.23  ;(rightward)
forrest@0:   res@pmLegendParallelPosF   =  0.73  ;(rightward)
forrest@0:   res@pmLegendParallelPosF   =  0.83  ;(rightward)
forrest@0: 
forrest@0: ; res@lgPerimOn             = False
forrest@0:   res@lgLabelFontHeightF     = 0.015
forrest@0:   res@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0: 
forrest@0: ;*******************************************************************
forrest@0: ; (A) time series plot: monthly ( 2 lines per plot)
forrest@0: ;*******************************************************************
forrest@0: 
forrest@0: ; x-axis in time series plot
forrest@0: 
forrest@0:   timeI = new((/ntime/),integer)
forrest@0:   timeF = new((/ntime/),float)
forrest@0:   timeI = ispan(1,ntime,1)
forrest@0:   timeF = year_start + (timeI-1)/12.
forrest@0:   timeF@long_name = "year" 
forrest@0: 
forrest@0:   plot_data = new((/2,ntime/),float)
forrest@0:   plot_data@long_name = "TgC/month"
forrest@0: 
forrest@0: ;----------------------------------------------
forrest@0: ; time series plot : per biome
forrest@0:  
forrest@0:   do m = 0, n_biome-1
forrest@0: 
forrest@0:      plot_name = "monthly_biome_"+ m
forrest@0: 
forrest@0:      wks = gsn_open_wks (plot_type,plot_name)   
forrest@0: 
forrest@0:      title = "Fire : "+ row_head(m)
forrest@0:      res@tiMainString = title
forrest@0: 
forrest@0:      plot_data(0,:) = yvalues_g(:,0,m)
forrest@0:      plot_data(1,:) = yvalues_g(:,1,m)
forrest@0:                                   
forrest@0:      plot = gsn_csm_xy(wks,timeF,plot_data,res)
forrest@0: 
forrest@0:      delete (wks)
forrest@0:      delete (plot)
forrest@0: 
forrest@0:      system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0:             "rm "+plot_name+"."+plot_type)
forrest@0:   end do
forrest@0: 
forrest@0: ;------------------------------------------
forrest@0: ; data for table : per biome
forrest@0: 
forrest@0: ; unit change from TgC/month to PgC/month
forrest@0:   unit_factor = 1.e-3
forrest@0: 
forrest@0:   score_max = 1.
forrest@0: 
forrest@0:   tmp_ob    = new((/ntime/),float)
forrest@0:   tmp_mod   = new((/ntime/),float)
forrest@0: 
forrest@0:   total_ob  = new((/n_biome/),float)
forrest@0:   total_mod = new((/n_biome/),float)
forrest@0:   Mscore2   = new((/n_biome/),float)
forrest@0: 
forrest@0:   do m = 0, n_biome-1
forrest@0: 
forrest@0:      tmp_ob  = yvalues_g(:,0,m) 
forrest@0:      tmp_mod = yvalues_g(:,1,m) 
forrest@0: 
forrest@0:      total_ob(m)  = avg(month_to_annual(tmp_ob, 0)) * unit_factor 
forrest@0:      total_mod(m) = avg(month_to_annual(tmp_mod,0)) * unit_factor
forrest@0:      
forrest@0:      cc_time = esccr(tmp_mod,tmp_ob,0)
forrest@0: 
forrest@0:      ratio = total_mod(m)/total_ob(m)
forrest@0: 
forrest@0:      good = ind(tmp_ob .ne. 0. .and. tmp_mod .ne. 0.)
forrest@0: 
forrest@0:      bias = sum( abs( tmp_mod(good)-tmp_ob(good) )/( abs(tmp_mod(good))+abs(tmp_ob(good)) ) )
forrest@0:      Mscore2(m) = (1.- (bias/dimsizes(good)))*score_max
forrest@0: 
forrest@0:      delete (good)
forrest@0:      
forrest@0:      text(m,0) = sprintf("%.2f",total_ob(m))
forrest@0:      text(m,1) = sprintf("%.2f",total_mod(m))
forrest@0:      text(m,2) = sprintf("%.2f",cc_time)
forrest@0:      text(m,3) = sprintf("%.2f",ratio)
forrest@0:      text(m,4) = sprintf("%.2f",Mscore2(m))
forrest@0:      text(m,5) = "<a href=./monthly_biome_"+m+".png>model_vs_ob</a>" 
forrest@0:   end do
forrest@0:  
forrest@0:   delete (tmp_ob)
forrest@0:   delete (tmp_mod)
forrest@0: 
forrest@0: ;--------------------------------------------
forrest@0: ; time series plot: all biome
forrest@0: 
forrest@0:      plot_name = "monthly_global"
forrest@0: 
forrest@0:      wks = gsn_open_wks (plot_type,plot_name)   
forrest@0: 
forrest@0:      title = "Fire : "+ row_head(n_biome)
forrest@0:      res@tiMainString = title
forrest@0: 
forrest@0:      do k = 0,ntime-1
forrest@0:         plot_data(0,k) = sum(yvalues_g(k,0,:))
forrest@0:         plot_data(1,k) = sum(yvalues_g(k,1,:))
forrest@0:      end do
forrest@0:                                   
forrest@0:      plot = gsn_csm_xy(wks,timeF,plot_data,res)
forrest@0: 
forrest@0:      delete (wks)
forrest@0:      delete (plot)
forrest@0: 
forrest@0:      system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0:             "rm "+plot_name+"."+plot_type)
forrest@0: 
forrest@0: ;------------------------------------------
forrest@0: ; data for table : global
forrest@0: 
forrest@0:   score_max = 1.
forrest@0: 
forrest@0:   tmp_ob  = ndtooned(yvalues_g(:,0,:))
forrest@0:   tmp_mod = ndtooned(yvalues_g(:,1,:))
forrest@0: 
forrest@0:   cc_time = esccr(tmp_mod,tmp_ob,0)
forrest@0: 
forrest@0:   ratio = sum(total_mod)/sum(total_ob) 
forrest@0: 
forrest@0:   good = ind(tmp_ob .ne. 0. .and. tmp_mod .ne. 0.)
forrest@0: 
forrest@0:   bias = sum( abs( tmp_mod(good)-tmp_ob(good) )/( abs(tmp_mod(good))+abs(tmp_ob(good)) ) )
forrest@0:   Mscore3 = (1.- (bias/dimsizes(good)))*score_max
forrest@0: 
forrest@0: ; print (Mscore3)
forrest@0: 
forrest@0:   delete (good) 
forrest@0: 
forrest@0:   text(nrow-1,0) = sprintf("%.2f",sum(total_ob))
forrest@0:   text(nrow-1,1) = sprintf("%.2f",sum(total_mod))
forrest@0:   text(nrow-1,2) = sprintf("%.2f",cc_time)
forrest@0:   text(nrow-1,3) = sprintf("%.2f",ratio)
forrest@0: ; text(nrow-1,4) = sprintf("%.2f",avg(Mscore2))
forrest@0:   text(nrow-1,4) = sprintf("%.2f",    Mscore3)
forrest@0:   text(nrow-1,5) = "<a href=./monthly_global.png>model_vs_ob</a>"
forrest@0: 
forrest@0: ;**************************************************
forrest@0: ; create html table
forrest@0: ;**************************************************
forrest@0: 
forrest@0:   header_text = "<H1>Fire Emissions from GFEDv2 (1997-2004) vs "+model_name+"</H1>" 
forrest@0: 
forrest@0:   header = (/"<HTML>" \
forrest@0:             ,"<HEAD>" \
forrest@0:             ,"<TITLE>CLAMP metrics</TITLE>" \
forrest@0:             ,"</HEAD>" \
forrest@0:             ,header_text \
forrest@0:             /) 
forrest@0:   footer = "</HTML>"
forrest@0: 
forrest@0:   table_header = (/ \
forrest@0:         "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
forrest@0:        ,"<tr>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >Biome Type</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
forrest@0:        ,"</tr>" \
forrest@0:        /)
forrest@0:   table_footer = "</table>"
forrest@0:   row_header = "<tr>"
forrest@0:   row_footer = "</tr>"
forrest@0: 
forrest@0:   lines = new(50000,string)
forrest@0:   nline = 0
forrest@0: 
forrest@0:   set_line(lines,nline,header)
forrest@0:   set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ;row of table
forrest@0: 
forrest@0:   do n = 0,nrow-1
forrest@0:      set_line(lines,nline,row_header)
forrest@0: 
forrest@0:      txt0  = row_head(n)
forrest@0:      txt1  = text(n,0)
forrest@0:      txt2  = text(n,1)
forrest@0:      txt3  = text(n,2)
forrest@0:      txt4  = text(n,3)
forrest@0:      txt5  = text(n,4)
forrest@0:      txt6  = text(n,5)
forrest@0: 
forrest@0:      set_line(lines,nline,"<th>"+txt0+"</th>")
forrest@0:      set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0:      set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0:      set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0:      set_line(lines,nline,"<th>"+txt4+"</th>")
forrest@0:      set_line(lines,nline,"<th>"+txt5+"</th>")
forrest@0:      set_line(lines,nline,"<th>"+txt6+"</th>")
forrest@0: 
forrest@0:      set_line(lines,nline,row_footer)
forrest@0:   end do
forrest@0: ;-----------------------------------------------
forrest@0:   set_line(lines,nline,table_footer)
forrest@0:   set_line(lines,nline,footer) 
forrest@0: 
forrest@0: ; Now write to an HTML file.
forrest@0: 
forrest@0:   output_html = "table_fire.html"
forrest@0: 
forrest@0:   idx = ind(.not.ismissing(lines))
forrest@0:   if(.not.any(ismissing(idx))) then
forrest@0:     asciiwrite(output_html,lines(idx))
forrest@0:   else
forrest@0:    print ("error?")
forrest@0:   end if
forrest@0: 
forrest@0:   delete (idx)
forrest@0: 
forrest@0: ;**************************************************************************************
forrest@0: ; update score
forrest@0: ;**************************************************************************************
forrest@0: 
forrest@0:   M_all  = Mscore1 + Mscore3
forrest@0:   M_fire = sprintf("%.2f", M_all) 
forrest@0: 
forrest@0:   if (isvar("compare")) then
forrest@0:      system("sed -e '1,/M_fire/s/M_fire/"+M_fire+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0:             "mv -f "+html_new2+" "+html_name2)
forrest@0:   end if
forrest@0: 
forrest@0:   system("sed s#M_fire#"+M_fire+"# "+html_name+" > "+html_new+";"+ \
forrest@0:          "mv -f "+html_new+" "+html_name)
forrest@0: 
forrest@0: ;***************************************************************************
forrest@0: ; get total score and write to file
forrest@0: ;***************************************************************************
forrest@0: 
forrest@0:   asciiwrite("M_save.fire", M_fire)
forrest@0: 
forrest@0:   delete (M_fire)
forrest@0: 
forrest@0: ;***************************************************************************
forrest@0: ; output plot and html
forrest@0: ;***************************************************************************
forrest@0:   output_dir = model_name+"/fire"
forrest@0: 
forrest@0:   system("mv *.png *.html " + output_dir) 
forrest@0: ;***************************************************************************
forrest@0: 
forrest@0: end
forrest@0: