! ================================================================================================================================
!  MODULE       : diffuco
!
!  CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
!  LICENCE      : IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF   This module calculates the limiting coefficients, both aerodynamic
!! and hydrological, for the turbulent heat fluxes.
!!
!!\n DESCRIPTION: The aerodynamic resistance R_a is used to limit
!! the transport of fluxes from the surface layer of vegetation to the point in the atmosphere at which
!! interaction with the LMDZ atmospheric circulation model takes place. The aerodynamic resistance is
!! calculated either within the module r_aerod (if the surface drag coefficient is provided by the LMDZ, and 
!! if the flag 'ldq_cdrag_from_gcm' is set to TRUE) or r_aero (if the surface drag coefficient must be calculated).\n
!!
!! Within ORCHIDEE, evapotranspiration is a function of the Evaporation Potential, but is modulated by a
!! series of resistances (canopy and aerodynamic) of the surface layer, here represented by beta.\n
!!
!! DESCRIPTION	:
!! \latexonly 
!!     \input{diffuco_intro.tex}
!! \endlatexonly
!! \n
!!
!! This module calculates the beta for several different scenarios: 
!! - diffuco_snow calculates the beta coefficient for sublimation by snow, 
!! - diffuco_inter calculates the beta coefficient for interception loss by each type of vegetation, 
!! - diffuco_bare calculates the beta coefficient for bare soil, 
!! - diffuco_trans_co2 calculates the beta coefficient for transpiration for each type of vegetation, using Farqhar's formula
!! - chemistry_bvoc calculates the beta coefficient for emissions of biogenic compounds \n
!!
!! Finally, the module diffuco_comb computes the combined $\alpha$ and $\beta$ coefficients for use 
!! elsewhere in the module. \n

!! RECENT CHANGE(S): Nathalie le 28 mars 2006 - sur proposition de Fred Hourdin, ajout
!! d'un potentiometre pour regler la resistance de la vegetation (rveg is now in pft_parameters)
!! October 2018: Removed diffuco_trans using Jarvis formula for calculation of beta coefficient
!!
!! REFERENCE(S)	: None
!!
!! SVN          :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/diffuco.f90 $
!! $Date: 2022-07-20 11:30:43 +0200 (mer. 20 juil. 2022) $
!! $Revision: 7709 $
!! \n
!_ ================================================================================================================================

MODULE diffuco

  ! modules used :
  USE constantes
  USE constantes_soil
  USE qsat_moisture
  USE sechiba_io_p
  USE ioipsl
  USE pft_parameters
  USE grid
  USE time, ONLY : one_day, dt_sechiba
  USE ioipsl_para 
  USE xios_orchidee
  USE chemistry, ONLY : chemistry_initialize, chemistry_bvoc, chemistry_clear
  IMPLICIT NONE

  ! public routines :
  PRIVATE
  PUBLIC :: diffuco_main, diffuco_initialize, diffuco_finalize, diffuco_clear

  INTERFACE Arrhenius_modified
     MODULE PROCEDURE Arrhenius_modified_0d, Arrhenius_modified_1d
  END INTERFACE

  !
  ! variables used inside diffuco module : declaration and initialisation
  !
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: wind                      !! Wind module (m s^{-1})
!$OMP THREADPRIVATE(wind)

CONTAINS


!!  =============================================================================================================================
!! SUBROUTINE:    diffuco_initialize
!!
!>\BRIEF	  Allocate module variables, read from restart file or initialize with default values
!!
!! DESCRIPTION:	  Allocate module variables, read from restart file or initialize with default values.
!!                Call chemistry_initialize for initialization of variables needed for the calculations of BVOCs.
!!
!! RECENT CHANGE(S): None
!!
!! REFERENCE(S): None
!! 
!! FLOWCHART: None
!! \n
!_ ==============================================================================================================================
  SUBROUTINE diffuco_initialize (kjit,    kjpindex, index,                  &
                                 rest_id, lalo,     neighbours, resolution, &
                                 rstruct, q_cdrag)
    
    !! 0. Variable and parameter declaration
    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number (-) 
    INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size (-)
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)   :: index            !! Indeces of the points on the map (-)
    INTEGER(i_std),INTENT (in)                         :: rest_id          !! _Restart_ file identifier (-)
    REAL(r_std),DIMENSION (kjpindex,2),   INTENT (in)  :: lalo             !! Geographical coordinates
    INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours  !! Vector of neighbours for each 
    REAL(r_std),DIMENSION (kjpindex,2), INTENT(in)     :: resolution       !! The size in km of each grid-box in X and Y
    
    !! 0.2 Output variables
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct          !! Structural resistance for the vegetation
    
    !! 0.3 Modified variables
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: q_cdrag          !! Surface drag coefficient  (-)
    
    !! 0.4 Local variables
    INTEGER                                            :: ilai
    INTEGER                                            :: jv
    INTEGER                                            :: ier
    CHARACTER(LEN=4)                                   :: laistring
    CHARACTER(LEN=80)                                  :: var_name        
    !_ ================================================================================================================================
    
    !! 1. Define flag ldq_cdrag_from_gcm. This flag determines if the cdrag should be taken from the GCM or be calculated. 
    !!    The default value is true if the q_cdrag variables was already initialized. This is the case when coupled to the LMDZ.

    !Config Key   = CDRAG_FROM_GCM
    !Config Desc  = Keep cdrag coefficient from gcm.
    !Config If    = OK_SECHIBA
    !Config Def   = y
    !Config Help  = Set to .TRUE. if you want q_cdrag coming from GCM (if q_cdrag on initialization is non zero).
    !Config         Keep cdrag coefficient from gcm for latent and sensible heat fluxes.
    !Config Units = [FLAG]
    IF ( ABS(MAXVAL(q_cdrag)) .LE. EPSILON(q_cdrag)) THEN
       ldq_cdrag_from_gcm = .FALSE.
    ELSE
       ldq_cdrag_from_gcm = .TRUE.
    ENDIF
    CALL getin_p('CDRAG_from_GCM', ldq_cdrag_from_gcm)
    IF (printlev>=2) WRITE(numout,*) "ldq_cdrag_from_gcm = ",ldq_cdrag_from_gcm

    !! 2. Allocate module variables
    ALLOCATE (wind(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'diffuco_initialize','Problem in allocate of variable wind','','')

    !! 3. Read variables from restart file
    IF (printlev>=3) WRITE (numout,*) 'Read DIFFUCO variables from restart file'

    CALL ioconf_setatt_p('UNITS', 's/m')
    CALL ioconf_setatt_p('LONG_NAME','Structural resistance')
    CALL restget_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, .TRUE., rstruct, "gather", nbp_glo, index_g)
    IF ( ALL(rstruct(:,:) == val_exp) ) THEN
       DO jv = 1, nvm
          rstruct(:,jv) = rstruct_const(jv)
       ENDDO
    ENDIF
    
    !! 4. Initialize chemistry module
    IF (printlev>=3) WRITE(numout,*) "ok_bvoc:",ok_bvoc
    IF ( ok_bvoc ) CALL chemistry_initialize(kjpindex, lalo, neighbours, resolution)
    
  END SUBROUTINE diffuco_initialize



!! ================================================================================================================================
!! SUBROUTINE    : diffuco_main
!!
!>\BRIEF	 The root subroutine for the module, which calls all other required
!! subroutines.
!! 
!! DESCRIPTION   : 

!! This is the main subroutine for the module. 
!! First it calculates the surface drag coefficient (via a call to diffuco_aero), using available parameters to determine
!! the stability of air in the surface layer by calculating the Richardson Nubmber. If a value for the 
!! surface drag coefficient is passed down from the atmospheric model and and if the flag 'ldq_cdrag_from_gcm' 
!! is set to TRUE, then the subroutine diffuco_aerod is called instead. This calculates the aerodynamic coefficient. \n
!!
!! Following this, an estimation of the saturated humidity at the surface is made (via a call
!! to qsatcalc in the module qsat_moisture). Following this the beta coefficients for sublimation (via 
!! diffuco_snow), interception (diffuco_inter), bare soil (diffuco_bare), and transpiration (via 
!! diffuco_trans_co2) are calculated in sequence. Finally 
!! the alpha and beta coefficients are combined (diffuco_comb). \n
!!
!! The surface drag coefficient is calculated for use within the module enerbil. It is required to to
!! calculate the aerodynamic coefficient for every flux. \n
!!
!! The various beta coefficients are used within the module enerbil for modifying the rate of evaporation, 
!! as appropriate for the surface. As explained in Chapter 2 of Guimberteau (2010), that module (enerbil) 
!! calculates the rate of evaporation essentially according to the expression $E = /beta E_{pot}$, where
!! E is the total evaporation and $E_{pot}$ the evaporation potential. If $\beta = 1$, there would be
!! essentially no resistance to evaporation, whereas should $\beta = 0$, there would be no evaporation and
!! the surface layer would be subject to some very stong hydrological stress. \n
!!
!! The following processes are calculated:
!! - call diffuco_aero for aerodynamic transfer coeficient
!! - call diffuco_snow for partial beta coefficient: sublimation
!! - call diffuco_inter for partial beta coefficient: interception for each type of vegetation
!! - call diffuco_bare for partial beta coefficient: bare soil
!! - call diffuco_trans_co2 for partial beta coefficient: transpiration for each type of vegetation, using Farqhar's formula
!! - call diffuco_comb for alpha and beta coefficient
!! - call chemistry_bvoc for alpha and beta coefficients for biogenic emissions
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): humrel, q_cdrag, vbeta, vbeta1, vbeta4,
!! vbeta2, vbeta3, rveget, cimean   
!!
!! REFERENCE(S) :				        
!! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
!! Circulation Model. Journal of Climate, 6, pp.248-273.
!! - de Rosnay, P, 1999. Représentation des interactions sol-plante-atmosphère dans le modèle de circulation générale
!! du LMD, 1999. PhD Thesis, Université Paris 6, available (25/01/12): 
!! http://www.ecmwf.int/staff/patricia_de_rosnay/publications.html#8
!! - Ducharne, A, 1997. Le cycle de l'eau: modélisation de l'hydrologie continentale, étude de ses interactions avec 
!! le climat, PhD Thesis, Université Paris 6
!! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation
!! sur le cycle de l'eau, PhD Thesis, available (25/01/12):
!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
!! - Lathière, J, 2005. Evolution des émissions de composés organiques et azotés par la biosphère continentale dans le 
!! modèle LMDz-INCA-ORCHIDEE, Université Paris 6
!!
!! FLOWCHART	:
!! \latexonly 
!!     \includegraphics[scale=0.5]{diffuco_main_flowchart.png}
!! \endlatexonly
!! \n
!_ ================================================================================================================================

  SUBROUTINE diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
     & zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, q_cdrag, qsurf, qair, pb, &
     & evap_bare_lim, evap_bare_lim_ns, evapot, evapot_corr, snow, flood_frac, flood_res, frac_nobio, snow_nobio, totfrac_nobio, &
     & swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
     & vbeta , vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, &
     & lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
     & hist_id, hist2_id)

  !! 0. Variable and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number (-) 
    INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size (-)
    INTEGER(i_std),INTENT (in)                         :: hist_id          !! _History_ file identifier (-)
    INTEGER(i_std),INTENT (in)                         :: hist2_id         !! _History_ file 2 identifier (-)
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)     :: index          !! Indeces of the points on the map (-)
    INTEGER(i_std),DIMENSION (kjpindex*(nlai+1)), INTENT (in) :: indexlai  !! Indeces of the points on the 3D map
    INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg       !! Indeces of the points on the 3D map (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: u                !! Eastward Lowest level wind speed (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: v                !! Northward Lowest level wind speed (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: zlev             !! Height of first layer (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: z0m              !! Surface roughness Length for momentum (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: z0h              !! Surface roughness Length for heat (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: roughheight      !! Effective height for roughness (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: temp_sol         !! Skin temperature (K)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: temp_air         !! Lowest level temperature (K)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: temp_growth      !! Growth temperature (°C) - Is equal to t2m_month
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: rau              !! Air Density (kg m^{-3})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: qsurf            !! Near surface air specific humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: qair             !! Lowest level air specific humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: snow             !! Snow mass (kg)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: flood_frac       !! Fraction of floodplains
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: flood_res        !! Reservoir in floodplains (estimation to avoid over-evaporation)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: pb               !! Surface level pressure (hPa)
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)      :: evap_bare_lim    !! Limit to the bare soil evaporation when the 
                                                                           !! 11-layer hydrology is used (-)
    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! Limit to the bare soil evaporation when the 
                                                                           !! 11-layer hydrology is used (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot           !! Soil Potential Evaporation (mm day^{-1}) 
                                                                           !! NdN This variable does not seem to be used at 
                                                                           !! all in diffuco
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot_corr      !! Soil Potential Evaporation
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio     !! Fraction of ice,lakes,cities,... (-)
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio     !! Snow on ice,lakes,cities,... (kg m^{-2})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: totfrac_nobio    !! Total fraction of ice+lakes+cities+... (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: swnet            !! Net surface short-wave flux (W m^{-2})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: swdown           !! Down-welling surface short-wave flux (W m^{-2})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: coszang          !! Cosine of the solar zenith angle (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ccanopy          !! CO2 concentration inside the canopy (ppm)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: veget            !! Fraction of vegetation type (-)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: veget_max        !! Max. fraction of vegetation type (LAI->infty)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: lai              !! Leaf area index (m^2 m^{-2})
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: qsintveg         !! Water on vegetation due to interception (kg m^{-2})
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: qsintmax         !! Maximum water on vegetation for interception 
                                                                           !! (kg m^{-2})
    REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! min+max+opt temps, vcmax, vjmax
                                                                           !! for photosynthesis (K ??)
    REAL(r_std),DIMENSION (kjpindex,2),   INTENT (in)  :: lalo               !! Geographical coordinates
    INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours    !! Vector of neighbours for each 
                                                                             !! grid point (1=N, 2=E, 3=S, 4=W)
    REAL(r_std),DIMENSION (kjpindex,2), INTENT(in)     :: resolution         !! The size in km of each grid-box in X and Y
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ptnlev1            !! 1st level of soil temperature (Kelvin)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: precip_rain        !! Rain precipitation expressed in mm/tstep
    REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT (in)  :: frac_age !! Age efficiency for isoprene emissions (from STOMATE)
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: tot_bare_soil      !! Total evaporating bare soil fraction 
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: frac_snow_veg      !! Snow cover fraction on vegeted area
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in):: frac_snow_nobio    !! Snow cover fraction on non-vegeted area

    !! 0.2 Output variables

    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vbeta            !! Total beta coefficient (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vbeta1           !! Beta for sublimation (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vbeta4           !! Beta for bare soil evaporation (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vbeta5           !! Beta for floodplains
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean           !! Mean stomatal conductance to CO2 (mol m-2 s-1) 
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2           !! Beta for interception loss (-)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3           !! Beta for transpiration (-)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot        !! Beta for potential transpiration
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget           !! Stomatal resistance for the whole canopy (s m^{-1})
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct          !! Structural resistance for the vegetation
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean           !! Mean leaf Ci (ppm)  
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out)  :: gpp              !! Assimilation ((gC m^{-2} dt_sechiba^{-1}), total area)  

    !! 0.3 Modified variables
 
    REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: humrel        !! Soil moisture stress (within range 0 to 1)
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)      :: q_cdrag       !! Surface drag coefficient  (-)

    !! 0.4 Local variables

    REAL(r_std),DIMENSION (kjpindex,nvm)     :: vbeta23            !! Beta for fraction of wetted foliage that will
                                                                   !! transpire once intercepted water has evaporated (-)
    REAL(r_std),DIMENSION (kjpindex)         :: raero              !! Aerodynamic resistance (s m^{-1})
    INTEGER(i_std)                           :: ilai
    CHARACTER(LEN=4)                         :: laistring
    CHARACTER(LEN=80)                        :: var_name           !! To store variables names for I/O
    REAL(r_std),DIMENSION(kjpindex)          :: qsatt              !! Surface saturated humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: cim                !! Intercellular CO2 over nlai 

    !! Based on ZunYin for correct transport
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: humrel_tmp         !! Soil moisture stress (within range 0 to 1)
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: vbeta3_tmp         !! Beta for transpiration (-)
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: rveget_tmp         !! Stomatal resistance for the whole canopy (s m^{-1})
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: rstruct_tmp        !! Structural resistance for the vegetation
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: cimean_tmp         !! Mean leaf Ci (ppm)
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: gsmean_tmp         !! Mean stomatal conductance to CO2 (umol m-2 s-1)
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: gpp_tmp            !! Assimilation ((gC m^{-2} s^{-1}), total area)
    REAL(r_std),DIMENSION (kjpindex,nvm)     :: cim_tmp            !! Intercellular CO2 over nlai

!_ ================================================================================================================================
    
    wind(:) = SQRT (u(:)*u(:) + v(:)*v(:))
  
  !! 1. Calculate the different coefficients

    IF (.NOT.ldq_cdrag_from_gcm) THEN
        ! Case 1a)
       CALL diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, &
                          qsurf, qair, snow, q_cdrag)
    ENDIF

    ! Case 1b)
    CALL diffuco_raerod (kjpindex, u, v, q_cdrag, raero)

  !! 2. Make an estimation of the saturated humidity at the surface

    CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)

  !! 3. Calculate the beta coefficient for sublimation
  
    CALL diffuco_snow (kjpindex, qair, qsatt, rau, u, v, q_cdrag, &
         snow, frac_nobio, totfrac_nobio, snow_nobio, frac_snow_veg, frac_snow_nobio, &
         vbeta1)


    CALL diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, &
         & flood_frac, flood_res, vbeta5)

  !! 4. Calculate the beta coefficient for interception

    CALL diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, &
       & qsintveg, qsintmax, rstruct, vbeta2, vbeta23) 


  !! 5. Calculate the beta coefficient for transpiration
    !! Dummy variables to compute vbeta3pot (humrel_temp = 1, no water stress)
    humrel_tmp(:,:) = un
    vbeta3_tmp = vbeta3
    rveget_tmp = rveget
    rstruct_tmp = rstruct
    cimean_tmp = cimean
    gsmean_tmp = gsmean
    gpp_tmp = gpp
    cim_tmp = cim

    CALL diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel, &
         assim_param, ccanopy, &
         veget, veget_max, lai, qsintveg, qsintmax, vbeta3, vbeta3pot, &
         rveget, rstruct, cimean, gsmean, gpp, &
         vbeta23, hist_id, indexveg, indexlai, index, kjit, cim)

   !! New resistance calc. with temporary (dummy) variables, just for vbeta3pot
    CALL diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel_tmp, &
         assim_param, ccanopy, &
         veget, veget_max, lai, qsintveg, qsintmax, vbeta3_tmp, vbeta3pot, &
         rveget_tmp, rstruct_tmp, cimean_tmp, gsmean_tmp, gpp_tmp, &
         vbeta23, hist_id, indexveg, indexlai, index, kjit, cim_tmp)

    !
    !biogenic emissions
    !
    IF ( ok_bvoc ) THEN
       CALL chemistry_bvoc (kjpindex, swdown, coszang, temp_air, &
            temp_sol, ptnlev1, precip_rain, humrel, veget_max, &
            lai, frac_age, lalo, ccanopy, cim, wind, snow, &
            veget, hist_id, hist2_id, kjit, index, &
            indexlai, indexveg)
    ENDIF
    !
    ! combination of coefficient : alpha and beta coefficient
    ! beta coefficient for bare soil
    !

    CALL diffuco_bare (kjpindex, tot_bare_soil, veget_max, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4)

  !! 6. Combine the alpha and beta coefficients

    ! Ajout qsintmax dans les arguments de la routine.... Nathalie / le 13-03-2006
    CALL diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, snow, &
       & veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3, vbeta4, &
       & evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax)

    CALL xios_orchidee_send_field("q_cdrag",q_cdrag)
    CALL xios_orchidee_send_field("raero",raero)
    CALL xios_orchidee_send_field("wind",wind)
    CALL xios_orchidee_send_field("qsatt",qsatt)
    CALL xios_orchidee_send_field("coszang",coszang)
    CALL xios_orchidee_send_field('cim', cim)

    IF ( .NOT. almaoutput ) THEN
       CALL histwrite_p(hist_id, 'raero', kjit, raero, kjpindex, index)
       CALL histwrite_p(hist_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
       CALL histwrite_p(hist_id, 'Wind', kjit, wind, kjpindex, index)
       CALL histwrite_p(hist_id, 'qsatt', kjit, qsatt, kjpindex, index)
       CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg)

       IF ( hist2_id > 0 ) THEN
          CALL histwrite_p(hist2_id, 'raero', kjit, raero, kjpindex, index)
          CALL histwrite_p(hist2_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
          CALL histwrite_p(hist2_id, 'Wind', kjit, wind, kjpindex, index)
          CALL histwrite_p(hist2_id, 'qsatt', kjit, qsatt, kjpindex, index)
       ENDIF
    ELSE
       CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg)
    ENDIF

    IF (printlev>=3) WRITE (numout,*) ' diffuco_main done '

  END SUBROUTINE diffuco_main

