Added simulation of soil carbon storage and fluxes. This formulation was

taken mostly from the LPJ model (Sitch, 2003), which in turn used a
Lloyd-Taylor model for the dependence of soil respiration on soil
temperature.  The dependence of soil respiration on soil moisture was
based on the formulation of Yi et al (2012) but modified to allow a small
respiration rate under saturated conditions.

At this point, we do not simulate the storage of carbon in living biomass.
Therefore, the flux of carbon into the soil (litterfall) is set equal to
the total NPP of the previous calendar year, spread evenly over the current
year.  As in the LPJ model, soil carbon is stored in 3 pools: litter (fast),
intermediate, and slow; with associated turnover times of 2.86 y, 33.3 y, and
1,000 y, respectively.  Litterfall enters the litter pool.  Carbon exits the
litter pool through respiration (RhLitter).  A fraction (fAir) of this
respired carbon is in the form of CO2 and is vented directly to the atmosphere
(RhLitter2Atm).  The remainder is sent to the intermediate and slow pools
in the proportions fInter and (1-fInter), respectively.  These pools also
respire carbon, which is assumed to be in the form of CO2 and vented directly
to the atmosphere.

There are several new output variables associated with this feature:
  OUT_RHET: Total heterotrophic respiration vented to the atmosphere
            (= RhLitter2Atm+RhInter+RhSlow)  [g C/m2d]
  OUT_NEE:  Net Ecosystem Exchange (= NPP-RHET) [g C/m2d]
  OUT_LITTERFALL: Flux of carbon from living biomass into litter pool [g C/m2d]
  OUT_CLITTER: Carbon density in the litter pool [g C/m2]
  OUT_CINTER: Carbon density in the intermediate pool [g C/m2]
  OUT_CSLOW: Carbon density in the slow pool [g C/m2]

This feature is part of the carbon cycle, controlled by the setting of the
CARBON option in the global parameter file.

Author: tbohn

src/LAKE.h

@@ -26,6 +26,7 @@
   2010-Dec-28 Added latitude to alblake() arglist.			TJB
   2011-Mar-01 Added rescale_snow_storage().  Added terms to argument
 	      list of initialize_lake().				TJB
+  2013-Jul-25 Added advect_carbon_storage().				TJB
 ******************************************************************************/
 
 //#ifndef LAKE_SET
@@ -70,6 +71,7 @@
 
 double adjflux(double, double, double ,double, double, double, double,
 	       double, double, double, double *, double *);
