!==================================================================================================================================
!  MODULE 	: hydraulic_arch
! 
!  CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
!
!  LICENCE      : IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF This module calculates the amount of water that is available for the leaves for transpiration,
!        accounting for effects of hydraulic architecture following (Hickler et al. 2006)      
!!
!!\n DESCRIPTION  : (1) Calculate soil water potential from soil water content as calculated in hydrol.f90
!!                  according to van Genuchten (1980). Soil water potential is weighted by the relative share of the
!!                  of theroots in each layer.
!!                  (2) Calculate pressure difference between leaves and soil (Whitehead 1998)
!!                  (3) Calculate root (Weatherly 1998), sapwood (Magnani et al. 2000, Tyree et al. 1994,
!!                  Sperry et al. 1998) and leave resistances (Sitch et al. 2003.
!!                  (4) Calculate the supply of water for transpiration by the leaves using Darcy's law (Whitehead 1998)
!! 
!! RECENT CHANGE(S):(1) Soil to root resistance has been added to the calculations of the soil water potential weighted 
!!                  by root density. The development was made by Emilie Joetzjer based on 
!!                  Williams et al., 2001 (tree physiology) and Duursma et al., 2012 (GMD). Also see Bonan et al. 2014 (GMD)
!! 
!! REFERENCE(S) : None
!!   
!! SVN     :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/trunk/ORCHIDEE/src_sechiba/hydraulic_arch.f90 $ 
!! $Date: 2023-11-22 20:28:54 +0100 (mer. 22 nov. 2023) $
!! $Revision: 8316 $
!! \n
!_ ================================================================================================================================
 
MODULE hydraulic_arch
 
  USE ioipsl
  USE xios_orchidee
 
  ! modules used :
  USE constantes
  USE constantes_soil
  USE constantes_soil_var
  USE pft_parameters
  USE dynamic_parameters
  USE ioipsl_para
  USE function_library, ONLY: wood_to_height_eff, cc_to_lai, & 
                              biomass_to_coupled_lai, biomass_to_lai, &
                              Arrhenius, Arrhenius_modified, &
                              Arrhenius_modified_nodim, get_printlev
  USE qsat_moisture
  USE mleb
  USE enerbil, ONLY: enerbil_main, enerbil_write

  IMPLICIT NONE

  ! Private and public routines
  PRIVATE
  PUBLIC hydraulic_arch_main, mleb_diffuco_trans_co2_short


  LOGICAL, SAVE                               :: l_first_sechiba = .TRUE.      !! Flag controlling the intialisation (true/false)
!$OMP THREADPRIVATE(l_first_sechiba)
  INTEGER(i_std), SAVE                        :: printlev_loc                  !! Local printlev in slowproc module
!$OMP THREADPRIVATE(printlev_loc)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: nvan                !! VanGenuchten coeficients n (unitless) ! RK:
!$OMP THREADPRIVATE(nvan)                                                 
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: avan                !! Van Genuchten coeficients a!@tex!$(mm^{-1})$!@endtex
!$OMP THREADPRIVATE(avan)                                                
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcr                 !! Residual volumetric water content !! !@tex !$(m^{3} !m^{-3})$!@endtex
!$OMP THREADPRIVATE(mcr)                                                 
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcs                 !! Saturated volumetric water content !@tex !$(m^{3} !m^{-3})$!@endtex
!$OMP THREADPRIVATE(mcs)

CONTAINS



!!  =============================================================================================================================
!! SUBROUTINE:    hydraulic_arch_main
!!
!>\BRIEF:	      call the relevant modules for the multi-layer energy budget and hydraulic stress
!!
!! DESCRIPTION:	  call the relevant modules for the multi-layer energy budget and hydraulic stress
!! This subroutine was preivously stored in sechiba.f90 and was named sechiba_mleb_hs
!! \n
!_ ==============================================================================================================================
  SUBROUTINE hydraulic_arch_main( &
       kjit,              kjpij,            kjpindex,     rest_id,        &
       rest_id_stom,      hist_id,          hist2_id,     hist_id_stom,   &
       hist_id_stom_IPCC, index,            indexveg,     ldrestart_read, &
       ldrestart_write,   njsc,             u,            v,              &
       qair,              precip_rain,      precip_snow,  lwdown,         &
       swnet,             swdown,           temp_air,                     &
       petAcoef,          peqAcoef,         petBcoef,     peqBcoef,       &
       pb,                tq_cdrag,         epot_air,     sinang,         &
       ccanopy,           emis,             soilcap,      soilflx,        &
       rau,               temp_growth,      vbeta23,      humrel,         &
       lai,               max_height_store, gs_distribution,  gs_diffuco_output, &
       gstot_component,   gstot_frac,       assim_param,   circ_class_biomass,   &
       circ_class_n,      ksoil,            lai_per_level, laieff_isotrop,  &
       Light_Abs_Tot,     Light_Tran_Tot,   qsintmax,      qsintveg,        &
       snowdz,            stempdiag,        swc,                            &
       frac_sub_snow,     frac_evap_intercept,             resist_intercept,&
       frac_evap_bare,    frac_evap_flood,  raero,                          &
       veget_max,         veget,            z_array_out,                    & 
       temp_sol_new,      tsol_rad,         temp_sol_add,  vevapflo,        &
       vevapnu,           vevapsno,         stressed,      unstressed,      &
       transpot,          vevapwet,         transpir,                       &
       e_frac,                                                              &
       vbeta2sum,         vbeta3sum,        fluxsens,      fluxlat,         &
       vevapp,            resist_equiv_evap_veg,                            &
       resist_equiv_tot,  resist_equiv_evap,                                &
       evapot,            evapot_corr,                                      &
       qsurf,             temp_sol,         pgflux,        u_speed,         &
       frac_evap_transp,  resist_transp,    vbeta3pot,     rveget, rstruct, &
       cimean,            gsmean,           gpp,           transpir_supply, &
       leaf_ci,           warnings,         JJ_out,        assimi_lev,      &
       cresist,           info_limitphoto,  vessel_loss,   root_profile,    &
       profile_vbeta3,    profile_rveget)

!! 0.1 Input variables    
    INTEGER(i_std), INTENT(in)                               :: kjit              !! Time step number (unitless)
    INTEGER(i_std), INTENT(in)                               :: kjpij             !! Total size of the un-compressed grid 
    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only 
    INTEGER(i_std), INTENT(in)                               :: rest_id           !! _Restart_ file identifier (unitless)
    INTEGER(i_std), INTENT(in)                               :: rest_id_stom      !! STOMATE's _Restart_ file identifier 
    INTEGER(i_std), INTENT(in)                               :: hist_id           !! _History_ file identifier (unitless)
    INTEGER(i_std), INTENT(in)                               :: hist2_id          !! _History_ file 2 identifier (unitless)
    INTEGER(i_std), INTENT(in)                               :: hist_id_stom      !! STOMATE's _History_ file identifier 
    INTEGER(i_std), INTENT(in)                               :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file identifier
    INTEGER(i_std), DIMENSION(kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map. 
    INTEGER(i_std), DIMENSION(kjpindex*nvm), INTENT(in)      :: indexveg          !! Indeces of the points on the 3D map
    LOGICAL, INTENT(in)                                      :: ldrestart_read    !! Logical for _restart_ file to read 
    LOGICAL, INTENT(in)                                      :: ldrestart_write   !! Logical for _restart_ file to write 
    INTEGER(i_std), DIMENSION(kjpindex), INTENT (in)         :: njsc              !! Index of the dominant soil textural class 
                                                                                  !! in the grid cell used for hydrology, (1-nscm, unitless)
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: u                 !! Lowest level wind speed in direction u 
                                                                                  !! @tex $(m.s^{-1})$ @endtex 
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: v                 !! Lowest level wind speed in direction v 
                                                                                  !! @tex $(m.s^{-1})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: qair              !! Lowest level specific humidity 
                                                                                  !! @tex $(kg kg^{-1})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: precip_rain       !! Rain precipitation 
                                                                                  !! @tex $(kg m^{-2})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: precip_snow       !! Snow precipitation 
                                                                                  !! @tex $(kg m^{-2})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: lwdown            !! Down-welling long-wave flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: swnet             !! Net surface short-wave flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: swdown            !! Down-welling surface short-wave flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: temp_air          !! Air temperature (K)
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: petAcoef          !! Coefficients A for T from the Planetary 
                                                                                  !! Boundary Layer
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: peqAcoef          !! Coefficients A for q from the Planetary 
                                                                                  !! Boundary Layer
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: petBcoef          !! Coefficients B for T from the Planetary 
                                                                                  !! Boundary Layer
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: peqBcoef          !! Coefficients B for q from the Planetary 
                                                                                  !! Boundary Layer
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: pb                !! Surface pressure (hPa)
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: tq_cdrag          !! Surface drag coefficient (-) 
                                                                                  !! @tex $(m.s^{-1})$ @endtex 
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: epot_air          !! Air potential energy (??J)
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: sinang            !! Sine of the solar angle (unitless)
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: ccanopy           !! CO2 concentration in the canopy (ppm)
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: emis              !! 
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: soilcap           !! Soil calorific capacity (J K^{-1])
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: soilflx           !! Soil flux (W m^{-2})
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: rau               !! Air density (kg m^{-3})
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: temp_growth       !! Growth temperature (°C) - Is equal to t2m_month
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in)        :: vbeta23           !! Beta for fraction of wetted foliage that will
                                                                                  !! transpire once intercepted water has evaporated (-)
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in)        :: humrel          !! Soil moisture availability [0-1]
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in)         :: lai               !!
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)         :: max_height_store  !!
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in)  ::assimi_lev
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot+1), INTENT(in)::info_limitphoto !! Stored for mleb_diffuco_
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in)  ::JJ_out
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in)  :: gs_distribution   !!
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in)  :: gs_diffuco_output !!
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in)  :: gstot_component   !!
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in)  :: gstot_frac        !!
    REAL(r_std), DIMENSION(kjpindex,nvm,npco2), INTENT(in)        :: assim_param       !!
    REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nparts,nelements), INTENT(in)  :: circ_class_biomass   
    REAL(r_std), DIMENSION(kjpindex,nvm,ncirc), INTENT(in)                   :: circ_class_n
    REAL(r_std), DIMENSION(kjpindex,nslm,nstm), INTENT (in)       :: ksoil            !! Hydraulic conductivity mm/d  
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in)  :: lai_per_level    !! This is the LAI per vertical level
    REAL(r_std), DIMENSION(kjpindex,nlevels_tot,nvm), INTENT(in)  :: laieff_isotrop   !! Effective LAI
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT (in) :: Light_Abs_Tot    !! Absorbed radiation per level for photosynthesis
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT (in) :: Light_Tran_Tot   !! Absorbed radiation per level for photosynthesis
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in)        :: qsintmax              !! Maximum water on vegetation
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in)        :: qsintveg              !! Water on vegetation due to interception 
    REAL(r_std), DIMENSION(kjpindex,nsnow),INTENT(in)        :: snowdz                !! Snow depth at each snow layer [m]
    REAL(r_std), DIMENSION(kjpindex,nslm),INTENT(in)         :: stempdiag             !! Soil temperature (K)
    REAL(r_std), DIMENSION(kjpindex,nslm,nstm), INTENT (in)  :: swc                   !! Volumetric soil water content
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: frac_sub_snow         !! Fraction of gridcell on which snow sublimation occurs (-)
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in)        :: frac_evap_intercept   !! Fraction of gridcell on which intercepted water evaporation occurs (-)  
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in)        :: resist_intercept      !! Sum of resistances applied to intercepted water evaporation (s m^{-1})   
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: frac_evap_bare        !! Fraction of gridcell on which bare soil evaporation occurs (-)
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: frac_evap_flood       !! Fraction of gridcell on which floodplains evaporation occurs (-)
    REAL(r_std), DIMENSION(kjpindex), INTENT (in)            :: raero                 !! Aerodynamic resistance
    REAL(r_std), DIMENSION(kjpindex, nvm), INTENT(in)        :: veget_max             !! Max. fraction of vegetation type (LAI -> infty, unitless)
    REAL(r_std), DIMENSION(kjpindex, nvm), INTENT(in)        :: veget                 !! Vegetation fraction for each pft 
    REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nlevels_tot), INTENT(inout) :: z_array_out !! An output of h_array, to use in sechiba
    REAL(r_std),DIMENSION (kjpindex,nvm,nslm), INTENT(in)    :: root_profile          !! Normalized root mass/length fraction in each soil layer 
                                                                                      !! (0-1, unitless)

!! 0.2 Output variables
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)  :: cresist                   !! coefficient for resistances (??)
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)        :: temp_sol_new              !! New ground temperature (K)
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)        :: tsol_rad                  !! Radiative surface temperature (W/m2) 
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)        :: temp_sol_add              !! Additional energy to melt snow for snow ablation case (K) 
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)        :: vevapflo                  !! Floodplains evaporation                        
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)        :: vevapnu                   !! Bare soil evaporation (mm day^{-1})            
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)        :: vevapsno                  !! Snow evaporation (mm day^{-1})               
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out)    :: stressed                  !! GPP when the PFT is stressed (g C m-2 dt-1)
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out)    :: unstressed                !! GPP when the PFT is unstressed (g C m-2 dt-1)
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out)    :: transpot                  !! Potential transpiration                      
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out)    :: vevapwet                  !! Interception (mm day^{-1})                     
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out)    :: transpir                  !! Transpiration (mm day^{-1})                  
    REAL(r_std), DIMENSION(kjpindex,nvm,nslm,nstm), INTENT(out) :: e_frac             !! Fraction of water transpired supplied by individual layers (no units)            
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out)    :: vessel_loss               !! Proportion of conductivity lost due to cavitation in the xylem (no unit)

!! 0.3 Modified variables    
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: vbeta2sum             !! sum of vbeta2 coefficients across all PFTs (-)
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: vbeta3sum             !! sum of vbeta3 coefficients across all PFTs (-)
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: fluxsens              !! Sensible heat flux 
                                                                                      !! @tex $(W m^{-2})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: fluxlat               !! Latent heat flux 
                                                                                      !! @tex $(W m^{-2})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: vevapp                !! Total of evaporation 
                                                                                      !! @tex $(kg m^{-2} days^{-1})$ @endtex
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: resist_equiv_evap_veg !! Equivalent resistance to the fluxes involved in the evaporation from the vegetation (bare soil evaporation + interception loss + transpiration) (s m^{-1})
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: resist_equiv_tot      !! Equivalent resistance to the fluxes involved in the total evaporative (all evaporations + snow sublimation) (s m^{-1})
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: resist_equiv_evap     !! Equivalent resistance to the fluxes involved in the evaporation (bare soil evaporation + interception loss + transpiration + floodplains evaporation) (s m^{-1})
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: evapot                !!
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: evapot_corr           !!
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: qsurf                 !! 
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: temp_sol              !!
    REAL(r_std), DIMENSION(kjpindex), INTENT(inout)          :: pgflux                !! Net energy into snowpack(W/m^2)
    REAL(r_std), DIMENSION(kjpindex,jnlvls),   INTENT(inout) :: u_speed               !! Canopy wind speed profile - see comment in mleb_finalize
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: frac_evap_transp      !! Fraction of grid cell where transpiration occurs (-)
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: resist_transp         !! Sum of resistances applied to intercepted water evaporation (s m^{-1})   
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: vbeta3pot             !!
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: rveget                !!
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: rstruct               !!
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: cimean                !!
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: gsmean                !!
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: gpp                   !! gC m^{-2}pixel
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)      :: transpir_supply       !! Supply of water for transpiration (mm timestep^{-1}) 
    REAL(r_std), DIMENSION(kjpindex,nvm,nlai), INTENT(inout) :: leaf_ci               !! intercellular CO2 concentration (ppm)
    REAL(r_std), DIMENSION(kjpindex,nvm,nwarns), INTENT(inout) :: warnings            !! A warning counter
    REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot),INTENT(in)   :: profile_vbeta3
    REAL(r_std), DIMENSION (kjpindex,nvm,nlevels_tot),INTENT(in)   :: profile_rveget


