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

; using observed biome class

; landfrac applied to area only.

;

; 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"

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

; read biome data: observed

  biome_name_ob = "MODIS LandCover"

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

  filo = "land_class_"+model_grid+".nc"

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

  classob = tofloat(fo->LAND_CLASS)               

  delete (fo)

; observed data has 20 land-type classes

  nclass_ob = 20

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

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

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

; take into account landfrac

 area = area * landfrac

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

; 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

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

; html table1 data

; column (not including header column)

  col_head  = (/"Observed Fire_Flux (PgC/yr)" \

               ,"Model Fire_Flux (PgC/yr)" \

               ,"Correlation Coefficient" \

               ,"Ratio model/observed" \

               ,"M_score" \

               ,"Timeseries plot" \

               /)

  ncol = dimsizes(col_head)

; row (not including header row)                   

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

; using observed biome class:  

  row_head  = (/"Evergreen Needleleaf Forests" \

               ,"Evergreen Broadleaf Forests" \

               ,"Deciduous Needleleaf Forest" \

               ,"Deciduous Broadleaf Forests" \

               ,"Mixed Forests" \                      

               ,"Closed Bushlands" \                   

               ,"Open Bushlands" \                     

               ,"Woody Savannas (S. Hem.)" \           

               ,"Savannas (S. Hem.)" \                 

               ,"Grasslands" \                         

               ,"Permanent Wetlands" \                 

               ,"Croplands" \                                           

               ,"Cropland/Natural Vegetation Mosaic" \             

               ,"Barren or Sparsely Vegetated" \                             

               ,"Woody Savannas (N. Hem.)" \           

               ,"Savannas (N. Hem.)" \

               ,"All Biome" \                

               /)

; using model biome class:  

; row_head  = (/"Not Vegetated" \

;              ,"Needleleaf Evergreen Temperate Tree" \

;              ,"Needleleaf Evergreen Boreal Tree" \

;;             ,"Needleleaf Deciduous Boreal Tree" \

;              ,"Broadleaf Evergreen Tropical Tree" \

;              ,"Broadleaf Evergreen Temperate Tree" \

;              ,"Broadleaf Deciduous Tropical Tree" \

;              ,"Broadleaf Deciduous Temperate Tree" \

;;             ,"Broadleaf Deciduous Boreal Tree" \

;;             ,"Broadleaf Evergreen Shrub" \

;              ,"Broadleaf Deciduous Temperate Shrub" \

;              ,"Broadleaf Deciduous Boreal Shrub" \

;              ,"C3 Arctic Grass" \

;              ,"C3 Non-Arctic Grass" \

;              ,"C4 Grass" \

;              ,"Corn" \

;;             ,"Wheat" \                      

;              ,"All Biome" \                

;              /)  

  nrow = dimsizes(row_head)                  

; arrays to be passed to table. 

  text = new ((/nrow, ncol/),string ) 

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

; (A) 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)

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

; compute correlation coef: space

  landmask_1d = ndtooned(landmask)

  data_mod_1d = ndtooned(data_mod_m)

  data_ob_1d  = ndtooned(data_ob_m )

  area_1d     = ndtooned(area)

  landfrac_1d = ndtooned(landfrac)

  good = ind(landmask_1d .gt. 0.)

  global_mod = sum(data_mod_1d(good)*area_1d(good)) * 1.e-15 * 12.

  global_ob  = sum(data_ob_1d(good) *area_1d(good)) * 1.e-15 * 12.

  print (global_mod)

  print (global_ob) 

  global_area= sum(area_1d)

  global_land= sum(area_1d(good))

  print (global_area)

  print (global_land)

  cc_space = esccr(data_mod_1d(good)*landfrac_1d(good),data_ob_1d(good)*landfrac_1d(good),0)

  delete (landmask_1d)

  delete (landfrac_1d)

; delete (area_1d)

  delete (data_mod_1d)

  delete (data_ob_1d)

  delete (good)

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