!!  =============================================================================================================================
!! SUBROUTINE: diffuco_finalize
!!
!>\BRIEF          Write to restart file
!!
!! DESCRIPTION:   This subroutine writes the module variables and variables calculated in diffuco
!!                to restart file
!!
!! RECENT CHANGE(S): None
!! REFERENCE(S): None
!! FLOWCHART: None
!! \n
!_ ==============================================================================================================================
  SUBROUTINE diffuco_finalize (kjit, kjpindex, rest_id, rstruct )

    !! 0. Variable and parameter declaration
    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number (-) 
    INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size (-)
    INTEGER(i_std),INTENT (in)                         :: rest_id          !! _Restart_ file identifier (-)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: rstruct          !! Structural resistance for the vegetation

    !! 0.4 Local variables
    INTEGER                                            :: ilai
    CHARACTER(LEN=4)                                   :: laistring
    CHARACTER(LEN=80)                                  :: var_name        

!_ ================================================================================================================================
    
  !! 1. Prepare the restart file for the next simulation
    IF (printlev>=3) WRITE (numout,*) 'Complete restart file with DIFFUCO variables '

    CALL restput_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, rstruct, 'scatter',  nbp_glo, index_g)
    
  END SUBROUTINE diffuco_finalize


!! ================================================================================================================================
!! SUBROUTINE		 			: diffuco_clear
!!
!>\BRIEF					Housekeeping module to deallocate the variables
!! rstruct and raero
!!
!! DESCRIPTION				        : Housekeeping module to deallocate the variables
!! rstruct and raero
!!
!! RECENT CHANGE(S)                             : None
!!
!! MAIN OUTPUT VARIABLE(S)	                : None
!!
!! REFERENCE(S)				        : None
!!
!! FLOWCHART                                    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE diffuco_clear()

    ! Deallocate and reset variables in chemistry module
    CALL chemistry_clear

  END SUBROUTINE diffuco_clear


!! ================================================================================================================================
!! SUBROUTINE	: diffuco_aero
!!
!>\BRIEF	This module first calculates the surface drag 
!! coefficient, for cases in which the surface drag coefficient is NOT provided by the coupled 
!! atmospheric model LMDZ or when the flag ldq_cdrag_from_gcm is set to FALSE 
!!
!! DESCRIPTION	: Computes the surface drag coefficient, for cases 
!! in which it is NOT provided by the coupled atmospheric model LMDZ. The module first uses the 
!! meteorolgical input to calculate the Richardson Number, which is an indicator of atmospheric 
!! stability in the surface layer. The formulation used to find this surface drag coefficient is 
!! dependent on the stability determined. \n
!!
!! Designation of wind speed
!! \latexonly 
!!     \input{diffucoaero1.tex}
!! \endlatexonly
!!
!! Calculation of geopotential. This is the definition of Geopotential height (e.g. Jacobson 
!! eqn.4.47, 2005). (required for calculation of the Richardson Number)
!! \latexonly 
!!     \input{diffucoaero2.tex}
!! \endlatexonly
!! 
!! \latexonly 
!!     \input{diffucoaero3.tex}
!! \endlatexonly
!!
!! Calculation of the virtual air temperature at the surface (required for calculation
!! of the Richardson Number)
!! \latexonly 
!!     \input{diffucoaero4.tex}
!! \endlatexonly
!!
!! Calculation of the virtual surface temperature (required for calculation of th
!! Richardson Number)
!! \latexonly 
!!     \input{diffucoaero5.tex}
!! \endlatexonly
!!
!! Calculation of the squared wind shear (required for calculation of the Richardson
!! Number)
!! \latexonly 
!!     \input{diffucoaero6.tex}
!! \endlatexonly
!! 
!! Calculation of the Richardson Number. The Richardson Number is defined as the ratio 
!! of potential to kinetic energy, or, in the context of atmospheric science, of the
!! generation of energy by wind shear against consumption
!! by static stability and is an indicator of flow stability (i.e. for when laminar flow 
!! becomes turbulent and vise versa). It is approximated using the expression below:
!! \latexonly 
!!     \input{diffucoaero7.tex}
!! \endlatexonly
!!
!! The Richardson Number hence calculated is subject to a minimum value:
!! \latexonly 
!!     \input{diffucoaero8.tex}
!! \endlatexonly
!! 
!! Computing the drag coefficient. We add the add the height of the vegetation to the 
!! level height to take into account that the level 'seen' by the vegetation is actually 
!! the top of the vegetation. Then we we can subtract the displacement height.
!! \latexonly 
!!     \input{diffucoaero9.tex}
!! \endlatexonly
!! 
!! For the stable case (i.e $R_i$ $\geq$ 0)
!! \latexonly 
!!     \input{diffucoaero10.tex}
!! \endlatexonly
!!
!! \latexonly 
!!     \input{diffucoaero11.tex}
!! \endlatexonly
!!          
!! For the unstable case (i.e. $R_i$ < 0)
!! \latexonly 
!!     \input{diffucoaero12.tex}
!! \endlatexonly
!!
!! \latexonly 
!!     \input{diffucoaero13.tex}
!! \endlatexonly
!!               
!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
!! To prevent this, a minimum limit to the drag coefficient is defined as:
!!
!! \latexonly 
!!     \input{diffucoaero14.tex}
!! \endlatexonly
!! 
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): q_cdrag
!!
!! REFERENCE(S)	: 
!! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
!! Circulation Model. Journal of Climate, 6, pp.248-273
!! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation
!! sur le cycle de l'eau, PhD Thesis, available from:
!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
!! - Jacobson M.Z., Fundamentals of Atmospheric Modeling (2nd Edition), published Cambridge 
!! University Press, ISBN 0-521-54865-9
!!
!! FLOWCHART	:
!! \latexonly 
!!     \includegraphics[scale=0.5]{diffuco_aero_flowchart.png}
!! \endlatexonly
!! \n
!_ ================================================================================================================================

  SUBROUTINE diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, &
                           qsurf, qair, snow, q_cdrag)

  !! 0. Variable and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                          :: kjpindex, kjit   !! Domain size
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: u                !! Eastward Lowest level wind speed (m s^{-1}) 
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: v                !! Northward Lowest level wind speed (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: zlev             !! Height of first atmospheric layer (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: z0h               !! Surface roughness Length for heat (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: z0m               !! Surface roughness Length for momentum (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: roughheight      !! Effective roughness height (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: temp_sol         !! Ground temperature (K)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: temp_air         !! Lowest level temperature (K)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: qsurf            !! near surface specific air humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: qair             !! Lowest level specific air humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: snow             !! Snow mass (kg)

    !! 0.2 Output variables
    
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: q_cdrag          !! Surface drag coefficient  (-)

    !! 0.3 Modified variables

    !! 0.4 Local variables

    INTEGER(i_std)                                      :: ji, jv
    REAL(r_std)                                         :: speed, zg, zdphi, ztvd, ztvs, zdu2
    REAL(r_std)                                         :: zri, cd_neut, zscf, cd_tmp
!_ ================================================================================================================================

  !! 1. Initialisation

    ! test if we have to work with q_cdrag or to calcul it
    DO ji=1,kjpindex
       
       !! 1a).1 Designation of wind speed
       !! \latexonly 
       !!     \input{diffucoaero1.tex}
       !! \endlatexonly
       speed = wind(ji)
    
       !! 1a).2 Calculation of geopotentiel
       !! This is the definition of Geopotential height (e.g. Jacobson eqn.4.47, 2005). (required
       !! for calculation of the Richardson Number)
       !! \latexonly 
       !!     \input{diffucoaero2.tex}
       !! \endlatexonly
       zg = zlev(ji) * cte_grav
      
       !! \latexonly 
       !!     \input{diffucoaero3.tex}
       !! \endlatexonly
       zdphi = zg/cp_air
       
       !! 1a).3 Calculation of the virtual air temperature at the surface 
       !! required for calculation of the Richardson Number
       !! \latexonly 
       !!     \input{diffucoaero4.tex}
       !! \endlatexonly
       ztvd = (temp_air(ji) + zdphi / (un + rvtmp2 * qair(ji))) * (un + retv * qair(ji)) 
       
       !! 1a).4 Calculation of the virtual surface temperature 
       !! required for calculation of the Richardson Number
       !! \latexonly 
       !!     \input{diffucoaero5.tex}
       !! \endlatexonly
       ztvs = temp_sol(ji) * (un + retv * qsurf(ji))
     
       !! 1a).5 Calculation of the squared wind shear 
       !! required for calculation of the Richardson Number
       !! \latexonly 
       !!     \input{diffucoaero6.tex}
       !! \endlatexonly
       zdu2 = MAX(cepdu2,speed**2)
       
       !! 1a).6 Calculation of the Richardson Number
       !!  The Richardson Number is defined as the ratio of potential to kinetic energy, or, in the 
       !!  context of atmospheric science, of the generation of energy by wind shear against consumption
       !!  by static stability and is an indicator of flow stability (i.e. for when laminar flow 
       !!  becomes turbulent and vise versa).\n
       !!  It is approximated using the expression below:
       !!  \latexonly 
       !!     \input{diffucoaero7.tex}
       !! \endlatexonly
       zri = zg * (ztvd - ztvs) / (zdu2 * ztvd)
      
       !! The Richardson Number hence calculated is subject to a minimum value:
       !! \latexonly 
       !!     \input{diffucoaero8.tex}
       !! \endlatexonly       
       zri = MAX(MIN(zri,5.),-5.)
       
       !! 1a).7 Computing the drag coefficient
       !!  We add the add the height of the vegetation to the level height to take into account
       !!  that the level 'seen' by the vegetation is actually the top of the vegetation. Then we 
       !!  we can subtract the displacement height.
       !! \latexonly 
       !!     \input{diffucoaero9.tex}
       !! \endlatexonly

       !! 7.0 Snow smoothering
       !! Snow induces low levels of turbulence.
       !! Sensible heat fluxes can therefore be reduced of ~1/3. Pomeroy et al., 1998
       cd_neut = ct_karman ** 2. / ( LOG( (zlev(ji) + roughheight(ji)) / z0m(ji) ) * LOG( (zlev(ji) + roughheight(ji)) / z0h(ji) ) )
       
       !! 1a).7.1 - for the stable case (i.e $R_i$ $\geq$ 0)
       IF (zri .GE. zero) THEN
          
          !! \latexonly 
          !!     \input{diffucoaero10.tex}
          !! \endlatexonly
          zscf = SQRT(un + cd * ABS(zri))
         
          !! \latexonly 
          !!     \input{diffucoaero11.tex}
          !! \endlatexonly          
          cd_tmp=cd_neut/(un + trois * cb * zri * zscf)
       ELSE
          
          !! 1a).7.2 - for the unstable case (i.e. $R_i$ < 0)
          !! \latexonly 
          !!     \input{diffucoaero12.tex}
          !! \endlatexonly
          zscf = un / (un + trois * cb * cc * cd_neut * SQRT(ABS(zri) * &
               & ((zlev(ji) + roughheight(ji)) / z0m(ji))))

          !! \latexonly 
          !!     \input{diffucoaero13.tex}
          !! \endlatexonly               
          cd_tmp=cd_neut * (un - trois * cb * zri * zscf)
       ENDIF
       
       !! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
       !! To prevent this, a minimum limit to the drag coefficient is defined as:
       
       !! \latexonly 
       !!     \input{diffucoaero14.tex}
       !! \endlatexonly
       !!
       q_cdrag(ji) = MAX(cd_tmp, min_qc/MAX(speed,min_wind))

       ! In some situations it might be useful to give an upper limit on the cdrag as well. 
       ! The line here should then be uncommented.
      !q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind))

    END DO

    IF (printlev>=3) WRITE (numout,*) ' not ldqcdrag_from_gcm : diffuco_aero done '

  END SUBROUTINE diffuco_aero