!! 0.4 Local variables
    ! variables introduced as a consequence of running enerbil twice (if needed, as a result of restriction of transpiration by
    !    the supply term). If running for the second time, we should start with the same input variables as the first run
    
    INTEGER(i_std)                           :: ji, jv
    REAL(r_std), DIMENSION (kjpindex)        :: evapot_save           !! Soil Potential Evaporation (mm day^{-1})
    REAL(r_std), DIMENSION (kjpindex)        :: evapot_corr_save      !! Soil Potential Evaporation Correction (mm day^{-1})
    REAL(r_std), DIMENSION (kjpindex)        :: diffevap              !! Difference betwence vevapp and composing fluxes (Kg/m^2/s)
    REAL(r_std), DIMENSION (kjpindex)        :: temp_sol_save         !! Soil temperature (K)
    REAL(r_std), DIMENSION (kjpindex)        :: qsurf_save            !! Surface specific humidity (kg kg^{-1}) 
    REAL(r_std), DIMENSION (kjpindex)        :: tsol_rad_save         !! Tsol_rad (W m^{-2})
    REAL(r_std), DIMENSION (kjpindex)        :: netrad_save           !! Net radiation flux (W m^{-2})
    REAL(r_std), DIMENSION (kjpindex)        :: vbeta_save            !! Resistance coefficient (-)
    REAL(r_std), DIMENSION (kjpindex)        :: vbeta2sum_save        !! sum of vbeta2 coefficients across all PFTs (-)
    REAL(r_std), DIMENSION (kjpindex)        :: vbeta3sum_save        !! sum of vbeta3 coefficients across all PFTs (-)
    REAL(r_std), DIMENSION (kjpindex,nvm)    :: vbeta3_save           !! Beta for Transpiration (unitless)
    REAL(r_std), DIMENSION (kjpindex,nvm)    :: vbeta3pot_save        !! Beta for Potential Transpiration
    REAL(r_std), DIMENSION (kjpindex)        :: pgflux_save           !! Net energy into snowpack(W/m^2)
    REAL(r_std), DIMENSION (kjpindex)        :: qsol_sat_new          !! New saturated surface air moisture (kg kg^{-1})
    REAL(r_std), DIMENSION (kjpindex)        :: qair_new
    REAL(r_std), DIMENSION (kjpindex)        :: vbeta3_save_test
    REAL(r_std), DIMENSION (kjpindex)        :: netrad                !! Net radiation (W m^{-2})
    REAL(r_std), DIMENSION (kjpindex)        :: lwabs                 !! LW radiation absorbed by the surface (W m^{-2})
    REAL(r_std), DIMENSION (kjpindex)        :: lwnet                 !! Net Long-wave radiation (W m^{-2})
    REAL(r_std), DIMENSION (kjpindex)        :: fluxsubli             !! Energy of sublimation (mm day^{-1})

    REAL(r_std), DIMENSION (nlevels_tot,kjpindex,nvm)   :: transpir_supply_column   !! (SECHIBA TOP LEVEL) Supply of water for transpiration
    !+++CHECK+++
    ! Variable dimensions are for a single pixel. Causing 1+1 problems
    ! when running a larger domain. Needs to be corrected when implementing
    ! a global use of the multi-layer energy budget.
    REAL(r_std), DIMENSION (nlevels_tot)                :: Light_Abs_Tot_mean    !! (SECHIBA TOP LEVEL) total light absorption for a given canopy level
    REAL(r_std), DIMENSION (nlevels_tot)                :: Light_Alb_Tot_mean    !! (SECHIBA TOP LEVEL) total albedo for a given level
    !+++++++++++

    LOGICAL                                             :: hydrol_flag           !! flag that 'trips' the energy budget for each grid square 
    LOGICAL, DIMENSION (kjpindex)                       :: hydrol_flag2          !! flag that 'trips' the energy budget for each grid square
    LOGICAL, DIMENSION (nbp_glo)                        :: hydrol_flag2_glo      !! hydrol_flag2 on the full global domain for all processors
    LOGICAL, DIMENSION (kjpindex,nvm)                   :: hydrol_flag3          !! flag that 'trips' the energy budget for each grid square and PFT

    LOGICAL, DIMENSION (kjpindex,nlevels_tot)           :: hydrol_flag2_column   !! flag that 'trips' the energy budget for each grid square, and canopy level
    LOGICAL, DIMENSION (kjpindex,nvm,nlevels_tot)       :: hydrol_flag3_column   !! flag that 'trips' the energy budget for each grid square, canopy level and PFT

    REAL(r_std), DIMENSION (kjpindex,nvm)              :: ratio_vbeta3      !! Reduction rate of vbeta3 when transpir_supply
                                                                            !! is smaller than demand

    REAL(r_std), DIMENSION(kjpindex,nvm)               :: cresist_save       !!temporal
    REAL(r_std), DIMENSION(kjpindex,nvm)               :: vbeta2, vbeta3       !!temporal
    REAL(r_std), DIMENSION(kjpindex)                   :: vbeta1, vbeta4, vbeta5, vbeta       !!temporal

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

    ! Define a set of flags, required to signal if we need to recalculate the
    ! energy budget depending on the hydraulic stress being a limiting factor to the
    ! latent heat flux.
    hydrol_flag = .FALSE.          ! trip flag per cycle (no dimension)
    hydrol_flag2(:) = .FALSE.      ! trip flag per grid square
    hydrol_flag3(:,:) = .FALSE.    ! trip flag per grid square and PFT

    ! Store vbeta3 for the calculation of ratio_vbeta3 
    !! Local variables used in mleb. Will be removed as soon as mleb is cleaned up.
    vbeta1(:) = frac_sub_snow(:)
    vbeta4(:) = frac_evap_bare(:)
    vbeta5(:) = frac_evap_flood(:)
    DO ji=1,kjpindex
       vbeta2(ji,:) = frac_evap_intercept(ji,:) * raero(ji) / resist_intercept(ji,:)
       vbeta3(ji,:) = frac_evap_transp(ji,:) * raero(ji) / resist_transp(ji,:)
       vbeta3_save(ji,:) = vbeta3(ji,:)
       vbeta(ji) = vbeta4(ji) + SUM(vbeta3(ji,:) + vbeta2(ji,:))
    ENDDO
    ratio_vbeta3(:,:) = un

    ! Store values of inout variables before calculating the energy budget
    ! so it can be recalculated with the exact same initial conditions when
    ! needed. 
    evapot_save(:) = evapot(:)
    evapot_corr_save(:) = evapot_corr(:)
    temp_sol_save(:) = temp_sol(:)
    qsurf_save(:) = qsurf(:)
    tsol_rad_save(:) = tsol_rad(:)
    pgflux_save(:) = pgflux(:)


    IF (.not.ok_mleb) THEN

       ! Calculate the energy budget for the first time. The main reason we need
       ! this calculation is to have an estimate of transpir.
       CALL enerbil_main (kjit, kjpindex, netrad, lwabs, lwnet, fluxsubli, index, &
            & indexveg, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, qair, &
            & peqAcoef, peqBcoef, pb, rau, resist_equiv_evap, resist_equiv_evap_veg, resist_equiv_tot, &
            & frac_sub_snow, frac_evap_intercept, resist_intercept, &
            & frac_evap_transp, resist_transp, vbeta3pot, frac_evap_bare, frac_evap_flood, raero, &
            & emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
            & vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, &
            & temp_sol_new, qsurf, evapot, evapot_corr, precip_rain, pgflux, &
            & snowdz, temp_sol_add, qair_new, qsol_sat_new)

    ELSE

       !+++CHECK+++
       ! Add the inout variables for mleb
       CALL ipslerr_p(3,'hydraulic_arch_main','The recalculation of the energybudget',&
            'needs to be checked when mleb rather than enerbil_main is used',&
            'Save the intent(inout) variables for mleb')
       ! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
       ! when running a larger domain. Needs to be corrected when implementing
       ! a global use of the multi-layer energy budget.
       !+++++++++++

       CALL mleb_main (kjit, kjpindex, ldrestart_read, ldrestart_write, &
            & index, indexveg, lwdown, swnet, swdown, epot_air, temp_air, &
            & u, v, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, pb, rau, &
            & vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
            & emis, soilflx, soilcap, tq_cdrag, humrel, &
            & fluxsens, fluxlat, vevapp, transpir, transpot, vevapnu, vevapwet, &
            & vevapsno, vevapflo, temp_sol, tsol_rad, temp_sol_new, qsurf, &
            & evapot, evapot_corr, rest_id, hist_id, hist2_id, & 
            & ok_mleb_history_file, &
            & transpir_supply, hydrol_flag, hydrol_flag2, hydrol_flag3, vbeta2sum, & 
            & vbeta3sum, veget_max, qsol_sat_new, qair_new, &
            & veget, &
            & Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
            & laieff_isotrop, lai_per_level,&
            & z_array_out, transpir_supply_column, u_speed, &
            & profile_vbeta3, profile_rveget, max_height_store, &
            & precip_rain,  pgflux, snowdz, temp_sol_add) 

    END IF

    ! Store the proxy for unstressed ecosystem functioning. This could be gpp, transpiration,
    ! evaporation. They all differ conceptually and numerically but the idea is the same.
    ! In an initial approach we tried to calculate the stress as a ratio of the initial 
    ! transpiration over water supplied through the plant. This approach resulted in
    ! inconsistencies likely because the many feedbacks in the calculation of transpir. For
    ! this reason it was decided to drop this approach and use the ratio between the initial 
    ! and adjusted gpp instead. This ratio will be used in stomate to adjust the allocation.
    ! Note that the feedback on gpp is calculated based on the ratio between ::transpir_supply 
    ! and transpir at the half hourly time step. In stomate we use the ratio between stressed
    ! and unstressed but at the daily time step. Night time values should be ignored in the
    ! daily calculation. If night time values are considered to unstressed this will skew
    ! the result.
    unstressed(:,:) = gpp(:,:)

    IF (ok_hydrol_arch) THEN

       CALL hydraulic_arch_calc (kjpindex, circ_class_biomass, circ_class_n, transpir, stempdiag, &
            temp_air, swc, transpir_supply, njsc, veget, &
            transpir_supply_column, ksoil, e_frac, vessel_loss, root_profile)

       DO ji = 1, kjpindex

          !don't restrict gs during night
          DO jv = 1, nvm
             IF (sinang(ji) .LT. min_sechiba) THEN
                ! Night time we will never recalculate the energy budget
                hydrol_flag2(ji) = .FALSE.
                hydrol_flag3(ji,jv) = .FALSE.
             ELSE
                ! The ratio between ::transpir_supply and ::transpir is used to determine
                ! whether we will recalculate gpp, transpir, ... (thus enerbil). This
                ! is done every half hour. The waterstress that will be passed to stomate
                ! to adjust the allocation will be calculated at the daily time step.
                ! It is important NOT to include night time values in the daily mean
                ! because more southern sides have shorter days during the growing season
                ! the number of half hours with day light biases the waterstress if it 
                ! is assumed that there is no waterstress at night.
                ! Note that if transpir_suppy = zero this is likely caused by a soil for 
                ! which all 11 layers are below wilting point. If this is the case it is
                ! correct to have transpir_supply = 0 and thus gpp = 0. This behaviour is
                ! not frequently observed. Even not in very dry regions. It suggests that
                ! the plant response to soil drying is too slow, implying transpir_supply
                ! is too high the days/week before the drought. Indeed the scheme
                ! to calculate the maximum transpir_supply is static, i.e. it always calculates 
                ! the maximum transpir_supply (and therefore stomata might close too late or too 
                ! little). 
                IF (veget_max(ji,jv) .GT. zero .AND. &
                     veget(ji,jv) .GT. zero) THEN  
                   ! comparing the transpiration in tree-stand scale    
                   IF ( transpir_supply(ji,jv) .LT. transpir(ji,jv) ) THEN 
                      IF(printlev_loc>=3) WRITE(numout,*) 'hydrol_flag2&3 set to TRUE, pft is: ', jv
                      hydrol_flag2(ji) = .TRUE.
                      hydrol_flag3(ji,jv) = .TRUE.
                   ELSE
                      IF(printlev_loc>=3) WRITE(numout,*) 'hydrol_flag2 set to FALSE, pft is: ', jv
                   END IF
                ELSE
                   ! There is no vegetation here, so we can skip it. The flag at the 
                   ! pixel level (::hydrol_flag2) should not change its value.
                   hydrol_flag3(ji,jv) = .FALSE.
                ENDIF
             END IF ! sinang(ji) .LT. min_sechiba
          END DO ! jv = 1, nvm
       END DO ! ji = 1, kjpindex

       ! The energy budget is only calculated again, if the above conditions are satisfied 
       ! for any of the grid squares, and the feedback of waterstress on stomatal closure 
       ! (ok_gs_feedback) is activated
       IF (ok_gs_feedback) THEN
          ! Gather hydrol_flag2 on the full global domain and set hydrol_flag true 
          ! if at least one grid cell contains hydrol_flag2()=true.
          CALL gather(hydrol_flag2, hydrol_flag2_glo)
          IF (is_root_prc) THEN
             hydrol_flag=.FALSE.
             DO ji = 1, nbp_glo
                IF  (hydrol_flag2_glo(ji)) THEN
                   hydrol_flag = .TRUE.
                   EXIT
                END IF
             END DO
          END IF
          ! Copy hydrol_flag to all processors
          CALL bcast(hydrol_flag)
       ELSE
          ! If the feedback of water stress on stomatal closure is not activated, 
          ! then hydrol_flag is set to false, and no energy re-calculations are conducted
          hydrol_flag=.FALSE.
       END IF

       IF (hydrol_flag) THEN
          DO ji = 1, kjpindex
             IF (hydrol_flag2(ji)) THEN
                DO jv = 1, nvm
                   IF (.NOT. hydrol_flag3(ji, jv) .OR. lai(ji,jv) .LT. 0.5 ) THEN 
                      ! this is for cases when the supply term is higher than the transpiration
                      ! i.e. the normal state of affairs. The condition for lai avoids that
                      ! we are getting crazy values when recalculating the gpp. Towards the
                      ! end of the growing season the lai is decreasing which results in a very
                      ! high R_shoot. If there is serious water stress at the time that the lai 
                      ! is very low, the recalculated gpp values go through the roof. We use
                      ! this arbitrary threshold to avoid this behaviour. This implies that small
                      ! canopies will not experience water stress in the model.
                   ELSE
                      ! To calculate the gs_effective (the effective stomatal conductance), and 
                      ! GPP first we need to calculate a constrained vbeta3, and work back from 
                      ! there, using a tweaked version of diffuco_trans_co2. We calculate the 
                      ! vbeta3 based on the ::transpir_supply term for each grid and pft. We should 
                      ! keep the spatial weighting fraction,::veget_max, for vebta3. Even if the 
                      ! calculation of ::transpir_supply is at the stand scale.   
                      IF (veget_max(ji,jv) .GT. zero .AND. veget(ji,jv) .GT. zero) THEN       
                         vbeta3(ji,jv) = transpir_supply(ji,jv)  / &  
                              &  ( dt_sechiba * (un - vbeta1(ji)) * (qsol_sat_new(ji) - qair_new(ji)) * &
                              &  ( rau(ji) * (MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji)))) * tq_cdrag(ji)))
                      ELSE
                         ! No vegetation present for this condition, so vbeta3 is set equal to zero. 
                         vbeta3(ji,jv) = zero
                      ENDIF

                      ! Calculate the ratio between the stressed and unstressed vbeta3. It will be
                      ! used to correct some extreme values.
                      ratio_vbeta3(ji,jv) = vbeta3(ji,jv)/vbeta3_save(ji,jv)

                      ! Recalculate leaf_ci, gsmean and gpp for the new vbeta3 that matches the 
                      ! water limitation
                      CALL mleb_diffuco_trans_co2_short(assimi_lev, assim_param, ccanopy, gstot_frac, ji, jv, &
                           kjpindex, Light_Abs_Tot, Light_Tran_Tot, lai, lai_per_level, pb, qair, tq_cdrag, &
                           qsintmax, qsintveg, ratio_vbeta3, rstruct, rveget, swdown, temp_air, temp_growth, &
                           u, v, vbeta23, vbeta3, veget_max, veget, info_limitphoto, JJ_out, cimean, cresist, &
                           leaf_ci, vbeta3pot, gsmean, gpp)

                      ! Clean some unexpected/extreme values (less than 1% of the cases) 
                      IF(gpp(ji,jv) .EQ. unstressed(ji,jv)) THEN
                         gpp(ji,jv) = ratio_vbeta3(ji,jv) * unstressed(ji,jv)
                      ELSE IF( gpp(ji,jv) .LE. zero) THEN
                         gpp(ji,jv) = zero
                      ELSE IF( gpp(ji,jv) .LE. 0.1*ratio_vbeta3(ji,jv)*unstressed(ji,jv)) THEN
                         gpp(ji,jv) = 0.1*ratio_vbeta3(ji,jv)*unstressed(ji,jv)
                      ELSE IF( gpp(ji,jv) .GT. 100 ) THEN
                         IF (err_act.GT.1) THEN
                            WRITE(numout,*)'GPP',gpp(ji,jv)
                            WRITE(numout,*)'vbeta3',vbeta3(ji,jv)
                            WRITE(numout,*)'vbeta3_save',vbeta3_save(ji,jv)
                            WRITE(numout,*)'transpir_supply',transpir_supply(ji,jv)
                            WRITE(numout,*)'vbeta1',vbeta1(ji)
                            WRITE(numout,*)'qsol_sat_new',qsol_sat_new(ji)
                            WRITE(numout,*)'qair_new',qair_new(ji)
                            WRITE(numout,*)'rau',rau(ji)
                            WRITE(numout,*)'tq_cdrag',tq_cdrag(ji)
                            WRITE(numout,*)'transpir',transpir(ji,jv)
                            WRITE(numout,*)'rstruct',rstruct(ji,jv)
                            WRITE(numout,*)'rveget',rveget(ji,jv)
                            WRITE(numout,*)'gsmean',gsmean(ji,jv)
                            CALL ipslerr_p(2,'hydraulic_arch_main', &
                              'The recalculation of GPP is too high(>100)!',&
                              'If this happens a lot, you should take a closer look','')
                         ENDIF
                      ENDIF
                   END IF ! (hydrol_flag3(ji, jv) .eq. .FALSE.)
                ENDDO ! jv = 1, nvm
             ELSE
                ! do nothing
             END IF ! (hydrol_flag2(ji) .eq. .TRUE.)  

             frac_evap_transp(ji,:)=frac_evap_transp(ji,:)*ratio_vbeta3(ji,jv)
             
             resist_equiv_evap(ji) = un / MAX(10.0**(-20),( (un - frac_sub_snow(ji)) * (un - frac_evap_flood(ji)) * &
                  ( frac_evap_bare(ji) /raero(ji) + SUM(frac_evap_intercept(ji,:)/resist_intercept(ji,:)) + &
                  SUM(frac_evap_transp(ji,:)/resist_transp(ji,:)) ) + frac_evap_flood(ji)/raero(ji) ))

             resist_equiv_evap_veg(ji) = un / MAX(10.0**(-20),( (un - frac_sub_snow(ji)) * (un - frac_evap_flood(ji)) * &
                  ( frac_evap_bare(ji) /raero(ji) + SUM(frac_evap_intercept(ji,:)/resist_intercept(ji,:)) + &
                  SUM(frac_evap_transp(ji,:)/resist_transp(ji,:)) ) ))

             resist_equiv_tot(ji) = un /  MAX(10.0**(-20),( frac_sub_snow(ji) * (un - frac_evap_flood(ji)) / raero(ji) + &
                  (un - frac_sub_snow(ji)) * (un - frac_evap_flood(ji)) * ( frac_evap_bare(ji) /raero(ji) + &
                  SUM(frac_evap_intercept(ji,:)/resist_intercept(ji,:)) + SUM(frac_evap_transp(ji,:)/resist_transp(ji,:)) ) + &
                  frac_evap_flood(ji)/raero(ji)))

             vbeta(ji)=vbeta4(ji)+SUM(vbeta2(ji,:))+SUM(vbeta3(ji,:))
             vbeta3sum(ji)=SUM(frac_evap_transp(ji,:)/resist_transp(ji,:))
             
          END DO ! ji = 1, kjpindex
          
          
          ! mleb_main and enerbill_main have to be called here again 
          ! (to recalculate the energy budget because of the
          ! supply term restriction.

          ! Before 2nd call to enerbil, recall initial values from before 
          ! mleb was called the first time.
          evapot(:) = evapot_save(:)
          evapot_corr(:) = evapot_corr_save(:)
          temp_sol(:) = temp_sol_save(:)
          qsurf(:) = qsurf_save(:)
          tsol_rad(:) = tsol_rad_save(:)
          pgflux(:) = pgflux_save(:)

          IF (.not. ok_mleb) THEN

             ! Second call calculates the energy budget with the water limited value of vbeta3
             CALL enerbil_main (kjit, kjpindex, netrad, lwabs, lwnet, fluxsubli, index, &
                  & indexveg, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, qair, &
                  & peqAcoef, peqBcoef, pb, rau, resist_equiv_evap, resist_equiv_evap_veg, resist_equiv_tot, &
                  & frac_sub_snow, frac_evap_intercept, resist_intercept, &
                  & frac_evap_transp, resist_transp, vbeta3pot, frac_evap_bare, frac_evap_flood, raero, &
                  & emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
                  & vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, &
                  & temp_sol_new, qsurf, evapot, evapot_corr, precip_rain, pgflux, &
                  & snowdz, temp_sol_add, qair_new, qsol_sat_new)

          ELSE

             !+++CHECK+++
             ! (1) Add the inout variables for mleb
             ! (2) Variable dimensions of xxx_Tot_mean are for a single pixel. Causing 1+1 problems
             ! when running a larger domain. Needs to be corrected when implementing
             ! a global use of the multi-layer energy budget.
             !+++++++++++
             CALL mleb_main (kjit, kjpindex, ldrestart_read, ldrestart_write, &
                  & index, indexveg, lwdown, swnet, swdown, epot_air, temp_air, &
                  & u, v, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, pb, rau, &
                  & vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
                  & emis, soilflx, soilcap, tq_cdrag, humrel, &
                  & fluxsens, fluxlat, vevapp, transpir, transpot, vevapnu, vevapwet, &
                  & vevapsno, vevapflo, temp_sol, tsol_rad, temp_sol_new, qsurf, &
                  & evapot, evapot_corr, rest_id, hist_id, hist2_id, & 
                  & ok_mleb_history_file, &
                  & transpir_supply, hydrol_flag, hydrol_flag2, &
                  & hydrol_flag3, vbeta2sum, vbeta3sum, veget_max, qsol_sat_new, qair_new, &
                  & veget, &
                  & Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
                  & laieff_isotrop, lai_per_level,z_array_out, transpir_supply_column, &
                  & u_speed, profile_vbeta3, profile_rveget, max_height_store, &
                  & precip_rain,  pgflux, snowdz, temp_sol_add)
              
          END IF

       ELSE
          ! do nothing
       END IF ! (hydrol_flag .eq. .TRUE.) 

    END IF   ! IF (ok_hydrol_arch) THEN

    ! Make sure there is always a stressed and an unstressed value to avoid numerical
    ! issues in stomate. Store the proxy for stressed ecosystem functioning. This value 
    ! is used in stomate to calculate the waterstress used in allocation. See comment
    ! where unstressed(:,:) is defined.
    stressed(:,:) = gpp(:,:)


    !! 2.3.1 Write enerbil data to the history files.
    IF (ok_mleb) THEN
       CALL mleb_write (kjit, kjpindex, index, lwdown, temp_sol, &
            temp_sol_new, evapot, evapot_corr, hist_id, hist2_id, &
            vevapp, vevapwet, transpir, vevapnu, vevapsno, vevapflo)
    ELSE
       CALL enerbil_write (kjit, kjpindex, hist_id, hist2_id, index, &
            & evapot, evapot_corr, fluxsubli, lwdown, lwnet, netrad, transpir, &
            & vevapflo, vevapnu, vevapp, vevapsno, vevapwet)
    ENDIF

END SUBROUTINE hydraulic_arch_main


!! ================================================================================================================================
!! SUBROUTINE 	: hydraulic_arch_calc
!!
!>\BRIEF        Calculate the amount of water that is available for the leaves for transpiration,
!               accounting for effects of hydraulic architecture following (Hickler et al. 2006)    
!!
!!\n DESCRIPTION : (1) Calculate soil water potential from soil water content as calculated in hydrol.f90 
!!   according to van Genuchten (1980). Soil water potential is weighted by the relative share of the
!!   roots in each layer.
!!   (2) Calculate pressure difference between leaves and soil (Whitehead 1998)
!!   (3) Calculate root (Weatherly 1998), sapwood (Magnani et al. 2000, Tyree et al. 1994, Sperry et al. 1998)
!!   and leave resistances (Sitch et al. 2003).
!!   (4) Calculate the supply of water for transpiration by the leaves using Darcy's law (Whitehead 1998)
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): :: transpir_supply   
!!
!! REFERENCE(S)	: +++COMPLETE+++
!!                Hickler et al. 2006
!!                van Genuchten (1980)
!!                Whitehead 1998
!!                Weatherly 1998
!!                Magnani et al. 2000
!!                Tyree et al. 1994
!!                Sperry et al. 1998
!!                Sitch et al. 2003
!!
!! FLOWCHART    : 
!! 
!_ ================================================================================================================================

  SUBROUTINE hydraulic_arch_calc (kjpindex, circ_class_biomass, circ_class_n, transpir, stempdiag, &
       temp_air, swc, transpir_supply, njsc, veget, &
       transpir_supply_column, ksoil, e_frac, vessel_loss, root_profile) 


    !! 0. Variable declaration

    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                      :: kjpindex            !! Domain size - terrestrial pixels only (unitless)
    INTEGER(i_std), DIMENSION(:), INTENT (in)       :: njsc                !! Index of the dominant soil textural class 
                                                                           !! in the grid cell used for hydrology 
                                                                           !! (1-nscm, unitless)
    REAL(r_std),DIMENSION(:,:),INTENT(in)           :: stempdiag           !! Soil temperature (K)
    REAL(r_std),DIMENSION (:), INTENT (in)          :: temp_air            !! Air temperature (K)   
    REAL(r_std),DIMENSION (:,:,:), INTENT (in)      :: swc                 !! Volumetric soil water content
    REAL(r_std),DIMENSION (:,:,:), INTENT (in)      :: circ_class_n        !! number of trees within a circumference class (tree m-2)
    REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in)   :: circ_class_biomass  !! Biomass components of the model tree  
                                                                           !! within a circumference class
                                                                           !! class @tex $(g C ind^{-1})$ @endtex  
    REAL(r_std),DIMENSION (:,:), INTENT (in)        :: veget               !! Fraction of the vegetation that covers the soil
                                                                           !! also project canopy cover - calculated by Pgap
    REAL(r_std),DIMENSION (:,:,:), INTENT (in)      :: ksoil               !! Hydraulic conductivity mm/d 
    REAL(r_std), DIMENSION(:,:), INTENT (in)        :: transpir            !! Unstressed transpiration mm/dt 
    REAL(r_std),DIMENSION (:,:,:), INTENT(in)       :: root_profile        !! Normalized root mass/length fraction in each soil layer 
                                                                           !! (0-1, unitless)

    !! 0.2 Output variables
    REAL(r_std), DIMENSION(:,:), INTENT (out)       :: transpir_supply     !! Supply of water for transpiration. Same unit as transpir. 
                                                                           !! @tex $(mm m^{-2} leaf dt^{-1})$ @endtex
    REAL(r_std), DIMENSION(:,:,:,:), INTENT(out)    :: e_frac              !! Fraction of water transpired supplied by individual layers (no units),
                                                                           !! dimension(kjpindex,nvm,nslm,nstm)
    REAL(r_std), DIMENSION(:,:), INTENT(out)        :: vessel_loss         !! Conductivity lost due to cavitation in the xylem (no unit).

    !! 0.3 Modified variables


    !! 0.4 Local variables
    INTEGER(i_std)                                  :: ier                 !! Error handling
    INTEGER(i_std)                                  :: ipts,ivm,icirc      !! indices
    INTEGER(i_std)                                  :: ibdl,istm,jst       !! indices
    INTEGER(i_std)                                  :: tmp_ncirc           !! indices
    REAL(r_std)                                     :: R_root              !! Hydraulic resistance of the roots (MPa s m-3)
    REAL(r_std)                                     :: R_leaf              !! Hydraulic resistance of the leaves (MPa s m-3)
    REAL(r_std)                                     :: R_sap               !! Hydraulic resistance of sapwood (MPa s m-3)
    REAL(r_std)                                     :: R_shoot             !! Hydraulic resistance of shoot corrected for temperature (??)
    REAL(r_std), DIMENSION(ncirc)                   :: height              !! Vegetation height per circ_class_biomass (m)
    REAL(r_std), DIMENSION(kjpindex,nslm,nstm)      :: psi_soil            !! Soil water potential per soil layer per soiltile higher 
                                                                           !! than -5 MPa (threshold for hygroscopic water) (MPa)
    REAL(r_std), DIMENSION(ncirc)                   :: circ_class_transpir_supply  !! Supply of water for transpiration per model tree 
                                                                                   !! in each circumference class @tex $(m s^{-1}tree^{-1})$ @endtex
    REAL(r_std)                                     :: delta_P             !! Pressure difference between leaves and soil
    REAL(r_std), DIMENSION(nslm+1)                  :: z_soil              !! Variable to store depth of the different
                                                                           !! soil layers (m)
    REAL(r_std), DIMENSION(nvm)                     :: k_sap_cav           !! sapwood conductivity decreased by cavitation (m2 s-1MPa)
    REAL(r_std), DIMENSION(nvm)                     :: rpc                 !! scaling factor for integral (unitless)
    REAL(r_std), DIMENSION(kjpindex,nslm,nstm)      :: psi_soiltile        !! soil water potentential per soil layer per soil tile (MPa) 
    REAL(r_std), DIMENSION (kjpindex,nslm,nstm)     :: swc_adj             !! Volumetric soil water content, adjusted when it is lower 
                                                                           !! than or equal to the residual soil water content
    REAL(r_std), DIMENSION(kjpindex,nvm)            :: psi_soilroot        !! Sum of soil water potential in each soil layer weighted 
                                                                           !! by the share of the roots in each layer (MPa)
    INTEGER(i_std)                                  :: i                   !! counter
    REAL(r_std), DIMENSION(nlevels_tot, ncirc)      :: cent_box            !! Height of the mid-point of each box in the column that is made up
                                                                           !!    of the Pgap points
    REAL(r_std), DIMENSION(nlevels_tot)             :: delta_P_column      !! Pressure difference between leaves and soil
    REAL(r_std), DIMENSION(nlevels_tot,ncirc)       :: circ_class_transpir_supply_column !! Supply of water for transpiration per model tree 
                                                                                         !! in each circumference class @tex $(m s^{-1}tree^{-1})$ @endtex 

    REAL(r_std), DIMENSION(nlevels_tot,kjpindex,nvm):: transpir_supply_column !! Supply of water for transpiration 
                                                                              !! @tex $(mm m{-2} ground dt^{-1})$ @endtex
    REAL(r_std), DIMENSION(kjpindex,nslm,nstm)      :: e_max                  !! Maximum amount of available water for transpiration 
                                                                              !! per layer (mmol m2 s-1)
    REAL(r_std), DIMENSION(kjpindex,nslm,nstm)      :: res_soil            !! Soil to root resistance  per layer (s m-2)
    REAL(r_std), DIMENSION(kjpindex,nslm,nstm)      :: ksoilUNIT           !! Hydraulic conductivity (from hydrol.90)

    REAL(r_std), DIMENSION(nslm)                    :: lr                  !! Root length (m)
    REAL(r_std), DIMENSION(nslm)                    :: rs                  !! Half the mean distance between roots (m) 

    REAL(r_std),DIMENSION(kjpindex,nvm,ncirc)           :: R_root_circ     !! Hydraulic resistance of the roots corrected for temperature 
                                                                           !! per circumference class (MPa s m-3))
    REAL(r_std),DIMENSION(kjpindex,nvm,ncirc)           :: R_shoot_circ    !! Hydraulic resistance of the shoot corrected for temperature 
                                                                           !! per circumference class (MPa s m-3)
    REAL(r_std),DIMENSION(kjpindex,nvm)                 :: R_root_tot      !! Hydraulic resistance of the roots per pft
    REAL(r_std),DIMENSION(kjpindex,nvm)                 :: R_shoot_tot     !! Hydraulic resistance of the shoot per pft
    REAL(r_std), DIMENSION(nslm)                        :: nroot_tmp       !! Temporary variable to calculate the nroot
    REAL(r_std), DIMENSION(kjpindex,nslm,nvm)           :: tmp             !! Temporary variable to write the root profile to xios

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

    !! Initialize local printlev
    IF (l_first_sechiba) THEN
       printlev_loc=get_printlev('hydrolarch')
    ENDIF

    ! initialise some variables for the transpir_supply_column loop
    i = 0
    cent_box = zero
    delta_P_column = zero
    circ_class_transpir_supply_column = zero
    transpir_supply_column = zero

    ! Save the hydraulic resitances
    R_root_circ(:,:,:) = zero
    R_shoot_circ(:,:,:) = zero
    R_root_tot(:,:) = zero
    R_shoot_tot(:,:) = zero

    ! To avoid that the plants will not experience water stress in the first 
    ! day of the simulation. Transpir_supply is calculated before 
    ! there is any biomass available.
    transpir_supply(:,:) = val_exp
    psi_soiltile(:,:,:) = -un
    e_frac(:,:,:,:) = zero

    ! Initialize vessel_loss so it always has a value.
    IF (hack_vessel_loss.LT.zero .OR. &
         hack_vessel_loss.GT.un) THEN
       ! Set to zero. This is what we want for real simulations
       vessel_loss(:,:) = zero
    ELSE
       ! Use a constant fix value instead. This option
       ! was added for debugging and checking the sanity
       ! of the code
       vessel_loss(:,:) = hack_vessel_loss
    ENDIF

    ! Distinguish between the two soil map possibilities for these variables
    ! First we need to make sure there are allocated.
    
    IF (.NOT. ALLOCATED(nvan)) THEN
       ALLOCATE (nvan(nscm),stat=ier)
       IF (ier /= 0) CALL ipslerr_p(3,'hydraulic_arch_calc','Pb in allocate of variable nvan','','')
    ENDIF
   
    IF (.NOT. ALLOCATED(avan)) THEN
       ALLOCATE (avan(nscm),stat=ier)
       IF (ier /= 0) CALL ipslerr_p(3,'hydraulic_arch_calc','Pb in allocate of variable avan','','')
    ENDIF
       
    IF (.NOT. ALLOCATED(mcr)) THEN
       ALLOCATE (mcr(nscm),stat=ier)
       IF (ier /= 0) CALL ipslerr_p(3,'hydraulic_arch_calc','Pb in allocate of variable mcr','','')
    ENDIF
       
    IF (.NOT. ALLOCATED(mcs)) THEN
       ALLOCATE (mcs(nscm),stat=ier)
       IF (ier /= 0) CALL ipslerr_p(3,'hydraulic_arch_calc','Pb in allocate of variable mcs','','')
    ENDIF


    nvan(:) = nvan_usda(:)
    avan(:) = avan_usda(:)
    mcr(:) = mcr_usda(:)
    mcs(:) = mcs_usda(:)

    !! 1. Calculate plant soil water potential

    ! Soil water potential for each soiltile in each layer from 
    ! soil water content according to Van Genuchten (1980)

    DO ipts = 1,kjpindex  ! (loopref: ha1)

       ! Debug
       IF (printlev_loc>=4) THEN
          WRITE(numout,*) 'soil texture',njsc
          WRITE(numout,*) 'swc', swc(ipts,:,:)  
       ENDIF
       !-

       !! 2.1 Calculate soil water potential
       !  Swc calculated in hydrol.f90 cannot be lower than or equal 
       !  to the residual soil water content for each soiltile.
       DO ibdl = 1, nslm  ! (loopref: ha2)

          DO istm = 1,nstm ! (loopref: ha3)

             IF (swc(ipts,ibdl,istm) .LE. mcr(njsc(ipts))) THEN

                !+++CHECK+++
                !Explain why you do this
                swc_adj(ipts,ibdl,istm) = mcr(njsc(ipts)) + .00001
                !+++++++++++

             ELSE

                swc_adj(ipts,ibdl,istm) = swc(ipts,ibdl,istm)

             ENDIF

          ENDDO  ! istm = 1,nstm (loopref: ha3)

          ! Soil water potential for each soiltile in each layer from 
          ! soil water content according to Van Genuchten (1980). Convert 
          ! unit avan from mm-1 to m-1. Convert units soil matric potential
          ! from m to MPa by multiplying with cte_grav and rho_h2o. 
          ! Soiltile = fraction of each soil column 
          ! (three soil columns=bare soil, tree, grass)
          DO istm = 1,nstm   ! (loopref: ha4)

             IF (swc_adj(ipts,ibdl,istm) .GT. min_stomate .AND. &
                  (swc_adj(ipts,ibdl,istm) - mcr(njsc(ipts))) .GT. min_stomate) THEN
                
               
                psi_soiltile(ipts,ibdl,istm) = &
                     - ( cte_grav * rho_h2o * kilo_to_unit * & 
                     ( 1/(avan(njsc(ipts))*kilo_to_unit) ) * &
                     ( ( (swc_adj(ipts,ibdl,istm) - mcr(njsc(ipts)) ) / &
                     ( mcs(njsc(ipts)) - mcr(njsc(ipts)) ) ) ** &
                     ( -1 / (1-1/nvan(njsc(ipts))) ) - 1) ** &
                     (1/nvan(njsc(ipts)) ) ) / mega_to_unit

             ELSE

                psi_soiltile(ipts,ibdl,istm) = zero

             ENDIF

          ENDDO   ! DO istm = 1,nstm (loopref: ha4)

          ! Threshold for soil water potential: it cannot be lower than the 
          ! soil water potential for hygroscopic water (= -5 MPa, Larcher 2001)
          psi_soil(ipts,ibdl,:) = MAX(psi_soiltile(ipts,ibdl,:),-5.)

         
       ENDDO  ! DO ibdl = 1, nslm  ! (loopref: ha2)

       ! following initialisation of swc
       l_first_sechiba = .FALSE.

       ! Debug
       IF (printlev_loc>=4) THEN
          WRITE(numout,*) 'psi_soiltile', psi_soiltile(ipts,:,:)
          WRITE(numout,*) 'psi_soil', psi_soil
       ENDIF
       !-

       !! 3. Calculate soil to root and plant resistances
   
       !! 3.1 Calculate root profile

       ! This is now being done in hydrol.f90. Given hydraulic_arch is called
       ! before hydrol.f90, the root_profile runs one time step behind the
       ! soil moisture profile. This does not appear as an important issue.
       ! In case the static root profile is used the root profile is fixed and
       ! there is no issue at all. In case a dynamic root profile is used,
       ! this 30 minute delay would occur. See hydrol_dynamic_root_profile
       ! for some details/thinking on the issue.

       !! 3.2 Calculate water transport from the soil-root zone up to the leaves
       ! Initialize psi_soilroot to avoid problems with sending unitialized values
       ! to xios
       psi_soilroot(ipts,1) = xios_default_val

       DO ivm = 2,nvm ! (loopref: ha6)

          !! Calculate PFT-specific characteristics
          IF (is_tree(ivm)) THEN

             ! Calculate for each circumference classes
             tmp_ncirc = ncirc

             ! Calculate the effective height (m) from biomass. The effective
             ! value is used here because the water needs to be transported
             ! through both the belowground and aboveground biomass. Using the 
             ! height (thus not the effective height) would only account for the 
             ! height due to the aboveground biomass.
             height(:) = wood_to_height_eff(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm,pipe_tune2(ipts,ivm))

          ELSE ! is_tree(ivm)

             ! There are no circumference classes for grasses
             ! and crops. Biomasses are stored in first circumference
             ! class
             tmp_ncirc = 1

             ! Calculate height (m) from biomass. Use the function 
             ! biomass_to_lai to account for dynamic sla
             height(1) = biomass_to_lai(&
                  circ_class_biomass(ipts,ivm,1,ileaf,icarbon),ivm) * lai_to_height(ivm)
          ENDIF

          !! 3.2.1 Calculate soil to root resistance
          !  Calculate soil water potential weighted for root density
          !  Note that the second soiltile in swc is for trees.
          !  This implementation relies on a strong assumptions. First,
          !  roots are assumed to follow a rigid distribution irrespective 
          !  of how the water is distributed in the soil. So even if on 
          !  average lots of water is available in deeper layer, the modeled 
          !  plants will grow most roots in the top layers. A detailed description  
          !  allows for an adjustment of the soil water potential (psi_soilroot)  
          !  due to soil to root resistance. Psi_soil is weighted by the fraction of 
          !  evapotranspiration supplied by each layer (E_frac), which depends on 
          !  the soil to root resistance per layer (res_soil) accounting for 
          !  different root properties and hydraulic conductivity. 

          !Select soil tile for this pft
          jst = pref_soil_veg(ivm) 

          ! We can only make these calculations if there is actual root
          ! biomass.
          IF(SUM(circ_class_biomass(ipts,ivm,:,iroot,icarbon)) .GE. min_stomate ) THEN

             DO ibdl = 1, nslm

                ! Calculate root length and one-half the mean distance between roots (rs)
                ! following Bonan et al. (2014) GMD. eq A23
                
                ! Calculate root length (m) using parameter specific root length (srl). 
                ! Convert gC/ind to g/ind (deux).
                lr(ibdl)=root_profile(ipts,ivm,ibdl)*deux* &
                     SUM(circ_class_biomass(ipts,ivm,:,iroot,icarbon)*&
                     circ_class_n(ipts,ivm,:)) * srl(ivm)

                ! Grasses and crops have only roots to one meter deep
                ! hence lr(ibdl) will be zero fo the deeper soil layers
                IF (lr(ibdl).GT.min_stomate) THEN
                   rs(ibdl) = (un/(pi*lr(ibdl)))**(undemi)
                ELSE
                   rs(ibdl) = zero
                END IF

                ! Convert units for the hydraulic conductivity from mm/d to m/s 
                ksoilUNIT(ipts,ibdl,jst)=ksoil(ipts,ibdl,jst)/(mm_m * one_day)

                ! The soil to root resistance in each layer (res_soil) is calculated before
                ! the contribution from each layer to the evapotransporation e_frac. 
                ! froot gives the mean radius of fine roots (m).
                ! Unit of res_soil is s m-2. Further
                ! unit conversion is not needed because e_max is normalized.
                IF (rs(ibdl).GT.min_stomate) THEN
                   res_soil(ipts,ibdl,jst)=(LOG(rs(ibdl)/r_froot(ivm)))/&
                        (deux*pi*lr(ibdl)*ksoilUNIT(ipts,ibdl,jst))
                ELSE
                   ! If there are no roots in a specific layer the resistance is
                   ! infinite.
                   res_soil(ipts,ibdl,jst)=large_value
                END IF
               
                ! Psi_root is a poorly measured parameter. Since psi_leaf_min is 
                ! frequently measured and a psi_root that is lower than psi_leaf 
                ! will not result in more more water transport to the leaves anyway, we 
                ! prefer to set psi_root equal to psi_leaf rather than choosing an 
                ! arbitrary value that is unrelated to to the prescribed psi_leaf.
                psi_root(ivm)=psi_leaf_min(ivm)
                e_max(ipts,ibdl,jst)= (psi_soil(ipts,ibdl,jst)-psi_root(ivm))/&
                     res_soil(ipts,ibdl,jst)

                ! If we had to divide by an infinite resistance (where infinite was
                ! represented by large_value (=10e+33), e_max should less than 
                ! min_sechiba and will be set to zero.
                IF (e_max(ipts,ibdl,jst) .LT. min_sechiba) e_max(ipts,ibdl,jst) = zero

             END DO

             ! If there is a contribution to evapotransporation, the fraction
             ! coming from each layer is calculated first, followed by calculations 
             ! of the weighted soilroot water potential (psi_soilroot). When there is no 
             ! contribution, psi_soilroot is set very low.  
             IF( SUM(e_max(ipts,:,jst)) .GT. min_sechiba) THEN
                
                !+++CHECK+++
                ! There is no need to have the nst dimension in e_frac because each PFT
                ! is linked to just a single soil tile. Could be simplified.
                DO ibdl = 1, nslm
                   e_frac(ipts,ivm,ibdl,jst)=e_max(ipts,ibdl,jst)/SUM(e_max(ipts,:,jst))
                ENDDO
                !+++++++++++

                IF (hack_e_frac) THEN
                   ! Simplify the calculation of psi_soilroot by not accounting
                   ! for root length. Overwrite e_frac with root_profile.
                   e_frac(ipts,ivm,:,jst) = root_profile(ipts,ivm,:)
                END IF

                ! Calculate psi_soilroot
                psi_soilroot(ipts,ivm) = SUM( psi_soil(ipts,:,jst)*e_frac(ipts,ivm,:,jst))
                
             ELSE
                
                ! The calculation of psi_soil is not realistic because psi_root > psi_soil
                ! Use a suboptimal solution by prescribing e_frac and psi_soilroot
                e_frac(ipts,ivm,:,jst)= un/nslm
                psi_soilroot(ipts,ivm)=-4
                
             END IF

             ! Debug  
             IF(printlev_loc .GE. 4 .AND. ivm .eq. test_pft) THEN
                WRITE(numout,*)'hydraulic_arch_calc in biomass statement, ivm, jst',ivm, jst
                WRITE(numout,*)'psi_soil', psi_soil(ipts,:,jst)
                WRITE(numout,*)'lr',lr(:)
                WRITE(numout,*)'rs', rs(:)
                WRITE(numout,*)'ksoilUNIT', ksoilUNIT(ipts,:,jst)
                WRITE(numout,*)'res_soil',res_soil(ipts,:,jst)
                WRITE(numout,*)'e_max',e_max(ipts,:,jst)
                WRITE(numout,*)'e_frac', e_frac(ipts,ivm,:,jst)
                WRITE(numout,*)'psi_soilroot',psi_soilroot(ipts,ivm)
             ENDIF

          ELSE ! IF(biomass .gt. min_stomate) THEN

             ! Only do the calculations for biomass  
             e_frac(ipts,ivm,:,jst)= zero
             psi_soilroot(ipts,ivm) = SUM(psi_soil(ipts,:,jst)*e_frac(ipts,ivm,:,jst))

          ENDIF ! (biomass .GT. min_stomate) THEN  
   

          !! 3.2.2 Calculate plant resistances
          !  Calculate plant resistance to water transport according to the principle 
          !  of an electric circuit. The code is largely based on Hickler et al. 2006.
          
          !  Set transpir_supply to zero whem transpir is zero (e.g. at night),
          !  so their daily means are comparable.  
          IF ( SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon)) .GT. min_sechiba .AND. &
               SUM(circ_class_biomass(ipts,ivm,:,iroot,icarbon)) .GT. min_sechiba .AND. &
               transpir(ipts,ivm) .GT. min_sechiba ) THEN
             
             DO icirc = 1,tmp_ncirc  ! (loopref ha7)  

                ! This ensures that the loop goes over all circ_classes for trees and only
                ! the first circ_class for grasses and crops
                IF (circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon) .GT. min_sechiba .AND. &
                     circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) .GT. min_sechiba ) THEN 

                   !! 3.2.2.1 Root resistance
                   !  Root resistance is in the first order related to root mass 
                   !  and expression was proposed by Weatherly, 1982). Convert
                   !  units of circ_class_biomss i.e. gC/ind to kg/ind
                   !  Unit of R_root is (MPa s m-3 ind-1)
                   R_root = un / (k_root(ivm) * 2 * &
                        circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) &
                        / kilo_to_unit)

                   !! 3.2.2.2 Sapwood resistance
                   !  Calculate k_sap_cav as a function of psi_soilroot and the pft specific parameter
                   !  psi_50 which representing its vulnerability to cavitation, following 
                   !  Tyree et al. (1994) and Sperry et al. (1998).
                   IF (is_tree(ivm)) THEN
                      
                      !! !+++++CHECK: psi_50 should be lower than 0 to have cavitation, so it seems we never
                      !! reduce k_sap due to cavitation for the moment!
                      IF( (psi_soilroot(ipts,ivm) == zero) .AND. (psi_50(ivm) .LE. zero) )THEN

                         ! This seems to cause a strange numerical issue, even though we should
                         ! just be taking the exponential of zero (i.e., one).  
                         k_sap_cav(ivm) = k_sap(ivm)

                      ELSE
                         !++++ TEMP +++++
                         IF((psi_soilroot(ipts,ivm)/psi_50(ivm)) .LT. zero)THEN
                            !++++++ CHECK +++++++++
                            ! If this is true, we generally try to take this ratio to
                            ! an odd power, which results in an imaginary number.  So
                            ! we use this simple solution for the moment.
                         !   WRITE(numout,*) 'psi_soilroot', psi_soilroot(ipts,ivm)
                         !   WRITE(numout,*) 'ipts', ipts
                         !   WRITE(numout,*) 'ivm', ivm
                         !   CALL ipslerr_p (3,'hydraulic_arch_calc', &
                         !        'this quantity should never be negative', &
                         !        'See the output file for more details.','')
                            ! A temporary fix.  We cannot have the stop statement anymore,
                            ! since capping psi_soilroot at 0.0 above results in many PFTs
                            ! dying.  This needs to be addressed.
                            k_sap_cav(ivm) = k_sap(ivm)
                         ELSE
                            k_sap_cav(ivm) = k_sap(ivm) * &
                                 EXP(-((psi_soilroot(ipts,ivm)/psi_50(ivm)) ** &
                                 c_cavitation(ivm)))
                         ENDIF

                      ENDIF

                      ! If k_sap_cav is too low, we can get a divide by zero
                      ! error in the next couple lines or vessel_loss could exceed
                      ! one. So this is a protection against that.
                      IF(k_sap_cav(ivm) .LE. min_sechiba)THEN

                         k_sap_cav(ivm) = min_sechiba

                      ENDIF

                      ! Caculation of the PLC (percentage of lost conductivity
                      ! in the xylem).
                      IF (ok_vessel_mortality .AND. &
                           (hack_vessel_loss.LT.zero .OR. &
                           hack_vessel_loss.GT.un)) THEN

                         ! The ratio between k_sap_cav and k_sap
                         ! gives the fraction  of functional vessels remaining. The
                         ! PLC will be used in module stomate_turnover.f90 to
                         ! compute mortality.
                         vessel_loss(ipts,ivm) = un - (k_sap_cav(ivm) / k_sap(ivm))
                         
                      END IF

                      !  Following the expression by Magnani et al 2000, sapwood resistance
                      !  is given as a function of sapwood area. Convert Sapwood mass to 
                      !  sapwood area.
                      !  Units of R_sap is (MPa s m-3 ind-1).
                      R_sap = (height(icirc)**2 * tree_ff(ivm)*pipe_density(ivm)) / (k_sap_cav(ivm) * &
                           ( circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon) + &
                           circ_class_biomass(ipts,ivm,icirc,isapbelow,icarbon)))

                   ELSE  ! IF (is_tree(ivm)) THEN

                      !+++CHECK+++
                      ! We have about the same leaf and root restistances and masses
                      ! as for trees but lack a resistance in the sapwood. This is resulting
                      ! in a huge supply because we lack part of the resistance. This was 
                      ! the motivation to create an R_sap for grasses and roots. This equation 
                      ! is made up and should be carefully checked against physiological
                      ! studies for grasses and crops on the issue.
                      R_sap = (height(icirc)**2 * pipe_density(ivm)) / (k_sap(ivm) *  &
                           circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon))
                      !+++++++++++

                   ENDIF  ! IF (is_tree(ivm)) THEN 

                   !! 3.2.2.3 Leaf resistance
                   !  Leaf resistance related to foliar projective area (Hickler 2006) or 
                   !  LAI (more consistent with pipe model)? First fpc= 1-e**(-0.5*LAI) was used 
                   !  instead of lai because it is assumed that only sunlit leaves contribute to transpiration.
                   !  Then it was replaced by veget which in itself is calculated in Pgap and thus accounts for 
                   !  LAI and canopy structure. This was still inconsistent with the calculation of the
                   !  atmospheric demand for transpiration. To improve consistency the coupled
                   !  lai is now used as done in the energy budget calculations.
                   !  R_leaf = 1/(k_leaf * fpc)