; compute M_global

  score_max = 1.

  Mscore1 = cc_space * cc_space * score_max

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

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

; 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_space)+")"

  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

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

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

  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

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

; (B) Time series : per biome

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

 data_n = 2

 dsizes = dimsizes(data_mod)

 nyear  = dsizes(0)

 nmonth = dsizes(1)

 ntime  = nyear * nmonth

 year_start = 1997

 year_end   = 2004

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

; Calculate "nice" bins for binning the data

; using ob biome class

  nclass = nclass_ob

; using model biome class

; nclass = nclass_mod

  range  = fspan(0,nclass,nclass+1)

; Use this range information to grab all the values in a

; particular range, and then take an average.

  nx = dimsizes(range) - 1

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

; put data into bins

; using observed biome class

  base  = ndtooned(classob)

; using model biome class

; base  = ndtooned(classmod)

; output

  area_bin = new((/nx/),float)

  yvalues  = new((/ntime,data_n,nx/),float)

; Loop through each range, using base.

  do i=0,nx-1

     if (i.ne.(nx-1)) then

        idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))

     else

        idx = ind(base.ge.range(i))

     end if

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

;    for area 

     if (.not.any(ismissing(idx))) then 

        area_bin(i) = sum(area_1d(idx))

     else

        area_bin(i) = area_bin@_FillValue

     end if

;#############################################################

;using observed biome class: 

;     set the following 4 classes to _FillValue:

;     Water Bodies(0), Urban and Build-Up(13),

;     Permenant Snow and Ice(15), Unclassified(17)

      if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then

         area_bin(i) = yvalues@_FillValue

      end if

;############################################################# 

;#############################################################

; using model biome class:

;     set the following 4 classes to _FillValue:

;     (3)Needleleaf Deciduous Boreal Tree,

;     (8)Broadleaf Deciduous Boreal Tree,

;     (9)Broadleaf Evergreen Shrub,

;     (16)Wheat

;    if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then

;       area_bin(i) = area_bin@_FillValue

;    end if

;#############################################################  

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

; for data_mod and data_ob

  do n = 0,data_n-1

     t = -1

     do m = 0,nyear-1

     do k = 0,nmonth-1

        t = t + 1 

        if (n.eq.0) then

           data = ndtooned(data_ob(m,k,:,:))

        end if

        if (n.eq.1) then

           data = ndtooned(data_mod(m,k,:,:))

        end if

;       Calculate average

        if (.not.any(ismissing(idx))) then

           yvalues(t,n,i) = sum(data(idx)*area_1d(idx)) 

        else

           yvalues(t,n,i) = yvalues@_FillValue

        end if

;#############################################################

;using observed biome class: 

;     set the following 4 classes to _FillValue:

;     Water Bodies(0), Urban and Build-Up(13),

;     Permenant Snow and Ice(15), Unclassified(17)

      if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then

         yvalues(t,n,i) = yvalues@_FillValue

      end if

;############################################################# 

;#############################################################

; using model biome class:

;     set the following 4 classes to _FillValue:

;     (3)Needleleaf Deciduous Boreal Tree,

;     (8)Broadleaf Deciduous Boreal Tree,

;     (9)Broadleaf Evergreen Shrub,

;     (16)Wheat

;       if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then

;          yvalues(t,n,i) = yvalues@_FillValue

;       end if

;#############################################################  

     end do

     end do

     delete(data)

  end do 

    delete(idx)

  end do

  delete (base)

  delete (data_mod)

  delete (data_ob)

  global_bin = sum(area_bin)

  print (global_bin)

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

; get area_good

  good = ind(.not.ismissing(area_bin))

  area_g = area_bin(good)  

  n_biome = dimsizes(good)

  global_good = sum(area_g)

  print (global_good)

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

; data for tseries plot

  yvalues_g = new((/ntime,data_n,n_biome/),float)

  yvalues_g@units = "TgC/month"

; change unit to Tg C/month