!! ================================================================================================================================
!! SUBROUTINE    : diffuco_snow
!!
!>\BRIEF         This subroutine computes the beta coefficient for snow sublimation.
!!
!! DESCRIPTION   : This routine computes beta coefficient for snow sublimation, which
!! integrates the snow on both vegetation and other surface types (e.g. ice, lakes,
!! cities etc.) \n
!!
!! A critical depth of snow (snowcri) is defined to calculate the fraction of each grid-cell
!! that is covered with snow (snow/snowcri) while the remaining part is snow-free.
!! We also carry out a first calculation of sublimation (subtest) to lower down the beta
!! coefficient if necessary (if subtest > snow). This is a predictor-corrector test. 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::vbeta1 
!!
!! REFERENCE(S) :
!! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
!! Circulation Model. Journal of Climate, 6, pp. 248-273
!! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation
!! sur le cycle de l'eau, PhD Thesis, available from:
!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================
  
SUBROUTINE diffuco_snow (kjpindex, qair, qsatt, rau, u, v,q_cdrag, &
       & snow, frac_nobio, totfrac_nobio, snow_nobio, frac_snow_veg, frac_snow_nobio, &
       vbeta1)

  !! 0. Variable and parameter declaration
    
    !! 0.1 Input variables
 
    INTEGER(i_std), INTENT(in)                           :: kjpindex       !! Domain size (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: qair           !! Lowest level specific air humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: qsatt          !! Surface saturated humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: rau            !! Air density (kg m^{-3})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: u              !! Eastward Lowest level wind speed (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: v              !! Northward Lowest level wind speed (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: q_cdrag        !! Surface drag coefficient  (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: snow           !! Snow mass (kg m^{-2})
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio     !! Fraction of ice, lakes, cities etc. (-)
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio     !! Snow on ice, lakes, cities etc. (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: totfrac_nobio  !! Total fraction of ice, lakes, cities etc. (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: frac_snow_veg  !! Snow cover fraction on vegeted area
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in)  :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
    
    !! 0.2 Output variables

    REAL(r_std),DIMENSION (kjpindex), INTENT (out)       :: vbeta1         !! Beta for sublimation (dimensionless ratio) 
    
    !! 0.3 Modified variables

    !! 0.4 Local variables

    REAL(r_std)                                          :: subtest        !! Sublimation for test (kg m^{-2})
    REAL(r_std)                                          :: zrapp          !! Modified factor (ratio)
    REAL(r_std)                                          :: speed          !! Wind speed (m s^{-1})
    REAL(r_std)                                          :: vbeta1_add     !! Beta for sublimation (ratio)
    INTEGER(i_std)                                       :: ji, jv         !! Indices (-)
!_ ================================================================================================================================

  !! 1. Calculate beta coefficient for snow sublimation on the vegetation\n

    DO ji=1,kjpindex  ! Loop over # pixels - domain size

       ! Fraction of mesh that can sublimate snow
       vbeta1(ji) = (un - totfrac_nobio(ji)) * frac_snow_veg(ji)

       ! Limitation of sublimation in case of snow amounts smaller than the atmospheric demand. 
       speed = MAX(min_wind, wind(ji))

       subtest = dt_sechiba * vbeta1(ji) * speed * q_cdrag(ji) * rau(ji) * &
               & ( qsatt(ji) - qair(ji) )

       IF ( subtest .GT. min_sechiba ) THEN
          zrapp = snow(ji) / subtest
          IF ( zrapp .LT. un ) THEN
             vbeta1(ji) = vbeta1(ji) * zrapp
          ENDIF
       ENDIF

    END DO ! Loop over # pixels - domain size

  !! 2. Add the beta coefficients calculated from other surfaces types (snow on ice,lakes, cities...)

    DO jv = 1, nnobio ! Loop over # other surface types
!!$      !
!!$      IF ( jv .EQ. iice ) THEN
!!$        !
!!$        !  Land ice is of course a particular case
!!$        !
!!$        DO ji=1,kjpindex
!!$          vbeta1(ji) = vbeta1(ji) + frac_nobio(ji,jv)
!!$        ENDDO
!!$        !
!!$      ELSE
        !
        DO ji=1,kjpindex ! Loop over # pixels - domain size

           vbeta1_add = frac_nobio(ji,jv) * frac_snow_nobio(ji, jv)

           ! Limitation of sublimation in case of snow amounts smaller than
           ! the atmospheric demand. 
           speed = MAX(min_wind, wind(ji))
            
            !!     Limitation of sublimation by the snow accumulated on the ground 
            !!     A first approximation is obtained with the old values of
            !!     qair and qsol_sat: function of temp-sol and pb. (see call of qsatcalc)
           subtest = dt_sechiba * vbeta1_add * speed * q_cdrag(ji) * rau(ji) * &
                & ( qsatt(ji) - qair(ji) )

           IF ( subtest .GT. min_sechiba ) THEN
              zrapp = snow_nobio(ji,jv) / subtest
              IF ( zrapp .LT. un ) THEN
                 vbeta1_add = vbeta1_add * zrapp
              ENDIF
           ENDIF

           vbeta1(ji) = vbeta1(ji) + vbeta1_add

        ENDDO ! Loop over # pixels - domain size

!!$      ENDIF
      
    ENDDO ! Loop over # other surface types

    IF (printlev>=3) WRITE (numout,*) ' diffuco_snow done '

  END SUBROUTINE diffuco_snow


!! ================================================================================================================================
!! SUBROUTINE		 			: diffuco_flood 
!!
!>\BRIEF				       	This routine computes partial beta coefficient : floodplains
!!
!! DESCRIPTION				        : 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S)	                : vbeta5
!!
!! REFERENCE(S)				        : None
!!
!! FLOWCHART                                    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, &
       & flood_frac, flood_res, vbeta5)

    ! interface description
    ! input scalar 
    INTEGER(i_std), INTENT(in)                               :: kjpindex   !! Domain size
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair       !! Lowest level specific humidity
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qsatt      !! Surface saturated humidity
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: rau        !! Density
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u          !! Lowest level wind speed 
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: v          !! Lowest level wind speed
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q_cdrag    !! Surface drag coefficient  (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: flood_res  !! water mass in flood reservoir
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: flood_frac !! fraction of floodplains
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: evapot     !! Potential evaporation
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: evapot_corr!! Potential evaporation2
    ! output fields
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vbeta5     !! Beta for floodplains

    ! local declaration
    REAL(r_std)                                              :: subtest, zrapp, speed
    INTEGER(i_std)                                           :: ji, jv

!_ ================================================================================================================================
    !
    ! beta coefficient for sublimation for floodplains
    !
    DO ji=1,kjpindex
       !
       IF (evapot(ji) .GT. min_sechiba) THEN
          vbeta5(ji) = flood_frac(ji) *evapot_corr(ji)/evapot(ji)
       ELSE
          vbeta5(ji) = flood_frac(ji)
       ENDIF
       !
       ! -- Limitation of evaporation in case of water amounts smaller than
       !    the atmospheric demand. 
       
       !
       speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))
       !
       subtest = dt_sechiba * vbeta5(ji) * speed * q_cdrag(ji) * rau(ji) * &
               & ( qsatt(ji) - qair(ji) )
       !  
       IF ( subtest .GT. min_sechiba ) THEN
          zrapp = flood_res(ji) / subtest
          IF ( zrapp .LT. un ) THEN
             vbeta5(ji) = vbeta5(ji) * zrapp
          ENDIF
       ENDIF
       !
    END DO

    IF (printlev>=3) WRITE (numout,*) ' diffuco_flood done '

  END SUBROUTINE diffuco_flood


!! ================================================================================================================================
!! SUBROUTINE    : diffuco_inter
!!
!>\BRIEF	 This routine computes the partial beta coefficient
!! for the interception for each type of vegetation
!!
!! DESCRIPTION   : We first calculate the dry and wet parts of each PFT (wet part = qsintveg/qsintmax).
!! The former is submitted to transpiration only (vbeta3 coefficient, calculated in 
!! diffuco_trans_co2), while the latter is first submitted to interception loss 
!! (vbeta2 coefficient) and then to transpiration once all the intercepted water has been evaporated 
!! (vbeta23 coefficient). Interception loss is also submitted to a predictor-corrector test, 
!! as for snow sublimation. \n
!!
!! \latexonly 
!!     \input{diffucointer1.tex}
!! \endlatexonly
!! Calculate the wet fraction of vegetation as  the ration between the intercepted water and the maximum water 
!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
!!
!! \latexonly 
!!     \input{diffucointer2.tex}
!! \endlatexonly
!!
!! Calculation of $\beta_3$, the canopy transpiration resistance
!! \latexonly 
!!     \input{diffucointer3.tex}
!! \endlatexonly            
!! 
!! We here determine the limitation of interception loss by the water stored on the leaf. 
!! A first approximation of interception loss is obtained using the old values of
!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
!! \latexonly 
!!     \input{diffucointer4.tex}
!! \endlatexonly
!!
!! \latexonly
!!     \input{diffucointer5.tex}
!! \endlatexonly
!!
!! \latexonly 
!!     \input{diffucointer6.tex}
!! \endlatexonly
!!
!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
!! 'zqsvegrap'.
!! \latexonly 
!!     \input{diffucointer7.tex}
!! \endlatexonly
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::vbeta2, ::vbeta23
!!
!! REFERENCE(S) :
!! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
!! Circulation Model. Journal of Climate, 6, pp. 248-273
!! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation
!! sur le cycle de l'eau, PhD Thesis, available from:
!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
!! - Perrier, A, 1975. Etude physique de l'évaporation dans les conditions naturelles. Annales 
!! Agronomiques, 26(1-18): pp. 105-123, pp. 229-243
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, &
     & qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
   
  !! 0 Variable and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                           :: kjpindex   !! Domain size (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: qair       !! Lowest level specific air humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: qsatt      !! Surface saturated humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: rau        !! Air Density (kg m^{-3})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: u          !! Eastward Lowest level wind speed (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: v          !! Northward Lowest level wind speed (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: q_cdrag    !! Surface drag coefficient  (-)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: humrel     !! Soil moisture stress (within range 0 to 1)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: veget      !! vegetation fraction for each type (fraction)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: qsintveg   !! Water on vegetation due to interception (kg m^{-2})
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: qsintmax   !! Maximum water on vegetation (kg m^{-2})
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: rstruct    !! architectural resistance (s m^{-1})
    
    !! 0.2 Output variables
    
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)   :: vbeta2     !! Beta for interception loss (-)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)   :: vbeta23    !! Beta for fraction of wetted foliage that will 
                                                                       !! transpire (-)

    !! 0.4 Local variables

    INTEGER(i_std)                                       :: ji, jv                               !! (-), (-)
    REAL(r_std)                                          :: zqsvegrap, ziltest, zrapp, speed     !!
!_ ================================================================================================================================

  !! 1. Initialize

    vbeta2(:,:) = zero
    vbeta23(:,:) = zero
   
  !! 2. The beta coefficient for interception by vegetation. 
    
    DO jv = 2,nvm

      DO ji=1,kjpindex

         IF (veget(ji,jv) .GT. min_sechiba .AND. qsintveg(ji,jv) .GT. zero ) THEN

            zqsvegrap = zero
            IF (qsintmax(ji,jv) .GT. min_sechiba ) THEN

            !! \latexonly 
            !!     \input{diffucointer1.tex}
            !! \endlatexonly
            !!
            !! We calculate the wet fraction of vegetation as  the ration between the intercepted water and the maximum water 
            !! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
            !!
                zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
            END IF

            !! \latexonly 
            !!     \input{diffucointer2.tex}
            !! \endlatexonly
            speed = MAX(min_wind, wind(ji))

            !! Calculation of $\beta_3$, the canopy transpiration resistance
            !! \latexonly 
            !!     \input{diffucointer3.tex}
            !! \endlatexonly
            vbeta2(ji,jv) = veget(ji,jv) * zqsvegrap * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv)))
            
            !! We here determine the limitation of interception loss by the water stored on the leaf. 
            !! A first approximation of interception loss is obtained using the old values of
            !! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
            !! \latexonly 
            !!     \input{diffucointer4.tex}
            !! \endlatexonly
            ziltest = dt_sechiba * vbeta2(ji,jv) * speed * q_cdrag(ji) * rau(ji) * &
               & ( qsatt(ji) - qair(ji) )

            IF ( ziltest .GT. min_sechiba ) THEN

                !! \latexonly 
                !!     \input{diffucointer5.tex}
                !! \endlatexonly
                zrapp = qsintveg(ji,jv) / ziltest
                IF ( zrapp .LT. un ) THEN
                   
                    !! \latexonly 
                    !!     \input{diffucointer6.tex}
                    !! \endlatexonly
                    !!
		    !! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
                    !! 'zqsvegrap'.
                   IF ( humrel(ji,jv) >= min_sechiba ) THEN
                      vbeta23(ji,jv) = MAX(vbeta2(ji,jv) - vbeta2(ji,jv) * zrapp, zero)
                   ELSE
                      ! We don't want transpiration when the soil cannot deliver it
                      vbeta23(ji,jv) = zero
                   ENDIF
                    
                    !! \latexonly 
                    !!     \input{diffucointer7.tex}
                    !! \endlatexonly
                    vbeta2(ji,jv) = vbeta2(ji,jv) * zrapp
                ENDIF
            ENDIF
        END IF
