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

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

;components

 component = (/"Leaf","Wood","Fine_Root","Litter","Coarse_Woody_Debris","Soil"/)

 n_comp = dimsizes(component)

 field_pool = (/"LEAFC","WOODC","FROOTC","LITTERC","CWDC","SOILC"/)

 field_flux = (/"LEAFC_ALLOC","WOODC_ALLOC","FROOTC_ALLOC","LITTERC_LOSS","CWDC_LOSS","SOILC_HR"/)

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

; model name and grid       

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

 model_grid = "T42"

 model_name = "i01.06cn"

;model_name = "i01.06casa"

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

; read data: ob       

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

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

; 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

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

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

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

; read data: model, each component

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

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

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

 do k = 0,n_comp-1

 pool  = fm->$field_pool(k)$

 flux  = fm->$field_flux(k)$

;Units for these variables are:

;pool: g C/m^2

;flux: g C/m^2/s

 nsec_per_year = 60*60*24*365

 flux = flux *  nsec_per_year 

; casa only

; all the plant pools (leaf, wood, and fine root) and

; coarse woody debris (cwd) and litter pools for

; CASA need to be divided by 1200.  The soil flux

; and turnover time are fine and do not need to be adjusted.

; if (k .ne. n_comp-1) then

;     flux = flux/1200.

; end if    

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

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

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

; using observed biome class  

; nclass      = nclass_ob

; using model biome class

  nclass      = nclass_mod

  nclassn     = nclass + 1

  range       = fspan(0,nclassn-1,nclassn)

; print (range)

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

; particular range, and then take an average.

  nr           = dimsizes(range)

  nx           = nr-1

  xvalues      = new((/2,nx/),float)

  xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.

  dx           = xvalues(0,1) - xvalues(0,0)       ; range width

  dx4          = dx/4                              ; 1/4 of the range

  xvalues(1,:) = xvalues(0,:) - dx/5.

;==============================

; put data into bins

;==============================

; using observed biome class

; base_1D  = ndtooned(classob)

; using model biome class

  base_1D  = ndtooned(classmod)

  data1_1D = ndtooned(pool)

  data2_1D = ndtooned(flux)

; output

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

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

  do nd=0,1

;   See if we are doing data1 (nd=0) or data2 (nd=1).

    base = base_1D

    if(nd.eq.0) then

      data = data1_1D

    else

      data = data2_1D

    end if

; Loop through each range, using base.

    do i=0,nr-2

      if (i.ne.(nr-2)) 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

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

        count(nd,i)   = dimsizes(idx)

      else

        yvalues(nd,i) = yvalues@_FillValue

        count(nd,i)   = 0

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

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

         count(nd,i)   = 0

      end if

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

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

; Clean up for next time in loop.

      delete(idx)

    end do

    delete(data)

  end do

;============================

;compute turnover time

;============================

 u       = yvalues(0,:)

 v       = yvalues(1,:)

 u_count = count(0,:)

 v_count = count(1,:)

;print (dimsizes(u))

;print (dimsizes(v))

 good = ind(.not.ismissing(u) .and. .not.ismissing(v))

;print (good)

 uu       = u(good)

 vv       = v(good)

;print (dimsizes(uu))

;print (dimsizes(vv))

 uu_count = u_count(good)

 vv_count = v_count(good)  

 n_biome = dimsizes(uu)

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

 t_biome = uu/vv

;t_biome_avg = avg(t_biome)

 t_biome_avg = sum(uu*uu_count)/sum(vv*vv_count)

;print (t_biome)

;print (t_biome_avg)

;===========================

; for html table - biome

;===========================

  output_html = "table_"+component(k)+".html"

; column (not including header column)

  col_head = (/component(k)+" Flux",component(k)+" Pool",component(k)+" Turnover Time"/)

  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. 

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

 do i=0,nrow-2

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

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

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

 end do

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

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

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

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

; html table

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

  header_text = "<H1>"+component(k)+" Turnover Time:  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 Class</th>" \

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

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

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

     txt1  = row_head(n)

     txt2  = text4(n,0)

     txt3  = text4(n,1)

     txt4  = text4(n,2)

     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,row_footer)

  end do

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

  set_line(lines,nline,table_footer)

  set_line(lines,nline,footer) 

; Now write to an HTML file.

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

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

    asciiwrite(output_html,lines(idx))

  else

   print ("error?")

  end if

  delete (idx)

end do

end