! =================================================================================================================================
! MODULE       : stomate_litter
!
! CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
! LICENCE      : IPSL (2006)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF       Update litter and lignine content after litter fall and 
!! calculating litter decomposition.      
!!
!!\n DESCRIPTION: None
!!
!! RECENT CHANGE(S): None
!!
!! REFERENCE(S)	: None
!!
!! SVN          :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_litter.f90 $
!! $Date: 2024-03-07 18:38:03 +0100 (jeu. 07 mars 2024) $
!! $Revision: 8462 $
!! \n
!_ ================================================================================================================================

MODULE stomate_litter

  ! modules used:
  USE xios_orchidee
  USE ioipsl_para
  USE stomate_data
  USE constantes
  USE constantes_soil
  USE pft_parameters

  IMPLICIT NONE

  ! private & public routines

  PRIVATE
  PUBLIC littercalc,littercalc_clear, deadleaf

  LOGICAL, SAVE                        :: firstcall_litter = .TRUE.       !! first call
!$OMP THREADPRIVATE(firstcall_litter)

CONTAINS

!! ================================================================================================================================
!! SUBROUTINE   : littercalc_calc
!!
!!\BRIEF        Set the flag ::firstcall_litter to .TRUE. and as such activate section
!! 1.1 of the subroutine littercalc (see below).
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S) : None
!!
!! MAIN OUTPUT VARIABLE(S) : None
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE littercalc_clear
    firstcall_litter =.TRUE.
  END SUBROUTINE littercalc_clear


!! ================================================================================================================================
!! SUBROUTINE   : littercalc
!!
!!\BRIEF        Calculation of the litter decomposition and therefore of the 
!! heterotrophic respiration from litter following Parton et al. (1987).
!!
!! DESCRIPTION  : The littercal routine splits the litter in 4 pools: 
!! aboveground metaboblic, aboveground structural, belowground metabolic and 
!! belowground structural. the fraction (F) of plant material going to metabolic 
!! and structural is defined following Parton et al. (1987)
!! \latexonly
!! \input{littercalc1.tex}
!! \endlatexonly
!! \n
!! where L is the lignin content of the plant carbon pools considered and CN 
!! its CN ratio. L and CN are fixed parameters for each plant carbon pools,
!! therefore it is the ratio between each plant carbon pool within a PFT, which
!! controlled the part of the total litter, that will be considered as
!! recalcitrant (i.e. structural litter) or labile (i.e. metabolic litter).\n  
!! 
!! The routine calculates the fraction of aboveground litter which is metabolic
!! or structural (the litterpart variable) which is then used in lpj_fire.f90.\n 
!! 
!! In the section 2, the routine calculate the new plant material entering the
!! litter pools by phenological death of plants organs (corresponding to the
!! variable turnover) and by fire, herbivory and others non phenological causes
!! (variable bm_to_litter). This calculation is first done for each PFT and then
!! the values calculated for each PFT are added up. Following the same approach
!! the lignin content of the total structural litter is calculated and will be
!! then used as a factor control of the decomposition of the structural litter
!! (lignin_struc) in the section 5.1.2. A test is performed to avoid that we add
!! more lignin than structural litter. Finally, the variable litterpart is
!! updated.\n
!! 
!! In the section 3 and 4 the temperature and the moisture controlling the
!! decomposition are calculated for above and belowground. For aboveground
!! litter, air temperature and litter moisture are calculated in sechiba and used 
!! directly. For belowground, soil temperature and moisture are also calculated 
!! in sechiba but are modulated as a function of the soil depth. The modulation 
!! is a multiplying factor exponentially distributed between 0 (in depth) and 1
!! in surface.\n  
!! 
!! Then, in the section 5, the routine calculates the structural litter decomposition 
!! (C) following first order kinetics following Parton et al. (1987).
!! \latexonly
!! \input{littercalc2.tex}
!! \endlatexonly
!! \n
!! with k the decomposition rate of the structural litter. 
!! k corresponds to
!! \latexonly
!! \input{littercalc3.tex}
!! \endlatexonly
!! \n
!! with littertau the turnover rate, T a function of the temperature and M a function of
!! the moisture described below.\n
!!  
!! Then, the fraction of dead leaves (DL) composed by aboveground structural litter is
!! calculated as following
!! \latexonly
!! \input{littercalc4.tex}
!! \endlatexonly
!! \n
!! with k the decomposition rate of the structural litter previously
!! described.\n
!!
!! In the section 5.1, the fraction of decomposed structural litter
!! incorporated to the soil (Input) and its associated heterotrophic respiration are
!! calculated. For structural litter, the C decomposed could go in the active
!! soil carbon pool or in the slow carbon, as described in 
!! stomate_soilcarbon.f90.\n
!! \latexonly
!! \input{littercalc5.tex}
!! \endlatexonly
!! \n
!! with f a parameter describing the fraction of structural litter incorporated
!! into the considered soil carbon pool, C the amount of litter decomposed and L 
!! the amount of lignin in the litter. The litter decomposed which is not
!! incorporated into the soil is respired.\n
!!
!! In the section 5.2, the fraction of decomposed metabolic litter
!! incorporated to the soil and its associated heterotrophic respiration are
!! calculated with the same approaches presented for 5.1 but no control factor
!! depending on the lignin content are used.\n
!! 
!! In the section 6 the dead leaf cover is calculated through a call to the 
!! deadleaf subroutine presented below.\n
!!
!! In the section 7, if the flag SPINUP_ANALYTIC is set to true, we fill MatrixA
!! and VectorB following Lardy(2011).
!!
!! MAIN OUTPUT VARIABLES: ::deadleaf_cover, ::resp_hetero_litter, ::soilcarbon_input, 
!! ::control_temp, ::control_moist
!!
!! REFERENCES:
!! - Parton, WJ, Schimel, DS, Cole, CV, and Ojima, DS. 1987. Analysis
!! of factors controlling soil organic matter levels in Great Plains
!! grasslands. Soil Science Society of America journal (USA)
!! (51):1173-1179.
!! - 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){littercalcflow.jpg}
!! \endlatexonly
!! \n
!_ ================================================================================================================================

  SUBROUTINE littercalc (npts, &
       turnover, bm_to_litter, &
       veget_cov_max, tsurf, tsoil, soilhum, litterhum, &
       soilhumSAT, litterhumSAT, &
       litterpart, litter, dead_leaves, lignin_struc, &
       deadleaf_cover, resp_hetero_litter, &
       soilcarbon_input, control_temp, control_moist, &
       MatrixA, VectorB, bulk, clay, carbon)

    !! 0. Variable and parameter declaration
    
    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                                  :: npts               !! Domain size - number of grid pixels
    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: turnover         !! Turnover rates of plant biomass 
                                                                                      !! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex 
    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(in) :: bm_to_litter     !! Conversion of biomass to litter 
                                                                                      !! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex 
    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                  :: veget_cov_max      !! PFT "Maximal" coverage fraction of a PFT 
                                                                                      !! defined in the input vegetation map 
                                                                                      !! @tex $(m^2 m^{-2})$ @endtex 
    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: tsurf              !! Temperature (K) at the surface
    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)               :: tsoil              !! Soil temperature (K)
    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)               :: soilhum            !! Daily soil humidity of each soil layer 
                                                                                      !! (unitless)
    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)               :: soilhumSAT         !! Daily soil humidity of each soil layer 
                                                                                      !! (unitless)
    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: litterhum          !! Daily litter humidity (unitless)
    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: litterhumSAT       !! Daily litter humidity (unitless)
    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: clay               !! Clay fraction (unitless, 0-1) 
    REAL(r_std), DIMENSION(npts),INTENT(in)                     :: bulk               !! Bulk density (kg/m**3)  
    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in)          :: carbon             !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f

    !! 0.2 Output variables
    
    REAL(r_std), DIMENSION(npts), INTENT(out)                   :: deadleaf_cover     !! Fraction of soil covered by dead leaves 
                                                                                      !! over all PFTs (0-1, unitless)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)               :: resp_hetero_litter !! Litter heterotrophic respiration. The unit
                                                                                      !! is given by m^2 of ground.  
                                                                                      !! @tex $(gC dt_sechiba one\_day^{-1}) m^{-2})$ @endtex 
    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(out)         :: soilcarbon_input   !! Quantity of carbon going into carbon pools
                                                                                      !! from litter decomposition. The unit is  
                                                                                      !! given by m^2 of ground 
                                                                                      !! @tex $(gC m^{-2} dt\_slow^{-1})$ @endtex 
    REAL(r_std), DIMENSION(npts,nlevs), INTENT(out)             :: control_temp       !! Temperature control of heterotrophic 
                                                                                      !! respiration, above and below (0-1, 
                                                                                      !! unitless)
    REAL(r_std), DIMENSION(npts,nlevs), INTENT(out)             :: control_moist      !! Moisture control of heterotrophic 
                                                                                      !! respiration (0.25-1, unitless)
    REAL(r_std), DIMENSION(npts,nvm,nbpools,nbpools), INTENT(out) :: MatrixA          !! Matrix containing the fluxes between the
                                                                                      !! carbon pools per sechiba time step 
                                                                                      !! @tex $(gC.m^2.day^{-1})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nbpools), INTENT(out)         :: VectorB          !! Vector containing the litter increase per
                                                                                      !! sechiba time step
                                                                                      !! @tex $(gC m^{-2})$ @endtex
   
    !! 0.3 Modified variables
    
    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)       :: litterpart         !! Fraction of litter above the ground 
                                                                                      !! belonging to the different PFTs (0-1, 
                                                                                      !! unitless)
    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) :: litter   !! Metabolic and structural litter,above and
                                                                                      !! below ground. The unit is given by m^2 of 
                                                                                      !! ground @tex $(gC m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)       :: dead_leaves        !! Dead leaves per ground unit area, per PFT,
                                                                                      !! metabolic and structural in 
                                                                                      !! @tex $(gC m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout)       :: lignin_struc       !! Ratio Lignin content in structural litter,
                                                                                      !! above and below ground, (0-1, unitless)
    
    !! 0.4 Local variables
 
    REAL(r_std)                                                 :: dt                 !! Number of sechiba(fast processes) time-step per day 
    REAL(r_std), SAVE, DIMENSION(nparts,nlitt)                  :: litterfrac         !! The fraction of leaves, wood, etc. that 
                                                                                      !! goes into metabolic and structural 
                                                                                      !! litterpools (0-1, unitless)
!$OMP THREADPRIVATE(litterfrac)
    REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)                :: z_soil             !! Soil levels (m)
!$OMP THREADPRIVATE(z_soil)
    REAL(r_std), DIMENSION(npts)                                :: rpc                !! Integration constant for vertical root 
                                                                                      !! profiles (unitless)
    REAL(r_std), SAVE, DIMENSION(nlitt)                         :: litter_tau         !! Turnover time in litter pools (days)
!$OMP THREADPRIVATE(litter_tau)
    REAL(r_std), SAVE, DIMENSION(nlitt,ncarb,nlevs)             :: frac_soil          !! Fraction of litter that goes into soil 
                                                                                      !! (litter -> carbon, above and below). The
                                                                                      !! remaining part goes to the atmosphere
!$OMP THREADPRIVATE(frac_soil)
    REAL(r_std), DIMENSION(npts)                                :: tsoil_decomp       !! Temperature used for decompostition in 
                                                                                      !! soil (K)
    REAL(r_std), DIMENSION(npts)                                :: soilhum_decomp     !! Humidity used for decompostition in soil
                                                                                      !! (unitless)
    REAL(r_std), DIMENSION(npts)                                :: soilhumSAT_decomp  !! Humidity used for decompostition in soil
                                                                                      !! (unitless) for Moyano et al 2012 
    REAL(r_std), DIMENSION(npts)                                :: fd                 !! Fraction of structural or metabolic litter
                                                                                      !! decomposed (unitless)
    REAL(r_std), DIMENSION(npts,nelements)                      :: qd                 !! Quantity of structural or metabolic litter
                                                                                      !! decomposed @tex $(gC m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nlevs)                      :: old_struc          !! Old structural litter, above and below 
                                                                                      !! @tex $(gC m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nlitt,nlevs,nelements)      :: litter_inc_PFT     !! Increase of litter, per PFT, metabolic and
                                                                                      !! structural, above and below ground. The 
                                                                                      !! unit is given by m^2 of ground.  
                                                                                      !! @tex $(gC m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements)      :: litter_inc         !! Increase of metabolic and structural 
                                                                                      !! litter, above and below ground. The unit 
                                                                                      !! is given by m^2 of ground. 
                                                                                      !! @tex $(gC m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nlevs)                      :: lignin_struc_inc   !! Lignin increase in structural litter, 
                                                                                      !! above and below ground. The unit is given 
                                                                                      !! by m^2 of ground. 
                                                                                      !! @tex $(gC m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nlitt,nelements)            :: litter_pft         !! Metabolic and structural litter above the 
                                                                                      !! ground per PFT
    REAL(r_std), DIMENSION(npts)                                :: zdiff_min          !! Intermediate field for looking for minimum
                                                                                      !! of what? this is not used in the code. 
                                                                                      !! [??CHECK] could we delete it?
    CHARACTER(LEN=10), DIMENSION(nlitt)                         :: litter_str         !! Messages to write output information about
                                                                                      !! the litter
    CHARACTER(LEN=22), DIMENSION(nparts)                        :: part_str           !! Messages to write output information about
                                                                                      !! the plant
    CHARACTER(LEN=7), DIMENSION(ncarb)                          :: carbon_str         !! Messages to write output information about
                                                                                      !! the soil carbon
    CHARACTER(LEN=5), DIMENSION(nlevs)                          :: level_str          !! Messages to write output information about
                                                                                      !! the level (aboveground or belowground litter)
    INTEGER(i_std)                                              :: i,j,k,l,m          !! Indices (unitless)
    INTEGER(i_std)                                              :: ier                !! Error handling