!        ! Autre formulation possible pour l'evaporation permettant une transpiration sur tout le feuillage
!        !commenter si formulation Nathalie sinon Tristan
!        speed = MAX(min_wind, wind(ji))
!        
!        vbeta23(ji,jv) = MAX(zero, veget(ji,jv) * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv))) - vbeta2(ji,jv))

      END DO

    END DO

    IF (printlev>=3) WRITE (numout,*) ' diffuco_inter done '

  END SUBROUTINE diffuco_inter


!! ==============================================================================================================================
!! SUBROUTINE      : diffuco_bare
!!
!>\BRIEF	   This routine computes the partial beta coefficient corresponding to
!! bare soil
!!
!! DESCRIPTION	   : Bare soil evaporation is submitted to a maximum possible flow (evap_bare_lim)
!! 
!! Calculation of wind speed
!! \latexonly 
!!     \input{diffucobare1.tex}
!! \endlatexonly
!!             
!! The calculation of $\beta_4$
!! \latexonly 
!!     \input{diffucobare2.tex}
!! \endlatexonly
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::vbeta4
!!
!! REFERENCE(S)	 :
!! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
!! Circulation Model. Journal of Climate, 6, pp.248-273
!! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation
!! sur le cycle de l'eau, PhD Thesis, available from:
!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE diffuco_bare (kjpindex, tot_bare_soil, veget_max, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4)

    !! 0. Variable and parameter declaration

    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                         :: kjpindex       !! Domain size (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: tot_bare_soil  !! Total evaporating bare soil fraction
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: veget_max      !! Max. fraction of vegetation type (LAI->infty) 
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: vbeta2         !! Beta for Interception 
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: vbeta3         !! Beta for Transpiration 

    !! 0.2 Output variables
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vbeta4         !! Beta for bare soil evaporation (-)
    
    !! 0.3 Modified variables 
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: evap_bare_lim  !! limiting factor for bare soil evaporation
                                                                         !! when the 11-layer hydrology is used (-)
    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! limiting factor for bare soil evaporation 
                                                                         !! when the 11-layer hydrology is used (-)       
    !! 0.4 Local variables
    INTEGER(i_std)                                     :: ji
    REAL(r_std), DIMENSION(kjpindex)                   :: vegtot
    
!_ ================================================================================================================================

  !! 1. Calculation of the soil resistance and the beta (beta_4) for bare soil

    ! To use the new version of hydrol_split_soil and ensure the water conservation, we must have throughout sechiba at all time:
    !        a) evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)*soiltile(ji,:)*vegtot(ji))
    !        b) all the terms (vbeta4, evap_bare_lim, evap_bare_lim_ns) =0 if evap_bare_lim(ji) LE min_sechiba
    ! This must also be kept true in diffuco_comb
    
    DO ji = 1, kjpindex
       
       ! The limitation by 1-beta2-beta3 is due to the fact that evaporation under vegetation is possible
       !! \latexonly 
       !!     \input{diffucobare3.tex}
       !! \endlatexonly
       
       IF ( (evap_bare_lim(ji) .GT. min_sechiba) .AND. &  
            ! in this case we can't have vegtot LE min_sechina, cf hydrol_soil
            (un - SUM(vbeta2(ji,:)+vbeta3(ji,:)) .GT. min_sechiba) ) THEN 
          ! eventually, none of the left-hand term is close to zero

          vegtot(ji) = SUM(veget_max(ji,:)) 
          IF (evap_bare_lim(ji) < (un - SUM(vbeta2(ji,:)+vbeta3(ji,:)))) THEN 
             ! Standard case 
             vbeta4(ji) = evap_bare_lim(ji) 
          ELSE 
             vbeta4(ji) = un - SUM(vbeta2(ji,:)+vbeta3(ji,:)) 
             
             ! We now have to redefine evap_bare_lim & evap_bare_lim_ns 
             IF (evap_bare_lim(ji) .GT. min_sechiba) THEN 
                evap_bare_lim_ns(ji,:) = evap_bare_lim_ns(ji,:) * vbeta4(ji) / evap_bare_lim(ji) 
             ELSE ! we must re-invent evap_bare_lim_ns => uniform across soiltiles              
                evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji) 
             ENDIF
	 
             evap_bare_lim(ji) = vbeta4(ji) 
             ! consistent with evap_bare_lim(ji) = 
             ! SUM(evap_bare_lim_ns(ji,:)*soiltile(ji,:)*vegtot(ji)) 
             ! as SUM(soiltile(ji,:)) = 1 
          END IF
          
       ELSE ! instead of having very small vbeta4, we set everything to zero
          vbeta4(ji) = zero
          evap_bare_lim(ji) = zero
          evap_bare_lim_ns(ji,:) = zero
       ENDIF
       
    END DO
    
    IF (printlev>=3) WRITE (numout,*) ' diffuco_bare done '
    
  END SUBROUTINE diffuco_bare


