clamp: all/09x.carbon_sink.ncl@0c6405ab2ff4

     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) = "--"

   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) = "--"

   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

   608  data_n = 2

   610  dsizes = dimsizes(NPP3)

   611  nyear  = dsizes(0)

   612  nmonth = dsizes(1)

   613  ntime  = nyear * nmonth

   615  year_start = 1979

   616  year_end   = 2004

   618 ;*******************************************************************

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

   620 ;********************************************************************

   622 ; using model biome class

   623   nclass = nclass_mod

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

   627 ; print (range)

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

   629 ; particular range, and then take an average.

   631   nx = dimsizes(range) - 1

   633 ;==============================

   634 ; put data into bins

   635 ;==============================

   637 ; using observed biome class

   638 ; base  = ndtooned(classob)

   639 ; using model biome class

   640   base  = ndtooned(classmod)

   642 ; output

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

   646 ; Loop through each range, using base.

   648   do i=0,nx-1

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

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

   652      else

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

   654      end if

   656   do n = 0,data_n-1

   658      t = -1

   659      do m = 0,nyear-1

   660      do k = 0,nmonth-1

   662         t = t + 1

   664         if (n.eq.0) then

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

   666         end if

   668         if (n.eq.1) then

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

   670         end if

   672 ;       if (n.eq.2) then

   673 ;          data = ndtooned(Fire(m,k,:,:))

   674 ;       end if

   676 ;       Calculate average

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

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

   680         else

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

   682         end if

   684 ;#############################################################

   685 ; using model biome class:

   686 ;     set the following 4 classes to _FillValue:

   687 ;     (3)Needleleaf Deciduous Boreal Tree,

   688 ;     (8)Broadleaf Deciduous Boreal Tree,

   689 ;     (9)Broadleaf Evergreen Shrub,

   690 ;     (16)Wheat

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

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

   694         end if

   695 ;#############################################################

   697      end do

   698      end do

   700      delete(data)

   701   end do

   703     delete(idx)

   704   end do

   706   delete (base)

   707   delete (NPP3)

   708   delete (NEE3)

   709 ; delete (Fire)

   711 ;----------------------------------------------------------------

   712 ; data for tseries plot

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

   716   yvalues_g@units = "TgC/month"

   718 ; change unit to Tg C/month

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

   720   factor_unit = 1.e-12

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

   724 ;*******************************************************************

   725 ; general settings for line plot

   726 ;*******************************************************************

   728 ; res

   729   res                   = True

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

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

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

   734   res@trXMinF   = year_start

   735   res@trXMaxF   = year_end + 1

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

   738 ; res@vpWidthF  = 0.8

   739   res@vpWidthF  = 0.75

   741 ; res@gsnMaximize = True

   743 ;*******************************************************************

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

   745 ;*******************************************************************

   747 ; component = (/"NPP","NEE","Fire"/)

   748   component = (/"NPP","NEE"/)

   750 ; for x-axis in xyplot

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

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

   754   timeI = ispan(1,ntime,1)

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

   756   timeF@long_name = "year"

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

   759   plot_data@long_name = "TgC/month"

   761   do n = 0, data_n-1

   762   do m = 0, n_biome-1

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

   766      wks = gsn_open_wks (plot_type,plot_name)

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

   769      res@tiMainString = title

   770      res@tiMainFontHeightF = 0.025

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

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

   776      delete (wks)

   777      delete (plot)

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

   780            "rm "+plot_name+"."+plot_type)

   781   end do

   782   end do

   784   do n = 0, data_n-1

   786      plot_name = component(n)+"_monthly_global"

   788      wks = gsn_open_wks (plot_type,plot_name)

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

   791      res@tiMainString = title

   792      res@tiMainFontHeightF = 0.025

   794      do k = 0,ntime-1

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

   796      end do

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

   800      delete (wks)

   801      delete (plot)

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

   804            "rm "+plot_name+"."+plot_type)

   805   end do

   807   delete (plot_data)

   808   delete (timeI)

   809   delete (timeF)

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

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

   813 ;*******************************************************************

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

   816   yvalues_g!0 = "time"

   817   yvalues_g!1 = "case"

   818   yvalues_g!2 = "record"

   820   yvalues_a = month_to_annual(yvalues_g,0)

   822   delete (yvalues_g)

   824 ; for x-axis in xyplot

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

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

   828   timeI = ispan(1,nyear,1)

   829   timeF = year_start + (timeI-1)

   830   timeF@long_name = "year"

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

   833   plot_data@long_name = "TgC/year"

   835   do n = 0, data_n-1

   836   do m = 0, n_biome-1

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

   840      wks = gsn_open_wks (plot_type,plot_name)

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

   843      res@tiMainString = title

   844      res@tiMainFontHeightF = 0.025

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

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

   850      delete (wks)

   851      delete (plot)

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

   854            "rm "+plot_name+"."+plot_type)

   855   end do

   856   end do

   858   do n = 0, data_n-1

   860      plot_name = component(n)+"_annual_global"

   862      wks = gsn_open_wks (plot_type,plot_name)

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

   865      res@tiMainString = title

   866      res@tiMainFontHeightF = 0.025

   868      do k = 0,nyear-1

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

   870      end do

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

   874      delete (wks)

   875      delete (plot)

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

   878            "rm "+plot_name+"."+plot_type)

   879   end do

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

   882 ; output plot and html

   883 ;****************************************

   884   output_dir = model_name+"/carbon_sink"

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

   887 ;****************************************

   889 end

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