+void advect_carbon_storage(double, double, lake_var_struct *, cell_data_struct *);
 void advect_soil_veg_storage(double, double, double, double *, soil_con_struct *, veg_con_struct *, cell_data_struct *, veg_var_struct *, lake_con_struct);
 void advect_snow_storage(double, double, double, snow_data_struct *);
 void alblake(double, double, double *, double *, float *, float *, double, double, 

src/Makefile

@@ -40,6 +40,7 @@
 # 2013-Jul-25 Added compute_coszen.c.						TJB
 # 2013-Jul-25 Added calc_Nscale_factors.c, canopy_assimilation.c, faparl.c,
 # 	      photosynth.c.							TJB
+# 2013-Jul-25 Added compute_soil_resp.c and soil_carbon_balance.c.		TJB
 #
 # $Id$
 #
@@ -86,7 +87,7 @@ OBJS =  CalcAerodynamic.o CalcBlowingSnow.o SnowPackEnergyBalance.o \
 	calc_water_energy_balance_errors.o canopy_assimilation.o canopy_evap.o \
 	check_files.o check_state_file.o close_files.o cmd_proc.o \
 	compress_files.o compute_coszen.o compute_dz.o compute_pot_evap.o \
-	compute_treeline.o compute_zwt.o correct_precip.o \
+	compute_soil_resp.o compute_treeline.o compute_zwt.o correct_precip.o \
 	display_current_settings.o dist_prec.o estimate_T1.o faparl.o free_dist_prcp.o \
 	free_vegcon.o frozen_soil.o full_energy.o func_atmos_energy_bal.o \
 	func_atmos_moist_bal.o func_canopy_energy_bal.o \
@@ -104,7 +105,7 @@ OBJS =  CalcAerodynamic.o CalcBlowingSnow.o SnowPackEnergyBalance.o \
 	read_snowband.o read_soilparam.o read_soilparam_arc.o read_veglib.o \
 	read_vegparam.o redistribute_during_storm.o root_brent.o runoff.o \
 	set_output_defaults.o snow_intercept.o snow_melt.o \
-	snow_utility.o soil_conduction.o \
+	snow_utility.o soil_carbon_balance.o soil_conduction.o \
 	soil_thermal_eqn.o solve_snow.o \
 	surface_fluxes.o svp.o vicNl.o vicerror.o \
 	write_data.o write_forcing_file.o write_header.o write_layer.o \

src/README.txt

@@ -145,6 +145,64 @@ Added simulation of photosynthesis.
 
 
 
+Added simulation of soil carbon storage and fluxes.
+
+	Files Affected:
+
+	compute_soil_resp.c (new)
+	full_energy.c
+	initialize_lake.c
+	initialize_soil.c
+	initialize_veg.c
+	LAKE.h
+	lakes.eb.c
+	Makefile
+	output_list_utils.c
+	put_data.c
+	read_initial_model_state.c
+	soil_carbon_balance.c (new)
+	surface_fluxes.c
+	vicNl_def.h
+	vicNl.h
+	write_model_state.c
+
+	Description:
+
+	Added simulation of soil carbon storage and fluxes.  This formulation was
+	taken mostly from the LPJ model (Sitch, 2003), which in turn used a
+	Lloyd-Taylor model for the dependence of soil respiration on soil
+	temperature.  The dependence of soil respiration on soil moisture was
+	based on the formulation of Yi et al (2012) but modified to allow a small
+	respiration rate under saturated conditions.
+
+	At this point, we do not simulate the storage of carbon in living biomass.
+	Therefore, the flux of carbon into the soil (litterfall) is set equal to
+	the total NPP of the previous calendar year, spread evenly over the current
+	year.  As in the LPJ model, soil carbon is stored in 3 pools: litter (fast),
+	intermediate, and slow; with associated turnover times of 2.86 y, 33.3 y, and
+	1,000 y, respectively.  Litterfall enters the litter pool.  Carbon exits the
+	litter pool through respiration (RhLitter).  A fraction (fAir) of this
+	respired carbon is in the form of CO2 and is vented directly to the atmosphere
+	(RhLitter2Atm).  The remainder is sent to the intermediate and slow pools
+	in the proportions fInter and (1-fInter), respectively.  These pools also
+	respire carbon, which is assumed to be in the form of CO2 and vented directly
+	to the atmosphere.
+
+	There are several new output variables associated with this feature:
+	  OUT_RHET: Total heterotrophic respiration vented to the atmosphere
+	            (= RhLitter2Atm+RhInter+RhSlow)  [g C/m2d]
+	  OUT_NEE:  Net Ecosystem Exchange (= NPP-RHET) [g C/m2d]
+	  OUT_LITTERFALL: Flux of carbon from living biomass into litter pool [g C/m2d]
+	  OUT_CLITTER: Carbon density in the litter pool [g C/m2]
+	  OUT_CINTER: Carbon density in the intermediate pool [g C/m2]
+	  OUT_CSLOW: Carbon density in the slow pool [g C/m2]
+
+	This feature is part of the carbon cycle, controlled by the setting of the
+	CARBON option in the global parameter file.
+
+
+
+
 Bug Fixes:
 ----------
 

src/compute_soil_resp.c

@@ -0,0 +1,141 @@
+/********************************************************************************
+  filename  : compute_soil_resp.c
+  purpose   : Calculate soil respiration (heterotrophic respiration, or Rh)
+  interface : - input :
+                - Nnodes:  number of points (vertically) at which to evaluate soil respiration
+                - dZ[]:    array of layer thicknesses [m] pertaining to nodes
+                - dZTot:   total depth [m] of all layers (total depth of soil containing carbon)
+                - dt:      timestep length [h]
+                - T[]:     array of node temperatures [K]
+                - w[]:     array of node moisture fractions [fraction]
+                - CLitter: carbon storage in litter pool [gC/m2]
+                - CInter:  carbon storage in intermediate pool [gC/m2]
+                - CSlow:   carbon storage in slow pool [gC/m2]
+              Note: the litter pool is concentrated at the soil surface, while the intermediate and
+                    slow pools are both assumed to be distributed uniformly throughout soil column.
+
+              - constants:
+                - E0_LT:     Lloyd-Taylor E0 parameter [K]
+                - T0_LT:     Lloyd-Taylor E0 parameter [K]
+                - wminFM:    minimum soil moisture (fraction) at which soil respiration can occur
+                - wmaxFM:    maximum soil moisture (fraction) at which soil respiration can occur
+                - woptFM:    soil moisture (fraction) at which maximum soil respiration occurs
+                - Rhsat:     ratio of soil respiration rate under saturated conditions (w=wmaxFM) to
+                             that under optimal conditions (w=woptFM)
+                - Rfactor:   scaling factor to account for other (non-moisture) sources of inhibition
+                             of respiration
+                - tauLitter: Litter pool turnover time [y]
+                - tauInter:  Intermediate pool turnover time [y]
+                - tauSlow:   Slow pool turnover time [y]
+
+              - output:
+                - RhLitter:     Soil respiration from litter pool [gC/m2]
+                - RhInterTot:   Soil respiration from intermediate pool [gC/m2]
+                - RhSlowTot:    Soil respiration from slow pool [gC/m2]
+
+  programmer: Ted Bohn
+  date      : July 25, 2013
+  changes   :
+  references: 
+********************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vicNl.h>
+
+static char vcid[] = "$Id: $";
+
+void compute_soil_resp(int Nnodes,
+                       double *dZ,
+                       double dZTot,
+                       double dt,
+                       double *T,
+                       double *w,
+                       double CLitter,
+                       double CInter,
+                       double CSlow,
+                       double *RhLitter,
+                       double *RhInterTot,
+                       double *RhSlowTot)
+{
+  extern option_struct options;
+  int i;
+  double Tref;
+  double *TK;
+  double fTLitter;
+  double *fTSoil;
+  double fMLitter;
+  double *fMSoil;
+  double *CInterNode;
+  double *CSlowNode;
+  double *RhInter;
+  double *RhSlow;
+
+  /* Allocate temp arrays */
+  TK = (double*)calloc(Nnodes,sizeof(double));
+  fTSoil = (double*)calloc(Nnodes,sizeof(double));
+  fMSoil = (double*)calloc(Nnodes,sizeof(double));
+  CInterNode = (double*)calloc(Nnodes,sizeof(double));
+  CSlowNode = (double*)calloc(Nnodes,sizeof(double));
+  RhInter = (double*)calloc(Nnodes,sizeof(double));
+  RhSlow = (double*)calloc(Nnodes,sizeof(double));
+
+  /* Compute Lloyd-Taylor temperature dependence */
+  Tref = 10+KELVIN; /* reference temperature of 10 C */
+  for (i=0; i<Nnodes; i++) {
+    TK[i] = T[i]+KELVIN;
+    if (TK[i] < T0_LT) TK[i] = T0_LT;
+  }
+  fTLitter = exp( E0_LT * ( 1/(Tref-T0_LT) - 1/(TK[0]-T0_LT) ) );
+  for (i=0; i<Nnodes; i++) {
+    fTSoil[i] = exp( E0_LT * ( 1/(Tref-T0_LT) - 1/(TK[i]-T0_LT) ) );
+  }
+
+  /* Compute moisture dependence */
+  for (i=0; i<Nnodes; i++) {
+    if (w[i] < wminFM) w[i] = wminFM;
+    if (w[i] > wmaxFM) w[i] = wmaxFM;
+  }
+  if (w[0] <= woptFM)
+   fMLitter = (w[0]-wminFM)*(w[0]-wmaxFM)/((w[0]-wminFM)*(w[0]-wmaxFM)-(w[0]-woptFM)*(w[0]-woptFM));
+  else
+   fMLitter = Rhsat + (1-Rhsat)*(w[0]-wminFM)*(w[0]-wmaxFM)/((w[0]-wminFM)*(w[0]-wmaxFM)-(w[0]-woptFM)*(w[0]-woptFM));
+  if (fMLitter > 1.0) fMLitter = 1.0;
+  if (fMLitter < 0.0) fMLitter = 0.0;
+  for (i=0; i<Nnodes; i++) {
+    if (w[i] <= woptFM)
+     fMSoil[i] = (w[i]-wminFM)*(w[i]-wmaxFM)/((w[i]-wminFM)*(w[i]-wmaxFM)-(w[i]-woptFM)*(w[i]-woptFM));
+    else
+     fMSoil[i] = Rhsat + (1-Rhsat)*(w[i]-wminFM)*(w[i]-wmaxFM)/((w[i]-wminFM)*(w[i]-wmaxFM)-(w[i]-woptFM)*(w[i]-woptFM));
+    if (fMSoil[i] > 1.0) fMSoil[i] = 1.0;
+    if (fMSoil[i] < 0.0) fMSoil[i] = 0.0;
+  }
+
+  /* Compute C per node */
+  for (i=0; i<Nnodes; i++) {
+    CInterNode[i] = CInter * dZ[i]/dZTot;
+    CSlowNode[i] = CSlow * dZ[i]/dZTot;
+  }
+
+  /* Compute Rh for various pools, nodes; C fluxes in [gC/m2d] */
+  *RhLitter = Rfactor*(fTLitter*fMLitter/(tauLitter*365.25*24/dt))*CLitter;
+  *RhInterTot = 0;
+  *RhSlowTot = 0;
+  for (i=0; i<Nnodes; i++) {
+    RhInter[i] = Rfactor*(fTSoil[i]*fMSoil[i]/(tauInter*365.25*24/dt))*CInterNode[i];
+    RhSlow[i] = Rfactor*(fTSoil[i]*fMSoil[i]/(tauSlow*365.25*24/dt))*CSlowNode[i];
+    *RhInterTot += RhInter[i];
+    *RhSlowTot += RhSlow[i];
+  }
+
+  free((char *)TK);
+  free((char *)fTSoil);
+  free((char *)fMSoil);
+  free((char *)CInterNode);
+  free((char *)CSlowNode);
+  free((char *)RhInter);
+  free((char *)RhSlow);
+
+}
+

