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: 
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: 
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_l= "/fis/cgd/cseg/people/jeff/surface_data/" 
forrest@0:  film_l = "lnd_T42.nc"
forrest@0:  fm_l   = addfile (dirm_l+film_l,"r")
forrest@0:   
forrest@0:  landfrac = fm_l->landfrac
forrest@0:  area     = fm_l->area
forrest@0: 
forrest@0: ; change area from km**2 to m**2
forrest@0:   area = area * 1.e6             
forrest@0: ;---------------------------------------------------
forrest@0: ; read data: model, group 1
forrest@0: 
forrest@0:  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:  film = model_name1 + "_1980-2004_ANN_climo.nc"
forrest@0:  fm   = addfile (dirm+film,"r")
forrest@0:  
forrest@0:  NPP1   = fm->NPP
forrest@0: 
forrest@0:  leafc  = fm->LEAFC
forrest@0:  woodc  = fm->WOODC
forrest@0:  frootc = fm->FROOTC
forrest@0:  VegC   = leafc
forrest@0:  VegC   = leafc + woodc + frootc
forrest@0: 
forrest@0:  litterc = fm->LITTERC
forrest@0:  cwdc    = fm->CWDC
forrest@0:  LiCwC   = litterc
forrest@0:  LiCwC   = litterc + cwdc
forrest@0: 
forrest@0:  SoilC   = fm->SOILC
forrest@0: 
forrest@0:  delete (fm)
forrest@0: ;--------------------------------------------------- 
forrest@0: ; read data: model, group 2
forrest@0: 
forrest@0:  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:  film = model_name2 + "_1990-2004_ANN_climo.nc"
forrest@0:  fm   = addfile (dirm+film,"r")
forrest@0: 
forrest@0:  NPP2   = fm->NPP 
forrest@0:  NEE2   = fm->NEE 
forrest@0: 
forrest@0: ;---------------------------------------------------
forrest@0: ; Units for these variables are:
forrest@0: 
forrest@0: ;NPP1: g C/m^2/s
forrest@0: ;NPP2: g C/m^2/s
forrest@0: ;NEE2: g C/m^2/s
forrest@0: 
forrest@0: ;VegC:  g C/m^2
forrest@0: ;LiCwC: g C/m^2
forrest@0: ;SoilC: g C/m^2
forrest@0: 
forrest@0:  nsec_per_year = 60*60*24*365
forrest@0: 
forrest@0: ; change unit to g C/m^2/year
forrest@0:   
forrest@0:  NPP1 = NPP1 *  nsec_per_year
forrest@0:  NPP2 = NPP2 *  nsec_per_year
forrest@0:  NEE2 = NEE2 *  nsec_per_year 
forrest@0: 
forrest@0: ;---------------------------------------------------
forrest@0: ; take into account landfrac
forrest@0: 
forrest@0:  area(:,:)    = area(:,:)    * landfrac(:,:)
forrest@0:  NPP1(0,:,:)  = NPP1(0,:,:)  * landfrac(:,:)
forrest@0:  VegC(0,:,:)  = VegC(0,:,:)  * landfrac(:,:)
forrest@0:  LiCwC(0,:,:) = LiCwC(0,:,:) * landfrac(:,:)
forrest@0:  SoilC(0,:,:) = SoilC(0,:,:) * landfrac(:,:)
forrest@0:  NPP2(0,:,:)  = NPP2(0,:,:)  * landfrac(:,:)
forrest@0:  NEE2(0,:,:)  = NEE2(0,:,:)  * landfrac(:,:)
forrest@0: 
forrest@0:  data_n = 7
forrest@0:                 
forrest@0: ;*******************************************************************
forrest@0: ; Calculate "nice" bins for binning the data in equally spaced ranges
forrest@0: ;********************************************************************
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: ; print (range)
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: 
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:   yvalues = new((/data_n,nx/),float)
forrest@0:   count   = new((/data_n,nx/),float)
forrest@0: 
forrest@0:   do n = 0,data_n-1
forrest@0: 
forrest@0:     if(n.eq.0) then
forrest@0:       data = ndtooned(area)
forrest@0:     end if
forrest@0: 
forrest@0:     if(n.eq.1) then
forrest@0:       data = ndtooned(NPP1)
forrest@0:     end if
forrest@0: 
forrest@0:     if(n.eq.2) then
forrest@0:       data = ndtooned(VegC)
forrest@0:     end if
forrest@0: 
forrest@0:     if(n.eq.3) then
forrest@0:       data = ndtooned(LiCwC)
forrest@0:     end if
forrest@0: 
forrest@0:     if(n.eq.4) then
forrest@0:       data = ndtooned(SoilC)
forrest@0:     end if
forrest@0: 
forrest@0:     if(n.eq.5) then
forrest@0:       data = ndtooned(NPP2)
forrest@0:     end if
forrest@0: 
forrest@0:     if(n.eq.6) then
forrest@0:       data = ndtooned(NEE2)
forrest@0:     end if
forrest@0: 
forrest@0: ; Loop through each range, using base.
forrest@0: 
forrest@0:     do i=0,nx-1
forrest@0:       if (i.ne.(nx-1)) then
forrest@0: ;        print("")
forrest@0: ;        print("In range ["+range(i)+","+range(i+1)+")")
forrest@0:          idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
forrest@0:       else
forrest@0: ;        print("")
forrest@0: ;        print("In range ["+range(i)+",)")
forrest@0:          idx = ind(base.ge.range(i))
forrest@0:       end if
forrest@0: 
forrest@0: ;     Calculate average
forrest@0:  
forrest@0:       if(.not.any(ismissing(idx))) then
forrest@0:         if (n.eq.0) then 
forrest@0:            yvalues(n,i) = sum(data(idx))
forrest@0:         else 
forrest@0:            yvalues(n,i) = avg(data(idx))
forrest@0:         end if
forrest@0: 
forrest@0:         count(n,i)   = dimsizes(idx)
forrest@0:       else
forrest@0:         yvalues(n,i) = yvalues@_FillValue
forrest@0:         count(n,i)   = 0
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(n,i) = yvalues@_FillValue
forrest@0:          count(n,i)   = 0
forrest@0:       end if
forrest@0: ;#############################################################  
forrest@0: 
forrest@0: ;     print(n + ": " + count + " points, avg = " + yvalues(n,i))
forrest@0: 
forrest@0:       delete(idx)
forrest@0:     end do 
forrest@0: 
forrest@0:     delete(data)
forrest@0:   end do
forrest@0: 
forrest@0:   delete (base)
forrest@0:   delete (area)
forrest@0:   delete (NPP1)
forrest@0:   delete (VegC)
forrest@0:   delete (LiCwC)
forrest@0:   delete (SoilC)
forrest@0:   delete (NPP2)
forrest@0:   delete (NEE2)
forrest@0: 
forrest@0: ;----------------------------------------------------------------
forrest@0: ; data for table1
forrest@0: 
forrest@0:  good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))
forrest@0: ;print (good)
forrest@0: 
forrest@0:  w = yvalues(0,:)
forrest@0:  area_g = w(good)
forrest@0: 
forrest@0:  w = yvalues(1,:)
forrest@0:  NPP1_g = w(good)
forrest@0: 
forrest@0:  w = yvalues(2,:)
forrest@0:  VegC_g = w(good)
forrest@0: 
forrest@0:  w = yvalues(3,:)
forrest@0:  LiCwC_g = w(good)
forrest@0: 
forrest@0:  w = yvalues(4,:)
forrest@0:  SoilC_g = w(good)
forrest@0: 
forrest@0:  w = yvalues(5,:)
forrest@0:  NPP2_g = w(good)
forrest@0: 
forrest@0:  w = yvalues(6,:)
forrest@0:  NEE2_g = w(good) 
forrest@0: 
forrest@0:  n_biome = dimsizes(NPP1_g)
forrest@0: 
forrest@0:  NPP1_t    = new((/n_biome/),float)
forrest@0:  VegC_t    = new((/n_biome/),float)
forrest@0:  LiCwC_t   = new((/n_biome/),float)
forrest@0:  SoilC_t   = new((/n_biome/),float)
forrest@0:  NEE2_t    = new((/n_biome/),float)
forrest@0:  NPP_ratio = new((/n_biome/),float)
forrest@0: 
forrest@0:  NPP_ratio = NPP2_g/NPP1_g
forrest@0: 
forrest@0: ;-----------------------------------------------------------------
forrest@0: ; data for table2
forrest@0: 
forrest@0: ; change unit from g to Pg (Peta gram)
forrest@0:  factor_unit = 1.e-15
forrest@0: 
forrest@0:  NPP1_t    = NPP1_g  * area_g * factor_unit
forrest@0:  VegC_t    = VegC_g  * area_g * factor_unit
forrest@0:  LiCwC_t   = LiCwC_g * area_g * factor_unit
forrest@0:  SoilC_t   = SoilC_g * area_g * factor_unit
forrest@0:  NEE2_t    = NEE2_g  * area_g * factor_unit
forrest@0: 
forrest@0:  print (NPP1_t)
forrest@0:  
forrest@0: ;-------------------------------------------------------------
forrest@0: ; html table1 data
forrest@0: 
forrest@0: ; column (not including header column)
forrest@0: 
forrest@0:   col_head  = (/"Area (1.e12m2)" \
forrest@0:                ,"NPP (gC/m2/yr)" \
forrest@0:                ,"VegC (gC/m2)" \
forrest@0:                ,"Litter+CWD (gC/m2)" \
forrest@0:                ,"SoilC (gC/m2)" \
forrest@0:                ,"NPP ratio" \
forrest@0:                ,"NEE (gC/m2/yr)" \
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 Biome" \                
forrest@0:                /)  
forrest@0:   nrow = dimsizes(row_head)                  
forrest@0: 
forrest@0: ; arrays to be passed to table. 
forrest@0:   text4 = new ((/nrow, ncol/),string )
forrest@0:  
forrest@0:  do i=0,nrow-2
forrest@0:   text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12)
forrest@0:   text4(i,1) = sprintf("%.1f",NPP1_g(i))
forrest@0:   text4(i,2) = sprintf("%.1f",VegC_g(i))
forrest@0:   text4(i,3) = sprintf("%.1f",LiCwC_g(i))
forrest@0:   text4(i,4) = sprintf("%.1f",SoilC_g(i))
forrest@0:   text4(i,5) = sprintf("%.1f",NPP_ratio(i))
forrest@0:   text4(i,6) = sprintf("%.1f",NEE2_g(i))
forrest@0:  end do
forrest@0:   text4(nrow-1,0) = "-"
forrest@0:   text4(nrow-1,1) = "-"
forrest@0:   text4(nrow-1,2) = "-"
forrest@0: 
forrest@0: ;-------------------------------------------------------
forrest@0: ; create html table1
forrest@0: 
forrest@0:   header_text = "<H1>NEE and Carbon Stocks and Fluxes:  Model "+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:        ,"   <th bgcolor=DDDDDD >"+col_head(6)+"</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: ;----------------------------
forrest@0: ;row of table
forrest@0: 
forrest@0:   do n = 0,nrow-2
forrest@0:      set_line(lines,nline,row_header)
forrest@0: 
forrest@0:      txt0  = row_head(n)
forrest@0:      txt1  = text4(n,0)
forrest@0:      txt2  = text4(n,1)
forrest@0:      txt3  = text4(n,2)
forrest@0:      txt4  = text4(n,3)
forrest@0:      txt5  = text4(n,4)
forrest@0:      txt6  = text4(n,5)
forrest@0:      txt7  = text4(n,6)
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:      set_line(lines,nline,"<th>"+txt7+"</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_carbon_sink1.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 (col_head)
forrest@0:   delete (row_head)
forrest@0:   delete (text4)
forrest@0: 
forrest@0: end
forrest@0: