clamp: carbon_sink/20x.table+tseries.ncl@4be95183fbcd

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

     2 ; for casa (no fire)

     3 ; using model biome vlass

4 ;

     5 ; required command line input parameters:

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

7 ;

     8 ; histogram normalized by rain and compute correleration

     9 ;**************************************************************

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

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

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

    13 ;**************************************************************

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

    15 begin

    16 ; add line to ascci/html file

    18   nnewlines = dimsizes(newlines)

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

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

    21     return

    22   end if

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

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

    25   nline = nline + nnewlines

    26   return

    27 end

    28 ;**************************************************************

    29 ; Main code.

    30 begin

    32  plot_type     = "ps"

    33  plot_type_new = "png"

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

    36 ; model name and grid

    37 ;************************************************

    39  model_grid = "T42"

    41  model_name  = "casa"

    42  model_name1 = "i01.06casa"

    43  model_name2 = "i01.10casa"

    45 ;---------------------------------------------------

    46 ; get biome data: model

    48   biome_name_mod = "Model PFT Class"

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

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

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

    55   classmod = fm->CLASS_PFT

    57   delete (fm)

    59 ; model data has 17 land-type classes

    61   nclass_mod = 17

    63 ;--------------------------------------------------

    64 ; get model data: landfrac and area

    66  dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/"

    67  film_l = "lnd_T42.nc"

    68  fm_l   = addfile (dirm_l+film_l,"r")

    70  landfrac = fm_l->landfrac

    71  area     = fm_l->area

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

    74   area = area * 1.e6

    76 ;---------------------------------------------------

    77 ; take into account landfrac

    79   area = area * landfrac

    81   delete (landfrac)

    83 ;---------------------------------------------------

    84 ; read data: model, group 1

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

    87  film = model_name1 + "_1980-2004_ANN_climo.nc"

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

    90  NPP1   = fm->NPP

    92  leafc  = fm->LEAFC

    93  woodc  = fm->WOODC

    94  frootc = fm->FROOTC

    95  VegC   = leafc

    96  VegC   = leafc + woodc + frootc

    98  litterc = fm->LITTERC

    99  cwdc    = fm->CWDC

   100  LiCwC   = litterc

   101  LiCwC   = litterc + cwdc

   103  SoilC   = fm->SOILC

   105  delete (fm)

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

   107 ; read data: model, group 2

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

   110  film = model_name2 + "_1990-2004_ANN_climo.nc"

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

   113  NPP2   = fm->NPP

   114  NEE2   = fm->NEE

   116 ;---------------------------------------------------

   117 ; Units for these variables are:

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

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

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

   123 ;VegC:  g C/m^2

   124 ;LiCwC: g C/m^2

   125 ;SoilC: g C/m^2

   127  nsec_per_year = 60*60*24*365

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

   131  NPP1 = NPP1 *  nsec_per_year

   132  NPP2 = NPP2 *  nsec_per_year

   133  NEE2 = NEE2 *  nsec_per_year

   135  data_n = 7

   137 ;*******************************************************************

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

   139 ;********************************************************************

   141 ; using model biome class

   142   nclass      = nclass_mod

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

   146 ; print (range)

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

   148 ; particular range, and then take an average.

   150   nx           = dimsizes(range) - 1

   152 ;==============================

   153 ; put data into bins

   154 ;==============================

   156 ; using observed biome class

   157 ; base  = ndtooned(classob)

   158 ; using model biome class

   159   base  = ndtooned(classmod)

   161   area_1d = ndtooned(area)

   163 ; output

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

   166   yvalues_t = new((/data_n,nx/),float)

   168 ; Loop through each range, using base.

   170   do i=0,nx-1

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

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

   174      else

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

   176      end if

   178   do n = 0,data_n-1

   180      if (n.eq.0) then

   181         data = ndtooned(area)

   182      end if

   184      if (n.eq.1) then

   185         data = ndtooned(NPP1)

   186      end if

   188      if (n.eq.2) then

   189         data = ndtooned(VegC)

   190      end if

   192      if (n.eq.3) then

   193         data = ndtooned(LiCwC)

   194      end if

   196      if (n.eq.4) then

   197         data = ndtooned(SoilC)

   198      end if

   200      if (n.eq.5) then

   201         data = ndtooned(NPP2)

   202      end if

   204      if (n.eq.6) then

   205         data = ndtooned(NEE2)

   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)

   252 ;----------------------------------------------------------------

   253 ; data for table1

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

   256 ;print (good)

   258  area_g  = yvalues(0,good)

   259  NPP1_g  = yvalues(1,good)

   260  VegC_g  = yvalues(2,good)

   261  LiCwC_g = yvalues(3,good)

   262  SoilC_g = yvalues(4,good)

   263  NPP2_g  = yvalues(5,good)

   264  NEE2_g  = yvalues(6,good)

   266  n_biome = dimsizes(NPP1_g)

   268  NPP_ratio = NPP2_g/NPP1_g

   270 ;-----------------------------------------------------------------

   271 ; data for table2

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

   274  factor_unit = 1.e-15

   276  NPP1_t  = yvalues_t(1,good) * factor_unit

   277  VegC_t  = yvalues_t(2,good) * factor_unit

   278  LiCwC_t = yvalues_t(3,good) * factor_unit

   279  SoilC_t = yvalues_t(4,good) * factor_unit

   280  NEE2_t  = yvalues_t(6,good) * factor_unit

   282  delete (yvalues)

   283  delete (yvalues_t)

   285 ;-------------------------------------------------------------

   286 ; html table1 data

   288 ; column (not including header column)

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

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

   292                ,"VegC (gC/m2)" \

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

   294                ,"SoilC (gC/m2)" \

   295                ,"NPP_ratio" \

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

   297                /)

   299   ncol = dimsizes(col_head)

   301 ; row (not including header row)

   303 ; using model biome class:

   304   row_head  = (/"Not Vegetated" \

   305                ,"Needleleaf Evergreen Temperate Tree" \

   306                ,"Needleleaf Evergreen Boreal Tree" \

   307 ;              ,"Needleleaf Deciduous Boreal Tree" \

   308                ,"Broadleaf Evergreen Tropical Tree" \

   309                ,"Broadleaf Evergreen Temperate Tree" \

   310                ,"Broadleaf Deciduous Tropical Tree" \

   311                ,"Broadleaf Deciduous Temperate Tree" \

   312 ;              ,"Broadleaf Deciduous Boreal Tree" \

   313 ;              ,"Broadleaf Evergreen Shrub" \

   314                ,"Broadleaf Deciduous Temperate Shrub" \

   315                ,"Broadleaf Deciduous Boreal Shrub" \

   316                ,"C3 Arctic Grass" \

   317                ,"C3 Non-Arctic Grass" \

   318                ,"C4 Grass" \

   319                ,"Corn" \

   320 ;              ,"Wheat" \

   321                ,"All Biome" \

   322                /)

   323   nrow = dimsizes(row_head)

   325 ; arrays to be passed to table.

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

   328  do i=0,nrow-2

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

   330   text4(i,1) = sprintf("%.1f",NPP1_g(i))

   331   text4(i,2) = sprintf("%.1f",VegC_g(i))

   332   text4(i,3) = sprintf("%.1f",LiCwC_g(i))

   333   text4(i,4) = sprintf("%.1f",SoilC_g(i))

   334   text4(i,5) = sprintf("%.2f",NPP_ratio(i))

   335   text4(i,6) = sprintf("%.1f",NEE2_g(i))

   336  end do

   338 ;-------------------------------------------------------

   339 ; create html table1

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

   343   header = (/"<HTML>" \

   344             ,"<HEAD>" \

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

   346             ,"</HEAD>" \

   347             ,header_text \

   348             /)

   349   footer = "</HTML>"

   351   table_header = (/ \

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

   353        ,"<tr>" \

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

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

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

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

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

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

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

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

   362        ,"</tr>" \

   363        /)

   364   table_footer = "</table>"

   365   row_header = "<tr>"

   366   row_footer = "</tr>"

   368   lines = new(50000,string)

   369   nline = 0

   371   set_line(lines,nline,header)

   372   set_line(lines,nline,table_header)

   374 ;----------------------------

   375 ;row of table

   377   do n = 0,nrow-2

   378      set_line(lines,nline,row_header)

   380      txt0  = row_head(n)

   381      txt1  = text4(n,0)

   382      txt2  = text4(n,1)

   383      txt3  = text4(n,2)

   384      txt4  = text4(n,3)

   385      txt5  = text4(n,4)

   386      txt6  = text4(n,5)

   387      txt7  = text4(n,6)

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

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

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

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

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

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

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

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

   398      set_line(lines,nline,row_footer)

   399   end do

   400 ;----------------------------

   401   set_line(lines,nline,table_footer)

   402   set_line(lines,nline,footer)

   404 ; Now write to an HTML file.

   406   output_html = "table_carbon_sink1.html"

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

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

   410     asciiwrite(output_html,lines(idx))

   411   else

   412    print ("error?")

   413   end if

   415   delete (idx)

   417   delete (col_head)

   418   delete (row_head)

   419   delete (text4)

   420   delete (table_header)

   422 ;-----------------------------------------------------------------

   423 ; html table2 data

   425 ; column (not including header column)

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

   428                ,"VegC (PgC)" \

   429                ,"Litter+CWD (PgC)" \

   430                ,"SoilC (PgC)" \

   431                ,"NEE (PgC/yr)" \

   432                ,"NPP timeseries" \

   433                ,"NEE timeseries" \

   434                ,"Fire timeseries" \

   435                /)

   437   ncol = dimsizes(col_head)

   439 ; row (not including header row)

   441 ; using model biome class:

   442   row_head  = (/"Not Vegetated" \

   443                ,"Needleleaf Evergreen Temperate Tree" \

   444                ,"Needleleaf Evergreen Boreal Tree" \

   445 ;              ,"Needleleaf Deciduous Boreal Tree" \

   446                ,"Broadleaf Evergreen Tropical Tree" \

   447                ,"Broadleaf Evergreen Temperate Tree" \

   448                ,"Broadleaf Deciduous Tropical Tree" \

   449                ,"Broadleaf Deciduous Temperate Tree" \

   450 ;              ,"Broadleaf Deciduous Boreal Tree" \

   451 ;              ,"Broadleaf Evergreen Shrub" \

   452                ,"Broadleaf Deciduous Temperate Shrub" \

   453                ,"Broadleaf Deciduous Boreal Shrub" \

   454                ,"C3 Arctic Grass" \

   455                ,"C3 Non-Arctic Grass" \

   456                ,"C4 Grass" \

   457                ,"Corn" \

   458 ;              ,"Wheat" \

   459                ,"All Biome" \

   460                /)

   461   nrow = dimsizes(row_head)

   463 ; arrays to be passed to table.

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

   466  do i=0,nrow-2

   467   text4(i,0) = sprintf("%.1f",NPP1_t(i))

   468   text4(i,1) = sprintf("%.1f",VegC_t(i))

   469   text4(i,2) = sprintf("%.1f",LiCwC_t(i))

   470   text4(i,3) = sprintf("%.1f",SoilC_t(i))

   471   text4(i,4) = sprintf("%.1f",NEE2_t(i))

   472   text4(i,5) = "<a href=./NPP_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NPP_annual_biome_"+i+".png>annual_plot</a>"

   473   text4(i,6) = "<a href=./NEE_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NEE_annual_biome_"+i+".png>annual_plot</a>"

   474   text4(i,7) = "--"

   475  end do

   476   text4(nrow-1,0) = sprintf("%.1f",sum(NPP1_t))

   477   text4(nrow-1,1) = sprintf("%.1f",sum(VegC_t))

   478   text4(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t))

   479   text4(nrow-1,3) = sprintf("%.1f",sum(SoilC_t))

   480   text4(nrow-1,4) = sprintf("%.1f",sum(NEE2_t))

   481   text4(nrow-1,5) = "<a href=./NPP_monthly_global.png>monthly_plot</a> <br> <a href=./NPP_annual_global.png>annual_plot</a>"

   482   text4(nrow-1,6) = "<a href=./NEE_monthly_global.png>monthly_plot</a> <br> <a href=./NEE_annual_global.png>annual_plot</a>"

   483   text4(nrow-1,7) = "--"

   485 ;**************************************************

   486 ; create html table2

   487 ;**************************************************

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

   491   header = (/"<HTML>" \

   492             ,"<HEAD>" \

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

   494             ,"</HEAD>" \

   495             ,header_text \

   496             /)

   497   footer = "</HTML>"

   499   table_header = (/ \

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

   501        ,"<tr>" \

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

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

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

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

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

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

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

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

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

   511        ,"</tr>" \

   512        /)

   513   table_footer = "</table>"

   514   row_header = "<tr>"

   515   row_footer = "</tr>"

   517   lines = new(50000,string)

   518   nline = 0

   520   set_line(lines,nline,header)

   521   set_line(lines,nline,table_header)

   522 ;-----------------------------------------------

   523 ;row of table

   525   do n = 0,nrow-1

   526      set_line(lines,nline,row_header)

   528      txt0  = row_head(n)

   529      txt1  = text4(n,0)

   530      txt2  = text4(n,1)

   531      txt3  = text4(n,2)

   532      txt4  = text4(n,3)

   533      txt5  = text4(n,4)

   534      txt6  = text4(n,5)

   535      txt7  = text4(n,6)

   536      txt8  = text4(n,7)

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

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

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

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

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

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

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

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

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

   548      set_line(lines,nline,row_footer)

   549   end do

   550 ;-----------------------------------------------

   551   set_line(lines,nline,table_footer)

   552   set_line(lines,nline,footer)

   554 ; Now write to an HTML file.

   556   output_html = "table_carbon_sink2.html"

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

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

   560     asciiwrite(output_html,lines(idx))

   561   else

   562    print ("error?")

   563   end if

   565   delete (idx)

   567 ;---------------------------------------------------

   568 ; read data: model, time series

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

   571  film = model_name2 + "_Fire_C_1979-2004_monthly.nc"

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

   574  NPP3 = fm->NPP

   575  NEE3 = fm->NEE

   576 ;Fire = fm->COL_FIRE_CLOSS

   578  delete (fm)

   580 ; Units for these variables are:

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

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

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

   586  nsec_per_month = 60*60*24*30

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

   590  NPP3 = NPP3 * nsec_per_month

   591  NEE3 = NEE3 * nsec_per_month

   592 ;Fire = Fire * nsec_per_month

   594 ;data_n = 3

   595  data_n = 2

   597  dsizes = dimsizes(NPP3)

   598  nyear  = dsizes(0)

   599  nmonth = dsizes(1)

   600  ntime  = nyear * nmonth

   602  year_start = 1979

   603  year_end   = 2004

   605 ;*******************************************************************

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

   607 ;********************************************************************

   609 ; using model biome class

   610   nclass = nclass_mod

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

   614 ; print (range)

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

   616 ; particular range, and then take an average.

   618   nx = dimsizes(range) - 1

   620 ;==============================

   621 ; put data into bins

   622 ;==============================

   624 ; using observed biome class

   625 ; base  = ndtooned(classob)

   626 ; using model biome class

   627   base  = ndtooned(classmod)

   629 ; output

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

   633 ; Loop through each range, using base.

   635   do i=0,nx-1

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

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

   639      else

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

   641      end if

   643   do n = 0,data_n-1

   645      t = -1

   646      do m = 0,nyear-1

   647      do k = 0,nmonth-1

   649         t = t + 1

   651         if (n.eq.0) then

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

   653         end if

   655         if (n.eq.1) then

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

   657         end if

   659 ;       if (n.eq.2) then

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

   661 ;       end if

   663 ;       Calculate average

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

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

   667         else

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

   669         end if

   671 ;#############################################################

   672 ; using model biome class:

   673 ;     set the following 4 classes to _FillValue:

   674 ;     (3)Needleleaf Deciduous Boreal Tree,

   675 ;     (8)Broadleaf Deciduous Boreal Tree,

   676 ;     (9)Broadleaf Evergreen Shrub,

   677 ;     (16)Wheat

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

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

   681         end if

   682 ;#############################################################

   684      end do

   685      end do

   687      delete(data)

   688   end do

   690     delete(idx)

   691   end do

   693   delete (base)

   694   delete (NPP3)

   695   delete (NEE3)

   696 ; delete (Fire)

   698 ;----------------------------------------------------------------

   699 ; data for tseries plot

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

   703   yvalues_g@units = "TgC/month"

   705 ; change unit to Tg C/month

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

   707   factor_unit = 1.e-12

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

   711 ;*******************************************************************

   712 ; general settings for line plot

   713 ;*******************************************************************

   715 ; res

   716   res                   = True

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

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

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

   721   res@trXMinF   = year_start

   722   res@trXMaxF   = year_end + 1

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

   725 ; res@vpWidthF  = 0.8

   726   res@vpWidthF  = 0.75

   728 ; res@gsnMaximize = True

   730 ;*******************************************************************

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

   732 ;*******************************************************************

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

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

   737 ; for x-axis in xyplot

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

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

   741   timeI = ispan(1,ntime,1)

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

   743   timeF@long_name = "year"

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

   746   plot_data@long_name = "TgC/month"

   748   do n = 0, data_n-1

   749   do m = 0, n_biome-1

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

   753      wks = gsn_open_wks (plot_type,plot_name)

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

   756      res@tiMainString = title

   757      res@tiMainFontHeightF = 0.025

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

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

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

   764            "rm "+plot_name+"."+plot_type)

   766      clear (wks)

   767      delete (plot)

   768   end do

   769   end do

   771   do n = 0, data_n-1

   773      plot_name = component(n)+"_monthly_global"

   775      wks = gsn_open_wks (plot_type,plot_name)

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

   778      res@tiMainString = title

   779      res@tiMainFontHeightF = 0.025

   781      do k = 0,ntime-1

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

   783      end do

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

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

   788            "rm "+plot_name+"."+plot_type)

   790      clear (wks)

   791      delete (plot)

   792   end do

   794   delete (plot_data)

   795   delete (timeI)

   796   delete (timeF)

   798 ;*******************************************************************

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

   800 ;*******************************************************************

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

   803   yvalues_g!0 = "time"

   804   yvalues_g!1 = "case"

   805   yvalues_g!2 = "record"

   807   yvalues_a = month_to_annual(yvalues_g,0)

   809   delete (yvalues_g)

   811 ; for x-axis in xyplot

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

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

   815   timeI = ispan(1,nyear,1)

   816   timeF = year_start + (timeI-1)

   817   timeF@long_name = "year"

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

   820   plot_data@long_name = "TgC/year"

   822   do n = 0, data_n-1

   823   do m = 0, n_biome-1

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

   827      wks = gsn_open_wks (plot_type,plot_name)

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

   830      res@tiMainString = title

   831      res@tiMainFontHeightF = 0.025

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

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

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

   838            "rm "+plot_name+"."+plot_type)

   840      clear (wks)

   841      delete (plot)

   842   end do

   843   end do

   845   do n = 0, data_n-1

   847      plot_name = component(n)+"_annual_global"

   849      wks = gsn_open_wks (plot_type,plot_name)

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

   852      res@tiMainString = title

   853      res@tiMainFontHeightF = 0.025

   855      do k = 0,nyear-1

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

   857      end do

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

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

   862            "rm "+plot_name+"."+plot_type)

   864      clear (wks)

   865      delete (plot)

   866   end do

   868 end

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