; change unit from g to Tg (Tera gram)

  factor_unit = 1.e-12

  yvalues_g = yvalues(:,:,good) * factor_unit

  delete (good)

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

; general settings for line plot

  res                   = True               

  res@xyDashPatterns    = (/0,0/)          ; make lines solid

  res@xyLineThicknesses = (/2.0,2.0/)      ; make lines thicker

  res@xyLineColors      = (/"blue","red"/) ; line color

  res@trXMinF   = year_start

  res@trXMaxF   = year_end + 1

  res@vpHeightF = 0.4                 ; change aspect ratio of plot

; res@vpWidthF  = 0.8

  res@vpWidthF  = 0.75   

  res@tiMainFontHeightF = 0.025       ; size of title 

  res@tmXBFormat  = "f"               ; not to add trailing zeros

; res@gsnMaximize = True

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

; Add a boxed legend using the simple method

  res@pmLegendDisplayMode    = "Always"

; res@pmLegendWidthF         = 0.1

  res@pmLegendWidthF         = 0.08

  res@pmLegendHeightF        = 0.06

  res@pmLegendOrthogonalPosF = -1.17

; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

; res@pmLegendOrthogonalPosF = -0.30  ;(downward)

; res@pmLegendParallelPosF   =  0.18

  res@pmLegendParallelPosF   =  0.23  ;(rightward)

  res@pmLegendParallelPosF   =  0.73  ;(rightward)

  res@pmLegendParallelPosF   =  0.83  ;(rightward)

; res@lgPerimOn             = False

  res@lgLabelFontHeightF     = 0.015

  res@xyExplicitLegendLabels = (/"observed",model_name/)

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

; (A) time series plot: monthly ( 2 lines per plot)

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

; x-axis in time series plot

  timeI = new((/ntime/),integer)

  timeF = new((/ntime/),float)

  timeI = ispan(1,ntime,1)

  timeF = year_start + (timeI-1)/12.

  timeF@long_name = "year" 

  plot_data = new((/2,ntime/),float)

  plot_data@long_name = "TgC/month"

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

; time series plot : per biome

  do m = 0, n_biome-1

     plot_name = "monthly_biome_"+ m

     wks = gsn_open_wks (plot_type,plot_name)   

     title = "Fire : "+ row_head(m)

     res@tiMainString = title

     plot_data(0,:) = yvalues_g(:,0,m)

     plot_data(1,:) = yvalues_g(:,1,m)

     plot = gsn_csm_xy(wks,timeF,plot_data,res)

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

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

     clear (wks)  

     delete (plot)

  end do

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

; data for table : per biome

; unit change from TgC/month to PgC/month

  unit_factor = 1.e-3

  score_max = 1.

  tmp_ob    = new((/ntime/),float)

  tmp_mod   = new((/ntime/),float)

  total_ob  = new((/n_biome/),float)

  total_mod = new((/n_biome/),float)

  Mscore2   = new((/n_biome/),float)

  do m = 0, n_biome-1

     tmp_ob  = yvalues_g(:,0,m) 

     tmp_mod = yvalues_g(:,1,m) 

     total_ob(m)  = avg(month_to_annual(tmp_ob, 0)) * unit_factor 

     total_mod(m) = avg(month_to_annual(tmp_mod,0)) * unit_factor

     cc_time = esccr(tmp_mod,tmp_ob,0)

     ratio = total_mod(m)/total_ob(m)

     good = ind(tmp_ob .ne. 0. .and. tmp_mod .ne. 0.)

     bias = sum( abs( tmp_mod(good)-tmp_ob(good) )/( abs(tmp_mod(good))+abs(tmp_ob(good)) ) )

     Mscore2(m) = (1.- (bias/dimsizes(good)))*score_max

     delete (good)

     text(m,0) = sprintf("%.2f",total_ob(m))

     text(m,1) = sprintf("%.2f",total_mod(m))

     text(m,2) = sprintf("%.2f",cc_time)

     text(m,3) = sprintf("%.2f",ratio)

     text(m,4) = sprintf("%.2f",Mscore2(m))

     text(m,5) = "<a href=./monthly_biome_"+m+".png>model_vs_ob</a>" 

  end do

  delete (tmp_ob)

  delete (tmp_mod)

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