!_ ================================================================================================================================

    IF (printlev>=3) WRITE(numout,*) 'Entering littercalc'

  !! 1. Initialisations of the different fields during the first call of the routine 
    dt = dt_sechiba/one_day
    IF ( firstcall_litter ) THEN

       !! 1.1.3 litter fractions:
       !!   what fraction of leaves, wood, etc. goes into metabolic and structural litterpools
       DO k = 1, nparts

          litterfrac(k,imetabolic) = metabolic_ref_frac - metabolic_LN_ratio * LC(k) * CN(k)
          litterfrac(k,istructural) = un - litterfrac(k,imetabolic)

       ENDDO

       !! 1.1.4 residence times in litter pools (days)
       litter_tau(imetabolic) = tau_metabolic * one_year      ! .5 years
       litter_tau(istructural) = tau_struct * one_year     ! 3 years

       !! 1.1.5 decomposition flux fraction that goes into soil
       !       (litter -> carbon, above and below)
       !       1-frac_soil goes into atmosphere
       frac_soil(:,:,:) = zero

       ! structural litter: lignin fraction goes into slow pool + respiration,
       !                    rest into active pool + respiration
       frac_soil(istructural,iactive,iabove) = frac_soil_struct_aa
       frac_soil(istructural,iactive,ibelow) = frac_soil_struct_ab
       frac_soil(istructural,islow,iabove) = frac_soil_struct_sa
       frac_soil(istructural,islow,ibelow) = frac_soil_struct_sb

       ! metabolic litter: all goes into active pool + respiration.
       !   Nothing into slow or passive pool.
       frac_soil(imetabolic,iactive,iabove) = frac_soil_metab_aa
       frac_soil(imetabolic,iactive,ibelow) = frac_soil_metab_ab

       
       !! 1.2 soil levels
       ALLOCATE(z_soil(0:nslm), stat=ier)
       IF ( ier /= 0 ) CALL ipslerr_p(3,'littercalc','Pb in allocate of z_soil','','')

       z_soil(0) = zero
       z_soil(1:nslm) = diaglev(1:nslm)

       
       !! 1.3 messages
       litter_str(imetabolic) = 'metabolic'
       litter_str(istructural) = 'structural'

       carbon_str(iactive) = 'active'
       carbon_str(islow) = 'slow'
       carbon_str(ipassive) = 'passive'

       level_str(iabove) = 'above'
       level_str(ibelow) = 'below'

       part_str(ileaf) = 'leaves'
       part_str(isapabove) = 'sap above ground'
       part_str(isapbelow) = 'sap below ground'
       part_str(iheartabove) = 'heartwood above ground'
       part_str(iheartbelow) = 'heartwood below ground'
       part_str(iroot) = 'roots'
       part_str(ifruit) = 'fruits'
       part_str(icarbres) = 'carbohydrate reserve'

       IF (printlev >=2) THEN
          WRITE(numout,*) 'litter:'

          WRITE(numout,*) '   > C/N ratios: '
          DO k = 1, nparts
             WRITE(numout,*) '       ', part_str(k), ': ',CN(k)
          ENDDO

          WRITE(numout,*) '   > Lignine/C ratios: '
          DO k = 1, nparts
             WRITE(numout,*) '       ', part_str(k), ': ',LC(k)
          ENDDO

          WRITE(numout,*) '   > fraction of compartment that goes into litter: '
          DO k = 1, nparts
             DO m = 1, nlitt
                WRITE(numout,*) '       ', part_str(k), '-> ',litter_str(m), ':',litterfrac(k,m)
             ENDDO
          ENDDO

          WRITE(numout,*) '   > scaling depth for decomposition (m): ',z_decomp

          WRITE(numout,*) '   > minimal carbon residence time in litter pools (d):'
          DO m = 1, nlitt
             WRITE(numout,*) '       ',litter_str(m),':',litter_tau(m)
          ENDDO

          WRITE(numout,*) '   > litter decomposition flux fraction that really goes '
          WRITE(numout,*) '     into carbon pools (rest into the atmosphere):'
          DO m = 1, nlitt
             DO l = 1, nlevs
                DO k = 1, ncarb
                   WRITE(numout,*) '       ',litter_str(m),' ',level_str(l),' -> ',&
                        carbon_str(k),':', frac_soil(m,k,l)
                ENDDO
             ENDDO
          ENDDO
       END IF ! printlev >=2
       firstcall_litter = .FALSE.

    ENDIF

    
    !! 1.3 litter above the ground per PFT.
    DO j = 1, nvm ! Loop over # PFTs

       DO k = 1, nlitt !Loop over litter pool
          litter_pft(:,j,k,icarbon) = litterpart(:,j,k) * litter(:,k,j,iabove,icarbon)
       ENDDO

    ENDDO

    
    !! 1.4 set output to zero
    deadleaf_cover(:) = zero
    resp_hetero_litter(:,:) = zero
    soilcarbon_input(:,:,:) = zero

    
  !! 2. Add biomass to different litterpools (per m^2 of ground)
    
    !! 2.1 first, save old structural litter (needed for lignin fractions).
    !     above/below
    DO l = 1, nlevs !Loop over litter levels (above and below ground)
       DO m = 1,nvm !Loop over PFTs

          old_struc(:,m,l) = litter(:,istructural,m,l,icarbon)

       ENDDO
    ENDDO

    
    !! 2.2 update litter, dead leaves, and lignin content in structural litter
    litter_inc(:,:,:,:,:) = zero
    lignin_struc_inc(:,:,:) = zero

    DO j = 1,nvm !Loop over PFTs

       !! 2.2.1 litter
       DO k = 1, nlitt    !Loop over litter pools (metabolic and structural)

          !! 2.2.2 calculate litter increase (per m^2 of ground).
          !       Only a given fracion of fruit turnover is directly coverted into litter.
          !       Litter increase for each PFT, structural and metabolic, above/below
          litter_inc_PFT(:,j,k,iabove,icarbon) = &
               litterfrac(ileaf,k) * bm_to_litter(:,j,ileaf,icarbon) + &
               litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove,icarbon) + &
               litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove,icarbon) + &
               litterfrac(ifruit,k) * bm_to_litter(:,j,ifruit,icarbon) + &
               litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres,icarbon) + &
               litterfrac(ileaf,k) * turnover(:,j,ileaf,icarbon) + &
               litterfrac(isapabove,k) * turnover(:,j,isapabove,icarbon) + &
               litterfrac(iheartabove,k) * turnover(:,j,iheartabove,icarbon) + &
               litterfrac(ifruit,k) * turnover(:,j,ifruit,icarbon) + &
               litterfrac(icarbres,k) * turnover(:,j,icarbres,icarbon)

          litter_inc_PFT(:,j,k,ibelow,icarbon) = &
               litterfrac(isapbelow,k) * bm_to_litter(:,j,isapbelow,icarbon) + &
               litterfrac(iheartbelow,k) * bm_to_litter(:,j,iheartbelow,icarbon) + &
               litterfrac(iroot,k) * bm_to_litter(:,j,iroot,icarbon) + &
               litterfrac(isapbelow,k) * turnover(:,j,isapbelow,icarbon) + &
               litterfrac(iheartbelow,k) * turnover(:,j,iheartbelow,icarbon) + &
               litterfrac(iroot,k) * turnover(:,j,iroot,icarbon)

          ! litter increase, met/struct, above/below
          litter_inc(:,k,j,iabove,icarbon) = litter_inc(:,k,j,iabove,icarbon) + litter_inc_PFT(:,j,k,iabove,icarbon)
          litter_inc(:,k,j,ibelow,icarbon) = litter_inc(:,k,j,ibelow,icarbon) + litter_inc_PFT(:,j,k,ibelow,icarbon)

          !! 2.2.3 dead leaves, for soil cover.
          dead_leaves(:,j,k) = &
               dead_leaves(:,j,k) + &
               litterfrac(ileaf,k) * ( bm_to_litter(:,j,ileaf,icarbon) + turnover(:,j,ileaf,icarbon) )

          !! 2.2.4 lignin increase in structural litter
          IF ( k .EQ. istructural ) THEN

             lignin_struc_inc(:,j,iabove) = &
                  lignin_struc_inc(:,j,iabove) + &
                  LC(ileaf) * bm_to_litter(:,j,ileaf,icarbon) + &
                  LC(isapabove) * bm_to_litter(:,j,isapabove,icarbon) + &
                  LC(iheartabove) * bm_to_litter(:,j,iheartabove,icarbon) + &
                  LC(ifruit) * bm_to_litter(:,j,ifruit,icarbon) + &
                  LC(icarbres) * bm_to_litter(:,j,icarbres,icarbon) + &
                  LC(ileaf) * turnover(:,j,ileaf,icarbon) + &
                  LC(isapabove) * turnover(:,j,isapabove,icarbon) + &
                  LC(iheartabove) * turnover(:,j,iheartabove,icarbon) + &
                  LC(ifruit) * turnover(:,j,ifruit,icarbon) + &
                  LC(icarbres) * turnover(:,j,icarbres,icarbon)

             lignin_struc_inc(:,j,ibelow) = &
                  lignin_struc_inc(:,j,ibelow) + &
                  LC(isapbelow) * bm_to_litter(:,j,isapbelow,icarbon) + &
                  LC(iheartbelow) * bm_to_litter(:,j,iheartbelow,icarbon) + &
                  LC(iroot) * bm_to_litter(:,j,iroot,icarbon) + &
                  LC(isapbelow)*turnover(:,j,isapbelow,icarbon) + &
                  LC(iheartbelow)*turnover(:,j,iheartbelow,icarbon) + &
                  LC(iroot)*turnover(:,j,iroot,icarbon)

          ENDIF

       ENDDO
    ENDDO

    !! 2.2.5 add new litter (struct/met, above/below)
    litter(:,:,:,:,:) = litter(:,:,:,:,:) + litter_inc(:,:,:,:,:)

    !! 2.2.6 for security: can't add more lignin than structural litter (above/below)
    DO l = 1, nlevs !Loop over litter levels (above and below ground)
       DO m = 1,nvm !Lopp over PFTs

          lignin_struc_inc(:,m,l) = &
               MIN( lignin_struc_inc(:,m,l), litter_inc(:,istructural,m,l,icarbon) )

       ENDDO
    ENDDO

    !! 2.2.7 new lignin content: add old lignin and lignin increase, divide by 
    !!       total structural litter (above/below)
    DO l = 1, nlevs !Loop over litter levels (above and below ground)
       DO m = 1,nvm !Loop over PFTs
          WHERE( litter(:,istructural,m,l,icarbon) .GT. min_stomate )

       !MM : Soenke modif
       ! Best vectorization ?
