clamp: carbon_sink/20.table+tseries.ncl@0c6405ab2ff4

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

     2 ;using model biome vlass

3 ;

     4 ; required command line input parameters:

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

6 ;

     7 ; histogram normalized by rain and compute correleration

     8 ;**************************************************************

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

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

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

    12 ;**************************************************************

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

    14 begin

    15 ; add line to ascci/html file

    17   nnewlines = dimsizes(newlines)

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

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

    20     return

    21   end if

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

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

    24   nline = nline + nnewlines

    25   return

    26 end

    27 ;**************************************************************

    28 ; Main code.

    29 begin

    31  plot_type     = "ps"

    32  plot_type_new = "png"

    34 ;************************************************

    35 ; model name and grid

    36 ;************************************************

    38  model_grid = "T42"

    40  model_name  = "cn"

    41  model_name1 = "i01.06cn"

    42  model_name2 = "i01.10cn"

    44 ;---------------------------------------------------

    45 ; get biome data: model

    47   biome_name_mod = "Model PFT Class"

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

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

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

    54   classmod = fm->CLASS_PFT

    56   delete (fm)

    58 ; model data has 17 land-type classes

    60   nclass_mod = 17

    62 ;--------------------------------------------------

    63 ; get model data: landfrac and area

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

    66  film_l = "lnd_T42.nc"

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

    69  landfrac = fm_l->landfrac

    70  area     = fm_l->area

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

    73   area = area * 1.e6

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

    76 ; take into account landfrac

    78   area = area * landfrac

    80   delete (landfrac)

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

    83 ; read data: model, group 1

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

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

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

    89  NPP1   = fm->NPP

    91  leafc  = fm->LEAFC

    92  woodc  = fm->WOODC

    93  frootc = fm->FROOTC

    94  VegC   = leafc

    95  VegC   = leafc + woodc + frootc

    97  litterc = fm->LITTERC

    98  cwdc    = fm->CWDC

    99  LiCwC   = litterc

   100  LiCwC   = litterc + cwdc

   102  SoilC   = fm->SOILC

   104  delete (fm)

   105 ;---------------------------------------------------

   106 ; read data: model, group 2

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

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

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

   112  NPP2   = fm->NPP

   113  NEE2   = fm->NEE

   115 ;---------------------------------------------------

   116 ; Units for these variables are:

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

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

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

   122 ;VegC:  g C/m^2

   123 ;LiCwC: g C/m^2

   124 ;SoilC: g C/m^2

   126  nsec_per_year = 60*60*24*365

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

   130  NPP1 = NPP1 *  nsec_per_year

   131  NPP2 = NPP2 *  nsec_per_year

   132  NEE2 = NEE2 *  nsec_per_year

   134  data_n = 7

   136 ;*******************************************************************

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

   138 ;********************************************************************

   140 ; using model biome class

   141   nclass      = nclass_mod

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

   145 ; print (range)

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

   147 ; particular range, and then take an average.

   149   nx           = dimsizes(range) - 1

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

   152 ; put data into bins

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

   155 ; using observed biome class

   156 ; base  = ndtooned(classob)

   157 ; using model biome class

   158   base  = ndtooned(classmod)

   160   area_1d = ndtooned(area)

   162 ; output

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

   165   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) = "<a href=./Fire_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./Fire_annual_biome_"+i+".png>annual_plot</a>"

   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) = "<a href=./Fire_monthly_global.png>monthly_plot</a> <br> <a href=./Fire_annual_global.png>annual_plot</a>"

   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

   596  dsizes = dimsizes(NPP3)

   597  nyear  = dsizes(0)

   598  nmonth = dsizes(1)

   599  ntime  = nyear * nmonth

   601  year_start = 1979

   602  year_end   = 2004

   604 ;*******************************************************************

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

   606 ;********************************************************************

   608 ; using model biome class

   609   nclass = nclass_mod

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

   613 ; print (range)

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

   615 ; particular range, and then take an average.

   617   nx = dimsizes(range) - 1

   619 ;==============================

   620 ; put data into bins

   621 ;==============================

   623 ; using observed biome class

   624 ; base  = ndtooned(classob)

   625 ; using model biome class

   626   base  = ndtooned(classmod)

   628 ; output

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

   632 ; Loop through each range, using base.

   634   do i=0,nx-1

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

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

   638      else

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

   640      end if

   642   do n = 0,data_n-1

   644      t = -1

   645      do m = 0,nyear-1

   646      do k = 0,nmonth-1

   648         t = t + 1

   650         if (n.eq.0) then

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

   652         end if

   654         if (n.eq.1) then

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

   656         end if

   658         if (n.eq.2) then

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

   660         end if

   662 ;       Calculate average

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

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

   666         else

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

   668         end if

   670 ;#############################################################

   671 ; using model biome class:

   672 ;     set the following 4 classes to _FillValue:

   673 ;     (3)Needleleaf Deciduous Boreal Tree,

   674 ;     (8)Broadleaf Deciduous Boreal Tree,

   675 ;     (9)Broadleaf Evergreen Shrub,

   676 ;     (16)Wheat

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

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

   680         end if

   681 ;#############################################################

   683      end do

   684      end do

   686      delete(data)

   687   end do

   689     delete(idx)

   690   end do

   692   delete (base)

   693   delete (NPP3)

   694   delete (NEE3)

   695   delete (Fire)

   697 ;----------------------------------------------------------------

   698 ; data for tseries plot

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

   702   yvalues_g@units = "TgC/month"

   704 ; change unit to Tg C/month

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

   706   factor_unit = 1.e-12

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

   710 ;*******************************************************************

   711 ; general settings for line plot

   712 ;*******************************************************************

   714 ; res

   715   res                   = True

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

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

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

   720   res@trXMinF   = year_start

   721   res@trXMaxF   = year_end + 1

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

   724 ; res@vpWidthF  = 0.8

   725   res@vpWidthF  = 0.75

   727 ; res@gsnMaximize = True

   729 ;*******************************************************************

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

   731 ;*******************************************************************

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

   735 ; for x-axis in xyplot

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

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

   739   timeI = ispan(1,ntime,1)

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

   741   timeF@long_name = "year"

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

   744   plot_data@long_name = "TgC/month"

   746   do n = 0, data_n-1

   747   do m = 0, n_biome-1

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

   751      wks = gsn_open_wks (plot_type,plot_name)

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

   754      res@tiMainString = title

   755      res@tiMainFontHeightF = 0.025

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

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

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

   762            "rm "+plot_name+"."+plot_type)

   764      clear (wks)

   765      delete (plot)

   766   end do

   767   end do

   769   do n = 0, data_n-1

   771      plot_name = component(n)+"_monthly_global"

   773      wks = gsn_open_wks (plot_type,plot_name)

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

   776      res@tiMainString = title

   777      res@tiMainFontHeightF = 0.025

   779      do k = 0,ntime-1

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

   781      end do

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

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

   786            "rm "+plot_name+"."+plot_type)

   788      clear (wks)

   789      delete (plot)

   790   end do

   792   delete (plot_data)

   793   delete (timeI)

   794   delete (timeF)

   796 ;*******************************************************************

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

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

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

   801   yvalues_g!0 = "time"

   802   yvalues_g!1 = "case"

   803   yvalues_g!2 = "record"

   805   yvalues_a = month_to_annual(yvalues_g,0)

   807   delete (yvalues_g)

   809 ; for x-axis in xyplot

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

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

   813   timeI = ispan(1,nyear,1)

   814   timeF = year_start + (timeI-1)

   815   timeF@long_name = "year"

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

   818   plot_data@long_name = "TgC/year"

   820   do n = 0, data_n-1

   821   do m = 0, n_biome-1

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

   825      wks = gsn_open_wks (plot_type,plot_name)

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

   828      res@tiMainString = title

   829      res@tiMainFontHeightF = 0.025

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

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

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

   836            "rm "+plot_name+"."+plot_type)

   838      clear (wks)

   839      delete (plot)

   840   end do

   841   end do

   843   do n = 0, data_n-1

   845      plot_name = component(n)+"_annual_global"

   847      wks = gsn_open_wks (plot_type,plot_name)

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

   850      res@tiMainString = title

   851      res@tiMainFontHeightF = 0.025

   853      do k = 0,nyear-1

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

   855      end do

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

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

   860            "rm "+plot_name+"."+plot_type)

   862      clear (wks)

   863      delete (plot)

   864   end do

   866 end

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