! ================================================================================================================================= ! MODULE : stomate_soilcarbon ! ! CONTACT : orchidee-help _at_ listes.ipsl.fr ! ! LICENCE : IPSL (2006) ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! !>\BRIEF Calculate soil dynamics largely following the Century model !! !!\n DESCRIPTION: None !! !! RECENT CHANGE(S): None !! !! SVN : !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_soilcarbon.f90 $ !! $Date: 2018-01-24 08:45:20 +0100 (Wed, 24 Jan 2018) $ !! $Revision: 4906 $ !! \n !_ ================================================================================================================================ MODULE stomate_soilcarbon ! modules used: USE ioipsl_para USE stomate_data USE constantes USE xios_orchidee IMPLICIT NONE ! private & public routines PRIVATE PUBLIC soilcarbon,soilcarbon_clear ! Variables shared by all subroutines in this module LOGICAL, SAVE :: firstcall_soilcarbon = .TRUE. !! Is this the first call? (true/false) !$OMP THREADPRIVATE(firstcall_soilcarbon) CONTAINS !! ================================================================================================================================ !! SUBROUTINE : soilcarbon_clear !! !>\BRIEF Set the flag ::firstcall_soilcarbon to .TRUE. and as such activate sections 1.1.2 and 1.2 of the subroutine soilcarbon !! (see below). !! !_ ================================================================================================================================ SUBROUTINE soilcarbon_clear firstcall_soilcarbon=.TRUE. ENDSUBROUTINE soilcarbon_clear !! ================================================================================================================================ !! SUBROUTINE : soilcarbon !! !>\BRIEF Computes the soil respiration and carbon stocks, essentially !! following Parton et al. (1987). !! !! DESCRIPTION : The soil is divided into 3 carbon pools, with different !! characteristic turnover times : active (1-5 years), slow (20-40 years) !! and passive (200-1500 years).\n !! There are three types of carbon transferred in the soil:\n !! - carbon input in active and slow pools from litter decomposition,\n !! - carbon fluxes between the three pools,\n !! - carbon losses from the pools to the atmosphere, i.e., soil respiration.\n !! !! The subroutine performs the following tasks:\n !! !! Section 1.\n !! The flux fractions (f) between carbon pools are defined based on Parton et !! al. (1987). The fractions are constants, except for the flux fraction from !! the active pool to the slow pool, which depends on the clay content,\n !! \latexonly !! \input{soilcarbon_eq1.tex} !! \endlatexonly\n !! In addition, to each pool is assigned a constant turnover time.\n !! !! Section 2.\n !! The carbon input, calculated in the stomate_litter module, is added to the !! carbon stock of the different pools.\n !! !! Section 3.\n !! First, the outgoing carbon flux of each pool is calculated. It is !! proportional to the product of the carbon stock and the ratio between the !! iteration time step and the residence time:\n !! \latexonly !! \input{soilcarbon_eq2.tex} !! \endlatexonly !! ,\n !! Note that in the case of crops, the additional multiplicative factor !! integrates the faster decomposition due to tillage (following Gervois et !! al. (2008)). !! In addition, the flux from the active pool depends on the clay content:\n !! \latexonly !! \input{soilcarbon_eq3.tex} !! \endlatexonly !! ,\n !! Each pool is then cut from the carbon amount corresponding to each outgoing !! flux:\n !! \latexonly !! \input{soilcarbon_eq4.tex} !! \endlatexonly\n !! Second, the flux fractions lost to the atmosphere is calculated in each pool !! by subtracting from 1 the pool-to-pool flux fractions. The soil respiration !! is then the summed contribution of all the pools,\n !! \latexonly !! \input{soilcarbon_eq5.tex} !! \endlatexonly\n !! Finally, each carbon pool accumulates the contribution of the other pools: !! \latexonly !! \input{soilcarbon_eq6.tex} !! \endlatexonly !! !! Section 4.\n !! If the flag SPINUP_ANALYTIC is set to true, the matrix A is updated following !! Lardy (2011). !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUTS VARIABLE(S): carbon, resp_hetero_soil !! !! REFERENCE(S) : !! - Parton, W.J., D.S. Schimel, C.V. Cole, and D.S. Ojima. 1987. Analysis of !! factors controlling soil organic matter levels in Great Plains grasslands. !! Soil Sci. Soc. Am. J., 51, 1173-1179. !! - Gervois, S., P. Ciais, N. de Noblet-Ducoudre, N. Brisson, N. Vuichard, !! and N. Viovy (2008), Carbon and water balance of European croplands !! throughout the 20th century, Global Biogeochem. Cycles, 22, GB2022, !! doi:10.1029/2007GB003018. !! - Lardy, R, et al., A new method to determine soil organic carbon equilibrium, !! Environmental Modelling & Software (2011), doi:10.1016|j.envsoft.2011.05.016 !! !! FLOWCHART : !! \latexonly !! \includegraphics[scale=0.5]{soilcarbon_flowchart.jpg} !! \endlatexonly !! \n !_ ================================================================================================================================ SUBROUTINE soilcarbon (npts, clay, & soilcarbon_input, control_temp, control_moist, veget_cov_max, & carbon, resp_hetero_soil, MatrixA) !! 0. Variable and parameter declaration !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless) REAL(r_std), DIMENSION(npts), INTENT(in) :: clay !! Clay fraction (unitless, 0-1) REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in) :: soilcarbon_input !! Amount of carbon going into the carbon pools from litter decomposition \f$(gC m^{-2} day^{-1})$\f REAL(r_std), DIMENSION(npts,nlevs), INTENT(in) :: control_temp !! Temperature control of heterotrophic respiration (unitless: 0->1) REAL(r_std), DIMENSION(npts,nlevs), INTENT(in) :: control_moist !! Moisture control of heterotrophic respiration (unitless: 0.25->1) REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_cov_max !! Fractional coverage: maximum share of the pixel taken by a pft !! 0.2 Output variables REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: resp_hetero_soil !! Soil heterotrophic respiration \f$(gC m^{-2} (dt_sechiba one_day^{-1})^{-1})$\f !! 0.3 Modified variables REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout) :: carbon !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(inout) :: MatrixA !! Matrix containing the fluxes between the carbon pools !! per sechiba time step !! @tex $(gC.m^2.day^{-1})$ @endtex !! 0.4 Local variables REAL(r_std) :: dt !! Time step \f$(dt_sechiba one_day^{-1})$\f REAL(r_std), SAVE, DIMENSION(ncarb) :: carbon_tau !! Residence time in carbon pools (days) !$OMP THREADPRIVATE(carbon_tau) REAL(r_std), DIMENSION(npts,ncarb,ncarb) :: frac_carb !! Flux fractions between carbon pools !! (second index=origin, third index=destination) !! (unitless, 0-1) REAL(r_std), DIMENSION(npts,ncarb) :: frac_resp !! Flux fractions from carbon pools to the atmosphere (respiration) (unitless, 0-1) REAL(r_std), DIMENSION(npts,ncarb,nelements) :: fluxtot !! Total flux out of carbon pools \f$(gC m^{2})$\f REAL(r_std), DIMENSION(npts,ncarb,ncarb,nelements) :: flux !! Fluxes between carbon pools \f$(gC m^{2})$\f CHARACTER(LEN=7), DIMENSION(ncarb) :: carbon_str !! Name of the carbon pools for informative outputs (unitless) INTEGER(i_std) :: k,kk,m,j,ij !! Indices (unitless) REAL(r_std), DIMENSION(npts,nvm,ncarb) :: decomp_rate_soilcarbon !! Decomposition rate of the soil carbon pools (s) REAL(r_std), DIMENSION(npts,ncarb) :: tsoilpools !! Diagnostic for soil carbon turnover rate by pool (1/s) !_ ================================================================================================================================ !! printlev is the level of diagnostic information, 0 (none) to 4 (full) IF (printlev>=3) WRITE(numout,*) 'Entering soilcarbon' !! 1. Initializations dt = dt_sechiba/one_day !! 1.1 Get soil "constants" !! 1.1.1 Flux fractions between carbon pools: depend on clay content, recalculated each time ! From active pool: depends on clay content frac_carb(:,iactive,iactive) = zero frac_carb(:,iactive,ipassive) = frac_carb_ap frac_carb(:,iactive,islow) = un - (metabolic_ref_frac - active_to_pass_clay_frac*clay(:)) - frac_carb(:,iactive,ipassive) ! 1.1.1.2 from slow pool frac_carb(:,islow,islow) = zero frac_carb(:,islow,iactive) = frac_carb_sa frac_carb(:,islow,ipassive) = frac_carb_sp ! From passive pool frac_carb(:,ipassive,ipassive) = zero frac_carb(:,ipassive,iactive) = frac_carb_pa frac_carb(:,ipassive,islow) = frac_carb_ps IF ( firstcall_soilcarbon ) THEN !! 1.1.2 Residence times in carbon pools (days) carbon_tau(iactive) = carbon_tau_iactive * one_year ! 1.5 years. This is same as CENTURY. But, in Parton et al. (1987), it's weighted by moisture and temperature dependences. carbon_tau(islow) = carbon_tau_islow * one_year ! 25 years. This is same as CENTURY. But, in Parton et al. (1987), it's weighted by moisture and temperature dependences. carbon_tau(ipassive) = carbon_tau_ipassive * one_year ! 1000 years. This is same as CENTURY. But, in Parton et al. (1987), it's weighted by moisture and temperature dependences. !! 1.2 Messages : display the residence times carbon_str(iactive) = 'active' carbon_str(islow) = 'slow' carbon_str(ipassive) = 'passive' IF (printlev >= 2) THEN WRITE(numout,*) 'soilcarbon:' WRITE(numout,*) ' > minimal carbon residence time in carbon pools (d):' DO k = 1, ncarb ! Loop over carbon pools WRITE(numout,*) '(1, ::carbon_str(k)):',carbon_str(k),' : (1, ::carbon_tau(k)):',carbon_tau(k) ENDDO WRITE(numout,*) ' > flux fractions between carbon pools: depend on clay content' END IF firstcall_soilcarbon = .FALSE. ENDIF !! 1.3 Set soil respiration and decomposition rate to zero resp_hetero_soil(:,:) = zero decomp_rate_soilcarbon(:,:,:) = zero !! 2. Update the carbon stocks with the different soil carbon input carbon(:,:,:) = carbon(:,:,:) + soilcarbon_input(:,:,:) * dt !! 3. Fluxes between carbon reservoirs, and to the atmosphere (respiration) \n !! 3.1. Determine the respiration fraction : what's left after ! subtracting all the 'pool-to-pool' flux fractions ! Diagonal elements of frac_carb are zero ! VPP killer: ! frac_resp(:,:) = 1. - SUM( frac_carb(:,:,:), DIM=3 ) frac_resp(:,:) = un - frac_carb(:,:,iactive) - frac_carb(:,:,islow) - & frac_carb(:,:,ipassive) !! 3.2. Calculate fluxes DO m = 1, nvm ! Loop over # PFTs !! 3.2.1. Flux out of pools DO k = 1, ncarb ! Loop over carbon pools from which the flux comes ! Determine total flux out of pool ! S.L. Piao 2006/05/05 - for crop multiply tillage factor of decomposition ! Not crop IF ( natural(m) ) THEN fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * & control_moist(:,ibelow) * control_temp(:,ibelow) decomp_rate_soilcarbon(:,m,k)=dt/carbon_tau(k) * & control_moist(:,ibelow) * control_temp(:,ibelow) ! C3 crop ELSEIF ( (.NOT. natural(m)) .AND. (.NOT. is_c4(m)) ) THEN fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * & control_moist(:,ibelow) * control_temp(:,ibelow) * flux_tot_coeff(1) decomp_rate_soilcarbon(:,m,k)=dt/carbon_tau(k) * & control_moist(:,ibelow) * control_temp(:,ibelow)*flux_tot_coeff(1) ! C4 Crop ELSEIF ( (.NOT. natural(m)) .AND. is_c4(m) ) THEN fluxtot(:,k,icarbon) = dt/carbon_tau(k) * carbon(:,k,m) * & control_moist(:,ibelow) * control_temp(:,ibelow) * flux_tot_coeff(2) decomp_rate_soilcarbon(:,m,k)=dt/carbon_tau(k) * & control_moist(:,ibelow) *control_temp(:,ibelow)*flux_tot_coeff(2) ENDIF ! END - S.L. Piao 2006/05/05 - for crop multiply tillage factor of decomposition ! Carbon flux from active pools depends on clay content IF ( k .EQ. iactive ) THEN fluxtot(:,k,icarbon) = fluxtot(:,k,icarbon) * ( un - flux_tot_coeff(3) * clay(:) ) decomp_rate_soilcarbon(:,m,k)=decomp_rate_soilcarbon(:,m,k)* ( un - flux_tot_coeff(3) * clay(:) ) ENDIF ! Update the loss in each carbon pool carbon(:,k,m) = carbon(:,k,m) - fluxtot(:,k,icarbon) ! Fluxes towards the other pools (k -> kk) DO kk = 1, ncarb ! Loop over the carbon pools where the flux goes flux(:,k,kk,icarbon) = frac_carb(:,k,kk) * fluxtot(:,k,icarbon) ENDDO ENDDO ! End of loop over carbon pools !! 3.2.2 respiration !BE CAREFUL: Here resp_hetero_soil is divided by dt to have a value which corresponds to ! the sechiba time step but then in stomate.f90 resp_hetero_soil is multiplied by dt. ! Perhaps it could be simplified. Moreover, we must totally adapt the routines to the dt_sechiba/one_day ! time step and avoid some constructions that could create bug during future developments. ! ! VPP killer: ! resp_hetero_soil(:,m) = SUM( frac_resp(:,:) * fluxtot(:,:), DIM=2 ) / dt resp_hetero_soil(:,m) = & ( frac_resp(:,iactive) * fluxtot(:,iactive,icarbon) + & frac_resp(:,islow) * fluxtot(:,islow,icarbon) + & frac_resp(:,ipassive) * fluxtot(:,ipassive,icarbon) ) / dt !! 3.2.3 add fluxes to active, slow, and passive pools ! VPP killer: ! carbon(:,:,m) = carbon(:,:,m) + SUM( flux(:,:,:), DIM=2 ) DO k = 1, ncarb ! Loop over carbon pools carbon(:,k,m) = carbon(:,k,m) + & flux(:,iactive,k,icarbon) + flux(:,ipassive,k,icarbon) + flux(:,islow,k,icarbon) ENDDO ! Loop over carbon pools ENDDO ! End loop over PFTs !! 4. (Quasi-)Analytical Spin-up !! 4.1.1 Finish to fill MatrixA with fluxes between soil pools IF (spinup_analytic) THEN DO m = 2,nvm ! flux leaving the active pool MatrixA(:,m,iactive_pool,iactive_pool) = moins_un * & dt/carbon_tau(iactive) * & control_moist(:,ibelow) * control_temp(:,ibelow) * & ( 1. - flux_tot_coeff(3) * clay(:)) ! flux received by the active pool from the slow pool MatrixA(:,m,iactive_pool,islow_pool) = frac_carb(:,islow,iactive)*dt/carbon_tau(islow) * & control_moist(:,ibelow) * control_temp(:,ibelow) ! flux received by the active pool from the passive pool MatrixA(:,m,iactive_pool,ipassive_pool) = frac_carb(:,ipassive,iactive)*dt/carbon_tau(ipassive) * & control_moist(:,ibelow) * control_temp(:,ibelow) ! flux received by the slow pool from the active pool MatrixA(:,m,islow_pool,iactive_pool) = frac_carb(:,iactive,islow) *& dt/carbon_tau(iactive) * & control_moist(:,ibelow) * control_temp(:,ibelow) * & ( 1. - flux_tot_coeff(3) * clay(:) ) ! flux leaving the slow pool MatrixA(:,m,islow_pool,islow_pool) = moins_un * & dt/carbon_tau(islow) * & control_moist(:,ibelow) * control_temp(:,ibelow) ! flux received by the passive pool from the active pool MatrixA(:,m,ipassive_pool,iactive_pool) = frac_carb(:,iactive,ipassive)* & dt/carbon_tau(iactive) * & control_moist(:,ibelow) * control_temp(:,ibelow) *& ( 1. - flux_tot_coeff(3) * clay(:) ) ! flux received by the passive pool from the slow pool MatrixA(:,m,ipassive_pool,islow_pool) = frac_carb(:,islow,ipassive) * & dt/carbon_tau(islow) * & control_moist(:,ibelow) * control_temp(:,ibelow) ! flux leaving the passive pool MatrixA(:,m,ipassive_pool,ipassive_pool) = moins_un * & dt/carbon_tau(ipassive) * & control_moist(:,ibelow) * control_temp(:,ibelow) IF ( (.NOT. natural(m)) .AND. (.NOT. is_c4(m)) ) THEN ! C3crop ! flux leaving the active pool MatrixA(:,m,iactive_pool,iactive_pool) = MatrixA(:,m,iactive_pool,iactive_pool) * & flux_tot_coeff(1) ! flux received by the active pool from the slow pool MatrixA(:,m,iactive_pool,islow_pool)= MatrixA(:,m,iactive_pool,islow_pool) * & flux_tot_coeff(1) ! flux received by the active pool from the passive pool MatrixA(:,m,iactive_pool,ipassive_pool) = MatrixA(:,m,iactive_pool,ipassive_pool) * & flux_tot_coeff(1) ! flux received by the slow pool from the active pool MatrixA(:,m,islow_pool,iactive_pool) = MatrixA(:,m,islow_pool,iactive_pool) * & flux_tot_coeff(1) ! flux leaving the slow pool MatrixA(:,m,islow_pool,islow_pool) = MatrixA(:,m,islow_pool,islow_pool) * & flux_tot_coeff(1) ! flux received by the passive pool from the active pool MatrixA(:,m,ipassive_pool,iactive_pool) = MatrixA(:,m,ipassive_pool,iactive_pool) * & flux_tot_coeff(1) ! flux received by the passive pool from the slow pool MatrixA(:,m,ipassive_pool,islow_pool) = MatrixA(:,m,ipassive_pool,islow_pool) * & flux_tot_coeff(1) ! flux leaving the passive pool MatrixA(:,m,ipassive_pool,ipassive_pool) = MatrixA(:,m,ipassive_pool,ipassive_pool) *& flux_tot_coeff(1) ENDIF ! (.NOT. natural(m)) .AND. (.NOT. is_c4(m)) IF ( (.NOT. natural(m)) .AND. is_c4(m) ) THEN ! C4crop ! flux leaving the active pool MatrixA(:,m,iactive_pool,iactive_pool) = MatrixA(:,m,iactive_pool,iactive_pool) * & flux_tot_coeff(2) ! flux received by the active pool from the slow pool MatrixA(:,m,iactive_pool,islow_pool)= MatrixA(:,m,iactive_pool,islow_pool) * & flux_tot_coeff(2) ! flux received by the active pool from the passive pool MatrixA(:,m,iactive_pool,ipassive_pool) = MatrixA(:,m,iactive_pool,ipassive_pool) * & flux_tot_coeff(2) ! flux received by the slow pool from the active pool MatrixA(:,m,islow_pool,iactive_pool) = MatrixA(:,m,islow_pool,iactive_pool) * & flux_tot_coeff(2) ! flux leaving the slow pool MatrixA(:,m,islow_pool,islow_pool) = MatrixA(:,m,islow_pool,islow_pool) * & flux_tot_coeff(2) ! flux received by the passive pool from the active pool MatrixA(:,m,ipassive_pool,iactive_pool) = MatrixA(:,m,ipassive_pool,iactive_pool) * & flux_tot_coeff(2) ! flux received by the passive pool from the slow pool MatrixA(:,m,ipassive_pool,islow_pool) = MatrixA(:,m,ipassive_pool,islow_pool) * & flux_tot_coeff(2) ! flux leaving the passive pool MatrixA(:,m,ipassive_pool,ipassive_pool) = MatrixA(:,m,ipassive_pool,ipassive_pool) * & flux_tot_coeff(2) ENDIF ! (.NOT. natural(m)) .AND. is_c4(m) IF (printlev>=4) WRITE(numout,*)'Finish to fill MatrixA' ENDDO ! Loop over # PFTS ! 4.2 Add Identity for each submatrix(7,7) DO j = 1,nbpools MatrixA(:,:,j,j) = MatrixA(:,:,j,j) + un ENDDO ENDIF ! (spinup_analytic) ! Output diagnostics DO k = 1, ncarb ! Loop over carbon pools DO ij = 1, npts IF (SUM(decomp_rate_soilcarbon(ij,:,k)*veget_cov_max(ij,:)) > min_sechiba) THEN tsoilpools(ij,k) = 1./(SUM(decomp_rate_soilcarbon(ij,:,k)*veget_cov_max(ij,:))/dt_sechiba) ELSE tsoilpools(ij,k) = xios_default_val END IF END DO END DO CALL xios_orchidee_send_field("tSoilPools",tsoilpools) IF (printlev>=4) WRITE(numout,*) 'Leaving soilcarbon' END SUBROUTINE soilcarbon END MODULE stomate_soilcarbon