!!$       lignin_struc(:,:,:) = &
!!$            ( lignin_struc(:,:,:)*old_struc(:,:,:) + lignin_struc_inc(:,:,:) ) / &
!!$            litter(:,istructural,:,:,icarbon)

             lignin_struc(:,m,l) = lignin_struc(:,m,l) * old_struc(:,m,l)
             lignin_struc(:,m,l) = lignin_struc(:,m,l) + lignin_struc_inc(:,m,l)
             lignin_struc(:,m,l) = lignin_struc(:,m,l) / litter(:,istructural,m,l,icarbon)
          ELSEWHERE
             lignin_struc(:,m,l) = zero
          ENDWHERE
       ENDDO
    ENDDO

    
    !! 2.3 new litter fraction per PFT (for structural and metabolic litter, above
    !!       the ground).
    DO j = 1,nvm !Loop over PFTs

       WHERE ( litter(:,:,j,iabove,icarbon) .GT. min_stomate )

          litterpart(:,j,:) = &
               ( litter_pft(:,j,:,icarbon) + litter_inc_PFT(:,j,:,iabove,icarbon) ) / litter(:,:,j,iabove,icarbon)

       ELSEWHERE

          litterpart(:,j,:) = zero

       ENDWHERE

    ENDDO

    
  !! 3. Temperature control on decay: Factor between 0 and 1

    !! 3.1 above: surface temperature
    control_temp(:,iabove) = control_temp_func (npts, tsurf)

    
    !! 3.2 below: convolution of temperature and decomposer profiles
    !!            (exponential decomposer profile supposed)
   
    !! 3.2.1 rpc is an integration constant such that the integral of the root profile is 1.
    rpc(:) = un / ( un - EXP( -z_soil(nslm) / z_decomp ) )

    !! 3.2.2 integrate over the nslm levels
    tsoil_decomp(:) = zero

    DO l = 1, nslm

       tsoil_decomp(:) = &
            tsoil_decomp(:) + tsoil(:,l) * rpc(:) * &
            ( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )

    ENDDO

    control_temp(:,ibelow) = control_temp_func (npts, tsoil_decomp)

  !! 4. Moisture control. Factor between 0 and 1
    
    !! 4.1 above the ground: litter humidity

    IF (ok_moyano_soilhumsat)THEN !!ok_moyano_soilhumsat is true
       control_moist(:,iabove) = control_moist_func_moyano (npts,clay,bulk,carbon, litter ,litterhumSAT, veget_cov_max)
    ELSE !ok_moyano_soilhumsat is false by default
       control_moist(:,iabove) = control_moist_func (npts, litterhum)
    ENDIF

        CALL xios_orchidee_send_field("ControlMoistLitter",control_moist(:,iabove))
        CALL xios_orchidee_send_field("litterhumSAT",litterhumSAT(:))
    !
    !! 4.2 below: convolution of humidity and decomposer profiles
    !            (exponential decomposer profile supposed)

    !! 4.2.1 rpc is an integration constant such that the integral of the root profile is 1.
    rpc(:) = un / ( un - EXP( -z_soil(nslm) / z_decomp ) )

    !! 4.2.2 integrate over the nslm levels
    soilhum_decomp(:) = zero
    soilhumSAT_decomp(:) = zero

    IF (ok_moyano_soilhumsat)THEN !!ok_moyano_soilhumsat is true
      DO l = 1, nslm !Loop over soil levels
         soilhumSAT_decomp(:) = &
              soilhumSAT_decomp(:) + soilhumSAT(:,l) * rpc(:) * &
              ( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
      ENDDO
      control_moist(:,ibelow) = control_moist_func_moyano (npts,clay,bulk,carbon, litter, soilhumSAT_decomp,veget_cov_max)

     ELSE !ok_moyano_soilhumsat is false by default
      DO l = 1, nslm !Loop over soil levels

         soilhum_decomp(:) = &
              soilhum_decomp(:) + soilhum(:,l) * rpc(:) * &
              ( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
      ENDDO
      control_moist(:,ibelow) = control_moist_func (npts, soilhum_decomp)
    ENDIF

        CALL xios_orchidee_send_field("ControlMoistSoil",control_moist(:,ibelow))
        CALL xios_orchidee_send_field("soilhumSAT",soilhumSAT_decomp(:))

  !! 5. fluxes from litter to carbon pools and respiration

    DO l = 1, nlevs !Loop over litter levels (above and below ground)
       DO m = 1,nvm !Loop over PFTs

          !! 5.1 structural litter: goes into active and slow carbon pools + respiration

          !! 5.1.1 total quantity of structural litter which is decomposed
          fd(:) = dt/litter_tau(istructural) * &
               control_temp(:,l) * control_moist(:,l) * exp( -litter_struct_coef * lignin_struc(:,m,l) )

          DO k = 1,nelements
             
             qd(:,k) = litter(:,istructural,m,l,k) * fd(:)

          END DO

          litter(:,istructural,m,l,:) = litter(:,istructural,m,l,:) - qd(:,:)

          !! 5.1.2 decompose same fraction of structural part of dead leaves. Not exact
          !!       as lignine content is not the same as that of the total structural litter.
          ! to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
          ! we do this test to do this calcul only ones in 1,nlev loop
          if (l == iabove)  dead_leaves(:,m,istructural) = dead_leaves(:,m,istructural) * ( un - fd(:) )

          !! 5.1.3 non-lignin fraction of structural litter goes into
          !!       active carbon pool + respiration
          soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
               frac_soil(istructural,iactive,l) * qd(:,icarbon) * ( 1. - lignin_struc(:,m,l) ) / dt

      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
      ! the sechiba time step but then in stomate.f90 resp_hetero_litter 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.
          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
               ( 1. - frac_soil(istructural,iactive,l) ) * qd(:,icarbon) * &
               ( 1. - lignin_struc(:,m,l) ) / dt

          !! 5.1.4 lignin fraction of structural litter goes into
          !!       slow carbon pool + respiration
          soilcarbon_input(:,islow,m) = soilcarbon_input(:,islow,m) + &
               frac_soil(istructural,islow,l) * qd(:,icarbon) * lignin_struc(:,m,l) / dt

      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
      ! the sechiba time step but then in stomate.f90 resp_hetero_litter 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.
          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
               ( 1. - frac_soil(istructural,islow,l) ) * qd(:,icarbon) * lignin_struc(:,m,l) / dt

          
          !! 5.2 metabolic litter goes into active carbon pool + respiration
         
          !! 5.2.1 total quantity of metabolic litter that is decomposed
          fd(:) = dt/litter_tau(imetabolic) * control_temp(:,l) * control_moist(:,l)

          DO k = 1,nelements
          
             qd(:,k) = litter(:,imetabolic,m,l,k) * fd(:)

          END DO

          litter(:,imetabolic,m,l,:) = litter(:,imetabolic,m,l,:) - qd(:,:)

          !! 5.2.2 decompose same fraction of metabolic part of dead leaves.
          !  to avoid a multiple (for ibelow and iabove) modification of dead_leaves,
          !  we do this test to do this calcul only ones in 1,nlev loop
          if (l == iabove)  dead_leaves(:,m,imetabolic) = dead_leaves(:,m,imetabolic) * ( 1. - fd(:) )


          !! 5.2.3 put decomposed litter into carbon pool + respiration
          soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
               frac_soil(imetabolic,iactive,l) * qd(:,icarbon) / dt

      !BE CAREFUL: Here resp_hetero_litter is divided by dt to have a value which corresponds to
      ! the sechiba time step but then in stomate.f90 resp_hetero_litter 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.
          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
               ( 1. - frac_soil(imetabolic,iactive,l) ) * qd(:,icarbon) / dt

       ENDDO
    ENDDO

    
  !! 6. calculate fraction of total soil covered by dead leaves

    CALL deadleaf (npts, veget_cov_max, dead_leaves, deadleaf_cover)

  !! 7. (Quasi-)Analytical Spin-up : Start filling MatrixA

    IF (spinup_analytic) THEN

       MatrixA(:,:,:,:) = zero
       VectorB(:,:,:) = zero
       
       
       DO m = 1,nvm

          !- MatrixA : carbon fluxes leaving the litter
          
          MatrixA(:,m,istructural_above,istructural_above)= - dt/litter_tau(istructural) * &
               control_temp(:,iabove) * control_moist(:,iabove) * exp( -litter_struct_coef * lignin_struc(:,m,iabove) )
          
          MatrixA(:,m,istructural_below,istructural_below) = - dt/litter_tau(istructural) * &
               control_temp(:,ibelow) * control_moist(:,ibelow) * exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )
          
          MatrixA(:,m,imetabolic_above,imetabolic_above) = - dt/litter_tau(imetabolic) * & 
               control_temp(:,iabove) * control_moist(:,iabove)
          
          MatrixA(:,m,imetabolic_below,imetabolic_below) = - dt/litter_tau(imetabolic) * & 
               control_temp(:,ibelow) * control_moist(:,ibelow)
          
          
          !- MatrixA : carbon fluxes between the litter and the pools (the rest of the matrix is filled in stomate_soilcarbon.f90)
          
          MatrixA(:,m,iactive_pool,istructural_above) = frac_soil(istructural,iactive,iabove) * &
               dt/litter_tau(istructural) * &                    
               control_temp(:,iabove) * control_moist(:,iabove) * & 
               exp( -litter_struct_coef * lignin_struc(:,m,iabove) ) * &
               ( 1. - lignin_struc(:,m,iabove) ) 
          

          MatrixA(:,m,iactive_pool,istructural_below) = frac_soil(istructural,iactive,ibelow) * &
               dt/litter_tau(istructural) * &                     
               control_temp(:,ibelow) * control_moist(:,ibelow) * & 
               exp( -litter_struct_coef * lignin_struc(:,m,ibelow) ) * &
               ( 1. - lignin_struc(:,m,ibelow) ) 
          
          
          MatrixA(:,m,iactive_pool,imetabolic_above) =  frac_soil(imetabolic,iactive,iabove) * &
               dt/litter_tau(imetabolic) * control_temp(:,iabove) * control_moist(:,iabove) 
          
          
          MatrixA(:,m,iactive_pool,imetabolic_below) =  frac_soil(imetabolic,iactive,ibelow) * &
               dt/litter_tau(imetabolic) * control_temp(:,ibelow) * control_moist(:,ibelow)
          
          MatrixA(:,m,islow_pool,istructural_above) = frac_soil(istructural,islow,iabove) * &
               dt/litter_tau(istructural) * &                   
               control_temp(:,iabove) * control_moist(:,iabove) * &
               exp( -litter_struct_coef * lignin_struc(:,m,iabove) )* &
               lignin_struc(:,m,iabove)
          
          
          MatrixA(:,m,islow_pool,istructural_below) = frac_soil(istructural,islow,ibelow) * &
               dt/litter_tau(istructural) * &   
               control_temp(:,ibelow) * control_moist(:,ibelow) *  &
                  exp( -litter_struct_coef * lignin_struc(:,m,ibelow) )* &
                  lignin_struc(:,m,ibelow) 
          
          
          !- VectorB : carbon input -
          
          VectorB(:,m,istructural_above) = litter_inc_PFT(:,m,istructural,iabove,icarbon)
          VectorB(:,m,istructural_below) = litter_inc_PFT(:,m,istructural,ibelow,icarbon)
          VectorB(:,m,imetabolic_above) = litter_inc_PFT(:,m,imetabolic,iabove,icarbon)
          VectorB(:,m,imetabolic_below) = litter_inc_PFT(:,m,imetabolic,ibelow,icarbon)
          
          IF (printlev>=4) WRITE(numout,*) 'We filled MatrixA and VectorB' 
          
       ENDDO ! Loop over # PFTs
          
    ENDIF ! spinup analytic 

    IF (printlev>=4) WRITE(numout,*) 'Leaving littercalc'

  END SUBROUTINE littercalc


!! ==============================================================================================================================\n
!! SUBROUTINE   : deadleaf
!!
!>\BRIEF        This routine calculates the deadleafcover. 
!!
!! DESCRIPTION  : It first calculates the lai corresponding to the dead leaves (LAI) using 
!! the dead leaves carbon content (DL) the specific leaf area (sla) and the 
!! maximal coverage fraction of a PFT (vegetmax) using the following equations:
!! \latexonly
!! \input{deadleaf1.tex}
!! \endlatexonly
!! \n
!! Then, the dead leaf cover (DLC) is calculated as following:\n
!! \latexonly
!! \input{deadleaf2.tex}
!! \endlatexonly
!! \n
!! 
!! RECENT CHANGE(S) : None
!!
!! MAIN OUTPUT VARIABLE: ::deadleaf_cover
!! 
!! REFERENCE(S) : None
!!
!! FLOWCHART : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE deadleaf (npts, veget_cov_max, dead_leaves, deadleaf_cover)

  !! 0. Variable and parameter declaration
    
    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                          :: npts           !! Domain size - number of grid pixels (unitless)
    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in)  :: dead_leaves    !! Dead leaves per ground unit area, per PFT, 
                                                                          !! metabolic and structural  
                                                                          !! @tex $(gC m^{-2})$ @endtex
    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)          :: veget_cov_max  !! PFT "Maximal" coverage fraction of a PFT defined in 
                                                                          !! the input vegetation map 
                                                                          !! @tex $(m^2 m^{-2})$ @endtex 
    
    !! 0.2 Output variables
    
    REAL(r_std), DIMENSION(npts), INTENT(out)           :: deadleaf_cover !! Fraction of soil covered by dead leaves over all PFTs
                                                                          !! (0-1, unitless)

    !! 0.3 Modified variables

    !! 0.4 Local variables

    REAL(r_std), DIMENSION(npts)                        :: dead_lai       !! LAI of dead leaves @tex $(m^2 m^{-2})$ @endtex
    INTEGER(i_std)                                      :: j              !! Index (unitless)
