clamp: turnover/03.biome_model.ncl@4be95183fbcd

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

     2 ;using model biome vlass

3 ;

     4 ; required command line input parameters:

     5 ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl

6 ;

     7 ; histogram normalized by rain and compute correleration

     8 ;**************************************************************

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

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

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

    12 ;**************************************************************

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

    14 begin

    15 ; add line to ascci/html file

    17   nnewlines = dimsizes(newlines)

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

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

    20     return

    21   end if

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

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

    24   nline = nline + nnewlines

    25   return

    26 end

    27 ;**************************************************************

    28 ; Main code.

    29 begin

    31  plot_type     = "ps"

    32  plot_type_new = "png"

    34 ;components

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

    37  n_comp = dimsizes(component)

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

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

    42 ;************************************************

    43 ; model name and grid

    44 ;************************************************

    46  model_grid = "T42"

    48  model_name = "i01.06cn"

    49 ;model_name = "i01.06casa"

    51 ;************************************************

    52 ; read data: ob

    53 ;************************************************

    54 ;------------------------------------------------

    55 ; read biome data: observed

    57 ; biome_name_ob = "MODIS LandCover"

    59 ; diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"

    60 ; filo = "land_class_"+model_grid+".nc"

    62 ; fo = addfile(diro+filo,"r")

    64 ; classob = tofloat(fo->LAND_CLASS)

    66 ; delete (fo)

    68 ; observed data has 20 land-type classes

    70 ; nclass_ob = 20

    72 ;---------------------------------------------------

    73 ; read biome data: model

    75   biome_name_mod = "Model PFT Class"

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

    78   film = "class_pft_"+model_grid+".nc"

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

    82   classmod = fm->CLASS_PFT

    84   delete (fm)

    86 ; model data has 17 land-type classes

    88   nclass_mod = 17

    90 ;---------------------------------------------------

    91 ; read data: model, each component

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

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

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

    98  do k = 0,n_comp-1

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

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

   103 ;Units for these variables are:

   104 ;pool: g C/m^2

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

   107  nsec_per_year = 60*60*24*365

   109  flux = flux *  nsec_per_year

   111 ; casa only

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

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

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

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

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

   118 ;     flux = flux/1200.

   119 ; end if

   121 ;*******************************************************************

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

   123 ;********************************************************************

   125 ; using observed biome class

   126 ; nclass      = nclass_ob

   127 ; using model biome class

   128   nclass      = nclass_mod

   130   nclassn     = nclass + 1

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

   132 ; print (range)

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

   135 ; particular range, and then take an average.

   137   nr           = dimsizes(range)

   138   nx           = nr-1

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

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

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

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

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

   145 ;==============================

   146 ; put data into bins

   147 ;==============================

   149 ; using observed biome class

   150 ; base_1D  = ndtooned(classob)

   151 ; using model biome class

   152   base_1D  = ndtooned(classmod)

   154   data1_1D = ndtooned(pool)

   155   data2_1D = ndtooned(flux)

   157 ; output

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

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

   162   do nd=0,1

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

   166     base = base_1D

   168     if(nd.eq.0) then

   169       data = data1_1D

   170     else

   171       data = data2_1D

   172     end if

   174 ; Loop through each range, using base.

   176     do i=0,nr-2

   177       if (i.ne.(nr-2)) then

   178 ;        print("")

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

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

   181       else

   182 ;        print("")

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

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

   185       end if

   187 ;     Calculate average

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

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

   191         count(nd,i)   = dimsizes(idx)

   192       else

   193         yvalues(nd,i) = yvalues@_FillValue

   194         count(nd,i)   = 0

   195       end if

   197 ;#############################################################

   198 ;using observed biome class:

   199 ;     set the following 4 classes to _FillValue:

   200 ;     Water Bodies(0), Urban and Build-Up(13),

   201 ;     Permenant Snow and Ice(15), Unclassified(17)

   203 ;     if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then

   204 ;        yvalues(nd,i) = yvalues@_FillValue

   205 ;        count(nd,i)   = 0

   206 ;     end if

   207 ;#############################################################

   209 ;#############################################################

   210 ;using model biome class:

   211 ;     set the following 4 classes to _FillValue:

   212 ;     (3)Needleleaf Deciduous Boreal Tree,

   213 ;     (8)Broadleaf Deciduous Boreal Tree,

   214 ;     (9)Broadleaf Evergreen Shrub,

   215 ;     (16)Wheat

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

   218          yvalues(nd,i) = yvalues@_FillValue

   219          count(nd,i)   = 0

   220       end if

   221 ;#############################################################

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

   225 ; Clean up for next time in loop.

   227       delete(idx)

   228     end do

   230     delete(data)

   231   end do

   233 ;============================

   234 ;compute turnover time

   235 ;============================

   237  u       = yvalues(0,:)

   238  v       = yvalues(1,:)

   239  u_count = count(0,:)

   240  v_count = count(1,:)

   242 ;print (dimsizes(u))

   243 ;print (dimsizes(v))

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

   247 ;print (good)

   249  uu       = u(good)

   250  vv       = v(good)

   252 ;print (dimsizes(uu))

   253 ;print (dimsizes(vv))

   255  uu_count = u_count(good)

   256  vv_count = v_count(good)

   258  n_biome = dimsizes(uu)

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

   261  t_biome = uu/vv

   263 ;t_biome_avg = avg(t_biome)

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

   266 ;print (t_biome)

   267 ;print (t_biome_avg)

   269 ;===========================

   270 ; for html table - biome

   271 ;===========================

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

   275 ; column (not including header column)

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

   279   ncol = dimsizes(col_head)

   281 ; row (not including header row)

   283 ;----------------------------------------------------

   284 ; using observed biome class:

   285 ; row_head  = (/"Evergreen Needleleaf Forests" \

   286 ;              ,"Evergreen Broadleaf Forests" \

   287 ;              ,"Deciduous Needleleaf Forest" \

   288 ;              ,"Deciduous Broadleaf Forests" \

   289 ;              ,"Mixed Forests" \

   290 ;              ,"Closed Bushlands" \

   291 ;              ,"Open Bushlands" \

   292 ;              ,"Woody Savannas (S. Hem.)" \

   293 ;              ,"Savannas (S. Hem.)" \

   294 ;              ,"Grasslands" \

   295 ;              ,"Permanent Wetlands" \

   296 ;              ,"Croplands" \

   297 ;              ,"Cropland/Natural Vegetation Mosaic" \

   298 ;              ,"Barren or Sparsely Vegetated" \

   299 ;              ,"Woody Savannas (N. Hem.)" \

   300 ;              ,"Savannas (N. Hem.)" \

   301 ;              ,"All Biome" \

   302 ;              /)

   304 ;----------------------------------------------------

   305 ; using model biome class:

   306   row_head  = (/"Not Vegetated" \

   307                ,"Needleleaf Evergreen Temperate Tree" \

   308                ,"Needleleaf Evergreen Boreal Tree" \

   309 ;              ,"Needleleaf Deciduous Boreal Tree" \

   310                ,"Broadleaf Evergreen Tropical Tree" \

   311                ,"Broadleaf Evergreen Temperate Tree" \

   312                ,"Broadleaf Deciduous Tropical Tree" \

   313                ,"Broadleaf Deciduous Temperate Tree" \

   314 ;              ,"Broadleaf Deciduous Boreal Tree" \

   315 ;              ,"Broadleaf Evergreen Shrub" \

   316                ,"Broadleaf Deciduous Temperate Shrub" \

   317                ,"Broadleaf Deciduous Boreal Shrub" \

   318                ,"C3 Arctic Grass" \

   319                ,"C3 Non-Arctic Grass" \

   320                ,"C4 Grass" \

   321                ,"Corn" \

   322 ;              ,"Wheat" \

   323                ,"All Biome" \

   324                /)

   325   nrow = dimsizes(row_head)

   327 ; arrays to be passed to table.

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

   330  do i=0,nrow-2

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

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

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

   334  end do

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

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

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

   339 ;**************************************************

   340 ; html table

   341 ;**************************************************

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

   345   header = (/"<HTML>" \

   346             ,"<HEAD>" \

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

   348             ,"</HEAD>" \

   349             ,header_text \

   350             /)

   351   footer = "</HTML>"

   353   table_header = (/ \

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

   355        ,"<tr>" \

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

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

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

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

   360        ,"</tr>" \

   361        /)

   362   table_footer = "</table>"

   363   row_header = "<tr>"

   364   row_footer = "</tr>"

   366   lines = new(50000,string)

   367   nline = 0

   369   set_line(lines,nline,header)

   370   set_line(lines,nline,table_header)

   371 ;-----------------------------------------------

   372 ;row of table

   374   do n = 0,nrow-1

   375      set_line(lines,nline,row_header)

   377      txt1  = row_head(n)

   378      txt2  = text4(n,0)

   379      txt3  = text4(n,1)

   380      txt4  = text4(n,2)

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

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

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

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

   387      set_line(lines,nline,row_footer)

   388   end do

   389 ;-----------------------------------------------

   390   set_line(lines,nline,table_footer)

   391   set_line(lines,nline,footer)

   393 ; Now write to an HTML file.

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

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

   396     asciiwrite(output_html,lines(idx))

   397   else

   398    print ("error?")

   399   end if

   401   delete (idx)

   402 end do

   403 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--