clamp: beta/07.landfrac.ncl@4be95183fbcd

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

     2 ; take into account landfrac

     3 ; note: landfrac from lnd_T42.nc

     4 ;       <= lnd_diag_4.0 has correct landfrac

     5 ;          lnd_diag_3.1 has wrong   landfrac

6 ;

     7 ; using model biome

8 ;

     9 ; required command line input parameters:

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

    11 ;

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

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

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

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

    16 ;**************************************************************

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

    18 begin

    19 ; add line to ascci/html file

    21   nnewlines = dimsizes(newlines)

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

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

    24     return

    25   end if

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

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

    28   nline = nline + nnewlines

    29   return

    30 end

    31 ;**************************************************************

    32 ; Main code.

    33 begin

    35  plot_type     = "ps"

    36  plot_type_new = "png"

    38 ;************************************************

    39 ; read data: model

    40 ;************************************************

    41  co2_i = 283.1878

    42  co2_f = 364.1252

    44  model_grid = "T42"

    46 ;model_name_i = "i01.07cn"

    47 ;model_name_f = "i01.10cn"

    49  model_name_i = "i01.07casa"

    50  model_name_f = "i01.10casa"

    52  model_name = model_name_f

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

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

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

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

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

    61  xm     = fm_f->lon

    62  ym     = fm_f->lat

    64  npp_i  = fm_i->NPP

    65  npp_f  = fm_f->NPP

    67  delete (fm_i)

    68  delete (fm_f)

    70 ;Units for these variables are:

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

    73  nsec_per_year = 60*60*24*365

    75  npp_i = npp_i *  nsec_per_year

    76  npp_f = npp_f *  nsec_per_year

    78 ;--------------------------------------------------

    79 ;get landfrac data

    81  dirm= "/fis/cgd/cseg/people/jeff/surface_data/"

    82  film_l = "lnd_T42.nc"

    83  fm_l   = addfile (dirm+film_l,"r")

    85  landfrac = fm_l->landfrac

    87 ;npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:)

    88 ;npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:)

    90  npp_i = npp_i * conform(npp_i, landfrac, (/1,2/))

    91  npp_f = npp_f * conform(npp_f, landfrac, (/1,2/))

    93 ;===================================================

    94 ; read data: observed at stations

    95 ;===================================================

    97  station = (/"DukeFACE" \

    98             ,"AspenFACE" \

    99             ,"ORNL-FACE" \

   100             ,"POP-EUROFACE" \

   101             /)

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

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

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

   106 ;print (lon_ob)

   108  n_sta  = dimsizes(station)

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

   110  beta_4_ob = 0.60

   112 ;===================================================

   113 ; get model data at station

   114 ;===================================================

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

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

   120 ;print (npp_i_4)

   121 ;print (npp_f_4)

   123 ;============================

   124 ;compute beta_4

   125 ;============================

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

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

   131  beta_4_avg = avg(beta_4)

   133 ;print (beta_4)

   134 ;print (beta_4_avg)

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

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

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

   141  print (M_beta)

   143 ;=========================

   144 ; for html table - station

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

   147   output_html = "table_station.html"

   149 ; column (not including header column)

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

   153   ncol = dimsizes(col_head)

   155 ; row (not including header row)

   156   row_head = (/"DukeFACE" \

   157               ,"AspenFACE" \

   158               ,"ORNL-FACE" \

   159               ,"POP-EUROFACE" \

   160               ,"All Station" \

   161               /)

   162   nrow = dimsizes(row_head)

   164 ; arrays to be passed to table.

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

   167  do i=0,nrow-2

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

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

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

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

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

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

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

   175   text4(i,7) = "-"

   176  end do

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

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

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

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

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

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

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

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

   186 ;-----------

   187 ; html table

   188 ;-----------

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

   192   header = (/"<HTML>" \

   193             ,"<HEAD>" \

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

   195             ,"</HEAD>" \

   196             ,header_text \

   197             /)

   198   footer = "</HTML>"

   200   table_header = (/ \

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

   202        ,"<tr>" \

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

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

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

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

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

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

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

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

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

   212        ,"</tr>" \

   213        /)

   214   table_footer = "</table>"

   215   row_header = "<tr>"

   216   row_footer = "</tr>"

   218   lines = new(50000,string)

   219   nline = 0

   221   set_line(lines,nline,header)

   222   set_line(lines,nline,table_header)

   223 ;-----------------------------------------------

   224 ;row of table

   226   do n = 0,nrow-1

   227      set_line(lines,nline,row_header)

   229      txt1  = row_head(n)

   230      txt2  = text4(n,0)

   231      txt3  = text4(n,1)

   232      txt4  = text4(n,2)

   233      txt5  = text4(n,3)

   234      txt6  = text4(n,4)

   235      txt7  = text4(n,5)

   236      txt8  = text4(n,6)

   237      txt9  = text4(n,7)

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

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

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

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

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

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

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

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

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

   249      set_line(lines,nline,row_footer)

   250   end do

   251 ;-----------------------------------------------

   252   set_line(lines,nline,table_footer)

   253   set_line(lines,nline,footer)

   255 ; Now write to an HTML file.

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

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

   258     asciiwrite(output_html,lines(idx))

   259   else

   260    print ("error?")

   261   end if

   263   delete (col_head)

   264   delete (row_head)

   265   delete (text4)

   266   delete (table_header)

   267   delete (idx)

   269 ;------------------------------------------------

   270 ; read biome data: model

   272   biome_name_mod = "Model PFT Class"

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

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

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

   279   classmod = fm->CLASS_PFT

   281   delete (fm)

   283 ; model data has 17 land-type classes

   285   nclass_mod = 17

   287 ;------------------------------------------------

   288 ; read biome data: observed

   290   biome_name_ob = "MODIS LandCover"

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

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

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

   297   classob = tofloat(fo->LAND_CLASS)

   299   delete (fo)

   301 ; input observed data has 20 land-type classes

   303   nclass_ob = 20

   305 ;********************************************************************

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

   307 ;********************************************************************

   309 ; using observed biome class

   310 ; nclass      = nclass_ob

   311 ; using model biome class

   312   nclass      = nclass_mod

   314   nclassn     = nclass + 1

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

   316 ; print (range)

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

   319 ; particular range, and then take an average.

   321   nr           = dimsizes(range)

   322   nx           = nr-1

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

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

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

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

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

   329 ;==============================

   330 ; put data into bins

   331 ;==============================

   333 ; using observed biome class

   334 ; base_1D  = ndtooned(classob)

   335 ; using model biome class

   336   base_1D  = ndtooned(classmod)

   338   data1_1D = ndtooned(npp_i)

   339   data2_1D = ndtooned(npp_f)

   341 ; output

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

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

   346   do nd=0,1

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

   350     base = base_1D

   352     if(nd.eq.0) then

   353       data = data1_1D

   354     else

   355       data = data2_1D

   356     end if

   358 ; Loop through each range, using base.

   360     do i=0,nr-2

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

   362 ;        print("")

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

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

   365       else

   366 ;        print("")

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

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

   369       end if

   371 ;     Calculate average

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

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

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

   376       else

   377         yvalues(nd,i) = yvalues@_FillValue

   378         count(nd,i)   = 0

   379       end if

   381 ;#############################################################

   382 ;using observed biome class:

   383 ;

   384 ;     set the following 4 classes to _FillValue:

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

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

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

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

   390 ;        count(nd,i)   = 0

   391 ;     end if

   392 ;#############################################################

   394 ;#############################################################

   395 ;using model biome class:

   396 ;

   397 ;     set the following 4 classes to _FillValue:

   398 ;     (3)Needleleaf Deciduous Boreal Tree,

   399 ;     (8)Broadleaf Deciduous Boreal Tree,

   400 ;     (9)Broadleaf Evergreen Shrub,

   401 ;     (16)Wheat

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

   404          yvalues(nd,i) = yvalues@_FillValue

   405          count(nd,i)   = 0

   406       end if

   407 ;#############################################################

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

   411 ; Clean up for next time in loop.

   413       delete(idx)

   414     end do

   416     delete(data)

   417   end do

   419 ;============================

   420 ;compute beta

   421 ;============================

   423  u       = yvalues(0,:)

   424  v       = yvalues(1,:)

   425  u_count = count(0,:)

   426  v_count = count(1,:)

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

   430  uu       = u(good)

   431  vv       = v(good)

   432  uu_count = u_count(good)

   433  vv_count = v_count(good)

   435  n_biome = dimsizes(uu)

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

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

   440 ;beta_biome_avg = avg(beta_biome)

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

   443 ;print (beta_biome_avg)

   445 ;===========================

   446 ; for html table - biome

   447 ;===========================

   449   output_html = "table_biome.html"

   451 ; column (not including header column)

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

   455   ncol = dimsizes(col_head)

   457 ; row (not including header row)

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

   460 ; using observed biome class:

   461 ;

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

   463 ;              ,"Evergreen Broadleaf Forests" \

   464 ;              ,"Deciduous Needleleaf Forest" \

   465 ;              ,"Deciduous Broadleaf Forests" \

   466 ;              ,"Mixed Forests" \

   467 ;              ,"Closed Bushlands" \

   468 ;              ,"Open Bushlands" \

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

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

   471 ;              ,"Grasslands" \

   472 ;              ,"Permanent Wetlands" \

   473 ;              ,"Croplands" \

   474 ;              ,"Cropland/Natural Vegetation Mosaic" \

   475 ;              ,"Barren or Sparsely Vegetated" \

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

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

   478 ;              ,"All Biome" \

   479 ;              /)

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

   482 ; using model biome class:

   483 ;

   484   row_head  = (/"Not Vegetated" \

   485                ,"Needleleaf Evergreen Temperate Tree" \

   486                ,"Needleleaf Evergreen Boreal Tree" \

   487 ;              ,"Needleleaf Deciduous Boreal Tree" \

   488                ,"Broadleaf Evergreen Tropical Tree" \

   489                ,"Broadleaf Evergreen Temperate Tree" \

   490                ,"Broadleaf Deciduous Tropical Tree" \

   491                ,"Broadleaf Deciduous Temperate Tree" \

   492 ;              ,"Broadleaf Deciduous Boreal Tree" \

   493 ;              ,"Broadleaf Evergreen Shrub" \

   494                ,"Broadleaf Deciduous Temperate Shrub" \

   495                ,"Broadleaf Deciduous Boreal Shrub" \

   496                ,"C3 Arctic Grass" \

   497                ,"C3 Non-Arctic Grass" \

   498                ,"C4 Grass" \

   499                ,"Corn" \

   500 ;              ,"Wheat" \

   501                ,"All Biome" \

   502                /)

   504   nrow = dimsizes(row_head)

   506 ; arrays to be passed to table.

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

   509  do i=0,nrow-2

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

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

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

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

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

   515  end do

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

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

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

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

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

   522 ;**************************************************

   523 ; html table

   524 ;**************************************************

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

   528   header = (/"<HTML>" \

   529             ,"<HEAD>" \

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

   531             ,"</HEAD>" \

   532             ,header_text \

   533             /)

   534   footer = "</HTML>"

   536   table_header = (/ \

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

   538        ,"<tr>" \

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

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

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

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

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

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

   545        ,"</tr>" \

   546        /)

   547   table_footer = "</table>"

   548   row_header = "<tr>"

   549   row_footer = "</tr>"

   551   lines = new(50000,string)

   552   nline = 0

   554   set_line(lines,nline,header)

   555   set_line(lines,nline,table_header)

   556 ;-----------------------------------------------

   557 ;row of table

   559   do n = 0,nrow-1

   560      set_line(lines,nline,row_header)

   562      txt1  = row_head(n)

   563      txt2  = text4(n,0)

   564      txt3  = text4(n,1)

   565      txt4  = text4(n,2)

   566      txt5  = text4(n,3)

   567      txt6  = text4(n,4)

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

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

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

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

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

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

   576      set_line(lines,nline,row_footer)

   577   end do

   578 ;-----------------------------------------------

   579   set_line(lines,nline,table_footer)

   580   set_line(lines,nline,footer)

   582 ; Now write to an HTML file.

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

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

   585     asciiwrite(output_html,lines(idx))

   586   else

   587    print ("error?")

   588   end if

   590 end

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