clamp: beta/05.biome.ncl@4be95183fbcd

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

     2 ; required command line input parameters:

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

4 ;

     5 ; histogram normalized by rain and compute correleration

     6 ;**************************************************************

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

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

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

    10 ;**************************************************************

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

    12 begin

    13 ; add line to ascci/html file

    15   nnewlines = dimsizes(newlines)

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

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

    18     return

    19   end if

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

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

    22   nline = nline + nnewlines

    23   return

    24 end

    25 ;**************************************************************

    26 ; Main code.

    27 begin

    29  plot_type     = "ps"

    30  plot_type_new = "png"

    32 ;************************************************

    33 ; read data: model

    34 ;************************************************

    35  co2_i = 283.1878

    36  co2_f = 364.1252

    38  model_grid = "T42"

    40  model_name_i = "i01.07cn"

    41  model_name_f = "i01.10cn"

    43 ;model_name_i = "i01.07casa"

    44 ;model_name_f = "i01.10casa"

    46  model_name = model_name_f

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

    49  film_i = model_name_i + "_1990-2004_ANN_climo.nc"

    50  film_f = model_name_f + "_1990-2004_ANN_climo.nc"

    52  fm_i   = addfile (dirm+film_i,"r")

    53  fm_f   = addfile (dirm+film_f,"r")

    55  xm     = fm_f->lon

    56  ym     = fm_f->lat

    58  npp_i  = fm_i->NPP

    59  npp_f  = fm_f->NPP

    61 ;Units for these variables are:

    62 ;npp_i: g C/m^2/s

    64  nsec_per_year = 60*60*24*365

    66  npp_i = npp_i *  nsec_per_year

    67  npp_f = npp_f *  nsec_per_year

    69 ;===================================================

    70 ; read data: observed -station

    71 ;===================================================

    73  station = (/"DukeFACE" \

    74             ,"AspenFACE" \

    75             ,"ORNL-FACE" \

    76             ,"POP-EUROFACE" \

    77             /)

    79  lat_ob = (/ 35.58,  45.40,  35.54, 42.22/)

    80  lon_ob = (/-79.05, -89.37, -84.20, 11.48/)

    81  lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)

    82 ;print (lon_ob)

    84  n_sta  = dimsizes(station)

    85  beta_4_ob = new((/n_sta/),float)

    86  beta_4_ob = 0.60

    88 ;===================================================

    89 ; get model data at station

    90 ;===================================================

    92  npp_i_4  =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)

    94  npp_f_4  =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)

    96 ;print (npp_i_4)

    97 ;print (npp_f_4)

    98 ;============================

    99 ;compute beta_4

   100 ;============================

   101  beta_4 = new((/n_sta/),float)

   103  beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)

   105  beta_4_avg = avg(beta_4)

   107 ;print (beta_4)

   108 ;print (beta_4_avg)

   110 ;M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 3.

   112  bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))

   113  M_beta  = (1. - (bias/n_sta))*3.

   115  print (M_beta)

   117 ;=========================

   118 ; for html table - station

   119 ;=========================

   121   output_html = "table_station.html"

   123 ; column (not including header column)

   125   col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)

   127   ncol = dimsizes(col_head)

   129 ; row (not including header row)

   130   row_head = (/"DukeFACE" \

   131             ,"AspenFACE" \

   132             ,"ORNL-FACE" \

   133             ,"POP-EUROFACE" \

   134             ,"All Station" \

   135             /)

   136   nrow = dimsizes(row_head)

   138 ; arrays to be passed to table.

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

   141  do i=0,nrow-2

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

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

   144   text4(i,2) = sprintf("%.1f",co2_i)

   145   text4(i,3) = sprintf("%.1f",co2_f)

   146   text4(i,4) = sprintf("%.1f",npp_i_4(0,i))

   147   text4(i,5) = sprintf("%.1f",npp_f_4(0,i))

   148   text4(i,6) = sprintf("%.1f",beta_4(i))

   149   text4(i,7) = "-"

   150  end do

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

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

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

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

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

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

   157   text4(nrow-1,6) = sprintf("%.1f",beta_4_avg)

   158   text4(nrow-1,7) = sprintf("%.1f",avg(beta_4_ob))

   160 ;-----------

   161 ; html table

   162 ;-----------

   164   header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"

   166   header = (/"<HTML>" \

   167             ,"<HEAD>" \

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

   169             ,"</HEAD>" \

   170             ,header_text \

   171             /)

   172   footer = "</HTML>"

   174   table_header = (/ \

   175         "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \

   176        ,"<tr>" \

   177        ,"   <th bgcolor=DDDDDD >Station</th>" \

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

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

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

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

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

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

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

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

   186        ,"</tr>" \

   187        /)

   188   table_footer = "</table>"

   189   row_header = "<tr>"

   190   row_footer = "</tr>"

   192   lines = new(50000,string)

   193   nline = 0

   195   set_line(lines,nline,header)

   196   set_line(lines,nline,table_header)

   197 ;-----------------------------------------------

   198 ;row of table

   200   do n = 0,nrow-1

   201      set_line(lines,nline,row_header)

   203      txt1  = row_head(n)

   204      txt2  = text4(n,0)

   205      txt3  = text4(n,1)

   206      txt4  = text4(n,2)

   207      txt5  = text4(n,3)

   208      txt6  = text4(n,4)

   209      txt7  = text4(n,5)

   210      txt8  = text4(n,6)

   211      txt9  = text4(n,7)

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

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

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

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

   217      set_line(lines,nline,"<th>"+txt5+"</th>")

   218      set_line(lines,nline,"<th>"+txt6+"</th>")

   219      set_line(lines,nline,"<th>"+txt7+"</th>")

   220      set_line(lines,nline,"<th>"+txt8+"</th>")

   221      set_line(lines,nline,"<th>"+txt9+"</th>")

   223      set_line(lines,nline,row_footer)

   224   end do

   225 ;-----------------------------------------------

   226   set_line(lines,nline,table_footer)

   227   set_line(lines,nline,footer)

   229 ; Now write to an HTML file.

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

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

   232     asciiwrite(output_html,lines(idx))

   233   else

   234    print ("error?")

   235   end if

   237   delete (col_head)

   238   delete (row_head)

   239   delete (text4)

   240   delete (table_header)

   241   delete (idx)

   243 ;************************************************

   244 ; read data: observed-2

   245 ;************************************************

   247   ob_name = "MODIS MOD 15A2 2000-2005"

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

   250   filo  = "land_class_"+model_grid+".nc"

   252   fo = addfile(diro+filo,"r")

   254   classob = tofloat(fo->LAND_CLASS)

   256   nclass = 20

   258 ;*******************************************************************

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

   260 ;********************************************************************

   262   nclassn     = nclass + 1

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

   264 ; print (range)

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

   267 ; particular range, and then take an average.

   269   nr           = dimsizes(range)

   270   nx           = nr-1

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

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

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

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

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

   276 ; get data

   278   base_1D  = ndtooned(classob)

   279   data1_1D = ndtooned(npp_i)

   280   data2_1D = ndtooned(npp_f)

   282 ; output

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

   286   do nd=0,1

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

   290     base = base_1D

   292     if(nd.eq.0) then

   293       data = data1_1D

   294     else

   295       data = data2_1D

   296     end if

   298 ; Loop through each range, using base.

   300     do i=0,nr-2

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

   302 ;        print("")

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

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

   305       else

   306 ;        print("")

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

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

   309       end if

   311 ;     Calculate average

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

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

   315         count = dimsizes(idx)

   316       else

   317         yvalues(nd,i)    = yvalues@_FillValue

   318         count            = 0

   319       end if

   321 ;#############################################################

   322 ; set the following 4 classes to _FillValue:

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

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

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

   327          yvalues(nd,i)    = yvalues@_FillValue

   328          count            = 0

   329       end if

   330 ;#############################################################

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

   334 ; Clean up for next time in loop.

   336       delete(idx)

   337     end do

   339     delete(data)

   340   end do

   343 ;============================

   344 ;compute beta

   345 ;============================

   347  u = yvalues(0,:)

   348  v = yvalues(1,:)

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

   352  uu = u(good)

   353  vv = v(good)

   355  n_biome = dimsizes(uu)

   357  beta_biome = new((/n_biome/),float)

   359  beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)

   361  beta_biome_avg = avg(beta_biome)

   363 ;print (beta_biome_avg)

   365 ;===========================

   366 ; for html table - biome

   367 ;===========================

   369   output_html = "table_biome.html"

   371 ; column (not including header column)

   373   col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)

   375   ncol = dimsizes(col_head)

   377 ; row (not including header row)

   378 ; 3 classes removed: Water Bodies, Urban and Build-Up, Unclassified

   379 ; function "good" removed the last 2 classes

   380 ; text4(i,2) = sprintf("%.2f",uu(i+1)) remove the first class

   382   row_head  = (/"Evergreen Needleleaf Forests" \

   383                ,"Evergreen Broadleaf Forests" \

   384                ,"Deciduous Needleleaf Forest" \

   385                ,"Deciduous Broadleaf Forests" \

   386                ,"Mixed Forests" \

   387                ,"Closed Bushlands" \

   388                ,"Open Bushlands" \

   389                ,"Woody Savannas (S. Hem.)" \

   390                ,"Savannas (S. Hem.)" \

   391                ,"Grasslands" \

   392                ,"Permanent Wetlands" \

   393                ,"Croplands" \

   394                ,"Cropland/Natural Vegetation Mosaic" \

   395                ,"Barren or Sparsely Vegetated" \

   396                ,"Woody Savannas (N. Hem.)" \

   397                ,"Savannas (N. Hem.)" \

   398                ,"All Biome" \

   399                /)

   400   nrow = dimsizes(row_head)

   402 ; arrays to be passed to table.

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

   405  do i=0,nrow-2

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

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

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

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

   410   text4(i,4) = sprintf("%.1f",beta_biome(i))

   411  end do

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

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

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

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

   416   text4(nrow-1,4) = sprintf("%.1f",beta_biome_avg)

   418 ;**************************************************

   419 ; html table

   420 ;**************************************************

   422   header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"

   424   header = (/"<HTML>" \

   425             ,"<HEAD>" \

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

   427             ,"</HEAD>" \

   428             ,header_text \

   429             /)

   430   footer = "</HTML>"

   432   table_header = (/ \

   433         "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \

   434        ,"<tr>" \

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

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

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

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

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

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

   441        ,"</tr>" \

   442        /)

   443   table_footer = "</table>"

   444   row_header = "<tr>"

   445   row_footer = "</tr>"

   447   lines = new(50000,string)

   448   nline = 0

   450   set_line(lines,nline,header)

   451   set_line(lines,nline,table_header)

   452 ;-----------------------------------------------

   453 ;row of table

   455   do n = 0,nrow-1

   456      set_line(lines,nline,row_header)

   458      txt1  = row_head(n)

   459      txt2  = text4(n,0)

   460      txt3  = text4(n,1)

   461      txt4  = text4(n,2)

   462      txt5  = text4(n,3)

   463      txt6  = text4(n,4)

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

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

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

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

   469      set_line(lines,nline,"<th>"+txt5+"</th>")

   470      set_line(lines,nline,"<th>"+txt6+"</th>")

   472      set_line(lines,nline,row_footer)

   473   end do

   474 ;-----------------------------------------------

   475   set_line(lines,nline,table_footer)

   476   set_line(lines,nline,footer)

   478 ; Now write to an HTML file.

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

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

   481     asciiwrite(output_html,lines(idx))

   482   else

   483    print ("error?")

   484   end if

   486 end

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