!_ ================================================================================================================================
    
  !! 1. LAI of dead leaves
  
    dead_lai(:) = zero

    DO j = 1,nvm !Loop over PFTs
       dead_lai(:) = dead_lai(:) + ( dead_leaves(:,j,imetabolic) + dead_leaves(:,j,istructural) ) * sla(j) &
            * veget_cov_max(:,j)
    ENDDO

  !! 2. fraction of soil covered by dead leaves

    deadleaf_cover(:) = un - exp( - 0.5 * dead_lai(:) )

    IF (printlev>=4) WRITE(numout,*) 'Leaving deadleaf'

  END SUBROUTINE deadleaf


!! ================================================================================================================================
!! FUNCTION     : control_moist_func_moyano
!!
!>\BRIEF        Calculate moisture control for litter and soil C decomposition
!!
!! DESCRIPTION  : Calculate moisture control factor applied
!! to litter decomposition and to soil carbon decomposition in
!! stomate_soilcarbon.f90 using the following equation based on the Moyano et al., 2012 BG paper: \n
!! \latexonly
!! \input{control_moist_func1.tex}
!! \endlatexonly
!! \n
!! with M the moisture control factor and soilmoisutre, the soil moisture 
!! calculated in sechiba.
!! Then, the function is ranged between Moistcont_min and 1:\n
!! \latexonly
!! \input{control_moist_func2.tex}
!! \endlatexonly
!! \n
!! RECENT CHANGE(S) : None
!!
!! RETURN VALUE : ::moistfunc_result
!! 
!! REFERENCE(S) : None
!!
!! FLOWCHART : None
!! \n
!_ ================================================================================================================================
  
  FUNCTION control_moist_func_moyano (npts,clay,bulk,carbon,litter,moist_in,veget_cov_max) RESULT (moistfunc_result)

  !! 0. Variable and parameter declaration
    
    !! 0.1 Input variables
    INTEGER(i_std)                           :: k,m,n,i,l,e         !!indices              
    INTEGER(i_std), INTENT(in)               :: npts                !! Domain size - number of grid pixel (unitless)
    REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in            !! relative humidity (unitless)
    REAL(r_std), DIMENSION(npts), INTENT(in) :: clay                !! Clay fraction (unitless, 0-1) 
    REAL(r_std), DIMENSION(npts), INTENT(in) :: bulk                !! Bulk density (kg/m**3)  
    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in) :: carbon    !! Soil carbon pools: active, s     low, or passive, \f$(gC m^{2})$\f
    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) ::litter   !! Metabolic and structural litter,above and
                                                                    !! below ground. The unit is given by m^2 of ground
                                                                    !! @tex $(gC m^{-2})$ @endtex
    REAL(r_std),DIMENSION(npts,nvm),INTENT(in):: veget_cov_max      !! PFT "Maximal" coverage fraction of a PFT 
 
    !! 0.2 Output variables
   
    REAL(r_std), DIMENSION(npts)             :: moistfunc_result    !! Moisture control factor (0.25-1, unitless)

    !! 0.3 Modified variables

    !! 0.4 Local variables
    REAL(r_std), DIMENSION(101)              :: Mint                !! range [0,0.01,1] to compute PRSRmax
    REAL(r_std), DIMENSION(npts,101)         :: PRSR                !! Proportional Response of Soil Respiration  
    REAL(r_std), DIMENSION(npts,101)         :: SR                  !! Soil respiration compute using Moyano et al 2012 procedure  
    REAL(r_std)                              :: SRmax               !! maximum value of SR
    REAL(r_std), DIMENSION(npts,101)         :: SRnorm              !! SR normalized by the SRmax 
    REAL(r_std)                              :: ind                 !! i index of maximum SRnorm value
    REAL(r_std)                              :: SRmin               !! minimum value of SRnorm
    REAL(r_std), DIMENSION(npts,101)         :: SRsc                !! Soil respiration rescale between 0 and 1 following Moyano et al 2012 procedure  
    REAL(r_std)                              :: SRscmax             !! maximum value of SRsc
    REAL(r_std)                              :: indmc               !! i index of maximum mois_round values 
    REAL(r_std)                              :: moist_round         !! moist_in round to 2decimal points
    REAL(r_std), DIMENSION(npts)             :: carbon_gCgs         !! total Soil carbon pools: active +slow+ passive, (gC/gsoil)
    REAL(r_std), DIMENSION(npts)             :: carbon_gCm2s        !! total Soil carbon pools: active +slow+ passive     , (gC/gsoil)
    REAL(r_std), DIMENSION(npts)             :: clay_Moyano         !! Clay fraction borned within the values used by Moyano et al., (unitless, 0-1)