!!$                   R_leaf = un / (k_leaf(ivm) * &
!!$                        (un - EXP(-0.5 * circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon) * &
!!$                        sla(ivm))))
                   ! Units of R_leaf is (MPa s m-3), because biomass_to_couple_lai
                   ! takes the biomass of the leaves at the tree level, and
                   ! hence we get the leaf area per individual, which are in m2.
                   R_leaf = un / (k_leaf(ivm) * &
                        biomass_to_coupled_lai(circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon), &
                        veget(ipts,ivm),ivm))

                   !! 3.2.2.4 Temperature dependance of resistance
                   !  Decrease root resistance with soil temperature as a result of 
                   !  changes in water viscosity (Cochard et al 2000). grnd_80 is the 
                   !  index of the soil layer that is 80 cm deep. This 80 cm is arbitrary 
                   !  and was used because we already used the same variable in wind throw
                   !  so the variable was already calculated.
                   !  Decrease shoot resistance with air temperature (Cochard et al. 2000)
                   R_root_circ(ipts,ivm,icirc) = R_root / (a_viscosity(1) + a_viscosity(2) * &
                        (stempdiag(ipts,grnd_80) - zeroCelsius))               
                   R_shoot_circ(ipts,ivm,icirc) = (R_leaf + R_sap) / (a_viscosity(1) + a_viscosity(2) * &
                        (temp_air(ipts) - zeroCelsius))

                     ! (JR 150115) in this version we try to calculate delta_P against height
                     ! calculate the 'mid-point' of the P_gap points

                   IF (.NOT. ok_impose_can_structure ) THEN
                    !DO i = 1, nlevels_tot - 1
                    !    cent_box(i, icirc) = z_array_out3(1,ivm,1,i) + (z_array_out3(1,ivm,1,i+1) - z_array_out3(1,ivm,1,i)) / 2.
                    !END DO
                    !cent_box(nlevels_tot, icirc) = z_array_out3(1,ivm,1,nlevels_tot) + &
                    !          &         ( h_array_out3(1,ivm,1, nlevels_tot) - z_array_out3(1,ivm,1,nlevels_tot) ) / 2.0                     
                    ! DO i = 1, nlevels_tot
                    !     WRITE(numout, *) 'cent_box(i, icirc) is: ', cent_box(i, icirc)
                    ! END DO
                    !   DO i = 1, nlevels_tot
                    !       delta_P_column(i) = psi_leaf(ivm) - psi_soilroot(ipts,ivm) + &
                    !            (cent_box(i, icirc)*cte_grav*rho_h2o) / kilo_to_unit                     
                    !   END DO
                   END IF ! (ok_impose_can_structure ) THEN 

                   !! 3.2.3 Calculate pressure difference between leaves and soil 
                   !  Defined by Whitehead 1998. Convert units to MPa
                   ! Transpir_supply is to a large extend determined by delta_P. DeltaP 
                   ! in turn is sensitive to two prescribed parameters: psi_leaf_min 
                   ! and psi_root. Psi_leaf_min is rather well documented but using it 
                   ! implies that we always calculate the maximum transpir_supply which 
                   ! results in drying the soil too quickly in dry regions 
                   ! (CHECK:still true in cleaned version?). Psi_root can 
                   ! be tuned (currently set to psi_leaf, see later) because it is poorly observed. 
                   ! Some of the expected dynamics might not be captured, for example, when 
                   ! a drought is emerging, transpiration increases and then suddenly collapses 
                   ! (CHECK if this behaviour is still occurs). 
                   ! In reality it is expected that the soil would dry and that a drier soil would reduce 
                   ! the water the plant can transport. Owing to the static character of 
                   ! this scheme, it looks like its limits are reached. Further improving 
                   ! water stress would require a more dynamic scheme that calculates 
                   ! psi_leaf_min and psi_root as a function of the plant and soil water status and respectively. 
                   delta_P = psi_leaf_min(ivm) - psi_soilroot(ipts,ivm) + &
                        (height(icirc)*cte_grav*rho_h2o) / kilo_to_unit
           

                   !! 3.2.4 Calculate water supply for transpiration. 
                   !  First for a model tree in each circ_class (m/s/tree), then
                   !  for each pft.

                   !  Using Dracy's law (Slatyer 1967, Whitehead 1998)
                   circ_class_transpir_supply(icirc) = (-delta_P) / &
                        (R_root_circ(ipts,ivm,icirc) + R_shoot_circ(ipts,ivm,icirc))
                 

                   IF (.NOT. ok_impose_can_structure ) THEN 
                       DO i = 1, nlevels_tot
                           circ_class_transpir_supply_column(i, icirc) = (-delta_P_column(i)) / &
                                (R_root_circ(ipts,ivm,icirc) + R_shoot_circ(ipts,ivm,icirc))
                       END DO
                   END IF ! (ok_impose_can_structure ) THEN 

                ELSE  
                      
                   !set transpir_supply to zero if there is no leaf or root biomass in the circ_class 
                   circ_class_transpir_supply(icirc) = zero

             ENDIF

             !---TEMP--- 
             IF (printlev_loc>=4 .AND. ivm .EQ. test_pft) THEN
                WRITE(numout,*) 'root_profile', root_profile(ipts,ivm,:)
                WRITE(numout,*) 'k_sap_cavitation, biomass(iroot)', k_sap_cav(ivm), &
                     2 * circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) / 1000.
                WRITE(numout,*) 'sapwood area, ', ( circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon) + &
                     circ_class_biomass(ipts,ivm,icirc,isapbelow,icarbon)) / (tree_ff(ivm) * pipe_density(ivm)) / &
                     height(icirc)
                WRITE(numout,*) 'fpc, ',(1-EXP(-0.5*circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon)*sla(ivm)))
                WRITE(numout,*) 'R_root, R_sap, R_leaf', R_root_circ(ipts,ivm,icirc), R_sap, R_leaf
                WRITE(numout,*) 'R_shoot', R_shoot_circ(ipts,ivm,icirc) 
                WRITE(numout,*) 'delta_P', delta_P
             ENDIF
                
          ENDDO ! DO icirc = 1,tmp_ncirc (loopref ha7) 


          ! Calculating transpiration supply at different heights here:
          IF (.NOT. ok_impose_can_structure ) THEN 

              DO i = 1, nlevels_tot
                  transpir_supply_column(i, ipts, ivm) = &
                      (SUM(circ_class_transpir_supply_column(i, :) * circ_class_n(ipts,ivm,:))) * &
                      dt_sechiba * kilo_to_unit
                
              END DO ! i = 1, nlevels_tot



          ! Set transpir_supply to zero when negative.
          ! (When psisoilroot becomes very low, delta_P becomes positive and transpir_supply negative)
            DO i = 1, nlevels_tot
               IF (transpir_supply_column(i, ipts, ivm) .LT. zero) THEN
                   WRITE(numout,*) 'transpir_supply before set to zero when negative', transpir_supply(ipts,ivm), ivm
                   WRITE(numout,*) 'psisoilroot', psi_soilroot(ipts,ivm)
                   transpir_supply_column(i, ipts, ivm) = MAX(transpir_supply_column(i, ipts, ivm), zero)
               ELSE
                   !do nothing
               ENDIF   
            END DO ! i = 1, nlevels_tot

          END IF ! (ok_impose_can_structure ) THEN 


          ! Calculate transpir_supply per pft. Unit conversion, convert m/s/tree to mm/m^2/half_hour.
          ! Multiply with veget to have the same unit as transpir in sechiba, 
          ! transpir is calculated by leaf area not by land area.
          IF (is_tree(ivm)) THEN
             ! transpir_supply per PFT (sum of diameter classes) per timestep
             transpir_supply(ipts,ivm) = &
                  (SUM(circ_class_transpir_supply(:) * circ_class_n(ipts,ivm,:))) * &
                  dt_sechiba * kilo_to_unit * veget(ipts,ivm)
          ELSE
             transpir_supply(ipts,ivm) = &
                  (circ_class_transpir_supply(1) * circ_class_n(ipts,ivm,1)) * &
                  dt_sechiba * kilo_to_unit* veget(ipts,ivm)
          ENDIF

          !+++++++++++
          ! calculate the hydraulic resistances per PFT
          R_root_tot(ipts,ivm) = (SUM(R_root_circ(ipts,ivm,:)*circ_class_n(ipts,ivm,:))) 
          r_shoot_tot(ipts,ivm) = (SUM(R_shoot_circ(ipts,ivm,:)*circ_class_n(ipts,ivm,:)))

          ! Set transpir_supply to zero when negative.
          ! (When psisoilroot becomes very low, delta_P becomes positive and transpir_supply negative)
          IF (transpir_supply(ipts,ivm) .LT. zero) then
             IF(printlev_loc.GT.4 .and. ivm .eq. test_pft)THEN
                WRITE(numout,*) 'transpir_supply before set to zero when negative', transpir_supply(ipts,ivm), ivm
                WRITE(numout,*) 'ipts, ivm:', ipts,ivm
                WRITE(numout,*) 'transpir_supply: ', transpir_supply(ipts,ivm)
                WRITE(numout,*) 'psisoilroot: ', psi_soilroot(ipts,ivm)
             ENDIF
             transpir_supply(ipts,ivm) = MAX(transpir_supply(ipts,ivm), zero)
          ELSE
             !do nothing
          ENDIF   


       ! Hydraulic architecure is not defined i.e. no biomass
       ELSE

          transpir_supply(ipts,ivm) = zero
          psi_soilroot(ipts,ivm) = zero

