clamp: all/08.turnover.ncl@0c6405ab2ff4

     1 ;********************************************************

     2 ; hardwire: flux = flux/1200. (for casa only)

     3 ;**************************************************************

     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

     5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

     7 ;**************************************************************

     8 procedure set_line(lines:string,nline:integer,newlines:string)

     9 begin

    10 ; add line to ascci/html file

    12   nnewlines = dimsizes(newlines)

    13   if(nline+nnewlines-1.ge.dimsizes(lines))

    14     print("set_line: bad index, not setting anything.")

    15     return

    16   end if

    17   lines(nline:nline+nnewlines-1) = newlines

    18 ;  print ("lines = " + lines(nline:nline+nnewlines-1))

    19   nline = nline + nnewlines

    20   return

    21 end

    22 ;**************************************************************

    23 ; Main code.

    24 begin

    26  plot_type     = "ps"

    27  plot_type_new = "png"

    29 ;------------------------------------------------------

    30 ; edit table.html of current model for movel1_vs_model2

    32  if (isvar("compare")) then

    33     html_name2 = compare+"/table.html"

    34     html_new2  = html_name2 +".new"

    35  end if

    37 ;------------------------------------------------------

    38 ; edit table.html for current model

    40  html_name = model_name+"/table.html"

    41  html_new  = html_name +".new"

    43 ;---------------------------------------------------------------

    44 ;components

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

    47  n_comp = dimsizes(component)

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

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

    52 ;--------------------------------------------------

    53 ; get landfrac data

    55  film_l   = "lnd_"+ model_grid +".nc"

    56  fm_l     = addfile (dirs+film_l,"r")

    57  landfrac = fm_l->landfrac

    59  delete (fm_l)

    60 ;---------------------------------------------------

    61 ; read biome data: model

    63   biome_name_mod = "Model PFT Class"

    65   film_c   = "class_pft_"+ model_grid +".nc"

    66   fm_c     = addfile (dirs+film_c,"r")

    67   classmod = fm_c->CLASS_PFT

    69   delete (fm_c)

    71 ; model data has 17 land-type classes

    72   nclass_mod = 17

    74 ;********************************************************************

    75 ; use land-type class to bin the data in equally spaced ranges

    76 ;********************************************************************

    78 ; using model biome class

    79   nclass = nclass_mod

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

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

    84 ; particular range, and then take an average.

    86   nx = dimsizes(range) - 1

    88 ; for 2 data: pool and flux

    89   data_n = 2

    91 ; using model biome class

    93   base = ndtooned(classmod)

    95 ; output

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

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

   100 ;--------------------------------------------------

   101 ; read model data, each component:

   103  fm = addfile (dirm+film4,"r")

   105  do k = 0,n_comp-1

   107     pool  = fm->$field_pool(k)$

   108     flux  = fm->$field_flux(k)$

   110 ;   Units for these variables are:

   111 ;   pool: g C/m^2

   112 ;   flux: g C/m^2/s

   114     nsec_per_year = 60*60*24*365

   116     flux = flux *  nsec_per_year

   118     unit_p = "gC/m2"

   119     unit_f = "gC/m2/year"

   120     unit_t = "year"

   122 ;#############################################################

   123 ;   casa only

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

   125 ;   coarse woody debris (cwd) and litter pools for

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

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

   129     if (BGC .eq. "casa") then

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

   131           flux = flux/1200.

   132        end if

   133     end if

   134 ;##############################################################

   136 ;   take into account landfrac

   138     pool = pool * conform(pool,landfrac,(/1,2/))

   139     flux = flux * conform(flux,landfrac,(/1,2/))

   141 ; Loop through each range, using base

   143   do i=0,nx-1

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

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

   147      else

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

   149      end if

   151 ;    loop through each dataset

   153      do n = 0,data_n-1

   155         if (n .eq. 0) then

   156            data = ndtooned(pool)

   157         end if

   159         if (n .eq. 1) then

   160            data = ndtooned(flux)

   161         end if

   163 ;       Calculate average

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

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

   167            count(n,i)   = dimsizes(idx)

   168         else

   169            yvalues(n,i) = yvalues@_FillValue

   170            count(n,i)   = 0

   171         end if

   173 ;#############################################################

   174 ; using model biome class:

   175 ;

   176 ;     set the following 4 classes to _FillValue:

   177 ;     (3)Needleleaf Deciduous Boreal Tree,

   178 ;     (8)Broadleaf Deciduous Boreal Tree,

   179 ;     (9)Broadleaf Evergreen Shrub,

   180 ;     (16)Wheat

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

   183          yvalues(n,i) = yvalues@_FillValue

   184          count(n,i)   = 0

   185       end if

   186 ;#############################################################

   188       delete(data)

   189     end do                 ; n-loop

   191     delete(idx)

   192   end do                   ; i-loop

   194   delete (pool)

   195   delete (flux)

   197 ;============================

   198 ;compute turnover time

   199 ;============================

   201  u       = yvalues(0,:)

   202  v       = yvalues(1,:)

   203  u_count = count(0,:)

   204  v_count = count(1,:)

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

   208  uu       = u(good)

   209  vv       = v(good)

   210  uu_count = u_count(good)

   211  vv_count = v_count(good)

   213  n_biome = dimsizes(uu)

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

   216  t_biome = uu/vv

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

   220 ;===========================

   221 ; for html table - biome

   222 ;===========================

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

   226 ; column (not including header column)

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

   230   ncol = dimsizes(col_head)

   232 ; row (not including header row)

   234 ;----------------------------------------------------

   235 ; using model biome class:

   236   row_head  = (/"Not Vegetated" \

   237                ,"Needleleaf Evergreen Temperate Tree" \

   238                ,"Needleleaf Evergreen Boreal Tree" \

   239 ;              ,"Needleleaf Deciduous Boreal Tree" \

   240                ,"Broadleaf Evergreen Tropical Tree" \

   241                ,"Broadleaf Evergreen Temperate Tree" \

   242                ,"Broadleaf Deciduous Tropical Tree" \

   243                ,"Broadleaf Deciduous Temperate Tree" \

   244 ;              ,"Broadleaf Deciduous Boreal Tree" \

   245 ;              ,"Broadleaf Evergreen Shrub" \

   246                ,"Broadleaf Deciduous Temperate Shrub" \

   247                ,"Broadleaf Deciduous Boreal Shrub" \

   248                ,"C3 Arctic Grass" \

   249                ,"C3 Non-Arctic Grass" \

   250                ,"C4 Grass" \

   251                ,"Corn" \

   252 ;              ,"Wheat" \

   253                ,"All Biome" \

   254                /)

   255   nrow = dimsizes(row_head)

   257 ; arrays to be passed to table.

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

   260  do i=0,nrow-2

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

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

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

   264  end do

   265   text(nrow-1,0) = "-"

   266   text(nrow-1,1) = "-"

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

   269 ;**************************************************

   270 ; html table

   271 ;**************************************************

   273   header_text = "<H1>"+component(k)+" Turnover Time:  Model "+model_name+"</H1>"

   275   header = (/"<HTML>" \

   276             ,"<HEAD>" \

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

   278             ,"</HEAD>" \

   279             ,header_text \

   280             /)

   281   footer = "</HTML>"

   283   table_header = (/ \

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

   285        ,"<tr>" \

   286        ,"   <th bgcolor=DDDDDD >Biome Class</th>" \

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

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

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

   290        ,"</tr>" \

   291        /)

   292   table_footer = "</table>"

   293   row_header = "<tr>"

   294   row_footer = "</tr>"

   296   lines = new(50000,string)

   297   nline = 0

   299   set_line(lines,nline,header)

   300   set_line(lines,nline,table_header)

   301 ;-----------------------------------------------

   302 ; row of table

   304   do n = 0,nrow-1

   305      set_line(lines,nline,row_header)

   307      txt1  = row_head(n)

   308      txt2  = text(n,0)

   309      txt3  = text(n,1)

   310      txt4  = text(n,2)

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

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

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

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

   317      set_line(lines,nline,row_footer)

   318   end do

   319 ;-----------------------------------------------

   320   set_line(lines,nline,table_footer)

   321   set_line(lines,nline,footer)

   323 ; Now write to an HTML file

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

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

   327     asciiwrite(output_html,lines(idx))

   328   else

   329    print ("error?")

   330   end if

   332   delete (idx)

   334   delete (good)

   335   delete (t_biome)

   336   delete (text)

   338  end do          ; k-loop

   340  delete (fm)

   342 ;***************************************************************************

   343 ; output plot and html

   344 ;***************************************************************************

   345   output_dir = model_name+"/turnover"

   347   system("mv *.html " + output_dir)

   348 ;******************************

   350 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
permissions	-rw-r--r--