!_ ===============================================================================================================================

  !! In ok_moyano_soilhumsat total soil carbon per grid cells need to be computed and
  !!  converted from gC/m2soil to gC/gsoil
  carbon_gCm2s(:)=zero
  carbon_gCgs(:)=zero
  DO n = 1, npts ! loop over grids
   DO m = 1, nvm ! Loop over # PFTs
    DO k = 1, ncarb ! Loop over carbon pools 
     carbon_gCm2s(n)=carbon_gCm2s(n)+carbon(n,k,m)*veget_cov_max(n,m)
    ENDDO
   ENDDO
  ENDDO

  IF (ok_orga) THEN
     DO n = 1, npts ! loop over grids
       carbon_gCgs(n)=carbon_gCm2s(n) / (bulk(n) * 1E3 * soilheight)
       carbon_gCgs(n)= MAX(min_carbon_moyano,MIN(max_carbon_moyano,carbon_gCgs(n))) 
       clay_Moyano(n) = MAX(min_clay_moyano,MIN(max_clay_moyano,clay(n)))
    ENDDO

     !!1. compute Max(Prsr(M[0,0.01,1])); Prsr:Proportional Response of Soil
     !!Respiration (PRSR) and Soil respiration(SR): SR(M)=SRo x Prsr(M)/Max(Prsr(M[0,0.01,1]))
     !! M:soilhumSAT, SRo=1 (arbitrary)
     DO n = 1,npts
       Mint(:)= zero
       SRmax = zero
       DO i = 1,101
          Mint(i)=0.01*(i-1)
          IF (carbon_gCgs(n) .LT. limit_carbon_orga) THEN
              PRSR(n,i) = beta1 * Mint(i)  + beta2 * Mint(i)**2.0 + beta3 * Mint(i)**3.0 &
                         & + beta4 * clay_Moyano(n) + beta5 * clay_Moyano(n) * Mint(i) &
                         & + beta6 * carbon_gCgs(n) + intercept
          ELSE 
              PRSR(n,i) = beta1_orga * Mint(i)  + beta2_orga * Mint(i)**2.0 + &
                         & beta3_orga * Mint(i)**3.0 + intercept_orga
          ENDIF

         IF (i.LT.2) THEN
            SR(n,i) = PRSR(n,i)
            SRmax = MAX(SR(n,i),SRmax)
         ELSE
            SR(n,i) = PRSR(n,i) * SR(n,i-1)
            SRmax = MAX(SR(n,i),SRmax)
          ENDIF
       ENDDO
       SRnorm(n,:)= SRo * SR(n,:)/SRmax
     ENDDO
  ELSE
     DO n = 1, npts ! loop over grids
       carbon_gCgs(n)=carbon_gCm2s(n) / (bulk(n) * 1E3 * soilheight)
       carbon_gCgs(n)= MAX(min_carbon_moyano,MIN(limit_carbon_orga,carbon_gCgs(n)))
       clay_Moyano(n) = MAX(min_clay_moyano,MIN(max_clay_moyano,clay(n)))
    ENDDO

     !!1. compute Max(Prsr(M[0,0.01,1])); Prsr:Proportional Response of Soil
     !!Respiration (PRSR) and Soil respiration(SR): SR(M)=SRo x
     !Prsr(M)/Max(Prsr(M[0,0.01,1]))
     !! M:soilhumSAT, SRo=1 (arbitrary)
     DO n = 1,npts
       Mint(:)= zero
       SRmax = zero
       DO i = 1,101
          Mint(i)=0.01*(i-1)
           PRSR(n,i) = beta1 * Mint(i)  + beta2 * Mint(i)**2.0 + beta3 * Mint(i)**3.0 &
                      & + beta4 * clay_Moyano(n) + beta5 * clay_Moyano(n) * Mint(i) &
                      & + beta6 * carbon_gCgs(n) + intercept

          IF (i.LT.2) THEN
            SR(n,i) = PRSR(n,i)
            SRmax = MAX(SR(n,i),SRmax)
          ELSE
            SR(n,i) = PRSR(n,i) * SR(n,i-1)
            SRmax = MAX(SR(n,i),SRmax)
          ENDIF
       ENDDO
       SRnorm(n,:)= SRo * SR(n,:)/SRmax
     ENDDO
  ENDIF   

  !!2. Rescale SR values between 0 and 1 and defined SR for moist_in
    DO n = 1,npts
      ind = 1.0
      indmc = zero
      moist_round = zero
      SRscmax = zero
      SRsc(n,:) = zero
      SRmin = 1.0
         
      ind = MAXLOC(SRnorm(n,:),1)
      DO i = 1,ind
        SRmin = MIN(SRnorm(n,i),SRmin)
      ENDDO
      DO i = 1,101
        SRsc(n,i) = SRnorm(n,i)- SRmin
        IF (i.LE.ind)THEN
          SRscmax = MAX(SRsc(n,i),SRscmax)
        ENDIF
      ENDDO
      SRsc(n,:)=SRsc(n,:)/SRscmax
         
      moist_round=(NINT(moist_in(n)*100.))/100.
      indmc=INT(moist_round/0.01) + un
        
     moistfunc_result(n) = SRsc(n,indmc)
     moistfunc_result(n) = MAX( moistcontSAT_min, MIN( un, moistfunc_result(n) ))
    ENDDO

  END FUNCTION control_moist_func_moyano

