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

;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_l= "/fis/cgd/cseg/people/jeff/surface_data/" 

 film_l = "lnd_T42.nc"

 fm_l   = addfile (dirm_l+film_l,"r")

 landfrac = fm_l->landfrac

 area     = fm_l->area

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

  area = area * 1.e6             

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

; read data: model, group 1

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

 film = model_name1 + "_1980-2004_ANN_climo.nc"

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

 NPP1   = fm->NPP

 leafc  = fm->LEAFC

 woodc  = fm->WOODC

 frootc = fm->FROOTC

 VegC   = leafc

 VegC   = leafc + woodc + frootc

 litterc = fm->LITTERC

 cwdc    = fm->CWDC

 LiCwC   = litterc

 LiCwC   = litterc + cwdc

 SoilC   = fm->SOILC

 delete (fm)

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

; read data: model, group 2

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

 film = model_name2 + "_1990-2004_ANN_climo.nc"

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

 NPP2   = fm->NPP 

 NEE2   = fm->NEE 

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

; Units for these variables are:

;NPP1: g C/m^2/s

;NPP2: g C/m^2/s

;NEE2: g C/m^2/s

;VegC:  g C/m^2

;LiCwC: g C/m^2

;SoilC: g C/m^2

 nsec_per_year = 60*60*24*365

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

 NPP1 = NPP1 *  nsec_per_year

 NPP2 = NPP2 *  nsec_per_year

 NEE2 = NEE2 *  nsec_per_year 

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

; take into account landfrac

 area(:,:)    = area(:,:)    * landfrac(:,:)

 NPP1(0,:,:)  = NPP1(0,:,:)  * landfrac(:,:)

 VegC(0,:,:)  = VegC(0,:,:)  * landfrac(:,:)

 LiCwC(0,:,:) = LiCwC(0,:,:) * landfrac(:,:)

 SoilC(0,:,:) = SoilC(0,:,:) * landfrac(:,:)

 NPP2(0,:,:)  = NPP2(0,:,:)  * landfrac(:,:)

 NEE2(0,:,:)  = NEE2(0,:,:)  * landfrac(:,:)

 data_n = 7

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

; Calculate "nice" bins for binning the data in equally spaced ranges

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

; using model biome class

  nclass      = nclass_mod

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

; print (range)

; 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

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

  count   = new((/data_n,nx/),float)

  do n = 0,data_n-1

    if(n.eq.0) then

      data = ndtooned(area)

    end if

    if(n.eq.1) then

      data = ndtooned(NPP1)

    end if

    if(n.eq.2) then

      data = ndtooned(VegC)

    end if

    if(n.eq.3) then

      data = ndtooned(LiCwC)

    end if

    if(n.eq.4) then

      data = ndtooned(SoilC)

    end if

    if(n.eq.5) then

      data = ndtooned(NPP2)

    end if

    if(n.eq.6) then

      data = ndtooned(NEE2)

    end if

; Loop through each range, using base.

    do i=0,nx-1

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

;        print("")

;        print("In range ["+range(i)+","+range(i+1)+")")

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

      else

;        print("")

;        print("In range ["+range(i)+",)")

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

      end if

;     Calculate average

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

        if (n.eq.0) then 

           yvalues(n,i) = sum(data(idx))

        else 

           yvalues(n,i) = avg(data(idx))

        end if

        count(n,i)   = dimsizes(idx)

      else

        yvalues(n,i) = yvalues@_FillValue

        count(n,i)   = 0

      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(n,i) = yvalues@_FillValue

         count(n,i)   = 0

      end if

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

;     print(n + ": " + count + " points, avg = " + yvalues(n,i))

      delete(idx)

    end do 

    delete(data)

  end do

  delete (base)

  delete (area)

  delete (NPP1)

  delete (VegC)

  delete (LiCwC)

  delete (SoilC)

  delete (NPP2)

  delete (NEE2)

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

; data for table1

 good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))

;print (good)

 w = yvalues(0,:)

 area_g = w(good)

 w = yvalues(1,:)

 NPP1_g = w(good)

 w = yvalues(2,:)

 VegC_g = w(good)

 w = yvalues(3,:)

 LiCwC_g = w(good)

 w = yvalues(4,:)

 SoilC_g = w(good)

 w = yvalues(5,:)

 NPP2_g = w(good)

 w = yvalues(6,:)

 NEE2_g = w(good) 

 n_biome = dimsizes(NPP1_g)

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

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

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

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

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

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

 NPP_ratio = NPP2_g/NPP1_g

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

