clamp: beta/04.biome.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   nclass = 20

    31   plot_type     = "ps"

    32   plot_type_new = "png"

    33 ;************************************************

    34 ; read data: model

    35 ;************************************************

    36  co2_i = 283.1878

    37  co2_f = 364.1252

    39  model_grid = "T42"

    41 ;model_name_i = "i01.07cn"

    42 ;model_name_f = "i01.10cn"

    44  model_name_i = "i01.07casa"

    45  model_name_f = "i01.10casa"

    47  model_name = model_name_f

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

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

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

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

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

    56  npp_i  = fm_i->NPP

    57  npp_f  = fm_f->NPP

    59 ;************************************************

    60 ; read data: observed

    61 ;************************************************

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

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

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

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

    70   classob = tofloat(fo->LAND_CLASS)

    72 ;*******************************************************************

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

    74 ;********************************************************************

    75   nclassn     = nclass + 1

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

    77 ; print (range)

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

    80 ; particular range, and then take an average.

    82   nr           = dimsizes(range)

    83   nx           = nr-1

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

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

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

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

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

    90 ; get data

    92   DATA11_1D = ndtooned(classob)

    93   DATA12_1D = ndtooned(npp_i)

    94   DATA22_1D = ndtooned(npp_f)

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

    97   mn_yvalues   = new((/2,nx/),float)

    98   mx_yvalues   = new((/2,nx/),float)

   100   do nd=0,1

   102 ; See if we are doing model or observational data.

   104     if(nd.eq.0) then

   105       data_ob  = DATA11_1D

   106       data_mod = DATA12_1D

   107     else

   108       data_ob  = DATA11_1D

   109       data_mod = DATA22_1D

   110     end if

   112 ; Loop through each range and check for values.

   114     do i=0,nr-2

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

   116 ;        print("")

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

   118          idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))

   119       else

   120 ;        print("")

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

   122          idx = ind(data_ob.ge.range(i))

   123       end if

   125 ; Calculate average, and get min and max.

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

   128         yvalues(nd,i)    = avg(data_mod(idx))

   129         mn_yvalues(nd,i) = min(data_mod(idx))

   130         mx_yvalues(nd,i) = max(data_mod(idx))

   131         count = dimsizes(idx)

   132       else

   133         count            = 0

   134         yvalues(nd,i)    = yvalues@_FillValue

   135         mn_yvalues(nd,i) = yvalues@_FillValue

   136         mx_yvalues(nd,i) = yvalues@_FillValue

   137       end if

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

   140 ;     print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))

   142 ; Clean up for next time in loop.

   144       delete(idx)

   145     end do

   146     delete(data_ob)

   147     delete(data_mod)

   148   end do

   149 ;============================

   150 ;compute beta

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

   153  nsec_per_year = 60*60*24*365

   155  u = yvalues(0,:)

   156  v = yvalues(1,:)

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

   159  uu = u(good)* nsec_per_year

   160  vv = v(good)* nsec_per_year

   162  n_biome = dimsizes(uu)

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

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

   168  beta_biome_avg = avg(beta_biome)

   170  print (beta_biome_avg)

   171 ;*******************************************************************

   172 ; for html table

   173 ;*******************************************************************

   175 ; column (not including header column)

   177   col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta"/)

   179   ncol = dimsizes(col_head)

   181 ; row (not including header row)

   182   row_head  = (/"Water Bodies" \

   183                ,"Evergreen Needleleaf Forests" \

   184                ,"Evergreen Broadleaf Forests" \

   185                ,"Deciduous Needleleaf Forest" \

   186                ,"Deciduous Broadleaf Forests" \

   187                ,"Mixed Forests" \

   188                ,"Closed Bushlands" \

   189                ,"Open Bushlands" \

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

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

   192                ,"Grasslands" \

   193                ,"Permanent Wetlands" \

   194                ,"Croplands" \

   195                ,"Cropland/Natural Vegetation Mosaic" \

   196                ,"Permanent Snow and Ice" \

   197                ,"Barren or Sparsely Vegetated" \

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

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

   200                ,"All Biome" \

   201                /)

   202   nrow = dimsizes(row_head)

   204 ; arrays to be passed to table.

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

   207  do i=0,nrow-2

   208   text4(i,0) = sprintf("%.2f",co2_i)

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

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

   211   text4(i,3) = sprintf("%.2f",vv(i))

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

   213  end do

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

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

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

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

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

   220 ;**************************************************

   221 ; html table

   222 ;**************************************************

   223   output_html = "table_biome.html"

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

   227   header = (/"<HTML>" \

   228             ,"<HEAD>" \

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

   230             ,"</HEAD>" \

   231             ,header_text \

   232             /)

   233   footer = "</HTML>"

   235   table_header = (/ \

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

   237        ,"<tr>" \

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

   239        ,"   <th bgcolor=DDDDDD >CO2_i</th>" \

   240        ,"   <th bgcolor=DDDDDD >CO2_f</th>" \

   241        ,"   <th bgcolor=DDDDDD >NPP_i</th>" \

   242        ,"   <th bgcolor=DDDDDD >NPP_f</th>" \

   243        ,"   <th bgcolor=DDDDDD >Beta</th>" \

   244        ,"</tr>" \

   245        /)

   246   table_footer = "</table>"

   247   row_header = "<tr>"

   248   row_footer = "</tr>"

   250   lines = new(50000,string)

   251   nline = 0

   253   set_line(lines,nline,header)

   254   set_line(lines,nline,table_header)

   255 ;-----------------------------------------------

   256 ;row of table

   258   do n = 0,nrow-1

   259      set_line(lines,nline,row_header)

   261      txt1  = row_head(n)

   262      txt2  = text4(n,0)

   263      txt3  = text4(n,1)

   264      txt4  = text4(n,2)

   265      txt5  = text4(n,3)

   266      txt6  = text4(n,4)

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

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

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

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

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

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

   275      set_line(lines,nline,row_footer)

   276   end do

   277 ;-----------------------------------------------

   278   set_line(lines,nline,table_footer)

   279   set_line(lines,nline,footer)

   281 ; Now write to an HTML file.

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

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

   284     asciiwrite(output_html,lines(idx))

   285   else

   286    print ("error?")

   287   end if

   289 end

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