!! ================================================================================================================================
!! FUNCTION     : control_moist_func
!!
!>\BRIEF        Calculate moisture control for litter and soil C decomposition
!!
!! DESCRIPTION  : Calculate moisture control factor applied
!! to litter decomposition and to soil carbon decomposition in
!! stomate_soilcarbon.f90 using the following equation: \n
!! \latexonly
!! \input{control_moist_func1.tex}
!! \endlatexonly
!! \n
!! with M the moisture control factor and soilmoisutre, the soil moisture
!! calculated in sechiba.
!! Then, the function is ranged between Moistcont_min and 1:\n
!! \latexonly
!! \input{control_moist_func2.tex}
!! \endlatexonly
!! \n
!! RECENT CHANGE(S) : None
!!
!! RETURN VALUE : ::moistfunc_result
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART : None
!! \n
!_ ================================================================================================================================
     
  FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)

  !! 0. Variable and parameter declaration
       
  !! 0.1 Input variables
              
   INTEGER(i_std), INTENT(in)               :: npts                !! Domain size - number of grid pixel (unitless)
   REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in            !! relative humidity (unitless)
     
  !! 0.2 Output variables
        
  REAL(r_std), DIMENSION(npts)             :: moistfunc_result    !! Moisture control factor (0.25-1, unitless)
   
  !! 0.3 Modified variables
    
  !! 0.4 Local variables
    