; time series plot: all biome

     plot_name = "monthly_global"

     wks = gsn_open_wks (plot_type,plot_name)   

     title = "Fire : "+ row_head(n_biome)

     res@tiMainString = title

     do k = 0,ntime-1

        plot_data(0,k) = sum(yvalues_g(k,0,:))

        plot_data(1,k) = sum(yvalues_g(k,1,:))

     end do

     plot = gsn_csm_xy(wks,timeF,plot_data,res)

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

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

     clear (wks)  

     delete (plot)

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

; data for table : global

  score_max = 1.

  tmp_ob  = ndtooned(yvalues_g(:,0,:))

  tmp_mod = ndtooned(yvalues_g(:,1,:))

  cc_time = esccr(tmp_mod,tmp_ob,0)

  ratio = sum(total_mod)/sum(total_ob) 

  good = ind(tmp_ob .ne. 0. .and. tmp_mod .ne. 0.)

  bias = sum( abs( tmp_mod(good)-tmp_ob(good) )/( abs(tmp_mod(good))+abs(tmp_ob(good)) ) )

  Mscore3 = (1.- (bias/dimsizes(good)))*score_max

  print (Mscore3)

  delete (good)

  text(nrow-1,0) = sprintf("%.2f",sum(total_ob))

  text(nrow-1,1) = sprintf("%.2f",sum(total_mod))

  text(nrow-1,2) = sprintf("%.2f",cc_time)

  text(nrow-1,3) = sprintf("%.2f",ratio)

; text(nrow-1,4) = sprintf("%.2f",avg(Mscore2))

  text(nrow-1,4) = sprintf("%.2f",    Mscore3)

  text(nrow-1,5) = "<a href=./monthly_global.png>model_vs_ob</a>"

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

; create html table

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

  header_text = "<H1>Fire Emission (1997-2004):  Model "+model_name+"</H1>" 

  header = (/"<HTML>" \

            ,"<HEAD>" \

            ,"<TITLE>CLAMP metrics</TITLE>" \

            ,"</HEAD>" \

            ,header_text \

            /) 

  footer = "</HTML>"

  table_header = (/ \

        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \

       ,"<tr>" \

       ,"   <th bgcolor=DDDDDD >Biome Type</th>" \

       ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \

       ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \

       ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \

       ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \

       ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \

       ,"   <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \

       ,"</tr>" \

       /)

  table_footer = "</table>"

  row_header = "<tr>"

  row_footer = "</tr>"

  lines = new(50000,string)

  nline = 0

  set_line(lines,nline,header)

  set_line(lines,nline,table_header)

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

;row of table

  do n = 0,nrow-1

     set_line(lines,nline,row_header)

     txt0  = row_head(n)

     txt1  = text(n,0)

     txt2  = text(n,1)

     txt3  = text(n,2)

     txt4  = text(n,3)

     txt5  = text(n,4)

     txt6  = text(n,5)

     set_line(lines,nline,"<th>"+txt0+"</th>")

     set_line(lines,nline,"<th>"+txt1+"</th>")

     set_line(lines,nline,"<th>"+txt2+"</th>")

     set_line(lines,nline,"<th>"+txt3+"</th>")

     set_line(lines,nline,"<th>"+txt4+"</th>")

     set_line(lines,nline,"<th>"+txt5+"</th>")

     set_line(lines,nline,"<th>"+txt6+"</th>")

     set_line(lines,nline,row_footer)

  end do

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

  set_line(lines,nline,table_footer)

  set_line(lines,nline,footer) 

; Now write to an HTML file.

  output_html = "table_fire.html"

  idx = ind(.not.ismissing(lines))

  if(.not.any(ismissing(idx))) then

    asciiwrite(output_html,lines(idx))

  else

   print ("error?")

  end if

  delete (idx)

end