!                DO i = 1, nlevels_tot
!                    transpir_supply_column(i, ipts, :) = zero
!                END DO ! i = 1, nlevels_tot


       ENDIF ! biomass leaves and roots.gt.zero
          
    ENDDO ! ivm = 2,nvm (loopref: ha6)

 ENDDO ! kjpindexipts  ! DO ipts = 1,kjpindex (loopref: ha1) 

 ! Debug
 IF (printlev_loc>=4) THEN
    WRITE(numout,*) 'transpir_supply, ', transpir_supply(:,:)
 ENDIF
 !-

 ! Send output variables to xios
 CALL xios_orchidee_send_field("psi_soil",psi_soil)
 CALL xios_orchidee_send_field("psi_soilroot",psi_soilroot)
 CALL xios_orchidee_send_field("transpir_unstressed",transpir*one_day/dt_sechiba)
 CALL xios_orchidee_send_field("transpir_supply",transpir_supply*one_day/dt_sechiba)
 CALL xios_orchidee_send_field("R_root",R_root_tot)
 CALL xios_orchidee_send_field("R_shoot",R_shoot_tot)
 CALL xios_orchidee_send_field("VESSEL_LOSS",vessel_loss)

END SUBROUTINE hydraulic_arch_calc


!! ==============================================================================================================================
!! SUBROUTINE   : mleb_diffuco_trans_co2_short
!!
!>\BRIEF        This subroutine recalculate carbon assimilation and stomatal 
!! conductance following vbeta3 recalculation in sechiba.f90
!!
!! DESCRIPTION  :\n
!! *** General:\n 
!! cresist, sum_rveget_rstruct, gs, assimi, and ci are recalculated sequentially according to 
!! the original equations in diffuco_trans_co2. Because this routine places inside of the loop
!! for pixel and pft, the subroutine itself doesn't contain loops for pixel and pft but uses 
!! ji and jv as inputs. Every equation should be identical with diffuco_trans_co2, so any
!! updates and changes from diffuco_trans_co2 should be applied here, too.

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

