clamp: all/07.beta.ncl@0c6405ab2ff4

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

     2 ; hardwired: co2_i = 283.1878

     3 ;            co2_f = 364.1252

4 ;

     5 ;            beta_4_ob = 0.6

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

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

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

     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 ; edit table.html of current model for movel1_vs_model2

    35  if (isvar("compare")) then

    36     html_name2 = compare+"/table.html"

    37     html_new2  = html_name2 +".new"

    38  end if

    40 ;-------------------------------------------------------

    41 ; edit table.html for current model

    43  html_name = model_name+"/table.html"

    44  html_new  = html_name +".new"

    46 ;-------------------------------------------------------

    47 ; read model data

    49 ;###############################################################

    50 ; hardwired for model data

    52 ; these values correspond to the start and the end of model data

    53  co2_i = 283.1878

    54  co2_f = 364.1252

    56  film_i = film6

    57  film_f = film5

    58 ;###############################################################

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

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

    63  xm     = fm_f->lon

    64  ym     = fm_f->lat

    66  npp_i  = fm_i->NPP

    67  npp_f  = fm_f->NPP

    69  delete (fm_i)

    70  delete (fm_f)

    72 ;Units for these variables are:

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

    75  nsec_per_year = 60*60*24*365

    77  npp_i = npp_i *  nsec_per_year

    78  npp_f = npp_f *  nsec_per_year

    80  unit = "gC/m2/year"

    82 ;------------------------------

    83 ; get landfrac data

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

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

    87  landfrac = fm_l->landfrac

    89  npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:)

    90  npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:)

    92  delete (fm_l)

    93  delete (landfrac)

    95 ;-----------------------------

    96 ; read biome data: model

    98   biome_name_mod = "Model PFT Class"

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

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

   102   classmod = fm_c->CLASS_PFT

   104   delete (fm_c)

   106 ; model data has 17 land-type classes

   107   nclass_mod = 17

   109 ;---------------------------------------------------

   110 ; read data: observed at stations

   112  station = (/"DukeFACE" \

   113             ,"AspenFACE" \

   114             ,"ORNL-FACE" \

   115             ,"POP-EUROFACE" \

   116             /)

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

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

   120  lon_obx = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)

   122  n_sta  = dimsizes(station)

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

   125 ;###################################################

   126 ; this is a hardwired value

   127  beta_4_ob = 0.60

   128 ;###################################################

   129 ;---------------------------------------------------

   130 ; get model data at station

   132  npp_i_4  =linint2_points(xm,ym,npp_i,True,lon_obx,lat_ob,0)

   134  npp_f_4  =linint2_points(xm,ym,npp_f,True,lon_obx,lat_ob,0)

   136 ;---------------------------------------------------

   137 ;compute beta_4

   139  score_max = 3.

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

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

   145  beta_4_avg = avg(beta_4)

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

   148  Mbeta  = (1. - (bias/n_sta))*score_max

   149  M_beta = sprintf("%.2f", Mbeta)

   151 ;=========================

   152 ; for html table - station

   153 ;=========================

   155   output_html = "table_station.html"

   157 ; column (not including header column)

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

   161   ncol = dimsizes(col_head)

   163 ; row (not including header row)

   164   row_head = (/"DukeFACE" \

   165               ,"AspenFACE" \

   166               ,"ORNL-FACE" \

   167               ,"POP-EUROFACE" \

   168               ,"All Station" \

   169               /)

   170   nrow = dimsizes(row_head)

   172 ; arrays to be passed to table.

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

   175  do i=0,nrow-2

   176   text(i,0) = sprintf("%.1f",lat_ob(i))

   177   text(i,1) = sprintf("%.1f",lon_ob(i))

   178   text(i,2) = sprintf("%.1f",co2_i)

   179   text(i,3) = sprintf("%.1f",co2_f)

   180   text(i,4) = sprintf("%.1f",npp_i_4(0,i))

   181   text(i,5) = sprintf("%.1f",npp_f_4(0,i))

   182   text(i,6) = sprintf("%.2f",beta_4(i))

   183   text(i,7) = "-"

   184  end do

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

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

   187   text(nrow-1,2) = "-"

   188   text(nrow-1,3) = "-"

   189   text(nrow-1,4) = "-"

   190   text(nrow-1,5) = "-"

   191   text(nrow-1,6) = sprintf("%.2f",beta_4_avg)

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

   194 ;-----------

   195 ; html table

   196 ;-----------

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

   200   header = (/"<HTML>" \

   201             ,"<HEAD>" \

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

   203             ,"</HEAD>" \

   204             ,header_text \

   205             /)

   206   footer = "</HTML>"

   208   table_header = (/ \

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

   210        ,"<tr>" \

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

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

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

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

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

   216        ,"   <th bgcolor=DDDDDD >"+col_head(4)+"<br>("+unit+")</th>" \

   217        ,"   <th bgcolor=DDDDDD >"+col_head(5)+"<br>("+unit+")</th>" \

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

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

   220        ,"</tr>" \

   221        /)

   222   table_footer = "</table>"

   223   row_header = "<tr>"

   224   row_footer = "</tr>"

   226   lines = new(50000,string)

   227   nline = 0

   229   set_line(lines,nline,header)

   230   set_line(lines,nline,table_header)

   231 ;-----------------------------------------------

   232 ;row of table

   234   do n = 0,nrow-1

   235      set_line(lines,nline,row_header)

   237      txt1 = row_head(n)

   238      txt2 = text(n,0)

   239      txt3 = text(n,1)

   240      txt4 = text(n,2)

   241      txt5 = text(n,3)

   242      txt6 = text(n,4)

   243      txt7 = text(n,5)

   244      txt8 = text(n,6)

   245      txt9 = text(n,7)

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

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

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

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

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

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

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

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

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

   257      set_line(lines,nline,row_footer)

   258   end do

   259 ;-----------------------------------------------

   260   set_line(lines,nline,table_footer)

   261   set_line(lines,nline,footer)

   263 ; Now write to an HTML file.

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

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

   266     asciiwrite(output_html,lines(idx))

   267   else

   268    print ("error?")

   269   end if

   271   delete (col_head)

   272   delete (row_head)

   273   delete (text)

   274   delete (table_header)

   275   delete (idx)

   277 ;********************************************************************

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

   279 ;********************************************************************

   281 ; using model biome class

   282   nclass = nclass_mod

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

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

   287 ; particular range, and then take an average.

   289   nx = dimsizes(range) - 1

   291 ;==============================

   292 ; put data into bins

   293 ;==============================

   295 ; for model data and observed

   296   data_n = 2

   298 ; using model biome class

   300   base = ndtooned(classmod)

   302 ; output

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

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

   307 ; Loop through each range, using base

   309   do i=0,nx-1

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

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

   313      else

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

   315      end if

   317 ;    loop through each dataset

   319      do n = 0,data_n-1

   321         if (n .eq. 0) then

   322            data = ndtooned(npp_i)

   323         end if

   325         if (n .eq. 1) then

   326            data = ndtooned(npp_f)

   327         end if

   329 ;       Calculate average

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

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

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

   334         else

   335            yvalues(n,i) = yvalues@_FillValue

   336            count(n,i)   = 0

   337         end if

   339 ;#############################################################

   340 ; using model biome class:

   341 ;

   342 ;     set the following 4 classes to _FillValue:

   343 ;     (3)Needleleaf Deciduous Boreal Tree,

   344 ;     (8)Broadleaf Deciduous Boreal Tree,

   345 ;     (9)Broadleaf Evergreen Shrub,

   346 ;     (16)Wheat

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

   349          yvalues(n,i) = yvalues@_FillValue

   350          count(n,i)   = 0

   351       end if

   352 ;#############################################################

   354       delete(data)

   355     end do                 ; n-loop

   357     delete(idx)

   358   end do                   ; i-loop

   360   delete (base)

   361   delete (npp_i)

   362   delete (npp_f)

   364 ;============================

   365 ;compute beta

   366 ;============================

   368  u       = yvalues(0,:)

   369  v       = yvalues(1,:)

   370  u_count = count(0,:)

   371  v_count = count(1,:)

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

   375  uu       = u(good)

   376  vv       = v(good)

   377  uu_count = u_count(good)

   378  vv_count = v_count(good)

   380  n_biome = dimsizes(uu)

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

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

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

   387 ;===========================

   388 ; for html table - biome

   389 ;===========================

   391   output_html = "table_biome.html"

   393 ; column (not including header column)

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

   397   ncol = dimsizes(col_head)

   399 ; row (not including header row)

   401 ;----------------------------------------------------

   402 ; using model biome class:

   403 ;

   404   row_head  = (/"Not Vegetated" \

   405                ,"Needleleaf Evergreen Temperate Tree" \

   406                ,"Needleleaf Evergreen Boreal Tree" \

   407 ;              ,"Needleleaf Deciduous Boreal Tree" \

   408                ,"Broadleaf Evergreen Tropical Tree" \

   409                ,"Broadleaf Evergreen Temperate Tree" \

   410                ,"Broadleaf Deciduous Tropical Tree" \

   411                ,"Broadleaf Deciduous Temperate Tree" \

   412 ;              ,"Broadleaf Deciduous Boreal Tree" \

   413 ;              ,"Broadleaf Evergreen Shrub" \

   414                ,"Broadleaf Deciduous Temperate Shrub" \

   415                ,"Broadleaf Deciduous Boreal Shrub" \

   416                ,"C3 Arctic Grass" \

   417                ,"C3 Non-Arctic Grass" \

   418                ,"C4 Grass" \

   419                ,"Corn" \

   420 ;              ,"Wheat" \

   421                ,"All Biome" \

   422                /)

   424   nrow = dimsizes(row_head)

   426 ; arrays to be passed to table.

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

   429  do i=0,nrow-2

   430   text(i,0) = sprintf("%.1f",co2_i)

   431   text(i,1) = sprintf("%.1f",co2_f)

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

   433   text(i,3) = sprintf("%.1f",vv(i))

   434   text(i,4) = sprintf("%.2f",beta_biome(i))

   435  end do

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

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

   438   text(nrow-1,2) = "-"

   439   text(nrow-1,3) = "-"

   440   text(nrow-1,4) = sprintf("%.2f",beta_biome_avg)

   442 ;**************************************************

   443 ; html table

   444 ;**************************************************

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

   448   header = (/"<HTML>" \

   449             ,"<HEAD>" \

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

   451             ,"</HEAD>" \

   452             ,header_text \

   453             /)

   454   footer = "</HTML>"

   456   table_header = (/ \

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

   458        ,"<tr>" \

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

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

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

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

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

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

   465        ,"</tr>" \

   466        /)

   467   table_footer = "</table>"

   468   row_header = "<tr>"

   469   row_footer = "</tr>"

   471   lines = new(50000,string)

   472   nline = 0

   474   set_line(lines,nline,header)

   475   set_line(lines,nline,table_header)

   476 ;-----------------------------------------------

   477 ;row of table

   479   do n = 0,nrow-1

   480      set_line(lines,nline,row_header)

   482      txt1  = row_head(n)

   483      txt2  = text(n,0)

   484      txt3  = text(n,1)

   485      txt4  = text(n,2)

   486      txt5  = text(n,3)

   487      txt6  = text(n,4)

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

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

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

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

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

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

   496      set_line(lines,nline,row_footer)

   497   end do

   498 ;-----------------------------------------------

   499   set_line(lines,nline,table_footer)

   500   set_line(lines,nline,footer)

   502 ; Now write to an HTML file

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

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

   506     asciiwrite(output_html,lines(idx))

   507   else

   508    print ("error?")

   509   end if

   510   delete (idx)

   512 ;**************************************************************************************

   513 ; update score

   514 ;**************************************************************************************

   515   if (isvar("compare")) then

   516      system("sed -e '1,/M_beta/s/M_beta/"+M_beta+"/' "+html_name2+" > "+html_new2+";"+ \

   517             "mv -f "+html_new2+" "+html_name2)

   518   end if

   520   system("sed s#M_beta#"+M_beta+"# "+html_name+" > "+html_new+";"+ \

   521          "mv -f "+html_new+" "+html_name)

   523 ;***************************************************************************

   524 ; get total score and write to file

   525 ;***************************************************************************

   526   M_total = Mbeta

   528   asciiwrite("M_save.beta", M_total)

   530   delete (M_total)

   532 ;***************************************************************************

   533 ; output plot and html

   534 ;***************************************************************************

   535   output_dir = model_name+"/beta"

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

   538 ;***************************************************************************

   540 end

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