forrest@0: ;********************************************************  
forrest@0: ; hardwire: flux = flux/1200. (for casa only)
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: ;components
forrest@0: 
forrest@0:  component = (/"Leaf","Wood","Fine_Root","Litter","Coarse_Woody_Debris","Soil"/)
forrest@0:  n_comp = dimsizes(component)
forrest@0: 
forrest@0:  field_pool = (/"LEAFC","WOODC","FROOTC","LITTERC","CWDC","SOILC"/)
forrest@0:  field_flux = (/"LEAFC_ALLOC","WOODC_ALLOC","FROOTC_ALLOC","LITTERC_LOSS","CWDC_LOSS","SOILC_HR"/)
forrest@0: 
forrest@0: ;--------------------------------------------------
forrest@0: ; get landfrac data
forrest@0:  
forrest@0:  film_l   = "lnd_"+ model_grid +".nc"
forrest@0:  fm_l     = addfile (dirs+film_l,"r")  
forrest@0:  landfrac = fm_l->landfrac
forrest@0: 
forrest@0:  delete (fm_l)
forrest@0: ;---------------------------------------------------
forrest@0: ; read 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: ; use land-type class to bin 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: ; 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: ; for 2 data: pool and flux
forrest@0:   data_n = 2
forrest@0: 
forrest@0: ; using model biome class
forrest@0: 
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: ;--------------------------------------------------
forrest@0: ; read model data, each component:      
forrest@0: 
forrest@0:  fm = addfile (dirm+film4,"r")
forrest@0: 
forrest@0:  do k = 0,n_comp-1
forrest@0: 
forrest@0:     pool  = fm->$field_pool(k)$
forrest@0:     flux  = fm->$field_flux(k)$
forrest@0: 
forrest@0: ;   Units for these variables are:
forrest@0: ;   pool: g C/m^2
forrest@0: ;   flux: g C/m^2/s
forrest@0: 
forrest@0:     nsec_per_year = 60*60*24*365
forrest@0:   
forrest@0:     flux = flux *  nsec_per_year 
forrest@0: 
forrest@0:     unit_p = "gC/m2"
forrest@0:     unit_f = "gC/m2/year"
forrest@0:     unit_t = "year"
forrest@0: 
forrest@0: ;#############################################################
forrest@0: ;   casa only
forrest@0: ;   all the plant pools (leaf, wood, and fine root) and
forrest@0: ;   coarse woody debris (cwd) and litter pools for
forrest@0: ;   CASA need to be divided by 1200.  The soil flux
forrest@0: ;   and turnover time are fine and do not need to be adjusted.
forrest@0: 
forrest@0:     if (BGC .eq. "casa") then   
forrest@0:        if (k .ne. n_comp-1) then
forrest@0:           flux = flux/1200.
forrest@0:        end if    
forrest@0:     end if
forrest@0: ;##############################################################
forrest@0: 
forrest@0: ;   take into account landfrac
forrest@0: 
forrest@0:     pool = pool * conform(pool,landfrac,(/1,2/))
forrest@0:     flux = flux * conform(flux,landfrac,(/1,2/))
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: ;    loop through each dataset
forrest@0:  
forrest@0:      do n = 0,data_n-1
forrest@0: 
forrest@0:         if (n .eq. 0) then
forrest@0:            data = ndtooned(pool)
forrest@0:         end if
forrest@0: 
forrest@0:         if (n .eq. 1) then
forrest@0:            data = ndtooned(flux)
forrest@0:         end if
forrest@0: 
forrest@0: ;       Calculate average 
forrest@0: 
forrest@0:         if (.not.any(ismissing(idx))) then
forrest@0:            yvalues(n,i) = avg(data(idx))
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: ;
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:       delete(data)
forrest@0:     end do                 ; n-loop
forrest@0: 
forrest@0:     delete(idx)
forrest@0:   end do                   ; i-loop
forrest@0: 
forrest@0:   delete (pool)
forrest@0:   delete (flux)
forrest@0: 
forrest@0: ;============================
forrest@0: ;compute turnover time
forrest@0: ;============================
forrest@0: 
forrest@0:  u       = yvalues(0,:)
forrest@0:  v       = yvalues(1,:)
forrest@0:  u_count = count(0,:)
forrest@0:  v_count = count(1,:)
forrest@0: 
forrest@0:  good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0: 
forrest@0:  uu       = u(good)
forrest@0:  vv       = v(good)
forrest@0:  uu_count = u_count(good)
forrest@0:  vv_count = v_count(good)  
forrest@0: 
forrest@0:  n_biome = dimsizes(uu)
forrest@0:  t_biome = new((/n_biome/),float)
forrest@0: 
forrest@0:  t_biome = uu/vv
forrest@0: 
forrest@0:  t_biome_avg = sum(uu*uu_count)/sum(vv*vv_count)
forrest@0: 
forrest@0: ;===========================
forrest@0: ; for html table - biome
forrest@0: ;===========================
forrest@0: 
forrest@0:   output_html = "table_"+component(k)+".html"
forrest@0: 
forrest@0: ; column (not including header column)
forrest@0: 
forrest@0:   col_head = (/component(k)+" Flux",component(k)+" Pool",component(k)+" Turnover Time"/)
forrest@0: 
forrest@0:   ncol = dimsizes(col_head)
forrest@0: 
forrest@0: ; row (not including header row)                   
forrest@0: 
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:   text = new ((/nrow, ncol/),string )
forrest@0:  
forrest@0:  do i=0,nrow-2
forrest@0:   text(i,0) = sprintf("%.1f",vv(i))
forrest@0:   text(i,1) = sprintf("%.1f",uu(i))
forrest@0:   text(i,2) = sprintf("%.2f",t_biome(i))
forrest@0:  end do
forrest@0:   text(nrow-1,0) = "-"
forrest@0:   text(nrow-1,1) = "-"
forrest@0:   text(nrow-1,2) = sprintf("%.2f",t_biome_avg)
forrest@0: 
forrest@0: ;**************************************************
forrest@0: ; html table
forrest@0: ;**************************************************
forrest@0: 
forrest@0:   header_text = "<H1>"+component(k)+" Turnover Time:  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 Class</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(0)+"<br>("+unit_f+")</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(1)+"<br>("+unit_p+")</th>" \
forrest@0:        ,"   <th bgcolor=DDDDDD >"+col_head(2)+"<br>("+unit_t+")</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:      txt1  = row_head(n)
forrest@0:      txt2  = text(n,0)
forrest@0:      txt3  = text(n,1)
forrest@0:      txt4  = text(n,2)
forrest@0: 
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: 
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:   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:   delete (good)
forrest@0:   delete (t_biome)
forrest@0:   delete (text)
forrest@0: 
forrest@0:  end do          ; k-loop
forrest@0: 
forrest@0:  delete (fm)
forrest@0: 
forrest@0: ;***************************************************************************
forrest@0: ; output plot and html
forrest@0: ;***************************************************************************
forrest@0:   output_dir = model_name+"/turnover"
forrest@0: 
forrest@0:   system("mv *.html " + output_dir) 
forrest@0: ;******************************
forrest@0: 
forrest@0: end
forrest@0: