! Update litter and lignine content after litter fall.
! Calculate litter decomposition.
!
! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/stomate_litter.f90,v 1.9 2009/06/24 10:43:21 ssipsl Exp $
! IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
MODULE stomate_litter

  ! modules used:

  USE ioipsl
  USE stomate_data
  USE constantes
  USE pft_parameters

  IMPLICIT NONE

  ! private & public routines

  PRIVATE
  PUBLIC littercalc,littercalc_clear, deadleaf

  ! first call
  LOGICAL, SAVE                                                     :: firstcall = .TRUE.


CONTAINS

  SUBROUTINE littercalc_clear
    firstcall =.TRUE.
  END SUBROUTINE littercalc_clear


  SUBROUTINE littercalc (npts, dt, &
       turnover, bm_to_litter, &
       veget_max, tsurf, tsoil, soilhum, litterhum, &
       litterpart, litter, dead_leaves, lignin_struc, &
       deadleaf_cover, resp_hetero_litter, &
       soilcarbon_input, control_temp, control_moist)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                               :: npts
    ! time step in days
    REAL(r_std), INTENT(in)                                            :: dt
    ! Turnover rates (gC/(m**2 of ground)/day)
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)               :: turnover
    ! conversion of biomass to litter (gC/(m**2 of ground)) / day
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)               :: bm_to_litter
    ! veget_max
    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                  :: veget_max
    ! temperature (K) at the surface
    REAL(r_std), DIMENSION(npts), INTENT(in)                           :: tsurf
    ! soil temperature (K)
    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)                      :: tsoil
    ! daily soil humidity
    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)                      :: soilhum
    ! daily litter humidity
    REAL(r_std), DIMENSION(npts), INTENT(in)                           :: litterhum

    ! 0.2 modified fields

    ! fraction of litter above the ground belonging to different PFTs
    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)             :: litterpart
    ! metabolic and structural litter,above and below ground (gC/m**2 of ground)
    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs), INTENT(inout)  :: litter
    ! dead leaves on ground, per PFT, metabolic and structural,
    !   in gC/(m**2 of ground)
    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)             :: dead_leaves
    ! ratio Lignine/Carbon in structural litter, above and below ground, (gC/m**2)
    REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout)        :: lignin_struc

    ! 0.3 output

    ! fraction of soil covered by dead leaves
    REAL(r_std), DIMENSION(npts), INTENT(out)                          :: deadleaf_cover
    ! litter heterotrophic respiration (in gC/day/m**2 of ground)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                :: resp_hetero_litter
    ! quantity of carbon going into carbon pools from litter decomposition
    !   (gC/(m**2 of ground)/day)
    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(out)          :: soilcarbon_input
    ! temperature control of heterotrophic respiration, above and below
    REAL(r_std), DIMENSION(npts,nlevs), INTENT(out)                    :: control_temp
    ! moisture control of heterotrophic respiration
    REAL(r_std), DIMENSION(npts,nlevs), INTENT(out)                    :: control_moist

    ! 0.4 local

    ! what fraction of leaves, wood, etc. goes into metabolic and structural litterpools
    REAL(r_std), SAVE, DIMENSION(nparts,nlitt)                         :: litterfrac
    ! soil levels (m)
    REAL(r_std), SAVE, DIMENSION(0:nbdl)                               :: z_soil
    ! integration constant for vertical profiles
    REAL(r_std), DIMENSION(npts)                                       :: rpc
    ! residence time in litter pools (days)
    REAL(r_std), SAVE, DIMENSION(nlitt)                                :: litter_tau
    ! decomposition flux fraction that goes into soil (litter -> carbon, above and below)
    !   rest goes into atmosphere
    REAL(r_std), SAVE, DIMENSION(nlitt,ncarb,nlevs)                    :: frac_soil
    ! temperature used for decompostition in soil (K)
    REAL(r_std), DIMENSION(npts)                                       :: tsoil_decomp
    ! humidity used for decompostition in soil
    REAL(r_std), DIMENSION(npts)                                       :: soilhum_decomp
    ! fraction of structural or metabolic litter decomposed
    REAL(r_std), DIMENSION(npts)                                       :: fd
    ! quantity of structural or metabolic litter decomposed (gC/m**2)
    REAL(r_std), DIMENSION(npts)                                       :: qd
    ! old structural litter, above and below (gC/m**2)
    REAL(r_std), DIMENSION(npts,nvm,nlevs)                       :: old_struc
    ! increase of litter, per PFT, metabolic and structural,
    !   above and below ground (gC/m**2 of ground)
    REAL(r_std), DIMENSION(npts,nvm,nlitt,nlevs)                      :: litter_inc_PFT
    ! increase of metabolic and structural litter, above and below ground (gC/m**2 of ground)
    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs)                 :: litter_inc
    ! lignin increase in structural litter, above and below ground (gC/m**2 of ground)
    REAL(r_std), DIMENSION(npts,nvm,nlevs)                       :: lignin_struc_inc
    ! metabolic and structural litter above the ground per PFT
    REAL(r_std), DIMENSION(npts,nvm,nlitt)                            :: litter_pft
    ! intermediate array for looking for minimum
    REAL(r_std), DIMENSION(npts)                                       :: zdiff_min
    ! for messages
    CHARACTER(LEN=10), DIMENSION(nlitt)                                    :: litter_str
    CHARACTER(LEN=22), DIMENSION(nparts)                                   :: part_str
    CHARACTER(LEN=7), DIMENSION(ncarb)                                     :: carbon_str
    CHARACTER(LEN=5), DIMENSION(nlevs)                                     :: level_str
    ! Indices
    INTEGER(i_std)                                                    :: i,j,k,l,m

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering littercalc'

    !
    ! 1 Initialisations
    !

    IF ( firstcall ) THEN

       !
       ! 1.1 get soil "constants"
       !

       ! 1.1.1 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) = 1. - litterfrac(k,imetabolic)

       ENDDO

       ! 1.1.2 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.3 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
       !

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

       !
       ! 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'

       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

       firstcall = .FALSE.

    ENDIF

    !
    ! 1.3 litter above the ground per PFT.
    !

    DO j = 2, nvm

       DO k = 1, nlitt
          litter_pft(:,j,k) = litterpart(:,j,k) * litter(:,k,j,iabove)
       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
       DO m = 2,nvm

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

       ENDDO
    ENDDO

    !
    ! 2.2 update litter, dead leaves, and lignin content in structural litter
    !

    litter_inc(:,:,:,:) = zero
    lignin_struc_inc(:,:,:) = zero

    DO j = 2,nvm

       ! 2.2.1 litter

       DO k = 1, nlitt    ! 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) = &
               litterfrac(ileaf,k) * bm_to_litter(:,j,ileaf) + &
               litterfrac(isapabove,k) * bm_to_litter(:,j,isapabove) + &
               litterfrac(iheartabove,k) * bm_to_litter(:,j,iheartabove) + &
               litterfrac(ifruit,k) * bm_to_litter(:,j,ifruit) + &
               litterfrac(icarbres,k) * bm_to_litter(:,j,icarbres) + &
               litterfrac(ileaf,k) * turnover(:,j,ileaf) + &
               litterfrac(isapabove,k) * turnover(:,j,isapabove) + &
               litterfrac(iheartabove,k) * turnover(:,j,iheartabove) + &
               litterfrac(ifruit,k) * turnover(:,j,ifruit) + &
               litterfrac(icarbres,k) * turnover(:,j,icarbres)

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

          ! litter increase, met/struct, above/below

          litter_inc(:,k,j,iabove) = litter_inc(:,k,j,iabove) + litter_inc_PFT(:,j,k,iabove)
          litter_inc(:,k,j,ibelow) = litter_inc(:,k,j,ibelow) + litter_inc_PFT(:,j,k,ibelow)

          ! 2.2.3 dead leaves, for soil cover.

          dead_leaves(:,j,k) = &
               dead_leaves(:,j,k) + &
               litterfrac(ileaf,k) * ( bm_to_litter(:,j,ileaf) + turnover(:,j,ileaf) )

          ! 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) + &
                  LC(isapabove) * bm_to_litter(:,j,isapabove) + &
                  LC(iheartabove) * bm_to_litter(:,j,iheartabove) + &
                  LC(ifruit) * bm_to_litter(:,j,ifruit) + &
                  LC(icarbres) * bm_to_litter(:,j,icarbres) + &
                  LC(ileaf) * turnover(:,j,ileaf) + &
                  LC(isapabove) * turnover(:,j,isapabove) + &
                  LC(iheartabove) * turnover(:,j,iheartabove) + &
                  LC(ifruit) * turnover(:,j,ifruit) + &
                  LC(icarbres) * turnover(:,j,icarbres)

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

          ENDIF

       ENDDO
    ENDDO

    ! 3.2.5 add new litter (struct/met, above/below)

    litter(:,:,:,:) = litter(:,:,:,:) + litter_inc(:,:,:,:)

    ! 3.2.6 for security: can't add more lignin than structural litter (above/below)

    DO l = 1, nlevs
       DO m = 2,nvm

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

       ENDDO
    ENDDO

    ! 3.2.7 new lignin content: add old lignin and lignin increase, divide by 
    !       total structural litter (above/below)

    DO l = 1, nlevs
       DO m = 2,nvm
          WHERE( litter(:,istructural,m,l) .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)
          ELSEWHERE
             lignin_struc(:,m,l) = zero
          ENDWHERE
       ENDDO
    ENDDO

    !
    ! 3.3 new litter fraction per PFT (for structural and metabolic litter, above
    !       the ground).
    !

    DO j = 2,nvm

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

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

       ELSEWHERE

          litterpart(:,j,:) = zero

       ENDWHERE

    ENDDO

    !
    ! 4 Temperature control on decay: Factor between 0 and 1
    !

    !
    ! 4.1 above: surface temperature
    !

    control_temp(:,iabove) = control_temp_func (npts, tsurf)

    !
    ! 4.2 below: convolution of temperature 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(nbdl) / z_decomp ) )

    ! 4.2.2 integrate over the nbdl levels

    tsoil_decomp(:) = zero

    DO l = 1, nbdl

       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)

    !
    ! 5 Moisture control. Factor between 0 and 1
    !

    !
    ! 5.1 above the ground: litter humidity
    !

    control_moist(:,iabove) = control_moist_func (npts, litterhum)

    !
    ! 5.2 below: convolution of humidity and decomposer profiles
    !            (exponential decomposer profile supposed)
    !

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

    ! 5.2.2 integrate over the nbdl levels

    soilhum_decomp(:) = zero

    DO l = 1, nbdl

       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)

    !
    ! 6 fluxes from litter to carbon pools and respiration
    !

    DO l = 1, nlevs
       DO m = 2,nvm

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

          ! 6.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) )

          qd(:) = litter(:,istructural,m,l) * fd(:)

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

          ! 6.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(:) )

          ! 6.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(:) * ( 1. - lignin_struc(:,m,l) ) / dt

          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
               ( 1. - frac_soil(istructural,iactive,l) ) * qd(:) * &
               ( 1. - lignin_struc(:,m,l) ) / dt

          ! 6.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(:) * lignin_struc(:,m,l) / dt

          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
               ( 1. - frac_soil(istructural,islow,l) ) * qd(:) * lignin_struc(:,m,l) / dt

          !
          ! 6.2 metabolic litter goes into active carbon pool + respiration
          !

          ! 6.2.1 total quantity of metabolic litter that is decomposed

          fd(:) = dt/litter_tau(imetabolic) * control_temp(:,l) * control_moist(:,l)

          qd(:) = litter(:,imetabolic,m,l) * fd(:)

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

          ! 6.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(:) )


          ! 6.2.3 put decomposed litter into carbon pool + respiration

          soilcarbon_input(:,iactive,m) = soilcarbon_input(:,iactive,m) + &
               frac_soil(imetabolic,iactive,l) * qd(:) / dt

          resp_hetero_litter(:,m) = resp_hetero_litter(:,m) + &
               ( 1. - frac_soil(imetabolic,iactive,l) ) * qd(:) / dt

       ENDDO
    ENDDO

    !
    ! 7 calculate fraction of total soil covered by dead leaves
    !

    CALL deadleaf (npts, veget_max, dead_leaves, deadleaf_cover)

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving littercalc'

  END SUBROUTINE littercalc

  SUBROUTINE deadleaf (npts, veget_max, dead_leaves, deadleaf_cover)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                               :: npts
    ! dead leaves on ground, per PFT, metabolic and structural,
    !   in gC/(m**2 of ground)
    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(in)                :: dead_leaves
    !veget_max
    REAL(r_std),DIMENSION(npts,nvm),INTENT(in)                 :: veget_max
    ! 0.2. output
    ! fraction of soil covered by dead leaves
    REAL(r_std), DIMENSION(npts), INTENT(out)                          :: deadleaf_cover

    ! 0.3. local

    ! LAI of dead leaves
    REAL(r_std), DIMENSION(npts)                                       :: dead_lai
    ! Index
    INTEGER(i_std)                                                    :: j

    !
    ! 1 LAI of dead leaves
    !

    dead_lai(:) = zero

    DO j = 2,nvm
       dead_lai(:) = dead_lai(:) + ( dead_leaves(:,j,imetabolic) + dead_leaves(:,j,istructural) ) * sla(j) &
            * veget_max(:,j)
    ENDDO

    !
    ! 2 fraction of soil covered by dead leaves
    !

    deadleaf_cover(:) = 1. - exp( - 0.5 * dead_lai(:) )

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving deadleaf'

  END SUBROUTINE deadleaf

  FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                               :: npts
    ! relative humidity
    REAL(r_std), DIMENSION(npts), INTENT(in)                           :: moist_in

    ! 0.2 result

    ! moisture control factor
    REAL(r_std), DIMENSION(npts)                                       :: moistfunc_result

    moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
    moistfunc_result(:) = MAX( 0.25_r_std, MIN( un, moistfunc_result(:) ) )

  END FUNCTION control_moist_func

  FUNCTION control_temp_func (npts, temp_in) RESULT (tempfunc_result)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                               :: npts
    ! temperature (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                           :: temp_in

    ! 0.2 result

    ! temperature control factor
    REAL(r_std), DIMENSION(npts)                                       :: tempfunc_result


    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