!! ==============================================================================================================================
!! SUBROUTINE   : diffuco_trans_co2
!!
!>\BRIEF        This subroutine computes carbon assimilation and stomatal 
!! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987).
!!
!! DESCRIPTION  :\n
!! *** General:\n 
!! The equations are different depending on the photosynthesis mode (C3 versus C4).\n 
!! Assimilation and conductance are computed over 20 levels of LAI and then 
!! integrated at the canopy level.\n 
!! This routine also computes partial beta coefficient: transpiration for each 
!! type of vegetation.\n
!! There is a main loop on the PFTs, then inner loops on the points where 
!! assimilation has to be calculated.\n
!! This subroutine is called at each sechiba time step by sechiba_main.\n
!! *** Details:
!! - Integration at the canopy level\n
!! \latexonly
!! \input{diffuco_trans_co2_1.1.tex}
!! \endlatexonly
!! - Light''s extinction \n
!! The available light follows a simple Beer extinction law. 
!! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.\n
!! \latexonly
!! \input{diffuco_trans_co2_1.2.tex}
!! \endlatexonly
!! - Estimation of relative humidity of air (for calculation of the stomatal conductance)\n
!! \latexonly
!! \input{diffuco_trans_co2_1.3.tex}
!! \endlatexonly
!! - Calculation of the water limitation factor\n
!! \latexonly
!! \input{diffuco_trans_co2_2.1.tex}
!! \endlatexonly
!! - Calculation of temperature dependent parameters for C4 plants\n
!! \latexonly
!! \input{diffuco_trans_co2_2.2.tex}
!! \endlatexonly
!! - Calculation of temperature dependent parameters for C3 plants\n
!! \latexonly
!! \input{diffuco_trans_co2_2.3.tex}
!! \endlatexonly
!! - Vmax scaling\n 
!! Vmax is scaled into the canopy due to reduction of nitrogen 
!! (Johnson and Thornley,1984).\n
!! \latexonly
!! \input{diffuco_trans_co2_2.4.1.tex}
!! \endlatexonly
!! - Assimilation for C4 plants (Collatz et al., 1992)\n
!! \latexonly
!! \input{diffuco_trans_co2_2.4.2.tex}
!! \endlatexonly         
!! - Assimilation for C3 plants (Farqhuar et al., 1980)\n
!! \latexonly
!! \input{diffuco_trans_co2_2.4.3.tex}
!! \endlatexonly
!! - Estimation of the stomatal conductance (Ball et al., 1987)\n
!! \latexonly
!! \input{diffuco_trans_co2_2.4.4.tex}
!! \endlatexonly
!!
!! RECENT CHANGE(S):
!!
!!   2018/11 replaced t2m and q2m by temp_air and qair; prognostic variables on first atmospheric model layer
!!
!!   2006/06   N. de Noblet
!!                - addition of q2m and t2m as input parameters for the 
!!                calculation of Rveget
!!                - introduction of vbeta23
!!
!! MAIN OUTPUT VARIABLE(S): beta coefficients, resistances, CO2 intercellular 
!! concentration
!!
!! REFERENCE(S) :
!! - Ball, J., T. Woodrow, and J. Berry (1987), A model predicting stomatal 
!! conductance and its contribution to the control of photosynthesis under 
!! different environmental conditions, Prog. Photosynthesis, 4, 221– 224.
!! - Collatz, G., M. Ribas-Carbo, and J. Berry (1992), Coupled photosynthesis 
!! stomatal conductance model for leaves of C4 plants, Aust. J. Plant Physiol.,
!! 19, 519–538.
!! - Farquhar, G., S. von Caemmener, and J. Berry (1980), A biochemical model of 
!! photosynthesis CO2 fixation in leaves of C3 species, Planta, 149, 78–90.
!! - Johnson, I. R., and J. Thornley (1984), A model of instantaneous and daily
!! canopy photosynthesis, J Theor. Biol., 107, 531 545
!! - McMurtrie, R.E., Rook, D.A. and Kelliher, F.M., 1990. Modelling the yield of Pinus radiata on a
!! site limited by water and nitrogen. For. Ecol. Manage., 30: 381-413
!! - Bounoua, L., Hall, F. G., Sellers, P. J., Kumar, A., Collatz, G. J., Tucker, C. J., and Imhoff, M. L. (2010), Quantifying the 
!! negative feedback of vegetation to greenhouse warming: A modeling approach, Geophysical Research Letters, 37, Artn L23701, 
!! Doi 10.1029/2010gl045338
!! - Bounoua, L., Collatz, G. J., Sellers, P. J., Randall, D. A., Dazlich, D. A., Los, S. O., Berry, J. A., Fung, I., 
!! Tucker, C. J., Field, C. B., and Jensen, T. G. (1999), Interactions between vegetation and climate: Radiative and physiological 
!! effects of doubled atmospheric co2, Journal of Climate, 12, 309-324, Doi 10.1175/1520-0442(1999)012<0309:Ibvacr>2.0.Co;2
!! - Sellers, P. J., Bounoua, L., Collatz, G. J., Randall, D. A., Dazlich, D. A., Los, S. O., Berry, J. A., Fung, I., 
!! Tucker, C. J., Field, C. B., and Jensen, T. G. (1996), Comparison of radiative and physiological effects of doubled atmospheric
!! co2 on climate, Science, 271, 1402-1406, DOI 10.1126/science.271.5254.1402
!! - Lewis, J. D., Ward, J. K., and Tissue, D. T. (2010), Phosphorus supply drives nonlinear responses of cottonwood 
!! (populus deltoides) to increases in co2 concentration from glacial to future concentrations, New Phytologist, 187, 438-448, 
!! DOI 10.1111/j.1469-8137.2010.03307.x
!! - Kattge, J., Knorr, W., Raddatz, T., and Wirth, C. (2009), Quantifying photosynthetic capacity and its relationship to leaf 
!! nitrogen content for global-scale terrestrial biosphere models, Global Change Biology, 15, 976-991, 
!! DOI 10.1111/j.1365-2486.2008.01744.x
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel, &
                                assim_param, Ca, &
                                veget, veget_max, lai, qsintveg, qsintmax, vbeta3, vbeta3pot, rveget, rstruct, &
                                cimean, gsmean, gpp, &
                                vbeta23, hist_id, indexveg, indexlai, index, kjit, cim)

    !
    !! 0. Variable and parameter declaration
    !

    !
    !! 0.1 Input variables
    !
    INTEGER(i_std), INTENT(in)                               :: kjpindex         !! Domain size (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swdown           !! Downwelling short wave flux 
                                                                                 !! @tex ($W m^{-2}$) @endtex 
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: pb               !! Lowest level pressure (hPa)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qsurf            !! Near surface specific humidity 
                                                                                 !! @tex ($kg kg^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair             !! Specific humidity at first atmospheric model layer
                                                                                 !! @tex ($kg kg^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_air         !! Air temperature at first atmospheric model layer (K)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_growth      !! Growth temperature (°C) - Is equal to t2m_month
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: rau              !! air density @tex ($kg m^{-3}$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u                !! Lowest level wind speed 
                                                                                 !! @tex ($m s^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: v                !! Lowest level wind speed 
                                                                                 !! @tex ($m s^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q_cdrag          !! Surface drag coefficient  (-)
    REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in)  :: assim_param      !! min+max+opt temps (K), vcmax, vjmax for 
                                                                                 !! photosynthesis 
                                                                                 !! @tex ($\mu mol m^{-2} s^{-1}$) @endtex 
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: Ca               !! CO2 concentration inside the canopy
                                                                                 !! @tex ($\mu mol mol^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: humrel           !! Soil moisture stress (0-1,unitless) 
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: veget            !! Coverage fraction of vegetation for each PFT 
                                                                                 !! depending on LAI (0-1, unitless) 
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: veget_max        !! Maximum vegetation fraction of each PFT inside 
                                                                                 !! the grid box (0-1, unitless) 
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: lai              !! Leaf area index @tex ($m^2 m^{-2}$) @endtex
                                                                                 !! @tex ($m s^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: qsintveg         !! Water on vegetation due to interception 
                                                                                 !! @tex ($kg m^{-2}$) @endte
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: qsintmax         !! Maximum water on vegetation
                                                                                 !! @tex ($kg m^{-2}$) @endtex
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: vbeta23          !! Beta for fraction of wetted foliage that will 
                                                                                 !! transpire (unitless) 
    INTEGER(i_std),INTENT (in)                               :: hist_id          !! _History_ file identifier (-)    
    INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg       !! Indeces of the points on the 3D map (-)   
    INTEGER(i_std),DIMENSION (kjpindex*(nlai+1)), INTENT (in) :: indexlai  !! Indeces of the points on the 3D map
INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)   :: index            !! Indeces of the points on the map (-)
    INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number (-)        
    !
    !! 0.2 Output variables
    !
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: vbeta3           !! Beta for Transpiration (unitless)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: vbeta3pot        !! Beta for Potential Transpiration
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: rveget           !! stomatal resistance of vegetation 
                                                                                 !! @tex ($s m^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: rstruct          !! structural resistance @tex ($s m^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: cimean           !! mean intercellular CO2 concentration 
                                                                                 !! @tex ($\mu mol mol^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: gsmean           !! mean stomatal conductance to CO2 (umol m-2 s-1)
    REAL(r_Std),DIMENSION (kjpindex,nvm), INTENT (out)       :: gpp              !! Assimilation ((gC m^{-2} dt_sechiba^{-1}), total area)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: cim              !! Intercellular CO2 over nlai 
    !
    !! 0.3 Modified variables
    !
 
    !
    !! 0.4 Local variables
    !
    REAL(r_std),DIMENSION (kjpindex,nvm)  :: vcmax                               !! maximum rate of carboxylation 
                                                                                 !! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
    INTEGER(i_std)                        :: ji, jv, jl, limit_photo             !! indices (unitless)
    REAL(r_std), DIMENSION(kjpindex,nlai+1) :: info_limitphoto
    REAL(r_std), DIMENSION(kjpindex,nvm,nlai)  :: leaf_ci                        !! intercellular CO2 concentration (ppm)
    REAL(r_std), DIMENSION(kjpindex)      :: leaf_ci_lowest                      !! intercellular CO2 concentration at the lowest 
                                                                                 !! LAI level
                                                                                 !! @tex ($\mu mol mol^{-1}$) @endtex
    INTEGER(i_std), DIMENSION(kjpindex)   :: ilai                                !! counter for loops on LAI levels (unitless)
    REAL(r_std), DIMENSION(kjpindex)      :: zqsvegrap                           !! relative water quantity in the water 
                                                                                 !! interception reservoir (0-1,unitless) 
    REAL(r_std)                           :: speed                               !! wind speed @tex ($m s^{-1}$) @endtex
    ! Assimilation
    LOGICAL, DIMENSION(kjpindex)          :: assimilate                          !! where assimilation is to be calculated 
                                                                                 !! (unitless) 
    LOGICAL, DIMENSION(kjpindex)          :: calculate                           !! where assimilation is to be calculated for 
                                                                                 !! in the PFTs loop (unitless) 
    INTEGER(i_std)                        :: nic,inic,icinic                     !! counter/indices (unitless)
    INTEGER(i_std), DIMENSION(kjpindex)   :: index_calc                          !! index (unitless)
    INTEGER(i_std)                        :: nia,inia,nina,inina,iainia          !! counter/indices (unitless)
    INTEGER(i_std), DIMENSION(kjpindex)   :: index_assi,index_non_assi           !! indices (unitless)
    REAL(r_std), DIMENSION(kjpindex, nlai+1)      :: vc2                                 !! rate of carboxylation (at a specific LAI level) 
                                                                                 !! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex 
    REAL(r_std), DIMENSION(kjpindex, nlai+1)      :: vj2                                 !! rate of Rubisco regeneration (at a specific LAI 
                                                                                 !! level) @tex ($\mu mol e- m^{-2} s^{-1}$) @endtex 
    REAL(r_std), DIMENSION(kjpindex, nlai+1)      :: assimi                              !! assimilation (at a specific LAI level) 
                                                                                 !! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
                                                                                 !! (temporary variables)
    REAL(r_std), DIMENSION(kjpindex)             :: gstop                        !! stomatal conductance to H2O at topmost level 
                                                                                 !! @tex ($m s^{-1}$) @endtex
    REAL(r_std), DIMENSION(kjpindex,nlai+1)      :: gs                                  !! stomatal conductance to CO2 
                                                                                 !! @tex ($\mol m^{-2} s^{-1}$) @endtex
    REAL(r_std), DIMENSION(kjpindex,nlai)      :: templeafci  


    REAL(r_std), DIMENSION(kjpindex)      :: gamma_star                          !! CO2 compensation point (ppm)
                                                                                 !! @tex ($\mu mol mol^{-1}$) @endtex

    REAL(r_std), DIMENSION(kjpindex)      :: air_relhum                          !! air relative humidity at 2m 
                                                                                 !! @tex ($kg kg^{-1}$) @endtex
    REAL(r_std), DIMENSION(kjpindex)      :: VPD                                 !! Vapor Pressure Deficit (kPa)
    REAL(r_std), DIMENSION(kjpindex)      :: water_lim                           !! water limitation factor (0-1,unitless)

    REAL(r_std), DIMENSION(kjpindex)      :: gstot                               !! total stomatal conductance to H2O
                                                                                 !! Final unit is
                                                                                 !! @tex ($m s^{-1}$) @endtex
    REAL(r_std), DIMENSION(kjpindex)      :: assimtot                            !! total assimilation 
                                                                                 !! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
    REAL(r_std), DIMENSION(kjpindex)      :: Rdtot                               !! Total Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 m−2 s−1) 
    REAL(r_std), DIMENSION(kjpindex)      :: leaf_gs_top                         !! leaf stomatal conductance to H2O at topmost level 
                                                                                 !! @tex ($\mol H2O m^{-2} s^{-1}$) @endtex
    REAL(r_std), DIMENSION(nlai+1)        :: laitab                              !! tabulated LAI steps @tex ($m^2 m^{-2}$) @endtex
    REAL(r_std), DIMENSION(kjpindex)      :: qsatt                               !! surface saturated humidity at 2m (??) 
                                                                                 !! @tex ($g g^{-1}$) @endtex
    REAL(r_std), DIMENSION(nvm,nlai)      :: light                               !! fraction of light that gets through upper LAI    
                                                                                 !! levels (0-1,unitless)
    REAL(r_std), DIMENSION(kjpindex)      :: T_Vcmax                             !! Temperature dependance of Vcmax (unitless)
    REAL(r_std), DIMENSION(kjpindex)      :: S_Vcmax_acclim_temp                 !! Entropy term for Vcmax 
                                                                                 !! accounting for acclimation to temperature (J K-1 mol-1)
    REAL(r_std), DIMENSION(kjpindex)      :: T_Jmax                              !! Temperature dependance of Jmax
    REAL(r_std), DIMENSION(kjpindex)      :: S_Jmax_acclim_temp                  !! Entropy term for Jmax 
                                                                                 !! accounting for acclimation toxs temperature (J K-1 mol-1)
    REAL(r_std), DIMENSION(kjpindex)      :: T_gm                                !! Temperature dependance of gmw
    REAL(r_std), DIMENSION(kjpindex)      :: T_Rd                                !! Temperature dependance of Rd (unitless)
    REAL(r_std), DIMENSION(kjpindex)      :: T_Kmc                               !! Temperature dependance of KmC (unitless)
    REAL(r_std), DIMENSION(kjpindex)      :: T_KmO                               !! Temperature dependance of KmO (unitless)
    REAL(r_std), DIMENSION(kjpindex)      :: T_Sco                               !! Temperature dependance of Sco
    REAL(r_std), DIMENSION(kjpindex)      :: T_gamma_star                        !! Temperature dependance of gamma_star (unitless)    
    REAL(r_std), DIMENSION(kjpindex)      :: vc                                  !! Maximum rate of Rubisco activity-limited carboxylation (mumol CO2 m−2 s−1)
    REAL(r_std), DIMENSION(kjpindex)      :: vj                                  !! Maximum rate of e- transport under saturated light (mumol CO2 m−2 s−1)
    REAL(r_std), DIMENSION(kjpindex)      :: gm                                  !! Mesophyll diffusion conductance (molCO2 m−2 s−1 bar−1)
    REAL(r_std), DIMENSION(kjpindex)      :: g0var
    REAL(r_std), DIMENSION(kjpindex,nlai+1)      :: Rd                                  !! Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 m−2 s−1)
    REAL(r_std), DIMENSION(kjpindex)      :: Kmc                                 !! Michaelis–Menten constant of Rubisco for CO2 (mubar)
    REAL(r_std), DIMENSION(kjpindex)      :: KmO                                 !! Michaelis–Menten constant of Rubisco for O2 (mubar)
    REAL(r_std), DIMENSION(kjpindex)      :: Sco                                 !! Relative CO2 /O2 specificity factor for Rubisco (bar bar-1)
    REAL(r_std), DIMENSION(kjpindex)      :: gb_co2                              !! Boundary-layer conductance (molCO2 m−2 s−1 bar−1)
    REAL(r_std), DIMENSION(kjpindex)      :: gb_h2o                              !! Boundary-layer conductance (molH2O m−2 s−1 bar−1)
    REAL(r_std), DIMENSION(kjpindex)      :: fvpd                                !! Factor for describing the effect of leaf-to-air vapour difference on gs (-)
    REAL(r_std), DIMENSION(kjpindex)      :: low_gamma_star                      !! Half of the reciprocal of Sc/o (bar bar-1)
    REAL(r_std)                           :: N_Vcmax                             !! Nitrogen level dependance of Vcmacx and Jmax 
    REAL(r_std)                           :: fcyc                                !! Fraction of electrons at PSI that follow cyclic transport around PSI (-)
    REAL(r_std)                           :: z                                   !! A lumped parameter (see Yin et al. 2009) ( mol mol-1)                          
    REAL(r_std)                           :: Rm                                  !! Day respiration in the mesophyll (umol CO2 m−2 s−1)
    REAL(r_std)                           :: Cs_star                             !! Cs -based CO2 compensation point in the absence of Rd (ubar)
    REAL(r_std), DIMENSION(kjpindex)      :: Iabs                                !! Photon flux density absorbed by leaf photosynthetic pigments (umol photon m−2 s−1)
    REAL(r_std), DIMENSION(kjpindex)      :: Jmax                                !! Maximum value of J under saturated light (umol e− m−2 s−1)
    REAL(r_std), DIMENSION(kjpindex, nlai+1)      :: JJ                                  !! Rate of e− transport (umol e− m−2 s−1)
    REAL(r_std)                           :: J2                                  !! Rate of all e− transport through PSII (umol e− m−2 s−1)
    REAL(r_std)                           :: VpJ2                                !! e− transport-limited PEP carboxylation rate (umol CO2 m−2 s−1)
    REAL(r_std)                           :: A_1, A_3                            !! Lowest First and third roots of the analytical solution for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 m−2 s−1)
    REAL(r_std)                           :: A_1_tmp, A_3_tmp                            !! Temporary First and third roots of the analytical solution for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 m−2 s−1)
    REAL(r_std)                           :: Obs                                 !! Bundle-sheath oxygen partial pressure (ubar)
    REAL(r_std), DIMENSION(kjpindex, nlai+1)      :: Cc                                  !! Chloroplast CO2 partial pressure (ubar)
    REAL(r_std)                           :: ci_star                             !! Ci -based CO2 compensation point in the absence of Rd (ubar)        
    REAL(r_std)                           :: a,b,c,d,m,f,j,g,h,i,l,p,q,r         !! Variables used for solving the cubic equation (see Yin et al. (2009))
    REAL(r_std)                           :: QQ,UU,PSI,x1,x2,x3                      !! Variables used for solving the cubic equation (see Yin et al. (2009))
                                        
    REAL(r_std)                           :: cresist                             !! coefficient for resistances (??)
    REAL(r_std), DIMENSION(kjpindex)      :: laisum                              !! when calculating cim over nlai

! @defgroup Photosynthesis Photosynthesis
! @{   
    ! 1. Preliminary calculations\n
!_ ================================================================================================================================

    cim(:,:)=zero
    leaf_ci(:,:,:) = zero

    !
    ! 1.1 Calculate LAI steps\n
    ! The integration at the canopy level is done over nlai fixed levels.
    !! \latexonly
    !! \input{diffuco_trans_co2_1.1.tex}
    !! \endlatexonly
! @}
! @codeinc
    DO jl = 1, nlai+1
      laitab(jl) = laimax*(EXP(lai_level_depth*REAL(jl-1,r_std))-1.)/(EXP(lai_level_depth*REAL(nlai,r_std))-un)
    ENDDO
! @endcodeinc

! @addtogroup Photosynthesis
! @{   
    !
    ! 1.2 Calculate light fraction for each LAI step\n
    ! The available light follows a simple Beer extinction law. 
    ! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.
    !! \latexonly
    !! \input{diffuco_trans_co2_1.2.tex}
    !! \endlatexonly
! @}
! @codeinc
    DO jl = 1, nlai
      DO jv = 1, nvm
        light(jv,jl) = exp( -ext_coeff(jv)*laitab(jl) )
      ENDDO
    ENDDO 
! @endcodeinc
    !
    ! Photosynthesis parameters
    !

    ! Choice of downregulation. Note that downregulation_co2_new excludes
    IF (downregulation_co2_new) THEN
       ! Option used for CMIP6 version from 6.1.11 
       ! A minimum value is used for the CO2 concentration(Ca)
       ! For low CO2 concentration values(under downregulation_co2_minimum) the
       ! parametrization allows a behavior in the same way as for the preindustral period.
       DO jv= 1, nvm
           vcmax(:,jv) = assim_param(:,jv,ivcmax)*(un-downregulation_co2_coeff_new(jv) * &
                         (MAX(Ca(:),downregulation_co2_minimum)-downregulation_co2_baselevel) / &
                         (MAX(Ca(:),downregulation_co2_minimum)+20.))
       ENDDO
    ELSE IF (downregulation_co2) THEN
       ! Option used for CMIP6 version 6.1.0 up to 6.1.10
       DO jv= 1, nvm
          vcmax(:,jv) = assim_param(:,jv,ivcmax)*(un-downregulation_co2_coeff(jv) * &
                          log(Ca(:)/downregulation_co2_baselevel))
       ENDDO
    ELSE
       vcmax(:,:) = assim_param(:,:,ivcmax)
    ENDIF

!    DO jv = 1, nvm
!       vcmax(:,:) = Vcmax25(jv)
!    ENDDO

! @addtogroup Photosynthesis
! @{   
    !
    ! 1.3 Estimate relative humidity of air (for calculation of the stomatal conductance).\n
    !! \latexonly
    !! \input{diffuco_trans_co2_1.3.tex}
    !! \endlatexonly
! @}
    !

    CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
    air_relhum(:) = &
      ( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) ) / &
      ( qsatt(:)*pb(:) / (Tetens_1+qsatt(:)*Tetens_2 ) )


    VPD(:) = ( qsatt(:)*pb(:) / (Tetens_1+qsatt(:)*Tetens_2 ) ) &
         - ( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) )
    ! VPD is needed in kPa
    VPD(:) = VPD(:)/10.

    !
    ! 2. beta coefficient for vegetation transpiration
    !
    rstruct(:,1) = rstruct_const(1)
    rveget(:,:) = undef_sechiba
    !
    vbeta3(:,:) = zero
    vbeta3pot(:,:) = zero
    gsmean(:,:) = zero
    gpp(:,:) = zero
    !
    cimean(:,1) = Ca(:)
    !
    ! @addtogroup Photosynthesis
    ! @{   
    ! 2. Loop over vegetation types\n
    ! @} 
    !
    DO jv = 2,nvm
       gamma_star(:)=zero
       Kmo(:)=zero
       Kmc(:)=zero
       gm(:)=zero
       g0var(:) =zero

       Cc(:,:)=zero
       Vc2(:,:)=zero
       JJ(:,:)=zero
       info_limitphoto(:,:)=zero
       gs(:,:)=zero
       templeafci(:,:)=zero
       assimi(:,:)=zero
       Rd(:,:)=zero

      !
      ! @addtogroup Photosynthesis
      ! @{   
      !
      ! 2.1 Initializations\n
      !! \latexonly
      !! \input{diffuco_trans_co2_2.1.tex}
      !! \endlatexonly
      ! @}      
      !
      ! beta coefficient for vegetation transpiration
      !
      rstruct(:,jv) = rstruct_const(jv)
      cimean(:,jv) = Ca(:)
      !
      !! mask that contains points where there is photosynthesis
      !! For the sake of vectorisation [DISPENSABLE], computations are done only for convenient points.
      !! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not
      !! calculated (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation, 
      !! temperature and relative humidity).
      !! For the points where assimilation is not calculated, variables are initialized to specific values. 
      !! The assimilate(kjpindex) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation.
      !! The index_assi(kjpindex) array indexes the nia points with assimilation, whereas the index_no_assi(kjpindex) array indexes
      !! the nina points with no assimilation.
      nia=0
      nina=0
      !
      DO ji=1,kjpindex
         !
         IF ( ( lai(ji,jv) .GT. 0.01 ) .AND. &
              ( veget_max(ji,jv) .GT. min_sechiba ) ) THEN
            
            IF ( ( veget(ji,jv) .GT. min_sechiba ) .AND. &
                 ( swdown(ji) .GT. min_sechiba )   .AND. &
                 ( humrel(ji,jv) .GT. min_sechiba) .AND. &
                 ( temp_growth(ji) .GT. tphoto_min(jv) ) .AND. &
                 ( temp_growth(ji) .LT. tphoto_max(jv) ) ) THEN
               !
               assimilate(ji) = .TRUE.
               nia=nia+1
               index_assi(nia)=ji
               !
            ELSE
               !
               assimilate(ji) = .FALSE.
               nina=nina+1
               index_non_assi(nina)=ji
               !
            ENDIF
         ELSE
            !
            assimilate(ji) = .FALSE.
            nina=nina+1
            index_non_assi(nina)=ji
            !
         ENDIF
         !

      ENDDO
      !

      gstot(:) = zero
      gstop(:) = zero
      assimtot(:) = zero
      Rdtot(:)=zero
      leaf_gs_top(:) = zero
      !
      zqsvegrap(:) = zero
      WHERE (qsintmax(:,jv) .GT. min_sechiba)
      !! relative water quantity in the water interception reservoir
          zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv))
      ENDWHERE
      !
      !! Calculates the water limitation factor.
      water_lim(:) = humrel(:,jv)

      ! give a default value of ci for all pixel that do not assimilate
      DO jl=1,nlai
         DO inina=1,nina
            leaf_ci(index_non_assi(inina),jv,jl) = Ca(index_non_assi(inina)) 
         ENDDO
      ENDDO
      !
      ilai(:) = 1
      !
      ! Here is the calculation of assimilation and stomatal conductance
      ! based on the work of Farquahr, von Caemmerer and Berry (FvCB model) 
      ! as described in Yin et al. 2009
      ! Yin et al. developed a extended version of the FvCB model for C4 plants
      ! and proposed an analytical solution for both photosynthesis pathways (C3 and C4)
      ! Photosynthetic parameters used are those reported in Yin et al. 
      ! Except For Vcmax25, relationships between Vcmax25 and Jmax25 for which we use 
      ! Medlyn et al. (2002) and Kattge & Knorr (2007)
      ! Because these 2 references do not consider mesophyll conductance, we neglect this term
      ! in the formulations developed by Yin et al. 
      ! Consequently, gm (the mesophyll conductance) tends to the infinite
      ! This is of importance because as stated by Kattge & Knorr and Medlyn et al.,
      ! values of Vcmax and Jmax derived with different model parametrizations are not 
      ! directly comparable and the published values of Vcmax and Jmax had to be standardized
      ! to one consistent formulation and parametrization

      ! See eq. 6 of Yin et al. (2009)
      ! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001)
      T_KmC(:)        = Arrhenius(kjpindex,temp_air,298.,E_KmC(jv))
      T_KmO(:)        = Arrhenius(kjpindex,temp_air,298.,E_KmO(jv))
      T_Sco(:)        = Arrhenius(kjpindex,temp_air,298.,E_Sco(jv))
      T_gamma_star(:) = Arrhenius(kjpindex,temp_air,298.,E_gamma_star(jv))


      ! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001)
      T_Rd(:)         = Arrhenius(kjpindex,temp_air,298.,E_Rd(jv))


      ! For C3 plants, we assume that the Entropy term for Vcmax and Jmax 
      ! acclimates to temperature as shown by Kattge & Knorr (2007) - Eq. 9 and 10
      ! and that Jmax and Vcmax respond to temperature following a modified Arrhenius function
      ! (with a decrease of these parameters for high temperature) as in Medlyn et al. (2002) 
      ! and Kattge & Knorr (2007).
      ! In Yin et al. (2009), temperature dependance to Vcmax is based only on a Arrhenius function
      ! Concerning this apparent unconsistency, have a look to the section 'Limitation of 
      ! Photosynthesis by gm' of Bernacchi (2002) that may provide an explanation
      
      ! Growth temperature tested by Kattge & Knorr range from 11 to 35°C
      ! So, we limit the relationship between these lower and upper limits
      S_Jmax_acclim_temp(:) = aSJ(jv) + bSJ(jv) * MAX(11., MIN(temp_growth(:),35.))      
      T_Jmax(:)  = Arrhenius_modified(kjpindex,temp_air,298.,E_Jmax(jv),D_Jmax(jv),S_Jmax_acclim_temp)

      S_Vcmax_acclim_temp(:) = aSV(jv) + bSV(jv) * MAX(11., MIN(temp_growth(:),35.))   
      T_Vcmax(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_Vcmax(jv),D_Vcmax(jv),S_Vcmax_acclim_temp)
       

      
      vc(:) = vcmax(:,jv) * T_Vcmax(:)
      ! As shown by Kattge & Knorr (2007), we make use
      ! of Jmax25/Vcmax25 ratio (rJV) that acclimates to temperature for C3 plants
      ! rJV is written as a function of the growth temperature
      ! rJV = arJV + brJV * T_month 
      ! See eq. 10 of Kattge & Knorr (2007)
      ! and Table 3 for Values of arJV anf brJV 
      ! Growth temperature is monthly temperature (expressed in °C) - See first paragraph of
      ! section Methods/Data of Kattge & Knorr
      vj(:) = ( arJV(jv) + brJV(jv) *  MAX(11., MIN(temp_growth(:),35.)) ) * vcmax(:,jv) * T_Jmax(:)

      T_gm(:)    = Arrhenius_modified(kjpindex,temp_air,298.,E_gm(jv),D_gm(jv),S_gm(jv))
      gm(:) = gm25(jv) * T_gm(:) * MAX(1-stress_gm(jv), water_lim(:))

      g0var(:) = g0(jv)* MAX(1-stress_gs(jv), water_lim(:))
      ! @endcodeinc
      !
      KmC(:)=KmC25(jv)*T_KmC(:)
      KmO(:)=KmO25(jv)*T_KmO(:)
      Sco(:)=Sco25(jv)*T_sco(:)
      gamma_star(:) = gamma_star25(jv)*T_gamma_star(:)



      ! low_gamma_star is defined by Yin et al. (2009)
      ! as the half of the reciprocal of Sco - See Table 2
      low_gamma_star(:) = 0.5 / Sco(:)

      ! VPD expressed in kPa
      ! Note : MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:)))) is always between 0-1 not including 0 and 1
      fvpd(:) = 1. / ( 1. / MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:)))) - 1. ) &
                * MAX(1-stress_gs(jv), water_lim(:))

      ! leaf boundary layer conductance 
      ! conversion from a conductance in (m s-1) to (mol H2O m-2 s-1)
      ! from Pearcy et al. (1991, see below)
      gb_h2o(:) = gb_ref * 44.6 * (tp_00/temp_air(:)) * (pb(:)/pb_std) 

      ! conversion from (mol H2O m-2 s-1) to (mol CO2 m-2 s-1)
      gb_co2(:) = gb_h2o(:) / ratio_H2O_to_CO2

      !
      ! @addtogroup Photosynthesis
      ! @{   
      !
      ! 2.4 Loop over LAI discretized levels to estimate assimilation and conductance\n
      ! @}           
      !
      !! The calculate(kjpindex) array is of type logical to indicate wether we have to sum over this LAI fixed level (the LAI of
      !! the point for the PFT is lower or equal to the LAI level value). The number of such points is incremented in nic and the 
      !! corresponding point is indexed in the index_calc array.
      JJ(:,:)=zero
      vc2(:,:)=zero
      vj2(:,:)=zero
      Cc(:,:)=zero
      gs(:,:)=zero
      assimi(:,:)=zero
      Rd(:,:)=zero

      DO jl = 1, nlai
         !
         nic=0
         calculate(:) = .FALSE.
         !
         IF (nia .GT. 0) then
            DO inia=1,nia
               calculate(index_assi(inia)) = (laitab(jl) .LE. lai(index_assi(inia),jv) )
               IF ( calculate(index_assi(inia)) ) THEN
                  nic=nic+1
                  index_calc(nic)=index_assi(inia)
               ENDIF
            ENDDO
         ENDIF
         !
         ! @addtogroup Photosynthesis
         ! @{   
         !
         ! 2.4.1 Vmax is scaled into the canopy due to reduction of nitrogen 
         !! (Johnson and Thornley,1984).\n
         !! \latexonly
         !! \input{diffuco_trans_co2_2.4.1.tex}
         !! \endlatexonly
         ! @}           
         !
         N_Vcmax = ( un - .7_r_std * ( un - light(jv,jl) ) )
         !

         vc2(:,jl) = vc(:) * N_Vcmax * MAX(1-stress_vcmax(jv), water_lim(:))
         vj2(:,jl) = vj(:) * N_Vcmax * MAX(1-stress_vcmax(jv), water_lim(:))

         ! see Comment in legend of Fig. 6 of Yin et al. (2009)
         ! Rd25 is assumed to equal 0.01 Vcmax25 
         Rd(:,jl) = vcmax(:,jv) * N_Vcmax * 0.01 * T_Rd(:)  * MAX(1-stress_vcmax(jv), water_lim(:))

         Iabs(:)=swdown(:)*W_to_mol*RG_to_PAR*ext_coeff(jv)*light(jv,jl)
         
         ! eq. 4 of Yin et al (2009)
         Jmax(:)=vj2(:,jl)
         JJ(:,jl) = ( alpha_LL(jv) * Iabs(:) + Jmax(:) - sqrt((alpha_LL(jv) * Iabs(:) + Jmax(:) )**2. &
              - 4 * theta(jv) * Jmax(:) * alpha_LL(jv) * Iabs(:)) ) &
              / ( 2 * theta(jv))

         !
         IF ( is_c4(jv) )  THEN
            !
            ! @addtogroup Photosynthesis
            ! @{   
            !
            ! 2.4.2 Assimilation for C4 plants (Collatz et al., 1992)\n
            !! \latexonly
            !! \input{diffuco_trans_co2_2.4.2.tex}
            !! \endlatexonly
            ! @}           
            !
            !
            !
            IF (nic .GT. 0) THEN
               DO inic=1,nic

                  ! Analytical resolution of the Assimilation based Yin et al. (2009)
                  icinic=index_calc(inic)

                  ! Eq. 28 of Yin et al. (2009)
                  fcyc= 1. - ( 4.*(1.-fpsir(jv))*(1.+fQ(jv)) + 3.*h_protons(jv)*fpseudo(jv) ) / &
                       ( 3.*h_protons(jv) - 4.*(1.-fpsir(jv)))
                                    
                  ! See paragraph after eq. (20b) of Yin et al.
                  Rm=Rd(icinic,jl)/2.
                                
                  ! We assume that cs_star equals ci_star (see Comment in legend of Fig. 6 of Yin et al. (2009)
                  ! Equation 26 of Yin et al. (2009)
                  Cs_star = (gbs(jv) * low_gamma_star(icinic) * Oi - &
                       ( 1. + low_gamma_star(icinic) * alpha(jv) / 0.047) * Rd(icinic,jl) + Rm ) &
                       / ( gbs(jv) + kp(jv) ) 

                  ! eq. 11 of Yin et al (2009)
                  J2 = JJ(icinic,jl) / ( 1. - fpseudo(jv) / ( 1. - fcyc ) )

                  ! Equation right after eq. (20d) of Yin et al. (2009)
                  z = ( 2. + fQ(jv) - fcyc ) / ( h_protons(jv) * (1. - fcyc ))

                  VpJ2 = fpsir(jv) * J2 * z / 2.

                  A_3=9999.

                  ! See eq. right after eq. 18 of Yin et al. (2009)
                  DO limit_photo=1,2
                     ! Is Vc limiting the Assimilation
                     IF ( limit_photo .EQ. 1 ) THEN
                        a = 1. + kp(jv) / gbs(jv)
                        b = 0.
                        x1 = Vc2(icinic,jl)
                        x2 = KmC(icinic)/KmO(icinic)
                        x3 = KmC(icinic)
                        ! Is J limiting the Assimilation
                     ELSE
                        a = 1.
                        b = VpJ2
                        x1 = (1.- fpsir(jv)) * J2 * z / 3.
                        x2 = 7. * low_gamma_star(icinic) / 3.
                        x3 = 0.
                     ENDIF

                     m=fvpd(icinic)-g0var(icinic)/gb_co2(icinic)
                     d=g0var(icinic)*(Ca(icinic)-Cs_star) + fvpd(icinic)*Rd(icinic,jl)
                     f=(b-Rm-low_gamma_star(icinic)*Oi*gbs(jv))*x1*d + a*gbs(jv)*x1*Ca(icinic)*d
                     j=(b-Rm+gbs(jv)*x3 + x2*gbs(jv)*Oi)*m + (alpha(jv)*x2/0.047-1.)*d &
                          + a*gbs(jv)*(Ca(icinic)*m - d/gb_co2(icinic) - (Ca(icinic) - Cs_star ))
 
                     g=(b-Rm-low_gamma_star(icinic)*Oi*gbs(jv))*x1*m - (alpha(jv)*low_gamma_star(icinic)/0.047+1.)*x1*d &
                          + a*gbs(jv)*x1*(Ca(icinic)*m - d/gb_co2(icinic) - (Ca(icinic)-Cs_star ))
 
                     h=-((alpha(jv)*low_gamma_star(icinic)/0.047+1.)*x1*m + (a*gbs(jv)*x1*(m-1.))/gb_co2(icinic) )
                     i= ( b-Rm + gbs(jv)*x3 + x2*gbs(jv)*Oi )*d + a*gbs(jv)*Ca(icinic)*d
                     l= ( alpha(jv)*x2/0.047 - 1.)*m - (a*gbs(jv)*(m-1.))/gb_co2(icinic)
 
                     p = (j-(h-l*Rd(icinic,jl))) / l
                     q = (i+j*Rd(icinic,jl)-g) / l
                     r = -(f-i*Rd(icinic,jl)) / l 
 
                     ! See Yin et al. (2009) and  Baldocchi (1994)
                     QQ = ( (p**2._r_std) - 3._r_std * q) / 9._r_std
                     UU = ( 2._r_std* (p**3._r_std) - 9._r_std *p*q + 27._r_std *r) /54._r_std

                     IF ( (QQ .GE. 0._r_std) .AND. (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) ) THEN
                        PSI = ACOS(UU/(QQ**1.5_r_std))
                        A_3_tmp = -2._r_std * SQRT(QQ) * COS(( PSI + 4._r_std * PI)/3._r_std ) - p / 3._r_std
                        IF (( A_3_tmp .LT. A_3 )) THEN
                           A_3 = A_3_tmp
                           info_limitphoto(icinic,jl)=2.
                        ELSE
                        ! In case, J is not limiting the assimilation
                        ! we have to re-initialise a, b, x1, x2 and x3 values
                        ! in agreement with a Vc-limited assimilation 
                           a = 1. + kp(jv) / gbs(jv)
                           b = 0.
                           x1 = Vc2(icinic,jl)
                           x2 = KmC(icinic)/KmO(icinic)
                           x3 = KmC(icinic)
                           info_limitphoto(icinic,jl)=1.
                        ENDIF
                     ENDIF

                     IF ( ( A_3 .EQ. 9999. ) .OR. ( A_3 .LT. (-Rd(icinic,jl)) ) ) THEN
                        IF ( printlev>=4 ) THEN
                           WRITE(numout,*) 'We have a problem in diffuco_trans_co2 for A_3'
                           WRITE(numout,*) 'no real positive solution found for pft:',jv
                           WRITE(numout,*) 'temp_air:',temp_air(icinic)
                           WRITE(numout,*) 'vpd:',vpd(icinic)
                        END IF
                        A_3 = -Rd(icinic,jl)
                     ENDIF
                     assimi(icinic,jl) = A_3

                     IF ( ABS( assimi(icinic,jl) + Rd(icinic,jl) ) .LT. min_sechiba ) THEN
                        gs(icinic,jl) = g0var(icinic)
                        !leaf_ci keeps its initial value (Ca).
                     ELSE
                        ! Eq. 24 of Yin et al. (2009) 
                        Obs = ( alpha(jv) * assimi(icinic,jl) ) / ( 0.047 * gbs(jv) ) + Oi
                        ! Eq. 23 of Yin et al. (2009)
                        Cc(icinic,jl) = ( ( assimi(icinic,jl) + Rd(icinic,jl) ) * ( x2 * Obs + x3 ) + low_gamma_star(icinic) &
                             * Obs * x1 ) &
                             / MAX(min_sechiba, x1 - ( assimi(icinic,jl) + Rd(icinic,jl) ))
                        ! Eq. 22 of Yin et al. (2009)
                        leaf_ci(icinic,jv,jl) = ( Cc(icinic,jl) - ( b - assimi(icinic,jl) - Rm ) / gbs(jv) ) / a
                        ! Eq. 25 of Yin et al. (2009)
                        ! It should be Cs instead of Ca but it seems that 
                        ! other equations in Appendix C make use of Ca
                        gs(icinic,jl) = g0var(icinic) + ( assimi(icinic,jl) + Rd(icinic,jl) ) / &
                             ( Ca(icinic) - Cs_star ) * fvpd(icinic)             
                     ENDIF
                  ENDDO                  
               ENDDO
            ENDIF
         ELSE
            !
            ! @addtogroup Photosynthesis
            ! @{   
            !
            ! 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980)\n
            !! \latexonly
            !! \input{diffuco_trans_co2_2.4.3.tex}
            !! \endlatexonly
            ! @}           
            !
            !
            IF (nic .GT. 0) THEN
               DO inic=1,nic
                  icinic=index_calc(inic)
			
                  A_1=9999.

                  ! See eq. right after eq. 18 of Yin et al. (2009)
                  DO limit_photo=1,2
                     ! Is Vc limiting the Assimilation
                     IF ( limit_photo .EQ. 1 ) THEN
                        x1 = vc2(icinic,jl)
                        ! It should be O not Oi (comment from Vuichard)
                        x2 = KmC(icinic) * ( 1. + 2*gamma_star(icinic)*Sco(icinic) / KmO(icinic) )
                        ! Is J limiting the Assimilation
                     ELSE
                        x1 = JJ(icinic,jl)/4.
                        x2 = 2. * gamma_star(icinic)
                     ENDIF


                     ! See Appendix B of Yin et al. (2009)
                     a = g0var(icinic) * ( x2 + gamma_star(icinic) ) + &
                          ( g0var(icinic) / gm(icinic) + fvpd(icinic) ) * ( x1 - Rd(icinic,jl) )
                     b = Ca(icinic) * ( x1 - Rd(icinic,jl) ) - gamma_star(icinic) * x1 - Rd(icinic,jl) * x2
                     c = Ca(icinic) + x2 + ( 1./gm(icinic) + 1./gb_co2(icinic) ) * ( x1 - Rd(icinic,jl) ) 
                     d = x2 + gamma_star(icinic) + ( x1 - Rd(icinic,jl) ) / gm(icinic)
                     m = 1./gm(icinic) + ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * ( 1./gm(icinic) + 1./gb_co2(icinic) )  
   
                     p = - ( d + (x1 - Rd(icinic,jl) ) / gm(icinic) + a * ( 1./gm(icinic) + 1./gb_co2(icinic) ) + &
                          ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * c ) / m
   
                     q = ( d * ( x1 - Rd(icinic,jl) ) + a*c + ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * b ) / m
                     r = - a * b / m
   
                     ! See Yin et al. (2009) 
                     QQ = ( (p**2._r_std) - 3._r_std * q) / 9._r_std
                     UU = ( 2._r_std* (p**3._r_std) - 9._r_std *p*q + 27._r_std *r) /54._r_std
               
                     IF ( (QQ .GE. 0._r_std) .AND. (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) ) THEN
                        PSI = ACOS(UU/(QQ**1.5_r_std))
                        A_1_tmp = -2._r_std * SQRT(QQ) * COS( PSI / 3._r_std ) - p / 3._r_std
                        IF (( A_1_tmp .LT. A_1 )) THEN
                           A_1 = A_1_tmp
                           info_limitphoto(icinic,jl)=2.
                        ELSE
                        ! In case, J is not limiting the assimilation
                        ! we have to re-initialise x1 and x2 values
                        ! in agreement with a Vc-limited assimilation 
                           x1 = vc2(icinic,jl)
                           ! It should be O not Oi (comment from Vuichard)
                           x2 = KmC(icinic) * ( 1. + 2*gamma_star(icinic)*Sco(icinic) / KmO(icinic) )                           
                           info_limitphoto(icinic,jl)=1.
                        ENDIF
                     ENDIF
                  ENDDO
                  IF ( (A_1 .EQ. 9999.) .OR. ( A_1 .LT. (-Rd(icinic,jl)) ) ) THEN
                     IF ( printlev>=4 ) THEN
                        WRITE(numout,*) 'We have a problem in diffuco_trans_co2 for A_1'
                        WRITE(numout,*) 'no real positive solution found for pft:',jv
                        WRITE(numout,*) 'temp_air:',temp_air(icinic)
                        WRITE(numout,*) 'vpd:',vpd(icinic)
                     END IF
                     A_1 = -Rd(icinic,jl)
                  ENDIF
                  assimi(icinic,jl) = A_1

                  IF ( ABS( assimi(icinic,jl) + Rd(icinic,jl) ) .LT. min_sechiba ) THEN
                     gs(icinic,jl) = g0var(icinic)
                  ELSE
                     ! Eq. 18 of Yin et al. (2009)
                     Cc(icinic,jl) = ( gamma_star(icinic) * x1 + ( assimi(icinic,jl) + Rd(icinic,jl) ) * x2 )  &
                          / MAX( min_sechiba, x1 - ( assimi(icinic,jl) + Rd(icinic,jl) ) )
                     ! Eq. 17 of Yin et al. (2009)
                     leaf_ci(icinic,jv,jl) = Cc(icinic,jl) + assimi(icinic,jl) / gm(icinic) 
                     ! See eq. right after eq. 15 of Yin et al. (2009)
                     ci_star = gamma_star(icinic) - Rd(icinic,jl) / gm(icinic)
                     ! 
                     ! Eq. 15 of Yin et al. (2009)
                     gs(icinic,jl) = g0var(icinic) + ( assimi(icinic,jl) + Rd(icinic,jl) ) / ( leaf_ci(icinic,jv,jl) &
                          - ci_star ) * fvpd(icinic)
                  ENDIF
               ENDDO
            ENDIF
         ENDIF
         !
         IF (nic .GT. 0) THEN
            !
            DO inic=1,nic
               !
               ! @addtogroup Photosynthesis
               ! @{   
               !
               !! 2.4.4 Estimatation of the stomatal conductance (Ball et al., 1987).\n
               !! \latexonly
               !! \input{diffuco_trans_co2_2.4.4.tex}
               !! \endlatexonly
               ! @}           
               !
               icinic=index_calc(inic)
               !
               ! keep stomatal conductance of topmost level
               !
               IF ( jl .EQ. 1 ) THEN
                  leaf_gs_top(icinic) = gs(icinic,jl)
                  !
               ENDIF
               !
               ! @addtogroup Photosynthesis
               ! @{   
               !
               !! 2.4.5 Integration at the canopy level\n
               !! \latexonly
               !! \input{diffuco_trans_co2_2.4.5.tex}
               !! \endlatexonly
               ! @}           
               ! total assimilation and conductance
               assimtot(icinic) = assimtot(icinic) + &
                    assimi(icinic,jl) * (laitab(jl+1)-laitab(jl))
               Rdtot(icinic) = Rdtot(icinic) + &
                    Rd(icinic,jl) * (laitab(jl+1)-laitab(jl))
               gstot(icinic) = gstot(icinic) + &
                    gs(icinic,jl) * (laitab(jl+1)-laitab(jl))
               !
               ilai(icinic) = jl
               !
            ENDDO
            !
         ENDIF
      ENDDO  ! loop over LAI steps
 
      IF(jv==testpft) THEN
         templeafci(:,:)=leaf_ci(:,testpft,:)
         CALL histwrite_p(hist_id, 'Cc', kjit, Cc, kjpindex*(nlai+1), indexlai)
         CALL histwrite_p(hist_id, 'Vc', kjit, Vc2, kjpindex*(nlai+1), indexlai)
         CALL histwrite_p(hist_id, 'Vj', kjit, JJ, kjpindex*(nlai+1), indexlai)
         CALL histwrite_p(hist_id, 'limitphoto', kjit, info_limitphoto, kjpindex*(nlai+1), indexlai)
         CALL histwrite_p(hist_id, 'gammastar', kjit, gamma_star, kjpindex,index)
         CALL histwrite_p(hist_id, 'Kmo', kjit, Kmo, kjpindex,index)
         CALL histwrite_p(hist_id, 'Kmc', kjit, Kmc, kjpindex,index)
         CALL histwrite_p(hist_id, 'gm', kjit, gm, kjpindex, index)
         CALL histwrite_p(hist_id, 'gs', kjit, gs, kjpindex*(nlai+1), indexlai)
         CALL histwrite_p(hist_id, 'leafci', kjit, templeafci, kjpindex*(nlai), indexlai)
         CALL histwrite_p(hist_id, 'assimi', kjit, assimi, kjpindex*(nlai+1), indexlai)
         CALL histwrite_p(hist_id, 'Rd', kjit, Rd, kjpindex*(nlai+1), indexlai)
      ENDIF
      !! Calculated intercellular CO2 over nlai needed for the chemistry module
      cim(:,jv)=0.
      laisum(:)=0
      DO jl=1,nlai
         WHERE (laitab(jl) .LE. lai(:,jv) )
            cim(:,jv)= cim(:,jv)+leaf_ci(:,jv,jl)*(laitab(jl+1)-laitab(jl))
            laisum(:)=laisum(:)+ (laitab(jl+1)-laitab(jl))
         ENDWHERE
      ENDDO
      WHERE (laisum(:)>0)
         cim(:,jv)= cim(:,jv)/laisum(:)
      ENDWHERE


      !
      !! 2.5 Calculate resistances
      !
      IF (nia .GT. 0) THEN
         !
         DO inia=1,nia
            !
            iainia=index_assi(inia)

            !! Mean stomatal conductance for CO2 (mol m-2 s-1)
            gsmean(iainia,jv) = gstot(iainia)
            !
            ! cimean is the "mean ci" calculated in such a way that assimilation 
            ! calculated in enerbil is equivalent to assimtot
            !
            IF ( ABS(gsmean(iainia,jv)-g0var(iainia)*laisum(iainia)) .GT. min_sechiba) THEN
               cimean(iainia,jv) = (fvpd(iainia)*(assimtot(iainia)+Rdtot(iainia))) /&
                 (gsmean(iainia,jv)-g0var(iainia)*laisum(iainia)) + gamma_star(iainia) 
            ELSE
               cimean(iainia,jv) = gamma_star(iainia) 
            ENDIF
                 
            ! conversion from umol m-2 (PFT) s-1 to gC m-2 (mesh area) tstep-1
            gpp(iainia,jv) = assimtot(iainia)*12e-6*veget_max(iainia,jv)*dt_sechiba
            
            !
            ! conversion from mol/m^2/s to m/s
            !
            ! As in Pearcy, Schulze and Zimmermann
            ! Measurement of transpiration and leaf conductance
            ! Chapter 8 of Plant Physiological Ecology
            ! Field methods and instrumentation, 1991
            ! Editors:
            !
            !    Robert W. Pearcy,
            !    James R. Ehleringer,
            !    Harold A. Mooney,
            !    Philip W. Rundel
            !
            ! ISBN: 978-0-412-40730-7 (Print) 978-94-010-9013-1 (Online)

            gstot(iainia) =  mol_to_m_1 *(temp_air(iainia)/tp_00)*&
                 (pb_std/pb(iainia))*gstot(iainia)*ratio_H2O_to_CO2
            gstop(iainia) =  mol_to_m_1 * (temp_air(iainia)/tp_00)*&
                 (pb_std/pb(iainia))*leaf_gs_top(iainia)*ratio_H2O_to_CO2*&
                 laitab(ilai(iainia)+1)
            !
            rveget(iainia,jv) = un/gstop(iainia)

            !
            !
            ! rstruct is the difference between rtot (=1./gstot) and rveget
            !
            ! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif
            !rstruct(iainia,jv) = un/gstot(iainia) - &
            !     rveget(iainia,jv)
            rstruct(iainia,jv) = MAX( un/gstot(iainia) - &
                 rveget(iainia,jv), min_sechiba)
            !
            !
            !! wind is a global variable of the diffuco module.
            speed = MAX(min_wind, wind(iainia))
            !
            ! beta for transpiration
            !
            ! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin
            !! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode.
            !! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module.
            !vbeta3(iainia,jv) = veget_max(iainia,jv) * &
            !  (un - zqsvegrap(iainia)) * &
            !  (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
            !   rstruct(iainia,jv))))
            !! Global resistance of the canopy to evaporation
            cresist=(un / (un + speed * q_cdrag(iainia) * &
                 veget(iainia,jv)/veget_max(iainia,jv) * &
                 (rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv)))))

            IF ( humrel(iainia,jv) >= min_sechiba ) THEN
               vbeta3(iainia,jv) = veget(iainia,jv) * &
                    (un - zqsvegrap(iainia)) * cresist + &
                    MIN( vbeta23(iainia,jv), veget(iainia,jv) * &
                    zqsvegrap(iainia) * cresist )
            ELSE
               ! Because of a minimum conductance g0, vbeta3 cannot be zero even if humrel=0
               ! in the above equation.
               ! Here, we force transpiration to be zero when the soil cannot deliver it
               vbeta3(iainia,jv) = zero
            END IF

            ! vbeta3pot for computation of potential transpiration (needed for irrigation)
            vbeta3pot(iainia,jv) = MAX(zero, veget(iainia,jv) * cresist)
            !
            !
         ENDDO
         !
      ENDIF
      !
   END DO         ! loop over vegetation types
      
   IF (printlev>=3) WRITE (numout,*) ' diffuco_trans_co2 done '