!_ ================================================================================================================================
     
    moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
    moistfunc_result(:) = MAX( moistcont_min, MIN( un, moistfunc_result(:) ) )
     
  END FUNCTION control_moist_func

!! ================================================================================================================================
!! FUNCTION     : control_temp_func
!!
!>\BRIEF        Calculate temperature control for litter and soild C decomposition
!!
!! DESCRIPTION  : Calculate temperature control factor applied
!! to litter decomposition and to soil carbon decomposition in
!! stomate_soilcarbon.f90 using the following equation: \n
!! \latexonly
!! \input{control_temp_func1.tex}
!! \endlatexonly
!! \n
!! with T the temperature control factor, temp the temperature in Kelvin of 
!! the air (for aboveground litter) or of the soil (for belowground litter 
!! and soil)
!! Then, the function is limited in its maximal range to 1:\n
!! \latexonly
!! \input{control_temp_func2.tex}
!! \endlatexonly
!! \n
!! RECENT CHANGE(S) : None
!!
!! RETURN VALUE: ::tempfunc_result
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART : None
!! \n
!_ ================================================================================================================================

  FUNCTION control_temp_func (npts, temp_in) RESULT (tempfunc_result)

  !! 0. Variable and parameter declaration
    
    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                 :: npts            !! Domain size - number of land pixels (unitless)
    REAL(r_std), DIMENSION(npts), INTENT(in)   :: temp_in         !! Temperature (K)

    !! 0.2 Output variables
    REAL(r_std), DIMENSION(npts)               :: tempfunc_result !! Temperature control factor (0-1, unitless)

    !! 0.3 Modified variables

    !! 0.4 Local variables

!_ ================================================================================================================================

    tempfunc_result(:) = exp( soil_Q10 * ( temp_in(:) - (ZeroCelsius+tsoil_ref)) / Q10 )
    tempfunc_result(:) = MIN( un, tempfunc_result(:) )

  END FUNCTION control_temp_func

END MODULE stomate_litter
