clamp: beta/06.biome_ob.ncl@0c6405ab2ff4

     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)

    99 ;============================

   100 ;compute beta_4

   101 ;============================

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

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

   107  beta_4_avg = avg(beta_4)

   109 ;print (beta_4)

   110 ;print (beta_4_avg)

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

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

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

   117  print (M_beta)

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

   120 ; for html table - station

   121 ;=========================

   123   output_html = "table_station.html"

   125 ; column (not including header column)

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

   129   ncol = dimsizes(col_head)

   131 ; row (not including header row)

   132   row_head = (/"DukeFACE" \

   133             ,"AspenFACE" \

   134             ,"ORNL-FACE" \

   135             ,"POP-EUROFACE" \

   136             ,"All Station" \

   137             /)

   138   nrow = dimsizes(row_head)

   140 ; arrays to be passed to table.

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

   143  do i=0,nrow-2

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

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

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

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

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

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

   150   text4(i,6) = sprintf("%.2f",beta_4(i))

   151   text4(i,7) = "-"

   152  end do

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

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

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

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

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

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

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

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

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

   163 ; html table

   164 ;-----------

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

   168   header = (/"<HTML>" \

   169             ,"<HEAD>" \

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

   171             ,"</HEAD>" \

   172             ,header_text \

   173             /)

   174   footer = "</HTML>"

   176   table_header = (/ \

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

   178        ,"<tr>" \

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

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

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

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

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

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

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

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

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

   188        ,"</tr>" \

   189        /)

   190   table_footer = "</table>"

   191   row_header = "<tr>"

   192   row_footer = "</tr>"

   194   lines = new(50000,string)

   195   nline = 0

   197   set_line(lines,nline,header)

   198   set_line(lines,nline,table_header)

   199 ;-----------------------------------------------

   200 ;row of table

   202   do n = 0,nrow-1

   203      set_line(lines,nline,row_header)

   205      txt1  = row_head(n)

   206      txt2  = text4(n,0)

   207      txt3  = text4(n,1)

   208      txt4  = text4(n,2)

   209      txt5  = text4(n,3)

   210      txt6  = text4(n,4)

   211      txt7  = text4(n,5)

   212      txt8  = text4(n,6)

   213      txt9  = text4(n,7)

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

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

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

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

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

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

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

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

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

   225      set_line(lines,nline,row_footer)

   226   end do

   227 ;-----------------------------------------------

   228   set_line(lines,nline,table_footer)

   229   set_line(lines,nline,footer)

   231 ; Now write to an HTML file.

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

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

   234     asciiwrite(output_html,lines(idx))

   235   else

   236    print ("error?")

   237   end if

   239   delete (col_head)

   240   delete (row_head)

   241   delete (text4)

   242   delete (table_header)

   243   delete (idx)

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

   246 ; read data: observed-2

   247 ;************************************************

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

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

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

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

   256   classob = tofloat(fo->LAND_CLASS)

   258 ; observed data has 20 land-type classes

   259   nclass = 20

   261 ;*******************************************************************

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

   263 ;********************************************************************

   265   nclassn     = nclass + 1

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

   267 ; print (range)

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

   270 ; particular range, and then take an average.

   272   nr           = dimsizes(range)

   273   nx           = nr-1

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

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

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

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

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

   279 ; get data

   281   base_1D  = ndtooned(classob)

   282   data1_1D = ndtooned(npp_i)

   283   data2_1D = ndtooned(npp_f)

   285 ; output

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

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

   290   do nd=0,1

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

   294     base = base_1D

   296     if(nd.eq.0) then

   297       data = data1_1D

   298     else

   299       data = data2_1D

   300     end if

   302 ; Loop through each range, using base.

   304     do i=0,nr-2

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

   306 ;        print("")

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

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

   309       else

   310 ;        print("")

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

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

   313       end if

   315 ;     Calculate average

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

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

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

   320       else

   321         yvalues(nd,i) = yvalues@_FillValue

   322         count(nd,i)   = 0

   323       end if

   325 ;#############################################################

   326 ; set the following 4 classes to _FillValue:

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

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

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

   331          yvalues(nd,i) = yvalues@_FillValue

   332          count(nd,i)   = 0

   333       end if

   334 ;#############################################################

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

   338 ; Clean up for next time in loop.

   340       delete(idx)

   341     end do

   343     delete(data)

   344   end do

   346 ;============================

   347 ;compute beta

   348 ;============================

   350  u       = yvalues(0,:)

   351  v       = yvalues(1,:)

   352  u_count = count(0,:)

   353  v_count = count(1,:)

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

   357  uu       = u(good)

   358  vv       = v(good)

   359  uu_count = u_count(good)

   360  vv_count = v_count(good)

   362  n_biome = dimsizes(uu)

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

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

   367 ;beta_biome_avg = avg(beta_biome)

   368  beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)

   370  print (beta_biome_avg)

   372 ;===========================

   373 ; for html table - biome

   374 ;===========================

   376   output_html = "table_biome.html"

   378 ; column (not including header column)

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

   382   ncol = dimsizes(col_head)

   384 ; row (not including header row)

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

   386 ; function "good" removed the last 2 classes

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

   389   row_head  = (/"Evergreen Needleleaf Forests" \

   390                ,"Evergreen Broadleaf Forests" \

   391                ,"Deciduous Needleleaf Forest" \

   392                ,"Deciduous Broadleaf Forests" \

   393                ,"Mixed Forests" \

   394                ,"Closed Bushlands" \

   395                ,"Open Bushlands" \

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

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

   398                ,"Grasslands" \

   399                ,"Permanent Wetlands" \

   400                ,"Croplands" \

   401                ,"Cropland/Natural Vegetation Mosaic" \

   402                ,"Barren or Sparsely Vegetated" \

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

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

   405                ,"All Biome" \

   406                /)

   407   nrow = dimsizes(row_head)

   409 ; arrays to be passed to table.

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

   412  do i=0,nrow-2

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

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

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

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

   417   text4(i,4) = sprintf("%.2f",beta_biome(i))

   418  end do

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

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

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

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

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

   425 ;**************************************************

   426 ; html table

   427 ;**************************************************

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

   431   header = (/"<HTML>" \

   432             ,"<HEAD>" \

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

   434             ,"</HEAD>" \

   435             ,header_text \

   436             /)

   437   footer = "</HTML>"

   439   table_header = (/ \

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

   441        ,"<tr>" \

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

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

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

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

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

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

   448        ,"</tr>" \

   449        /)

   450   table_footer = "</table>"

   451   row_header = "<tr>"

   452   row_footer = "</tr>"

   454   lines = new(50000,string)

   455   nline = 0

   457   set_line(lines,nline,header)

   458   set_line(lines,nline,table_header)

   459 ;-----------------------------------------------

   460 ;row of table

   462   do n = 0,nrow-1

   463      set_line(lines,nline,row_header)

   465      txt1  = row_head(n)

   466      txt2  = text4(n,0)

   467      txt3  = text4(n,1)

   468      txt4  = text4(n,2)

   469      txt5  = text4(n,3)

   470      txt6  = text4(n,4)

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

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

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

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

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

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

   479      set_line(lines,nline,row_footer)

   480   end do

   481 ;-----------------------------------------------

   482   set_line(lines,nline,table_footer)

   483   set_line(lines,nline,footer)

   485 ; Now write to an HTML file.

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

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

   488     asciiwrite(output_html,lines(idx))

   489   else

   490    print ("error?")

   491   end if

   493 end

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