END SUBROUTINE diffuco_trans_co2


!! ================================================================================================================================
!! SUBROUTINE	   : diffuco_comb
!!
!>\BRIEF           This routine combines the previous partial beta 
!! coefficients and calculates the total alpha and complete beta coefficients.
!!
!! DESCRIPTION	   : Those integrated coefficients are used to calculate (in enerbil.f90) the total evapotranspiration 
!! from the grid-cell. \n
!!
!! In the case that air is more humid than surface, dew deposition can occur (negative latent heat flux). 
!! In this instance, for temperature above zero, all of the beta coefficients are set to 0, except for 
!! interception (vbeta2) and bare soil (vbeta4 with zero soil resistance). The amount of water that is 
!! intercepted by leaves is calculated based on the value of LAI of the surface. In the case of freezing 
!! temperatures, water is added to the snow reservoir, and so vbeta4 and vbeta2 are set to 0, and the 
!! total vbeta is set to 1.\n
!!
!! \latexonly 
!!     \input{diffucocomb1.tex}
!! \endlatexonly
!!
!! The beta and alpha coefficients are initially set to 1.
!! \latexonly 
!!     \input{diffucocomb2.tex}
!! \endlatexonly
!!
!! If snow is lower than the critical value:
!! \latexonly 
!!     \input{diffucocomb3.tex}
!! \endlatexonly
!! If in the presence of dew:
!! \latexonly 
!!     \input{diffucocomb4.tex}
!! \endlatexonly
!!
!! Determine where the water goes (soil, vegetation, or snow)
!! when air moisture exceeds saturation.
!! \latexonly 
!!     \input{diffucocomb5.tex}
!! \endlatexonly
!!
!! If it is not freezing dew is put into the interception reservoir and onto the bare soil. If it is freezing, 
!! water is put into the snow reservoir. 
!! Now modify vbetas where necessary: for soil and snow
!! \latexonly 
!!     \input{diffucocomb6.tex}
!! \endlatexonly
!!
!! and for vegetation
!! \latexonly 
!!     \input{diffucocomb7.tex}
!! \endlatexonly
!!
!! Then compute part of dew that can be intercepted by leafs.
!!
!! There will be no transpiration when air moisture is too high, under any circumstance
!! \latexonly 
!!     \input{diffucocomb8.tex}
!! \endlatexonly
!!
!! There will also be no interception loss on bare soil, under any circumstance.
!! \latexonly 
!!     \input{diffucocomb9.tex}
!! \endlatexonly
!!
!! The flowchart details the 'decision tree' which underlies the module. 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): vbeta1, vbeta4, humrel, vbeta2, vbeta3, vbeta
!!
!! REFERENCE(S) :
!! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
!! Circulation Model. Journal of Climate, 6, pp.248-273
!! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation
!! sur le cycle de l'eau, PhD Thesis, available from:
!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
!!
!! FLOWCHART    :
!! \latexonly 
!!     \includegraphics[scale=0.25]{diffuco_comb_flowchart.png}
!! \endlatexonly
!! \n
!_ ================================================================================================================================

  SUBROUTINE diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, &
       & snow, veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3 , vbeta4, &
       & evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax)    
    
    ! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006

  !! 0. Variable and parameter declaration
    
    !! 0.1 Input variables
    
    INTEGER(i_std), INTENT(in)                           :: kjpindex   !! Domain size (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: rau        !! Air Density (kg m^{-3})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: u          !! Eastward Lowest level wind speed (m s^{-1}) 
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: v          !! Nortward Lowest level wind speed (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: q_cdrag    !! Surface drag coefficient  (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: pb         !! Lowest level pressure (hPa)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: qair       !! Lowest level specific air humidity (kg kg^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: temp_sol   !! Skin temperature (K)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: temp_air   !! Lower air temperature (K)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: snow       !! Snow mass (kg)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: veget      !! Fraction of vegetation type (fraction)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: lai        !! Leaf area index (m^2 m^{-2})
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: qsintmax   !! Maximum water on vegetation (kg m^{-2})
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)        :: tot_bare_soil!! Total evaporating bare soil fraction 
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: veget_max   !! Max. fraction of vegetation type (LAI->infty)

    !! 0.2 Output variables
    
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)       :: vbeta      !! Total beta coefficient (-)

    !! 0.3 Modified variables 
    
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)     :: vbeta1     !! Beta for sublimation (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)     :: vbeta4     !! Beta for Bare soil evaporation (-) 
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel     !! Soil moisture stress (within range 0 to 1)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2     !! Beta for interception loss (-)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3     !! Beta for Transpiration (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)     :: evap_bare_lim  !! limiting factor for bare soil evaporation  
                                                                           !! when the 11-layer hydrology is used (-)
    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout):: evap_bare_lim_ns !! limiting factor for bare soil evaporation 
                                                                             !! when the 11-layer hydrology is used (-)    
    !! 0.4 Local variables
    
    INTEGER(i_std)                                       :: ji, jv
    REAL(r_std)                                          :: zevtest, zsoil_moist, zrapp
    REAL(r_std), DIMENSION(kjpindex)                     :: qsatt
    LOGICAL, DIMENSION(kjpindex)                         :: toveg, tosnow
    REAL(r_std)                                          :: coeff_dew_veg
    REAL(r_std), DIMENSION(kjpindex)                     :: vegtot