SUBROUTINE mleb_diffuco_trans_co2_short(assimi_lev, assim_param, Ca, gstot_frac, ji, jv, kjpindex, & 
                Light_Abs_Tot, Light_Tran_Tot, lai, lai_per_level, pb, qair, q_cdrag, &
                qsintmax, qsintveg, ratio_vbeta3, rstruct, rveget, swdown, temp_air, temp_growth, u, &
                v, vbeta23, vbeta3, veget_max, veget, info_limitphoto, JJ_out, cimean, cresist, &
                new_leaf_ci, vbeta3pot, gsmean, gpp)
    !
    !! 0. Variable and parameter declaration
    !
    !
    !! 0.1 Input variables
    !
    REAL(r_std), DIMENSION(:,:,:), INTENT(in)                     :: assimi_lev
    REAL(r_std), DIMENSION(:,:,:), INTENT(in)                     :: assim_param    !! min+max+opt temps (K), vcmax, vjmax 
                                                                                    !! for photosynthesis ($\mumolm^{-2}s^{-1}$)
    REAL(r_std), DIMENSION(kjpindex), INTENT(in)                  :: Ca             !! CO2 concentration inside the canopy
                                                                                    !! @tex ($\mu mol mol^{-1}$) @endtex
    REAL(r_std), DIMENSION(:,:,:), INTENT(in)                     :: gstot_frac     !! fraction of gs over nlevel_tot (-)
    INTEGER(i_std), INTENT(in)                                    :: ji             !! index for pixel
    INTEGER(i_std), INTENT(in)                                    :: jv             !! index for PFT
    INTEGER(i_std), INTENT(in)                                    :: kjpindex       !! Domain size (unitless)
    REAL(r_std), DIMENSION(:,:,:), INTENT(in)                     :: Light_Abs_Tot  !! Absorbed radiation per level for photosynthesis
    REAL(r_std), DIMENSION(:,:,:), INTENT(in)                     :: Light_Tran_Tot !! Transmitted radiation per level for photosynthesis
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in)              :: lai            !! Surface foliaire
    REAL(r_std), DIMENSION(:,:,:), INTENT(in)                     :: lai_per_level  !! This is the LAI per vertical level
                                                                                    !! @tex $(m^{2} m^{-2})$
    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)                  :: q_cdrag        !! Surface drag coefficient (-)
    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)              :: qsintveg       !! Water on vegetation due to interception 
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in)              :: ratio_vbeta3   !! Reduction rate of vbeta3 when transpir_supply 
                                                                                    !! is smaller than demand
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in)              :: rstruct        !! structural resistance @tex ($s m^{-1}$) @endtex
                                                                                    !! @tex ($\mu mol mol^{-1}$) @endtex
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in)              :: rveget         !! stomatal resistance of vegetation 
    REAL(r_std), DIMENSION(:), INTENT(in)                         :: swdown         !! Downwelling short wave flux 
                                                                                    !! @tex ($W m^{-2}$) @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)                  :: 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,nvm), INTENT(in)              :: vbeta23        !! Beta for fraction of wetted foliage that will 
                                                                                    !! transpire (unitless) 
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in)              :: vbeta3         !! Beta for Transpiration (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)              :: veget          !! Fraction of vegetation type (-)
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot+1), INTENT(in):: info_limitphoto!!Save information of limitation for photosynthesis
    REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in)  :: JJ_out

    !! 0.2 Output variables
    !
    !! 0.3 Modified variables
    !
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)        :: cimean       !! mean intercellular CO2 concentration 
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)        :: cresist      !! coefficient for resistances (-) 
    REAL(r_std), DIMENSION(kjpindex,nvm,nlai), INTENT(inout)   :: new_leaf_ci  !! intercellular CO2 concentration (ppm)
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)        :: vbeta3pot    !! Beta for Potential Transpiration
                                                                               !! @tex ($s m^{-1}$) @endtex
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout)        :: gsmean       !! mean stomatal conductance to CO2 (umol m-2 s-1)
    REAL(r_Std), DIMENSION(kjpindex,nvm), INTENT(inout)        :: gpp          !! Assimilation ((gC m^{-2} s^{-1}), total area)

    !! 0.4 Local variables
    !	
    REAL(r_std)                           :: assim_coeff_a, assim_coeff_b, &
                           assim_coeff_x1, assim_coeff_x2, assim_coeff_x3      !! Coefficient based on analytic solution in diffuco
    REAL(r_std), DIMENSION(nlevels_tot)   :: assimi                            !! assimilation (at a specific LAI level) 
    REAL(r_std)                           :: assimtot                          !! total assimilation 
    REAL(r_std)                           :: Cs_star                           !! Cs -based CO2 compensation point 
                                                                               !! in the absence of Rd (ubar)
    REAL(r_std)                           :: ci_star                           !! Ci -based CO2 compensation point in the absence of Rd (ubar)        
    REAL(r_std)                           :: fcyc                              !! Fraction of electrons at PSI that follow cyclic 
                                                                               !! transport around PSI (-)
    REAL(r_std)                           :: fvpd                              !! Factor for describing the effect of leaf-to-air vapour 
                                                                               !! difference on gs (-)
    REAL(r_std)                           :: gamma_star                        !! CO2 compensation point (ppm)
    REAL(r_std)                           :: gm                                !! Mesophyll diffusion conductance 
                                                                               !! (molCO2 mâ2 sâ1 barâ1) 
    REAL(r_std)                           :: gb                                !! Boundary-layer conductance (mol m^{-2} s^{-1} bar{-1})
    REAL(r_std), DIMENSION(nlevels_tot)   :: gs                                !! stomatal conductance to CO2 
                                                                               !! @tex ($\mol m^{-2} s^{-1}$) @endtex
    REAL(r_std)                           :: gstot                             !! total stomatal conductance to H2O
    REAL(r_std)                           :: g0var
    INTEGER(i_std)                        :: ilev                              !! index
    REAL(r_std)                           :: J2                                !! Rate of all e-transport through PSII 
                                                                               !! (umol e- m-2 s-1)
    REAL(r_std)                           :: jwind                             !! Wind module (m s^{-1})
    REAL(r_std)                           :: Kmc                               !! Michaelis–Menten constant of Rubisco for CO2 (mubar)
    REAL(r_std)                           :: KmO                               !! Michaelis–Menten constant of Rubisco for O2 (mubar)
    REAL(r_std)                           :: lai_total                         !! total LAI for column in question
    REAL(r_std)                           :: leafN_top                         !! Leaf N content - top of canopy (gN m-2[leaf])
    REAL(r_std), DIMENSION(nlevels_tot)   :: light_tran_to_level               !! Cumulative amount of light transmitted
                                                                               !! to a given level
    REAL(r_std)                           :: low_gamma_star                    !! Half of the reciprocal of Sc/o (bar bar-1)
    REAL(r_std)                           :: N_profile                         !! N-content reduction factor for a specific canopy depth (unitless) 
    REAL(r_std)                           :: N_Vcmax                           !! Nitrogen level dependance of Vcmacx and Jmax 
    REAL(r_std)                           :: ntot                              !! Total canopy (eg leaf) N (gN m-2[ground]) 
    REAL(r_std)                           :: nue                               !! Nitrogen use Efficiency with impact of leaf 
    REAL(r_std)                           :: Obs                               !! Bundle-sheath oxygen partial pressure (ubar)
    REAL(r_std), DIMENSION(kjpindex)      :: qsatt                             !! surface saturated humidity at 2m (??) 
    REAL(r_std), DIMENSION(nlevels_tot)   :: Rd                                !! Day respiration (respiratory CO2 release other than
                                                                               !! by photorespiration) (mumol CO2 m-2 s-1)
    REAL(r_std)                           :: Rm                                !! Day respiration in the mesophyll (umol CO2 m-2 s-1)
    REAL(r_std)                           :: Rdtot                             !! Total Day respiration (respiratory CO2 release other than 
                                                                               !! by photorespiration) (mumol CO2 m-2 s-1) 
    REAL(r_std)                           :: Sco                               !! Relative CO2 /O2 specificity factor for Rubisco (bar bar-1) 
    REAL(r_std)                           :: S_Vcmax_acclim_temp               !! Entropy term for Vcmax 
    REAL(r_std)                           :: speed                             !! wind speed @tex ($m s^{-1}$) @endtex
    REAL(r_std)                           :: sum_rveget_rstruct                !! recalculated rveget
                                                                               !! Final unit is @tex ($m s^{-1}$) @endtex
    REAL(r_std)                           :: T_gamma_star                      !! Temperature dependance of gamma_star (unitless)    
    REAL(r_std)                           :: T_gm                              !! Temperature dependance of gmw
    REAL(r_std)                           :: T_Kmc                             !! Temperature dependance of KmC (unitless)
    REAL(r_std)                           :: T_KmO                             !! Temperature dependance of KmO (unitless)
    REAL(r_std)                           :: T_Sco                             !! Temperature dependance of Sco
    REAL(r_std)                           :: T_Rd                              !! Temperature dependance of Rd (unitless)
    REAL(r_std)                           :: T_Vcmax                           !! Temperature dependance of Vcmax (unitless)
    REAL(r_std)                           :: VPD                               !! Vapor Pressure Deficit (kPa)
                                                                               !! age (umol CO2 (gN)-1 s-1)
    REAL(r_std)                           :: VpJ2                              !! e- transport-limited PEP carboxylation rate (umol CO2 m-2 s-1)
    REAL(r_std)                           :: vcmax                             !! maximum rate of carboxylation 
                                                                               !! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
    REAL(r_std)                           :: vc                                !! Maximum rate of Rubisco activity-limited carboxylation 
                                                                               !! (mumol CO2 m−2 s−1)
    REAL(r_std)                           :: vc2                               !! rate of carboxylation (at a specific LAI level) 
                                                                               !! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex 
    REAL(r_std)                           :: water_lim                         !! water limitation factor (0-1,unitless)
    REAL(r_std)                           :: zqsvegrap                         !! relative water quantity in the water 
                                                                               !! interception reservoir (0-1,unitless) 
    REAL(r_std)                           :: z                                 !! A lumped parameter (see Yin et al. 2009) ( mol mol-1) 
    REAL(r_std)                           :: Cc_loc                            !! Cellular CO2 concentration
!_ ================================================================================================================================
 
     ! 0. Initialization ++ CHECK dimension

     lai_total = zero
     assimi(:) = zero
     gs(:) = zero
     gstot = zero
     assimtot = zero
     Rd(:) = zero
     Rdtot=zero
     vbeta3pot(ji,jv) = zero
     gsmean(ji,jv) = zero
     gpp(ji,jv) = zero
     light_tran_to_level(:)=un
     cimean(ji,jv) = Ca(ji)
     vc2 = zero

     ! 1.1 Photosynthesis parameters
     !+++CHECK+++
     ! Why is total lai calculated as the sum of lai_per_layer rather than the more
     ! common approach by using ccc_to_lai ? lai(ji,jv) = cc_to_lai( circ_class_biomass(ji,jv,:, ileaf,&
     !      icarbon),circ_class_n(ji,jv,:),jv)?
     ! Is there a difference between both approaches? Probably best to test before changing this 
     ! code. I tried to test but the test requires running the multi-layer energy budget which
     ! crashed in stomate_soilcarbon.f90 see ticket 457.
     DO ilev=1,nlevels_tot
        !! total LAI for column in question
        lai_total = lai_total + lai_per_level(ji,jv,ilev)
     ENDDO

     ! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001)
     T_KmC        = Arrhenius(temp_air(ji),298.,E_KmC(jv))
     T_KmO        = Arrhenius(temp_air(ji),298.,E_KmO(jv))
     T_Sco        = Arrhenius(temp_air(ji),298.,E_Sco(jv))
     T_gamma_star = Arrhenius(temp_air(ji),298.,E_gamma_star(jv))
     ! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001)
     T_Rd = Arrhenius(temp_air(ji),298.,E_Rd(jv))
        
     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
     ! We derive its value from Equation 7 of Yin et al. (2009)
     ! assuming a constant O2 concentration (Oi)
     low_gamma_star = gamma_star/Oi

     ! 1.2 taking N into account for the calculations of Vcmax, as in
     ! diffuco_trans_co2.
     nue = assim_param(ji,jv,inue)
     ntot = assim_param(ji,jv,ileafn)
     IF (lai_total .GT. min_sechiba) THEN
             leafN_top = ntot * ext_coeff_N(jv) / ( 1.-exp(-ext_coeff_N(jv) * lai_total) )
     ELSE
             leafN_top = zero
     ENDIF
     vcmax = nue*leafN_top
     
     S_Vcmax_acclim_temp = aSV(jv) + bSV(jv) * MAX(11.,MIN(temp_growth(ji),35.))   
     T_Vcmax = Arrhenius_modified_nodim(temp_air(ji),298., &
                        E_Vcmax(jv),D_Vcmax(jv),S_Vcmax_acclim_temp)
     vc = vcmax * T_Vcmax

     ! 1.3 Estimate relative humidity of air
     CALL qsatcalc (kjpindex, temp_air, pb, qsatt)

     VPD = ( qsatt(ji)*pb(ji) / (Tetens_1+qsatt(ji)*Tetens_2 ) ) - &
              ( qair(ji) * pb(ji) / (Tetens_1+qair(ji)* Tetens_2) )
     ! VPD is needed in kPa
     ! We limit the impact of VPD in the range of [0:6] kPa
     ! This seems to be the typical range of values
     ! that the vapor pressure difference can take, although
     ! there is a paper by Mott and Parkhurst, Plant, Cell & Environment,
     ! Vol. 14, pp. 509--519 (1991) which gives a value of 15 kPa.
     VPD = MAX(zero,MIN(VPD/10.,6.))
     fvpd = 1. / ( 1. / (a1(jv) - b1(jv) * VPD) - 1. )
     
     gb = gb_ref * 44.6 * (tp_00/temp_air(ji)) * (pb(ji)/pb_std)
     ! conversion from (mol H2O m-2 s-1) to (mol CO2 m-2 s-1) 
     gb = gb / ratio_H2O_to_CO2

     ! 1.4 Intermediate variables needed for recalculation
     jwind = SQRT (u(ji)*u(ji) + v(ji)*v(ji))
     speed = MAX(min_wind, jwind)

     zqsvegrap = zero
     IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
         !! relative water quantity in the water interception reservoir
         zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
     ENDIF

     T_gm  = Arrhenius_modified_nodim(temp_air(ji),298.,E_gm(jv),D_gm(jv),S_gm(jv))
     ! water_lim is always 1 for this module, because it in ok_hydrol_arch already
     water_lim = un
     gm = gm25(jv) * T_gm * MAX(1-stress_gm(jv), water_lim)
     g0var = g0(jv)* MAX(1-stress_gs(jv), water_lim)

     ! give a default value of ci for all pixel that do not assimilate
     DO ilev=1,nlevels_tot
             new_leaf_ci(ji,jv,ilev) = Ca(ji)
     ENDDO

     ! 2. Recalculating assimilation
     
     ! 2.1 Recalculating resistance
     ! recalculation_1) creist
     ! If the leaf has a full film of water then the original equation for
     ! vbeta3 doesn't include cresist in it so the connection is lost.
     ! In that case, cresist is recalculated using a simple ratio of a reduction of
     ! vbeta3. Note that this approach follows the way of calculation in diffuco but 
     ! overlooks stomatal closure through water film.
     IF (zqsvegrap .EQ. un) THEN
         cresist(ji,jv) = cresist(ji,jv) *ratio_vbeta3(ji,jv)
     ELSE
         cresist(ji,jv)  = MAX(vbeta3(ji,jv) /veget(ji,jv) , &
                    (vbeta3(ji,jv) - vbeta23(ji,jv)) /(un - zqsvegrap)/veget(ji,jv))
     ENDIF

     ! recalculation_2) sum_rveget_rstruct
     IF (cresist(ji,jv) .GT. zero .AND. speed .GT. min_sechiba .AND. &
           q_cdrag(ji) .GT. min_sechiba .AND. rveg_pft(jv) .GT. min_sechiba) THEN
        
        sum_rveget_rstruct = (un/cresist(ji,jv) - un) / (speed * q_cdrag(ji) * &
             (veget(ji,jv)/veget_max(ji,jv)) * rveg_pft(jv))
     ELSE
        sum_rveget_rstruct = 9999.
     END IF

     vbeta3pot(ji,jv) = MAX(zero, veget_max(ji,jv) *cresist(ji,jv))

     ! recalculation_3) gstot
     IF (sum_rveget_rstruct .ne. 9999.) THEN
            gstot = un / sum_rveget_rstruct
     ELSE
            ! Need to convert the units of g0
            gstot = g0(jv) * mol_to_m_1 * (temp_air(ji)/tp_00) * &
                      (pb_std/pb(ji))*ratio_H2O_to_CO2
     END IF

     ! One equation normliznig gstot using
     ! lai_coupled is omitted because lai_coupled is the same with lai at this
     ! moment. It should be applied if the change in lai_coupled is made.

     ! we now do a unit conversion, from m/s to mol/m^2/s, the reverse
     !   of the calculation of section 2.5,

     gstot =  gstot / (mol_to_m_1 *(temp_air(ji)/tp_00)* &
                (pb_std/pb(ji))*ratio_H2O_to_CO2)

     DO ilev = nlevels_tot-1,1,-1
             light_tran_to_level(ilev)=light_tran_to_level(ilev+1)*Light_Tran_Tot(ji,jv,ilev+1)
     ENDDO

     DO ilev = 1, nlevels_tot
        IF (Light_Abs_Tot(ji,jv,ilev) .GT. min_sechiba) THEN
           IF(lai_per_level(ji,jv,ilev) .LE. min_sechiba)THEN
                ! This is not necessarily a problem.  LAI is updated in slowproc, whil
                ! absorbed light is updated in condveg, which is called after diffuco.
                ! So we might be one timestep off.  It generally should be at night,
                ! though, so there is no absorbed light.  If this persists, that
                ! is a real problem.
                IF(printlev_loc >= 3)THEN
                   WRITE(numout,*) 'WARNING: We have absorbed light but no LAI ' // &
                        'in this level. Skipping.'
                   WRITE(numout,'(A,3I8)') 'WARNING: ji,jv,ilev: ',ji,jv,ilev
                   WRITE(numout,'(A,2F20.10)') 'WARNING: Light_Abs_Tot(ji,jv,ilev), ' // &
                         'lai_per_level(ji,jv,ilev): ',&
                         Light_Abs_Tot(ji,jv,ilev), lai_per_level(ji,jv,ilev)
                ENDIF
                CYCLE
           ENDIF ! lai_per_level
           
           !
           ! nitrogen is scaled according to the transmitted light
           ! before the caclulation was (1 - exp(lai)) we replaced
           ! it by the transmitted light calculated by the 2stream
           ! model, maybe ( un - Light_Tran_Tot(iainia,jv,ilev) )
           ! could be replaced by the directly calculated absorbed
           ! light Light_Abs_Tot(iainia,jv,ilev).
           ! Notice that light in the old code was the cumulative light transmitted to
           ! the layer, while Light_Tran_Tot is the relative light transmitted
           ! through the layer.  Therefore, we use a different variable.
           !
           N_profile = ( un - .7_r_std * ( un - light_tran_to_level(ilev) ) )
           ! see Comment in legend of Fig. 6 of Yin et al. (2009)
           ! Rd25 is assumed to equal 0.01 Vcmax25         
           ! Effects of drought stress on respiration are less 
           ! well known: increases, decreases and no change have been found.
           ! For now there is no effect of drought stress on dark respiration 
           ! when hydrol_arch is used.
           vc2 = vc * N_profile * MAX(1-stress_vcmax(jv), water_lim)
           Rd(ilev) = vcmax * N_profile * 0.01 * T_Rd        

           IF (printlev_loc >= 4) THEN
                WRITE(numout,*) 'mleb_diffuco_trans_co2: recalculation starts'
           ENDIF

           ! 2.2 Assimilation for C4 plants 
           IF ( is_c4(jv) )  THEN
           
                ! Analytical resolution of the Assimilation based 
                ! Yin et al. (2009)
                ! 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(ilev)/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 * Oi - ( 1. + &
                   low_gamma_star * alpha(jv) / 0.047) * Rd(ilev) + Rm )/ &
                   ( gbs(jv) + kp(jv) )

                 ! eq. 11 of Yin et al (2009)
                 J2 = JJ_out(ji,jv,ilev) / ( 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.
                 ! See paragraph after eq. (20b) of Yin et al.
                 Rm=Rd(ilev)/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)
                 IF ( info_limitphoto(ji,jv,ilev) .EQ. 1 ) THEN
                    assim_coeff_a = 1. + kp(jv) / gbs(jv)
                    assim_coeff_b = 0.
                    assim_coeff_x1 = vc2
                    assim_coeff_x2 = KmC/KmO
                    assim_coeff_x3 = KmC
                 ! Is J limiting the Assimilation
                 ELSE IF (info_limitphoto(ji,jv,ilev) .EQ. 2) THEN
                    assim_coeff_a = 1.
                    assim_coeff_b = VpJ2
                    assim_coeff_x1 = (1.- fpsir(jv)) * J2 * z / 3.
                    assim_coeff_x2 = 7. * low_gamma_star / 3.
                    assim_coeff_x3 = 0.
                 ELSE
                    IF ( printlev_loc >= 4) THEN
                        WRITE(numout,*) 'limit_photo was not defined: no assimilation happened'
                        WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
                    ENDIF
                    assimi(ilev) = -Rd(ilev) 
                    CYCLE
                 ENDIF

                 ! recalculation_4) gs
                 gs(ilev) = gstot*gstot_frac(ji,jv,ilev) / lai_per_level(ji,jv,ilev)

                 IF ( gs(ilev) .LT. min_sechiba ) THEN
                     gs(ilev) = g0(jv)
                 ENDIF
                 ! IF (gs(iainia) .gt. gstot(iainia)) THEN
                 !    gs(iainia) = gstot(iainia)
                 ! END IF

                 IF (gs(ilev) .GE. g0(jv)) THEN
                    ! recalculation_5) assimi (C4) Equations copied from
                    ! diffuco_trans_co2 in diffuco.f90 and resuffled to isolate
                    ! assimi
                    assimi(ilev) = ( (gs(ilev)-g0(jv))*(Ca(ji)-Cs_star) - Rd(ilev)*fvpd ) / &
                                   ( fvpd +(gs(ilev)-g0(jv))/gb )
                 ELSE
                 
                    assimi(ilev) = -Rd(ilev)
                 ENDIF
                 
                 ! recalculation_6) Ci. Equations copied from diffuco_trans_co2 in diffuco.f90
                 Obs = ( alpha(jv) * assimi_lev(ji,jv,ilev) ) / ( 0.047 * gbs(jv) ) + Oi
                 ! Eq. 23 of Yin et al. (2009)
                 Cc_loc = (( assimi_lev(ji,jv,ilev) + Rd(ilev)) * &
                     ( assim_coeff_x2 * Obs + assim_coeff_x3 ) +&
                     low_gamma_star * Obs * assim_coeff_x1 )/ &
                     MAX(min_sechiba, assim_coeff_x1 - ( assimi_lev(ji,jv,ilev) + Rd(ilev) ))
                 ! Eq. 22 of Yin et al. (2009)
                 new_leaf_ci(ji,jv,ilev) = ( Cc_loc - ( assim_coeff_b - assimi(ilev) - Rm ) / &
                        gbs(jv) ) / assim_coeff_a
       
                 IF(printlev_loc >=4) THEN
                    WRITE(numout,*) 'after recalculation'
                    WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
                    WRITE(numout,*) 'cresist',cresist(ji,jv)
                    WRITE(numout,*) 'sum_rveget_rstruct',sum_rveget_rstruct
                    WRITE(numout,*) 'gs',gs(ilev)
                    WRITE(numout,*) 'assimi',assimi(ilev)
                    WRITE(numout,*) 'Rd', Rd(ilev)
                 ENDIF
 
                 IF (assimi(ilev) .lt. zero) THEN
                    assimi(ilev) = zero
                    !    WRITE(numout,*) 'assimi(iainia,jv) is negative ***************************************'
                    !    WRITE(numout,*) 'assimi(iainia,jv) is: ', assimi(iainia,jv)
                    !    WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
                    !    WRITE(numout,*) 'g0(jv) is: ', g0(jv)
                    !    WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
                    !    WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
                    !    WRITE(numout,*) 'Cs_star is: ', Cs_star
                    !    WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
                    !    WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
                    !    WRITE(numout,*) '*********************************************************'
                 ELSE IF (assimi(ilev) .gt. 500.) THEN
                    !    WRITE(numout,*) 'assimi(iainia,jv) is very large'
                    !    WRITE(numout,*) 'assimi(iainia,jv) is: ', assimi(iainia,jv)
                    !    WRITE(numout,*) 'temp_air(iainia) is: ', temp_air(iainia)
                    !    WRITE(numout,*) 'vbeta3(iainia, jv) is: ', vbeta3(iainia, jv)
                    !    WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
                    !    WRITE(numout,*) 'g0(jv) is: ', g0(jv)
                    !    WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
                    !    WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
                    !    WRITE(numout,*) 'Cs_star is: ', Cs_star
                    !    WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
                    !    WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
                    !    WRITE(numout,*) 'gstot(iainia) is: ', gstot(iainia)
                    !    WRITE(numout,*) 'sum_rveget_rstruct(iainia,jv) is: ', sum_rveget_rstruct(iainia,jv)
                    !    WRITE(numout,*) 'speed is: ', speed
                    !    WRITE(numout,*) 'q_cdrag(iainia) is: ',q_cdrag(iainia)
                    !    WRITE(numout,*) 'rveg_pft(jv) is: ', rveg_pft(jv)
                    !    WRITE(numout,*) 'gstot_frac(iainia,jv,ilev) is: ', gstot_frac(iainia,jv,ilev)
                    !    WRITE(numout,*) 'lai_per_level(iainia,jv,ilev) is: ', lai_per_level(iainia,jv,ilev)
                    !    WRITE(numout,*) '*********************************************************'            
                 END IF


           ELSE !is_c4
           
           ! 2.2 Assimilation for C3 plants 
                 ! Eq. 17 of Yin et al. (2009)
                 ! See eq. right after eq. 15 of Yin et al. (2009)
                 ci_star = gamma_star
                 
                 IF ( info_limitphoto(ji,jv,ilev) .EQ. 1 ) THEN
                    assim_coeff_x1 = vc2
                    ! It should be O not Oi (comment from Vuichard)
                    assim_coeff_x2 = KmC * ( 1. + 2*gamma_star*Sco/ KmO )
                 
                 ! Is J limiting the Assimilation
                 ELSE IF ( info_limitphoto(ji,jv,ilev) .EQ. 2) THEN
                    assim_coeff_x1 = JJ_out(ji,jv,ilev)/4.
                    assim_coeff_x2 = 2. * gamma_star
                 
                 ELSE
                    IF ( printlev_loc >= 4) THEN
                        WRITE(numout,*) 'limit_photo was not defined: no assimilation happened'
                        WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
                    ENDIF
                    assimi(ilev) = -Rd(ilev) 
                    CYCLE
                 ENDIF   
                 
                 ! recalculation_4) gs
                 gs(ilev) = gstot * gstot_frac(ji,jv,ilev) / lai_per_level(ji,jv,ilev)
                 ! IF (gs(iainia) .gt. gstot(iainia)) THEN
                 !    gs(iainia) = gstot(iainia)
                 ! END IF
                 IF ( gs(ilev) .LT. min_sechiba ) THEN
                     gs(ilev) = g0(jv)
                 ENDIF


                 IF(gs(ilev) .GE. g0(jv)) THEN
                                         
                 ! recalcaultion_5) assimi (C3)
                 ! Unlike assimi for c4 plants, assimi calculation for c3 plants contains Cc_loc
                 ! inside which needs 5 equations to combine including the cubic equation 
                 ! (Yin and Struik (2009), Eq 15 to 19), which can, in turn, make the 
                 ! recalculation quite complex. To avoid that, only Eq 15 and 16 were used
                 ! assuming that gm and Cc_loc will be constant. The result from the simplified
                 ! (below) equation is not 100% precise but precise enough to  accept 
                 ! (~0.3% lower) when it is checked in diffuco.
                 !
                 assimi(ilev) =(gb * gs(ilev) * (Ca(ji) * (gs(ilev) - g0(jv)) + &
                                  gamma_star * (g0(jv) - gs(ilev)) -Rd(ilev)*fvpd))/ &
                                  (gb * ( -g0(jv) +gs(ilev)*fvpd + gs(ilev)) + &
                                   gs(ilev)*(gs(ilev) - g0(jv)))
                 ELSE
                    assimi(ilev) = -Rd(ilev)
                    IF(printlev_loc >=4) THEN
                        WRITE(numout,*) 'gs < g0. assimi will be set to -Rd'
                        WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
                    ENDIF
                 ENDIF   


                 ! recalculation_6) Ci
                 ! Recalculated Ci is too low!! Decrease in assimi decreases Cc_loc way too much, which
                 ! results in too low leaf_ci. It can give a huge impact on C13, so
                 ! C13=y is prohibited for now. If assimi_lev is used instead of
                 ! assimi for Cc_loc calculation (keeping the Cc_loc without stress) the change in
                 ! leaf_ci is reasonable, but it was not applied because ci is not
                 ! written in restart file and it doesn't seem right.   

                 Cc_loc = (gamma_star * assim_coeff_x1 + (assimi(ilev)+Rd(ilev)) * assim_coeff_x2) / &
                        max(min_sechiba, assim_coeff_x1 - (assimi(ilev)+Rd(ilev)))
                 new_leaf_ci(ji,jv,ilev) = Cc_loc + assimi(ilev)/gm


                 IF(test_pft == jv .AND. printlev_loc >=4) THEN
                    WRITE(numout,*) 'after recalculation, C3'
                    WRITE(numout,*) 'ji,jv,ilev',ji,jv,ilev
                    WRITE(numout,*) 'cresist',cresist(ji,jv)
                    WRITE(numout,*) 'sum_rveget_rstruct',sum_rveget_rstruct
                    WRITE(numout,*) 'gs',gs(ilev)
                    WRITE(numout,*) 'assimi_recal',assimi(ilev)
                    WRITE(numout,*) 'Rd', Rd(ilev)
                    WRITE(numout,*) 'Cc_recal',Cc_loc
                    WRITE(numout,*) 'Ci_recal',new_leaf_ci(ji,jv,ilev)
                 ENDIF

                 IF (assimi(ilev) .lt. min_sechiba) THEN
                    assimi(ilev) = zero
                    !    WRITE(numout,*) 'assimi(iainia,jv) is negative ***************************************'
                    !    WRITE(numout,*) 'assimi(iainia,jv) is: ', assimi(iainia,jv)
                    !    WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
                    !    WRITE(numout,*) 'g0(jv) is: ', g0(jv)
                    !    WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
                    !    WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
                    !    WRITE(numout,*) 'Cs_star is: ', Cs_star
                    !    WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
                    !    WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
                    !    WRITE(numout,*) '*********************************************************'
                 ELSE IF (assimi(ilev) .gt. 500.) THEN
                    !    WRITE(numout,*) 'assimi(iainia,jv) is very large'
                    !    WRITE(numout,*) 'assimi(iainia,jv) is: ', assimi(iainia,jv)
                    !    WRITE(numout,*) 'temp_air(iainia) is: ', temp_air(iainia)
                    !    WRITE(numout,*) 'vbeta3(iainia, jv) is: ', vbeta3(iainia, jv)
                    !    WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
                    !    WRITE(numout,*) 'g0(jv) is: ', g0(jv)
                    !    WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
                    !    WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
                    !    WRITE(numout,*) 'Cs_star is: ', Cs_star
                    !    WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
                    !    WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
                    !    WRITE(numout,*) 'gstot(iainia) is: ', gstot(iainia)
                    !    WRITE(numout,*) 'sum_rveget_rstruct(iainia,jv) is: ', sum_rveget_rstruct(iainia,jv)
                    !    WRITE(numout,*) 'speed is: ', speed
                    !    WRITE(numout,*) 'q_cdrag(iainia) is: ',q_cdrag(iainia)
                    !    WRITE(numout,*) 'rveg_pft(jv) is: ', rveg_pft(jv)
                    !    WRITE(numout,*) 'gstot_frac(iainia,jv,ilev) is: ', gstot_frac(iainia,jv,ilev)
                    !    WRITE(numout,*) 'lai_per_level(iainia,jv,ilev) is: ', lai_per_level(iainia,jv,ilev)
                    !    WRITE(numout,*) '*********************************************************'                                
                 END IF
           ENDIF !is_c4
        
           IF(test_pft == jv .AND. printlev_loc>=4)THEN
                WRITE(numout,*) 'assimi ',ilev,ji,assimi(ilev)
                WRITE(numout,*) 'lai_per_level',lai_per_level(ji,jv,ilev)
           ENDIF
    
        ENDIF !IF Light
     ENDDO !DO ilev

     !Initialize gstot to calculate again with recalculated gs
     gstot = zero
     !! 2.3 Integration at the canopy level
     DO ilev = 1,nlevels_tot
        assimtot = assimtot + assimi(ilev) * lai_per_level(ji,jv,ilev)
        Rdtot = Rdtot + &
             Rd(ilev) * lai_per_level(ji,jv,ilev)
        gstot = gstot + &
             gs(ilev) * lai_per_level(ji,jv,ilev)
     ENDDO

     IF (err_act.GT.1 .AND. assimtot .gt. 20) THEN
        WRITE(numout,*) 'WARNING: assimtot is very high'
        WRITE(numout,*) 'assimtot is: ', assimtot
        DO ilev=1,nlevels_tot
           WRITE (numout,*) 'level: ', ilev, ' pft: ', jv , ' assimi_store: ', assimi(ilev), &
                            &' lai_per_level: ', lai_per_level(ji, jv, ilev)
        END DO
     END IF

     !! 2.4 Calculate resistances
     ! Copy from diffuco.f90

     !! Mean stomatal conductance for CO2 (mol m-2 s-1)
     gsmean(ji,jv) = gstot

     ! cimean is completely diagnostic and not
     ! used for anything.  If the analytical photosynthesis is unable to find
     ! a solution for the A_1 coefficients in any of the levels, assimtot will
     ! be equal to -Rd, which means we will be dividing by zero.  This loop
     ! is to prevent that, and since cimean is diagnostic is doesn't matter
     ! what the value is.
     IF ( ABS(gsmean(ji,jv)-g0var*lai_total) .GT. min_sechiba) THEN
        cimean(ji,jv) = (fvpd*(assimtot+Rdtot)) /&
                (gsmean(ji,jv)-g0var*lai_total) + gamma_star
     ELSE
        cimean(ji,jv) = gamma_star
     ENDIF

     gpp(ji,jv) = assimtot*12e-6*veget_max(ji,jv)*dt_sechiba

END SUBROUTINE mleb_diffuco_trans_co2_short


!! ADD HEADER

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

 SUBROUTINE jnetcdf_trans_co2(kjpindex, james_netcdf_flag, gs_distribution_in, gs_output_in, &
                               lai_per_level_in, light_tran_to_level_in, gs_diffuco_output_in, &
                               gstot_output_in) 
 
    USE netcdf

    IMPLICIT NONE

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                                         :: kjpindex         !! Domain size (unitless)

    CHARACTER(LEN=80)                                                  :: jamescdfname          
 
    REAL(r_std), DIMENSION (nvm, nlevels_tot), INTENT(in)              :: gs_distribution_in     
    REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in)               :: gs_output_in
    REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in)               :: gstot_output_in
    REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in)               :: lai_per_level_in
    REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in)               :: light_tran_to_level_in
    REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in)               :: gs_diffuco_output_in


    LOGICAL                                                            :: james_netcdf_flag    


    !! 0.2 Output variables

    !! 0.3 Modified variables

    !! 0.4 Local variables

    INTEGER(i_std), SAVE                     :: gsd_id
!$OMP THREADPRIVATE(gsd_id)
    INTEGER(i_std), SAVE                     :: gso_id
!$OMP THREADPRIVATE(gso_id)
    INTEGER(i_std), SAVE                     :: gsto_id