; data for table2

; change unit from g to Pg (Peta gram)

 factor_unit = 1.e-15

 NPP1_t    = NPP1_g  * area_g * factor_unit

 VegC_t    = VegC_g  * area_g * factor_unit

 LiCwC_t   = LiCwC_g * area_g * factor_unit

 SoilC_t   = SoilC_g * area_g * factor_unit

 NEE2_t    = NEE2_g  * area_g * factor_unit

 print (NPP1_t)

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

; html table1 data

; column (not including header column)

  col_head  = (/"Area (1.e12m2)" \

               ,"NPP (gC/m2/yr)" \

               ,"VegC (gC/m2)" \

               ,"Litter+CWD (gC/m2)" \

               ,"SoilC (gC/m2)" \

               ,"NPP_ratio" \

               ,"NEE (gC/m2/yr)" \

               /)

  ncol = dimsizes(col_head)

; row (not including header row)                   

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

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

 do i=0,nrow-2

  text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12)

  text4(i,1) = sprintf("%.1f",NPP1_g(i))

  text4(i,2) = sprintf("%.1f",VegC_g(i))

  text4(i,3) = sprintf("%.1f",LiCwC_g(i))

  text4(i,4) = sprintf("%.1f",SoilC_g(i))

  text4(i,5) = sprintf("%.1f",NPP_ratio(i))

  text4(i,6) = sprintf("%.1f",NEE2_g(i))

 end do

  text4(nrow-1,0) = "-"

  text4(nrow-1,1) = "-"

  text4(nrow-1,2) = "-"

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

; create html table1

  header_text = "<H1>NEE and Carbon Stocks and Fluxes:  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>" \

       ,"   <th bgcolor=DDDDDD >"+col_head(6)+"</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-2

     set_line(lines,nline,row_header)

     txt0  = row_head(n)

     txt1  = text4(n,0)

     txt2  = text4(n,1)

     txt3  = text4(n,2)

     txt4  = text4(n,3)

     txt5  = text4(n,4)

     txt6  = text4(n,5)

     txt7  = text4(n,6)

     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,"<th>"+txt7+"</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_carbon_sink1.html"

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

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

    asciiwrite(output_html,lines(idx))

  else

   print ("error?")

  end if

  delete (idx)

  delete (col_head)

  delete (row_head)

  delete (text4)

  delete (table_header)

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

; html table2 data

; column (not including header column)

  col_head  = (/"NPP (PgC/yr)" \

               ,"VegC (PgC)" \

               ,"Litter+CWD (PgC)" \

               ,"SoilC (PgC)" \

               ,"NEE (PgC/yr)" \

               ,"NPP timeseries" \

               ,"NEE timeseries" \

               ,"Fire timeseries" \

               /)

  ncol = dimsizes(col_head)

; row (not including header row)                   

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

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

 do i=0,nrow-2

  text4(i,0) = sprintf("%.1f",NPP1_t(i))

  text4(i,1) = sprintf("%.1f",VegC_t(i))

  text4(i,2) = sprintf("%.1f",LiCwC_t(i))

  text4(i,3) = sprintf("%.1f",SoilC_t(i))

  text4(i,4) = sprintf("%.1f",NEE2_t(i))

  text4(i,5) = "-"

  text4(i,6) = "-"

  text4(i,7) = "-"

 end do

  text4(nrow-1,0) = sprintf("%.1f",sum(NPP1_t))

  text4(nrow-1,1) = sprintf("%.1f",sum(VegC_t))

  text4(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t))

  text4(nrow-1,3) = sprintf("%.1f",sum(SoilC_t))

  text4(nrow-1,4) = sprintf("%.1f",sum(NEE2_t))

  text4(nrow-1,5) = "-"

  text4(nrow-1,6) = "-"

  text4(nrow-1,7) = "-"

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

; create html table1

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

  header_text = "<H1>NEE and Carbon Stocks and Fluxes (per biome):  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>" \

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

       ,"   <th bgcolor=DDDDDD >"+col_head(7)+"</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  = text4(n,0)

     txt2  = text4(n,1)

     txt3  = text4(n,2)

     txt4  = text4(n,3)

     txt5  = text4(n,4)

     txt6  = text4(n,5)

     txt7  = text4(n,6)

     txt8  = text4(n,7)

     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,"<th>"+txt7+"</th>")

     set_line(lines,nline,"<th>"+txt8+"</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_carbon_sink2.html"

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

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

    asciiwrite(output_html,lines(idx))

  else

   print ("error?")

  end if

end