!_ ================================================================================================================================
    
    !! 1 If we are in presence of dew
     
    CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)

    
    !! 1.1 Determine where the water goes 
    !! Determine where the water goes (soil, vegetation, or snow)
    !! when air moisture exceeds saturation.
    !! \latexonly 
    !!     \input{diffucocomb5.tex}
    !! \endlatexonly
    toveg(:) = .FALSE.
    tosnow(:) = .FALSE.
    DO ji = 1, kjpindex
      IF ( qsatt(ji) .LT. qair(ji) ) THEN
          IF (temp_air(ji) .GT. tp_00) THEN
              !! If it is not freezing dew is put into the 
              !! interception reservoir and onto the bare soil.
              toveg(ji) = .TRUE.
          ELSE
              !! If it is freezing water is put into the 
              !! snow reservoir.
              tosnow(ji) = .TRUE.
          ENDIF
      ENDIF
    END DO

    !! 1.2 Now modify vbetas where necessary.
    
    !! 1.2.1 Soil and snow
    !! \latexonly 
    !!     \input{diffucocomb6.tex}
    !! \endlatexonly

    ! We need to keep consistency between evap_bare_lim, evap_bare_lim_ns and vbeta4 (thus vevapnu)
    ! or we have a water conservation issue in hydrol_split_soil

    DO ji = 1, kjpindex

       IF ( toveg(ji) ) THEN

          vbeta1(ji) = zero
          vegtot(ji) = SUM(veget_max(ji,:))

          IF ( (tot_bare_soil(ji) .GT. min_sechiba) .AND. (vegtot(ji).GT. min_sechiba) ) THEN
             
             vbeta4(ji) = tot_bare_soil(ji)
             
             ! We now have to redefine evap_bare_lim & evap_bare_lim_ns
             IF (evap_bare_lim(ji) .GT. min_sechiba) THEN               
                evap_bare_lim_ns(ji,:) = evap_bare_lim_ns(ji,:) * vbeta4(ji) / evap_bare_lim(ji)
             ELSE ! we must re-invent evap_bare_lim_ns => uniform across soiltiles             
                evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji)               
             ENDIF

             evap_bare_lim(ji) = vbeta4(ji)
             ! consistent with evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)*soiltile(ji,:)*vegtot(ji))
             ! as SUM(soiltile(ji,:)) = 1

          ELSE          
             vbeta4(ji) = zero
             evap_bare_lim_ns(ji,:) = zero
             evap_bare_lim(ji) = zero
          ENDIF
       ENDIF

       IF ( tosnow(ji) ) THEN
          vbeta1(ji) = un
          vbeta4(ji) = zero
          evap_bare_lim_ns(ji,:) = zero
          evap_bare_lim(ji) = zero
       ENDIF

    ENDDO

    !! 1.2.2 Vegetation and interception loss
    !! \latexonly 
    !!     \input{diffucocomb7.tex}
    !! \endlatexonly
    DO jv = 1, nvm
      
      DO ji = 1, kjpindex
               
        IF ( toveg(ji) ) THEN
           IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
              
              ! Compute part of dew that can be intercepted by leafs.
              IF ( lai(ji,jv) .GT. min_sechiba) THEN
                IF (lai(ji,jv) .GT. 1.5) THEN
                   coeff_dew_veg= &
                         &   dew_veg_poly_coeff(6)*lai(ji,jv)**5 &
                         & - dew_veg_poly_coeff(5)*lai(ji,jv)**4 &
                         & + dew_veg_poly_coeff(4)*lai(ji,jv)**3 &
                         & - dew_veg_poly_coeff(3)*lai(ji,jv)**2 &
                         & + dew_veg_poly_coeff(2)*lai(ji,jv) &
                         & + dew_veg_poly_coeff(1)
                 ELSE
                    coeff_dew_veg=un
                 ENDIF
              ELSE
                 coeff_dew_veg=zero
              ENDIF
              IF (jv .EQ. 1) THEN
                 ! This line may not work with CWRR when frac_bare is distributed among three soiltiles
                 ! Fortunately, qsintmax(ji,1)=0 (LAI=0 in PFT1) so we never pass here
                 vbeta2(ji,jv) = coeff_dew_veg*tot_bare_soil(ji)
              ELSE
                 vbeta2(ji,jv) = coeff_dew_veg*veget(ji,jv)
              ENDIF
           ELSE
              vbeta2(ji,jv) = zero ! if qsintmax=0, vbeta2=0
           ENDIF
        ENDIF
        IF ( tosnow(ji) ) vbeta2(ji,jv) = zero
        
      ENDDO
      
    ENDDO

    !! 1.2.3 Vegetation and transpiration  
    !! There will be no transpiration when air moisture is too high, under any circumstance
    !! \latexonly 
    !!     \input{diffucocomb8.tex}
    !! \endlatexonly
    DO jv = 1, nvm
      DO ji = 1, kjpindex
        IF ( qsatt(ji) .LT. qair(ji) ) THEN
          vbeta3(ji,jv) = zero
          humrel(ji,jv) = zero
        ENDIF
      ENDDO
    ENDDO
    
    
    !! 1.2.4 Overrules 1.2.2
    !! There will also be no interception loss on bare soil, under any circumstance.
    !! \latexonly 
    !!     \input{diffucocomb9.tex}
    !! \endlatexonly
    DO ji = 1, kjpindex
       IF ( qsatt(ji) .LT. qair(ji) ) THEN
          vbeta2(ji,1) = zero
       ENDIF
    ENDDO

    !! 2  Now calculate vbeta in all cases (the equality needs to hold for enerbil to be consistent)

    DO ji = 1, kjpindex
          vbeta(ji) = vbeta4(ji) + SUM(vbeta2(ji,:)) + SUM(vbeta3(ji,:))

          IF (vbeta(ji) .LT. min_sechiba) THEN
             vbeta(ji) = zero
             vbeta4(ji) = zero
             vbeta2(ji,:)= zero
             vbeta3(ji,:)= zero
             evap_bare_lim_ns(ji,:) = zero
             evap_bare_lim(ji) = zero
          END IF
    ENDDO 

    CALL xios_orchidee_send_field("evap_bare_lim",evap_bare_lim) 
    CALL xios_orchidee_send_field("evap_bare_lim_ns",evap_bare_lim_ns)

    IF (printlev>=3) WRITE (numout,*) ' diffuco_comb done '

  END SUBROUTINE diffuco_comb