!$OMP THREADPRIVATE(gsto_id)
    INTEGER(i_std), SAVE                     :: lpl_id
!$OMP THREADPRIVATE(lpl_id)
    INTEGER(i_std), SAVE                     :: lttl_id
!$OMP THREADPRIVATE(lttl_id)
    INTEGER(i_std), SAVE                     :: gsdo_id
!$OMP THREADPRIVATE(gsdo_id)

    INTEGER(i_std)                           :: ii, jj
    INTEGER(i_std), SAVE                     :: mleb_step_count    
!$OMP THREADPRIVATE(mleb_step_count)

    INTEGER(i_std), SAVE                     :: iret, iret2, fid, fid2
!$OMP THREADPRIVATE(iret,iret2,fid,fid2)
    INTEGER(i_std)                           :: TimeDimID, jtid, ts_id, HeightDimID, i
    INTEGER(i_std)                           :: height_dim_id, height_dim_id2, time_dim_id
    INTEGER(i_std)                           :: pftno_dim_id
    INTEGER(i_std)                           :: dimids(3) 
    INTEGER(i_std), SAVE                     :: tap_id
!$OMP THREADPRIVATE(tap_id)
    INTEGER, DIMENSION(3)                    :: jstartnc, jcountnc    


    jamescdfname = '/home/users/jryder/james.ryder/outputcdf.nc'
 
    IF (james_netcdf_flag) THEN

        james_netcdf_flag = .FALSE.
        
        ! mleb_step_count = 0
 
        ! WRITE(numout,*) 'james_netcdf_flag call here'

        ! … define the parameters here
 
        iret = nf90_create(jamescdfname, nf90_write, fid)
        
  
        ! define the dimensions

        iret = nf90_def_dim(fid, "pftno", nvm, pftno_dim_id)        
        iret = nf90_def_dim(fid, "height", nlevels_tot, height_dim_id)
        iret = nf90_def_dim(fid, "time", nf90_unlimited, time_dim_id)

        dimids = (/ pftno_dim_id, height_dim_id, time_dim_id /)

 
        ! define the variables
        iret = nf90_def_var (fid, "gs_distribution_in", nf90_real4, dimids, gsd_id)
        iret = nf90_def_var (fid, "gs_output_in", nf90_real4, dimids, gso_id)
        iret = nf90_def_var (fid, "gstot_output_in", nf90_real4, dimids, gsto_id)
        iret = nf90_def_var (fid, "lai_per_level_in", nf90_real4, dimids, lpl_id)
        iret = nf90_def_var (fid, "light_tran_to_level_in", nf90_real4, dimids, lttl_id)
        iret = nf90_def_var (fid, "gs_diffuco_output_in", nf90_real4, dimids, gsdo_id)

 
        iret = nf90_enddef (fid)
         
    END IF !(james_netcdf_flag)


    ! open the netCDF file at each step
    iret = nf90_open(jamescdfname, nf90_Write, fid) 

    jstartnc = (/ 1, 1, 1 /)  ! vector of integers specifying the index in the variable from 
                              !   which the first of the data values will be read
    jcountnc = (/ nvm, nlevels_tot, 1 /)  ! vector of integers specifying the number of indices 
                                     !   selected along each dimension

    mleb_step_count = mleb_step_count + 1    ! this is the time step counter
    
        
    jstartnc(3) = mleb_step_count   ! (2) refers to the time dimension at present (it 
                                     !   must always come last)
 

    iret = nf90_put_var(fid, gsd_id, gs_distribution_in, start = jstartnc, count = jcountnc) 
    iret = nf90_put_var(fid, gso_id, gs_output_in, start = jstartnc, count = jcountnc) 
    iret = nf90_put_var(fid, gsto_id, gstot_output_in, start = jstartnc, count = jcountnc) 
    iret = nf90_put_var(fid, lpl_id, lai_per_level_in, start = jstartnc, count = jcountnc) 
    iret = nf90_put_var(fid, lttl_id, light_tran_to_level_in, start = jstartnc, count = jcountnc) 
    iret = nf90_put_var(fid, gsdo_id, gs_diffuco_output_in, start = jstartnc, count = jcountnc) 



    ! close the netCDF file and write to disk                   
    iret = nf90_close(fid)

 
 END SUBROUTINE jnetcdf_trans_co2

 
 ! ------------------------------------------------------------------------------------------------------------------------------------- 



END MODULE hydraulic_arch
