clamp: all/09.carbon_sink.ncl@4be95183fbcd

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

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

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

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

     5 ;**************************************************************

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

     7 begin

     8 ; add line to ascci/html file

    10   nnewlines = dimsizes(newlines)

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

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

    13     return

    14   end if

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

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

    17   nline = nline + nnewlines

    18   return

    19 end

    20 ;**************************************************************

    21 ; Main code.

    22 begin

    24  plot_type     = "ps"

    25  plot_type_new = "png"

    27 ;----------------------------------------------------------

    28 ; edit current model for movel1_vs_model2

    30  if (isvar("compare")) then

    31     html_name2 = compare+"/table.html"

    32     html_new2  = html_name2 +".new"

    33  end if

    35 ;----------------------------------------------------------

    36 ; edit table.html for current model

    38  html_name = model_name+"/table.html"

    39  html_new  = html_name +".new"

    41 ;----------------------------------------------------------

    42 ; get biome data: model

    44   biome_name_mod = "Model PFT Class"

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

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

    48   classmod = fm_c->CLASS_PFT

    50   delete (fm_c)

    52 ; model data has 17 land-type classes

    53   nclass_mod = 17

    55 ;------------------------------------

    56 ; get model landfrac and area

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

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

    60  landfrac = fm_l->landfrac

    61  area     = fm_l->area

    63  delete (fm_l)

    65 ; change area from km**2 to m**2

    66  area = area * 1.e6

    68 ;-------------------------------------

    69 ; take into account landfrac

    71 ; area = area * landfrac

    73 ; delete (landfrac)

    75 ;---------------------------------------------------

    76 ; read data: model, group 1

    78  fm      = addfile (dirm+film4,"r")

    80  NPP1    = fm->NPP

    82  leafc   = fm->LEAFC

    83  woodc   = fm->WOODC

    84  frootc  = fm->FROOTC

    85  VegC    = leafc

    86  VegC    = leafc + woodc + frootc

    88  litterc = fm->LITTERC

    89  cwdc    = fm->CWDC

    90  LiCwC   = litterc

    91  LiCwC   = litterc + cwdc

    93  SoilC   = fm->SOILC

    95  delete (fm)

    96 ;---------------------------------------------------

    97 ; read data: model, group 2

    99  fm     = addfile (dirm+film5,"r")

   101  NPP2   = fm->NPP

   102  NEE2   = fm->NEE

   103  GPP2   = fm->GPP

   105  delete (fm)

   106 ;---------------------------------------------------

   107 ; Units for these variables are:

   109 ;NPP1: g C/m^2/s

   110 ;NPP2: g C/m^2/s

   111 ;NEE2: g C/m^2/s

   112 ;GPP2: g C/m^2/s

   114 ;VegC:  g C/m^2

   115 ;LiCwC: g C/m^2

   116 ;SoilC: g C/m^2

   118  nsec_per_year = 60*60*24*365

   120 ; change unit to g C/m^2/year

   122  NPP1 = NPP1 *  nsec_per_year * conform(NPP1,landfrac,(/1,2/))

   123  NPP2 = NPP2 *  nsec_per_year * conform(NPP2,landfrac,(/1,2/))

   124  NEE2 = NEE2 *  nsec_per_year * conform(NEE2,landfrac,(/1,2/))

   125  GPP2 = GPP2 *  nsec_per_year * conform(GPP2,landfrac,(/1,2/))

   127  VegC  = VegC  * conform(VegC,landfrac,(/1,2/))

   128  LiCwC = LiCwC * conform(LiCwC,landfrac,(/1,2/))

   129  SoilC = SoilC * conform(SoilC,landfrac,(/1,2/))

   131  data_n = 8

   133 ;*******************************************************************

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

   135 ;********************************************************************

   137 ; using model biome class

   138   nclass = nclass_mod

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

   142 ; print (range)

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

   144 ; particular range, and then take an average.

   146   nx     = dimsizes(range) - 1

   148 ;==============================

   149 ; put data into bins

   150 ;==============================

   152 ; using observed biome class

   153 ; base  = ndtooned(classob)

   154 ; using model biome class

   155   base  = ndtooned(classmod)

   157   area_1d = ndtooned(area)

   159 ; output

   161   yvalues   = new((/data_n,nx/),float) ; (per m2)

   162   yvalues_t = new((/data_n,nx/),float) ; (per biome)

   164 ; Loop through each range, using base.

   166   do i=0,nx-1

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

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

   170      else

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

   172      end if

   174   do n = 0,data_n-1

   176      if (n.eq.0) then

   177         data = ndtooned(area)

   178      end if

   180      if (n.eq.1) then

   181         data = ndtooned(NPP1)

   182      end if

   184      if (n.eq.2) then

   185         data = ndtooned(VegC)

   186      end if

   188      if (n.eq.3) then

   189         data = ndtooned(LiCwC)

   190      end if

   192      if (n.eq.4) then

   193         data = ndtooned(SoilC)

   194      end if

   196      if (n.eq.5) then

   197         data = ndtooned(NPP2)

   198      end if

   200      if (n.eq.6) then

   201         data = ndtooned(NEE2)

   202      end if

   204      if (n.eq.7) then

   205         data = ndtooned(GPP2)

   206      end if

   208 ;    Calculate sum and average

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

   211         if (n.eq.0) then

   212            yvalues(n,i)   = sum(data(idx))

   213            yvalues_t(n,i) = sum(data(idx))

   214         else

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

   216            yvalues_t(n,i) = sum(data(idx)*area_1d(idx))

   217         end if

   218      else

   219         yvalues(n,i)   = yvalues@_FillValue

   220         yvalues_t(n,i) = yvalues@_FillValue

   221      end if

   223 ;#############################################################

   224 ; using model biome class:

   225 ;     set the following 4 classes to _FillValue:

   226 ;     (3)Needleleaf Deciduous Boreal Tree,

   227 ;     (8)Broadleaf Deciduous Boreal Tree,

   228 ;     (9)Broadleaf Evergreen Shrub,

   229 ;     (16)Wheat

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

   232          yvalues(n,i)   = yvalues@_FillValue

   233          yvalues_t(n,i) = yvalues@_FillValue

   234       end if

   235 ;#############################################################

   237      delete (data)

   238   end do

   240   delete (idx)

   241   end do

   243   delete (base)

   244   delete (area)

   245   delete (NPP1)

   246   delete (VegC)

   247   delete (LiCwC)

   248   delete (SoilC)

   249   delete (NPP2)

   250   delete (NEE2)

   251   delete (GPP2)

   253 ;----------------------------------------------------------------

   254 ; data for table1

   256  good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))

   257 ;print (good)

   259  area_g  = yvalues(0,good)

   260  NPP1_g  = yvalues(1,good)

   261  VegC_g  = yvalues(2,good)

   262  LiCwC_g = yvalues(3,good)

   263  SoilC_g = yvalues(4,good)

   264  NPP2_g  = yvalues(5,good)

   265  NEE2_g  = yvalues(6,good)

   266  GPP2_g  = yvalues(7,good)

   268  NPP_ratio = NPP2_g/NPP1_g

   270  n_biome = dimsizes(NPP1_g)

   272 ;-----------------------------------------------------------------

   273 ; data for table2

   275 ; change unit from g to Pg (Peta gram)

   276  factor_unit = 1.e-15

   278  NPP1_t  = yvalues_t(1,good) * factor_unit

   279  VegC_t  = yvalues_t(2,good) * factor_unit

   280  LiCwC_t = yvalues_t(3,good) * factor_unit

   281  SoilC_t = yvalues_t(4,good) * factor_unit

   282  NEE2_t  = yvalues_t(6,good) * factor_unit

   283  GPP2_t  = yvalues_t(7,good) * factor_unit

   285  delete (yvalues)

   286  delete (yvalues_t)

   288 ;-------------------------------------------------------------

   289 ; html table1 data

   291 ; column (not including header column)

   293   col_head  = (/"Area (1.e12m2)" \

   294                ,"NPP (gC/m2/yr)" \

   295                ,"VegC (gC/m2)" \

   296                ,"Litter+CWD (gC/m2)" \

   297                ,"SoilC (gC/m2)" \

   298                ,"NPP_ratio" \

   299                ,"NEE (gC/m2/yr)" \

   300                ,"GPP (gC/m2/yr)" \

   301                /)

   303   ncol = dimsizes(col_head)

   305 ; row (not including header row)

   307 ; using model biome class:

   308   row_head  = (/"Not Vegetated" \

   309                ,"Needleleaf Evergreen Temperate Tree" \

   310                ,"Needleleaf Evergreen Boreal Tree" \

   311 ;              ,"Needleleaf Deciduous Boreal Tree" \

   312                ,"Broadleaf Evergreen Tropical Tree" \

   313                ,"Broadleaf Evergreen Temperate Tree" \

   314                ,"Broadleaf Deciduous Tropical Tree" \

   315                ,"Broadleaf Deciduous Temperate Tree" \

   316 ;              ,"Broadleaf Deciduous Boreal Tree" \

   317 ;              ,"Broadleaf Evergreen Shrub" \

   318                ,"Broadleaf Deciduous Temperate Shrub" \

   319                ,"Broadleaf Deciduous Boreal Shrub" \

   320                ,"C3 Arctic Grass" \

   321                ,"C3 Non-Arctic Grass" \

   322                ,"C4 Grass" \

   323                ,"Corn" \

   324 ;              ,"Wheat" \

   325                ,"All Biome" \

   326                /)

   327   nrow = dimsizes(row_head)

   329 ; arrays to be passed to table.

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

   332  do i=0,nrow-2

   333   text(i,0) = sprintf("%.1f",area_g(i)*1.e-12)

   334   text(i,1) = sprintf("%.1f",NPP1_g(i))

   335   text(i,2) = sprintf("%.1f",VegC_g(i))

   336   text(i,3) = sprintf("%.1f",LiCwC_g(i))

   337   text(i,4) = sprintf("%.1f",SoilC_g(i))

   338   text(i,5) = sprintf("%.2f",NPP_ratio(i))

   339   text(i,6) = sprintf("%.1f",NEE2_g(i))

   340   text(i,7) = sprintf("%.1f",GPP2_g(i))

   341  end do

   343 ;-------------------------------------------------------

   344 ; create html table1

   346   header_text = "<H1>NEE and Carbon Stocks and Fluxes:  Model "+model_name+"</H1>"

   348   header = (/"<HTML>" \

   349             ,"<HEAD>" \

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

   351             ,"</HEAD>" \

   352             ,header_text \

   353             /)

   354   footer = "</HTML>"

   356   table_header = (/ \

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

   358        ,"<tr>" \

   359        ,"   <th bgcolor=DDDDDD >Biome Type</th>" \

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

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

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

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

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

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

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

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

   368        ,"</tr>" \

   369        /)

   370   table_footer = "</table>"

   371   row_header = "<tr>"

   372   row_footer = "</tr>"

   374   lines = new(50000,string)

   375   nline = 0

   377   set_line(lines,nline,header)

   378   set_line(lines,nline,table_header)

   380 ;----------------------------

   381 ;row of table

   383   do n = 0,nrow-2

   384      set_line(lines,nline,row_header)

   386      txt0  = row_head(n)

   387      txt1  = text(n,0)

   388      txt2  = text(n,1)

   389      txt3  = text(n,2)

   390      txt4  = text(n,3)

   391      txt5  = text(n,4)

   392      txt6  = text(n,5)

   393      txt7  = text(n,6)

   394      txt8  = text(n,7)

   396      set_line(lines,nline,"<th>"+txt0+"</th>")

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

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

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

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

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

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

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

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

   406      set_line(lines,nline,row_footer)

   407   end do

   408 ;----------------------------

   409   set_line(lines,nline,table_footer)

   410   set_line(lines,nline,footer)

   412 ; Now write to an HTML file.

   414   output_html = "table_per_m2.html"

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

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

   418     asciiwrite(output_html,lines(idx))

   419   else

   420    print ("error?")

   421   end if

   423   delete (idx)

   425   delete (col_head)

   426   delete (row_head)

   427   delete (text)

   428   delete (table_header)

   430 ;-----------------------------------------------------------------

   431 ; html table2 data

   433 ; column (not including header column)

   435   col_head  = (/"NPP (PgC/yr)" \

   436                ,"VegC (PgC)" \

   437                ,"Litter+CWD (PgC)" \

   438                ,"SoilC (PgC)" \

   439                ,"NEE (PgC/yr)" \

   440                ,"GPP (PgC/yr)" \

   441                ,"NPP timeseries" \

   442                ,"NEE timeseries" \

   443                ,"Fire timeseries" \

   444                /)

   446   ncol = dimsizes(col_head)

   448 ; row (not including header row)

   450 ; using model biome class:

   451   row_head  = (/"Not Vegetated" \

   452                ,"Needleleaf Evergreen Temperate Tree" \

   453                ,"Needleleaf Evergreen Boreal Tree" \

   454 ;              ,"Needleleaf Deciduous Boreal Tree" \

   455                ,"Broadleaf Evergreen Tropical Tree" \

   456                ,"Broadleaf Evergreen Temperate Tree" \

   457                ,"Broadleaf Deciduous Tropical Tree" \

   458                ,"Broadleaf Deciduous Temperate Tree" \

   459 ;              ,"Broadleaf Deciduous Boreal Tree" \

   460 ;              ,"Broadleaf Evergreen Shrub" \

   461                ,"Broadleaf Deciduous Temperate Shrub" \

   462                ,"Broadleaf Deciduous Boreal Shrub" \

   463                ,"C3 Arctic Grass" \

   464                ,"C3 Non-Arctic Grass" \

   465                ,"C4 Grass" \

   466                ,"Corn" \

   467 ;              ,"Wheat" \

   468                ,"All Biome" \

   469                /)

   470   nrow = dimsizes(row_head)

   472 ; arrays to be passed to table.

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

   475  do i=0,nrow-2

   476   text(i,0) = sprintf("%.1f",NPP1_t(i))

   477   text(i,1) = sprintf("%.1f",VegC_t(i))

   478   text(i,2) = sprintf("%.1f",LiCwC_t(i))

   479   text(i,3) = sprintf("%.1f",SoilC_t(i))

   480   text(i,4) = sprintf("%.1f",NEE2_t(i))

   481   text(i,5) = sprintf("%.1f",GPP2_t(i))

   482   text(i,6) = "<a href=./NPP_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NPP_annual_biome_"+i+".png>annual_plot</a>"

   483   text(i,7) = "<a href=./NEE_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NEE_annual_biome_"+i+".png>annual_plot</a>"

   484   text(i,8) = "<a href=./Fire_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./Fire_annual_biome_"+i+".png>annual_plot</a>"

   485  end do

   487   text(nrow-1,0) = sprintf("%.1f",sum(NPP1_t))

   488   text(nrow-1,1) = sprintf("%.1f",sum(VegC_t))

   489   text(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t))

   490   text(nrow-1,3) = sprintf("%.1f",sum(SoilC_t))

   491   text(nrow-1,4) = sprintf("%.1f",sum(NEE2_t))

   492   text(nrow-1,5) = sprintf("%.1f",sum(GPP2_t))

   493   text(nrow-1,6) = "<a href=./NPP_monthly_global.png>monthly_plot</a> <br> <a href=./NPP_annual_global.png>annual_plot</a>"

   494   text(nrow-1,7) = "<a href=./NEE_monthly_global.png>monthly_plot</a> <br> <a href=./NEE_annual_global.png>annual_plot</a>"

   495   text(nrow-1,8) = "<a href=./Fire_monthly_global.png>monthly_plot</a> <br> <a href=./Fire_annual_global.png>annual_plot</a>"

   497 ;**************************************************

   498 ; create html table2

   499 ;**************************************************

   501   header_text = "<H1>NEE and Carbon Stocks and Fluxes (per biome):  Model "+model_name+"</H1>"

   503   header = (/"<HTML>" \

   504             ,"<HEAD>" \

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

   506             ,"</HEAD>" \

   507             ,header_text \

   508             /)

   509   footer = "</HTML>"

   511   table_header = (/ \

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

   513        ,"<tr>" \

   514        ,"   <th bgcolor=DDDDDD >Biome Type</th>" \

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

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

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

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

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

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

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

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

   523        ,"   <th bgcolor=DDDDDD >"+col_head(8)+"</th>" \

   524        ,"</tr>" \

   525        /)

   526   table_footer = "</table>"

   527   row_header = "<tr>"

   528   row_footer = "</tr>"

   530   lines = new(50000,string)

   531   nline = 0

   533   set_line(lines,nline,header)

   534   set_line(lines,nline,table_header)

   535 ;-----------------------------------------------

   536 ;row of table

   538   do n = 0,nrow-1

   539      set_line(lines,nline,row_header)

   541      txt0  = row_head(n)

   542      txt1  = text(n,0)

   543      txt2  = text(n,1)

   544      txt3  = text(n,2)

   545      txt4  = text(n,3)

   546      txt5  = text(n,4)

   547      txt6  = text(n,5)

   548      txt7  = text(n,6)

   549      txt8  = text(n,7)

   550      txt9  = text(n,8)

   552      set_line(lines,nline,"<th>"+txt0+"</th>")

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

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

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

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

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

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

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

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

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

   563      set_line(lines,nline,row_footer)

   564   end do

   565 ;-----------------------------------------------

   566   set_line(lines,nline,table_footer)

   567   set_line(lines,nline,footer)

   569 ; Now write to an HTML file.

   571   output_html = "table_per_biome.html"

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

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

   575     asciiwrite(output_html,lines(idx))

   576   else

   577    print ("error?")

   578   end if

   580   delete (idx)

   582 ;---------------------------------------------------

   583 ; read model data, time series:

   585  fm   = addfile (dirm+film7,"r")

   587  NPP3 = fm->NPP

   588  NEE3 = fm->NEE

   589  Fire = fm->COL_FIRE_CLOSS

   591  delete (fm)

   593 ; Units for these variables are:

   595 ;NPP3: g C/m^2/s

   596 ;NEE3: g C/m^2/s

   597 ;Fire: g C/m^2/s

   599  nsec_per_month = 60*60*24*30

   601 ; change unit to g C/m^2/month

   603  NPP3 = NPP3 * nsec_per_month * conform(NPP3,landfrac,(/2,3/))

   604  NEE3 = NEE3 * nsec_per_month * conform(NEE3,landfrac,(/2,3/))

   605  Fire = Fire * nsec_per_month * conform(Fire,landfrac,(/2,3/))

   607  data_n = 3

   609  dsizes = dimsizes(NPP3)

   610  nyear  = dsizes(0)

   611  nmonth = dsizes(1)

   612  ntime  = nyear * nmonth

   614  year_start = 1979

   615  year_end   = 2004

   617 ;*******************************************************************

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

   619 ;********************************************************************

   621 ; using model biome class

   622   nclass = nclass_mod

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

   626 ; print (range)

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

   628 ; particular range, and then take an average.

   630   nx = dimsizes(range) - 1

   632 ;==============================

   633 ; put data into bins

   634 ;==============================

   636 ; using observed biome class

   637 ; base  = ndtooned(classob)

   638 ; using model biome class

   639   base  = ndtooned(classmod)

   641 ; output

   643   yvalues = new((/ntime,data_n,nx/),float)

   645 ; Loop through each range, using base.

   647   do i=0,nx-1

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

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

   651      else

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

   653      end if

   655   do n = 0,data_n-1

   657      t = -1

   658      do m = 0,nyear-1

   659      do k = 0,nmonth-1

   661         t = t + 1

   663         if (n.eq.0) then

   664            data = ndtooned(NPP3(m,k,:,:))

   665         end if

   667         if (n.eq.1) then

   668            data = ndtooned(NEE3(m,k,:,:))

   669         end if

   671         if (n.eq.2) then

   672            data = ndtooned(Fire(m,k,:,:))

   673         end if

   675 ;       Calculate average

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

   678            yvalues(t,n,i) = sum(data(idx)*area_1d(idx))

   679         else

   680            yvalues(t,n,i) = yvalues@_FillValue

   681         end if

   683 ;#############################################################

   684 ; using model biome class:

   685 ;     set the following 4 classes to _FillValue:

   686 ;     (3)Needleleaf Deciduous Boreal Tree,

   687 ;     (8)Broadleaf Deciduous Boreal Tree,

   688 ;     (9)Broadleaf Evergreen Shrub,

   689 ;     (16)Wheat

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

   692            yvalues(t,n,i) = yvalues@_FillValue

   693         end if

   694 ;#############################################################

   696      end do

   697      end do

   699      delete(data)

   700   end do

   702     delete(idx)

   703   end do

   705   delete (base)

   706   delete (NPP3)

   707   delete (NEE3)

   708   delete (Fire)

   710 ;----------------------------------------------------------------

   711 ; data for tseries plot

   713   yvalues_g = new((/ntime,data_n,n_biome/),float)

   715   yvalues_g@units = "TgC/month"

   717 ; change unit to Tg C/month

   718 ; change unit from g to Tg (Tera gram)

   719   factor_unit = 1.e-12

   721   yvalues_g = yvalues(:,:,good) * factor_unit

   723 ;*******************************************************************

   724 ; general settings for line plot

   725 ;*******************************************************************

   727 ; res

   728   res                   = True

   729   res@xyDashPatterns    = (/0/)                ; make lines solid

   730   res@xyLineThicknesses = (/2.0/)          ; make lines thicker

   731   res@xyLineColors      = (/"blue"/) ; line color

   733   res@trXMinF   = year_start

   734   res@trXMaxF   = year_end + 1

   736   res@vpHeightF = 0.4                 ; change aspect ratio of plot

   737 ; res@vpWidthF  = 0.8

   738   res@vpWidthF  = 0.75

   740 ; res@gsnMaximize = True

   742 ;*******************************************************************

   743 ; (A) 1 component in each biome: monthly

   744 ;*******************************************************************

   746   component = (/"NPP","NEE","Fire"/)

   748 ; for x-axis in xyplot

   750   timeI = new((/ntime/),integer)

   751   timeF = new((/ntime/),float)

   752   timeI = ispan(1,ntime,1)

   753   timeF = year_start + (timeI-1)/12.

   754   timeF@long_name = "year"

   756   plot_data = new((/ntime/),float)

   757   plot_data@long_name = "TgC/month"

   759   do n = 0, data_n-1

   760   do m = 0, n_biome-1

   762      plot_name = component(n)+"_monthly_biome_"+ m

   764      wks = gsn_open_wks (plot_type,plot_name)

   766      title = component(n)+ ": "+ row_head(m)

   767      res@tiMainString = title

   768      res@tiMainFontHeightF = 0.025

   770      plot_data(:) = yvalues_g(:,n,m)

   772      plot=gsn_csm_xy(wks,timeF,plot_data,res)

   774      delete (wks)

   775      delete (plot)

   777     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   778            "rm "+plot_name+"."+plot_type)

   779   end do

   780   end do

   782   do n = 0, data_n-1

   784      plot_name = component(n)+"_monthly_global"

   786      wks = gsn_open_wks (plot_type,plot_name)

   788      title = component(n)+ ": Global"

   789      res@tiMainString = title

   790      res@tiMainFontHeightF = 0.025

   792      do k = 0,ntime-1

   793         plot_data(k) = sum(yvalues_g(k,n,:))

   794      end do

   796      plot=gsn_csm_xy(wks,timeF,plot_data,res)

   798      delete (wks)

   799      delete (plot)

   801     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   802            "rm "+plot_name+"."+plot_type)

   803   end do

   805   delete (plot_data)

   806   delete (timeI)

   807   delete (timeF)

   809 ;*******************************************************************

   810 ; (B) 1 component in each biome: annually

   811 ;*******************************************************************

   813   yvalues_a = new((/nyear,data_n,n_biome/),float)

   814   yvalues_g!0 = "time"

   815   yvalues_g!1 = "case"

   816   yvalues_g!2 = "record"

   818   yvalues_a = month_to_annual(yvalues_g,0)

   820   delete (yvalues_g)

   822 ; for x-axis in xyplot

   824   timeI = new((/nyear/),integer)

   825   timeF = new((/nyear/),float)

   826   timeI = ispan(1,nyear,1)

   827   timeF = year_start + (timeI-1)

   828   timeF@long_name = "year"

   830   plot_data = new((/nyear/),float)

   831   plot_data@long_name = "TgC/year"

   833   do n = 0, data_n-1

   834   do m = 0, n_biome-1

   836      plot_name = component(n)+"_annual_biome_"+ m

   838      wks = gsn_open_wks (plot_type,plot_name)

   840      title = component(n)+ ": "+ row_head(m)

   841      res@tiMainString = title

   842      res@tiMainFontHeightF = 0.025

   844      plot_data(:) = yvalues_a(:,n,m)

   846      plot=gsn_csm_xy(wks,timeF,plot_data,res)

   848      delete (wks)

   849      delete (plot)

   851     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   852            "rm "+plot_name+"."+plot_type)

   853   end do

   854   end do

   856   do n = 0, data_n-1

   858      plot_name = component(n)+"_annual_global"

   860      wks = gsn_open_wks (plot_type,plot_name)

   862      title = component(n)+ ": Global"

   863      res@tiMainString = title

   864      res@tiMainFontHeightF = 0.025

   866      do k = 0,nyear-1

   867         plot_data(k) = sum(yvalues_a(k,n,:))

   868      end do

   870      plot=gsn_csm_xy(wks,timeF,plot_data,res)

   872      delete (wks)

   873      delete (plot)

   875     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   876            "rm "+plot_name+"."+plot_type)

   877   end do

   879 ;****************************************

   880 ; output plot and html

   881 ;****************************************

   882   output_dir = model_name+"/carbon_sink"

   884   system("mv *.png *.html " + output_dir)

   885 ;****************************************

   887 end

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