!! ================================================================================================================================
!! SUBROUTINE	: diffuco_raerod
!!
!>\BRIEF	Computes the aerodynamic resistance, for cases in which the
!! surface drag coefficient is provided by the coupled atmospheric model LMDZ and  when the flag
!! 'ldq_cdrag_from_gcm' is set to TRUE
!!
!! DESCRIPTION	: Simply computes the aerodynamic resistance, for cases in which the
!! surface drag coefficient is provided by the coupled atmospheric model LMDZ. If the surface drag coefficient
!! is not provided by the LMDZ or signalled by the flag 'ldq_cdrag_from_gcm' set to FALSE, then the subroutine
!! diffuco_aero is called instead of this one.
!!
!! Calculation of the aerodynamic resistance, for diganostic purposes. First calculate wind speed:
!! \latexonly 
!!     \input{diffucoaerod1.tex}
!! \endlatexonly       
!!
!! next calculate ::raero
!! \latexonly 
!!     \input{diffucoaerod2.tex}
!! \endlatexonly
!! 
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::raero
!!
!! REFERENCE(S)	:
!! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
!! Circulation Model. Journal of Climate, 6, pp.248-273
!! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influence de l'irrigation
!! sur le cycle de l'eau, PhD Thesis, available from:
!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
!!
!! FLOWCHART    :  None
!! \n
!_ ================================================================================================================================

  SUBROUTINE diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
    
    IMPLICIT NONE
    
  !! 0. Variable and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                     :: kjpindex     !! Domain size (-)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: u            !! Eastward Lowest level wind velocity (m s^{-1}) 
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: v            !! Northward Lowest level wind velocity (m s^{-1})
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: q_cdrag      !! Surface drag coefficient  (-)
    
    !! 0.2 Output variables 
    
    REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: raero        !! Aerodynamic resistance (s m^{-1})
     
    !! 0.3 Modified variables

    !! 0.4 Local variables
    
    INTEGER(i_std)                                 :: ji           !! (-)
    REAL(r_std)                                    :: speed        !! (m s^{-1})
!_ ================================================================================================================================
   
  !! 1. Simple calculation of the aerodynamic resistance, for diganostic purposes.

    DO ji=1,kjpindex

       !! \latexonly 
       !!     \input{diffucoaerod1.tex}
       !! \endlatexonly       
       speed = MAX(min_wind, wind(ji))

       !! \latexonly 
       !!     \input{diffucoaerod2.tex}
       !! \endlatexonly
       raero(ji) = un / (q_cdrag(ji)*speed)
       
    ENDDO
  
  END SUBROUTINE diffuco_raerod


  FUNCTION Arrhenius (kjpindex,temp,ref_temp,energy_act) RESULT ( val_arrhenius )
    !! 0.1 Input variables

    INTEGER(i_std),INTENT(in)                     :: kjpindex          !! Domain size (-)
    REAL(r_std),DIMENSION(kjpindex),INTENT(in)    :: temp              !! Temperature (K)
    REAL(r_std), INTENT(in)                       :: ref_temp          !! Temperature of reference (K)
    REAL(r_std),INTENT(in)                        :: energy_act        !! Activation Energy (J mol-1)
    
    !! 0.2 Result

    REAL(r_std), DIMENSION(kjpindex)              :: val_arrhenius     !! Temperature dependance based on
                                                                       !! a Arrhenius function (-)
    
    val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*RR*(temp(:))))
  END FUNCTION Arrhenius

  FUNCTION Arrhenius_modified_1d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
    !! 0.1 Input variables

    INTEGER(i_std),INTENT(in)                     :: kjpindex          !! Domain size (-)
    REAL(r_std),DIMENSION(kjpindex),INTENT(in)    :: temp              !! Temperature (K)
    REAL(r_std), INTENT(in)                       :: ref_temp          !! Temperature of reference (K)
    REAL(r_std),INTENT(in)                        :: energy_act        !! Activation Energy (J mol-1)
    REAL(r_std),INTENT(in)                        :: energy_deact      !! Deactivation Energy (J mol-1)
    REAL(r_std),DIMENSION(kjpindex),INTENT(in)    :: entropy           !! Entropy term (J K-1 mol-1)
        
    !! 0.2 Result

    REAL(r_std), DIMENSION(kjpindex)              :: val_arrhenius     !! Temperature dependance based on
                                                                       !! a Arrhenius function (-)
    
    val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*RR*(temp(:))))  &
         * (1. + EXP( (ref_temp * entropy(:) - energy_deact) / (ref_temp * RR ))) &
         / (1. + EXP( (temp(:) * entropy(:) - energy_deact) / ( RR*temp(:))))
         
  END FUNCTION Arrhenius_modified_1d

  FUNCTION Arrhenius_modified_0d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
    !! 0.1 Input variables

    INTEGER(i_std),INTENT(in)                     :: kjpindex          !! Domain size (-)
    REAL(r_std),DIMENSION(kjpindex),INTENT(in)    :: temp              !! Temperature (K)
    REAL(r_std), INTENT(in)                       :: ref_temp          !! Temperature of reference (K)
    REAL(r_std),INTENT(in)                        :: energy_act        !! Activation Energy (J mol-1)
    REAL(r_std),INTENT(in)                        :: energy_deact      !! Deactivation Energy (J mol-1)
    REAL(r_std),INTENT(in)                        :: entropy           !! Entropy term (J K-1 mol-1)
        
    !! 0.2 Result

    REAL(r_std), DIMENSION(kjpindex)              :: val_arrhenius     !! Temperature dependance based on
                                                                       !! a Arrhenius function (-)
    
    val_arrhenius(:)=EXP(((temp(:)-ref_temp)*energy_act)/(ref_temp*RR*(temp(:))))  &
         * (1. + EXP( (ref_temp * entropy - energy_deact) / (ref_temp * RR ))) &
         / (1. + EXP( (temp(:) * entropy - energy_deact) / ( RR*temp(:))))
         
  END FUNCTION Arrhenius_modified_0d


END MODULE diffuco