src/full_energy.c

@@ -107,6 +107,7 @@ int  full_energy(char                 NEWCELL,
   2012-Aug-28 Added accumulation of rain and snow over lake to grid cell
 	      totals.								TJB
   2013-Jul-25 Added photosynthesis terms.					TJB
+  2013-Jul-25 Added soil carbon terms.						TJB
 
 **********************************************************************/
 {
@@ -407,11 +408,7 @@ int  full_energy(char                 NEWCELL,
               veg_var[dist][iveg][band].rsLayer[cidx] = HUGE_RESIST;
             }
             veg_var[dist][iveg][band].aPAR = 0;
-          }
-        }
-        if (dmy[rec].hour == 0) {
-          for (dist=0; dist<Ndist; dist++) {
-	    for(band=0; band<Nbands; band++) {
+            if (dmy->hour == 0) {
               calc_Nscale_factors(veg_lib[veg_class].NscaleFlag,
                                   veg_con[iveg].CanopLayerBnd,
                                   veg_lib[veg_class].LAI[dmy[rec].month-1],
@@ -421,6 +418,10 @@ int  full_energy(char                 NEWCELL,
                                   dmy[rec],
                                   veg_var[dist][iveg][band].NscaleFactor);
             }
+            if (dmy[rec].month == 1 && dmy[rec].day == 1) {
+              veg_var[dist][iveg][band].AnnualNPPPrev = veg_var[dist][iveg][band].AnnualNPP;
+              veg_var[dist][iveg][band].AnnualNPP = 0;
+            }
           }
         }
       }

src/initialize_lake.c

@@ -74,6 +74,7 @@ int initialize_lake (lake_var_struct   *lake,
   2012-Feb-08 Renamed depth_full_snow_cover to max_snow_distrib_slope
 	      and clarified the descriptions of the SPATIAL_SNOW
 	      option.							TJB
+  2013-Jul-25 Added soil carbon terms.					TJB
 **********************************************************************/
 {
   extern option_struct options;
@@ -331,6 +332,16 @@ int initialize_lake (lake_var_struct   *lake,
       lake->soil.pot_evap[i]       = 0.0;
     }
   }
+  if (options.CARBON) {
+    lake->soil.RhLitter = 0.0;
+    lake->soil.RhLitter2Atm = 0.0;
+    lake->soil.RhInter = 0.0;
+    lake->soil.RhSlow = 0.0;
+    lake->soil.RhTot = 0.0;
+    lake->soil.CLitter = 0.0;
+    lake->soil.CInter = 0.0;
+    lake->soil.CSlow = 0.0;
+  }
 
   return(0);
 

src/initialize_soil.c

@@ -26,6 +26,7 @@ void initialize_soil (cell_data_struct **cell,
   2009-Dec-11 Removed min_liq and options.MIN_LIQ.			TJB
   2011-Mar-01 Now initializes more cell data structure terms, including
 	      asat and zwt.						TJB
+  2013-Jul-25 Added soil carbon terms.					TJB
 **********************************************************************/
 {
   extern option_struct options;
@@ -53,6 +54,9 @@ void initialize_soil (cell_data_struct **cell,
       }
       compute_runoff_and_asat(soil_con, tmp_moist, 0, &(cell[veg][band].asat), &tmp_runoff);
       wrap_compute_zwt(soil_con, &(cell[veg][band]));
+      cell[veg][band].CLitter = 0;
+      cell[veg][band].CInter = 0;
+      cell[veg][band].CSlow = 0;
     }
   }
 

src/initialize_veg.c

@@ -21,6 +21,7 @@ void initialize_veg(veg_var_struct      **veg_var,
            types to include the wetland vegetation fraction when the 
            lake model is active.                                  LCB
   2013-Jul-25 Added photosynthesis terms.				TJB
+  2013-Jul-25 Added AnnualNPP.						TJB
 
 **********************************************************************/
 {
@@ -45,6 +46,8 @@ void initialize_veg(veg_var_struct      **veg_var,
         veg_var[i][j].Rgrowth = 0.0;
         veg_var[i][j].Raut = 0.0;
         veg_var[i][j].NPP = 0.0;
+        veg_var[i][j].AnnualNPP = 0.0;
+        veg_var[i][j].AnnualNPPPrev = 0.0;
         for ( k = 0 ; k < options.Ncanopy ; k++ ) {
           veg_var[i][j].NscaleFactor[k] = 0.0;
           veg_var[i][j].aPARLayer[k] = 0.0;

src/lakes.eb.c

@@ -1931,6 +1931,7 @@ int water_balance (lake_var_struct *lake, lake_con_struct lake_con, int dt, dist
   2011-Mar-07 Fixed bug in computation of lake->soil.runoff, baseflow, etc .	TJB
   2011-Mar-31 Fixed typo in declaration of frost_fract.				TJB
   2011-Sep-22 Added logic to handle lake snow cover extent.			TJB
+  2013-Jul-25 Added soil carbon terms.						TJB
 **********************************************************************/
 {
   extern option_struct   options;
@@ -2385,6 +2386,10 @@ int water_balance (lake_var_struct *lake, lake_con_struct lake_con, int dt, dist
     }
   }
 
+  if (options.CARBON) {
+    advect_carbon_storage(lakefrac, newfraction, lake, &(cell[WET][iveg][band]));
+  }
+
   free((char*)delta_moist);
   free((char*)moist);
 
@@ -2766,3 +2771,39 @@ void rescale_snow_energy_fluxes(double oldfrac,
 
 }
 
+void advect_carbon_storage(double lakefrac,
+                           double newfraction,
+                           lake_var_struct  *lake,
+                           cell_data_struct *cell)
+/**********************************************************************
+  advect_carbon_storage	Ted Bohn	2013
+
+  Function to update carbon storage in the lake and wetland soil columns
+  to account for changes in wetland area.
+
+  Modifications:
+**********************************************************************/
+{
+
+  extern option_struct options;
+  int i,k;
+
+  if (newfraction > lakefrac) { // lake grew, wetland shrank
+
+    if (newfraction < SMALL) newfraction = SMALL;
+    lake->soil.CLitter = (lakefrac*lake->soil.CLitter + (newfraction-lakefrac)*cell->CLitter)/newfraction;
+    lake->soil.CInter = (lakefrac*lake->soil.CInter + (newfraction-lakefrac)*cell->CInter)/newfraction;
+    lake->soil.CSlow = (lakefrac*lake->soil.CSlow + (newfraction-lakefrac)*cell->CSlow)/newfraction;
+
+  }
+
+  else if (newfraction < lakefrac) { // lake shrank, wetland grew
+
+    if ((1-newfraction) < SMALL) newfraction = 1 - SMALL;
+    cell->CLitter = ((lakefrac-newfraction)*lake->soil.CLitter + (1-lakefrac)*cell->CLitter)/(1-newfraction);
+    cell->CInter = ((lakefrac-newfraction)*lake->soil.CInter + (1-lakefrac)*cell->CInter)/(1-newfraction);
+    cell->CSlow = ((lakefrac-newfraction)*lake->soil.CSlow + (1-lakefrac)*cell->CSlow)/(1-newfraction);
+
+  }
+
+}