!  ==============================================================================================================================\n
!  MODULE 	: sechiba
! 
!  CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
!
!  LICENCE      : IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF        Structures the calculation of atmospheric and hydrological 
!! variables by calling diffuco_main, enerbil_main, hydrolc_main (or hydrol_main),
!! enerbil_fusion, condveg_main and thermosoil_main. Note that sechiba_main
!! calls slowproc_main and thus indirectly calculates the biogeochemical
!! processes as well.
!!
!!\n DESCRIPTION  : :: shumdiag, :: litterhumdiag and :: stempdiag are not 
!! saved in the restart file because at the first time step because they 
!! are recalculated. However, they must be saved as they are in slowproc 
!! which is called before the modules which calculate them.
!! 
!! RECENT CHANGE(S): None 
!! 
!! REFERENCE(S) : None
!!   
!! SVN     :
!! $HeadURL$ 
!! $Date$
!! $Revision$
!! \n
!_ ================================================================================================================================
 
MODULE sechiba
 
  USE ioipsl
  USE xios_orchidee
  USE constantes
  USE constantes_soil
  USE pft_parameters
  USE diffuco
  USE condveg
  USE enerbil
  USE hydrol                                                                     !! 
  USE hydrolc
  USE thermosoil
  USE thermosoilc
  USE sechiba_io
  USE slowproc
  USE routing
  use ioipsl_para


  IMPLICIT NONE

  PRIVATE
  PUBLIC sechiba_main, sechiba_initialize, sechiba_clear

  INTEGER(i_std), SAVE                             :: printlev_loc   !! local printlev for this module
!$OMP THREADPRIVATE(printlev_loc)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexveg       !! indexing array for the 3D fields of vegetation
!$OMP THREADPRIVATE(indexveg)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlai       !! indexing array for the 3D fields of vegetation
!$OMP THREADPRIVATE(indexlai)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnobio     !! indexing array for the 3D fields of other surfaces (ice,
                                                                     !! lakes, ...)
!$OMP THREADPRIVATE(indexnobio)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsoil      !! indexing array for the 3D fields of soil types (kjpindex*nstm)
!$OMP THREADPRIVATE(indexsoil)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexgrnd      !! indexing array for the 3D ground heat profiles (kjpindex*ngrnd)
!$OMP THREADPRIVATE(indexgrnd)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlayer     !! indexing array for the 3D fields of soil layers in CWRR (kjpindex*nslm)
!$OMP THREADPRIVATE(indexlayer)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnbdl      !! indexing array for the 3D fields of diagnostic soil layers (kjpindex*nbdl)
!$OMP THREADPRIVATE(indexnbdl)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexalb       !! indexing array for the 2 fields of albedo
!$OMP THREADPRIVATE(indexalb)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsnow      !! indexing array for the 3D fields snow layers
!$OMP THREADPRIVATE(indexsnow)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: veget          !! Fraction of vegetation type (unitless, 0-1)       
!$OMP THREADPRIVATE(veget)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: veget_max      !! Max. fraction of vegetation type (LAI -> infty, unitless)
!$OMP THREADPRIVATE(veget_max)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tot_bare_soil  !! Total evaporating bare soil fraction 
!$OMP THREADPRIVATE(tot_bare_soil)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: height         !! Vegetation Height (m)
!$OMP THREADPRIVATE(height)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: totfrac_nobio  !! Total fraction of continental ice+lakes+cities+...
                                                                     !! (unitless, 0-1)
!$OMP THREADPRIVATE(totfrac_nobio)
  REAL(r_std),ALLOCATABLE, SAVE,  DIMENSION (:,:)  :: wat_flux0      !! Water flux in the first soil layers exported for soil C calculations
!$OMP THREADPRIVATE(wat_flux0)
  REAL(r_std),ALLOCATABLE, SAVE,  DIMENSION (:,:,:) :: wat_flux       !! Water flux in the soil layers exported for soil C calculations
!$OMP THREADPRIVATE(wat_flux)
  REAL(r_std),ALLOCATABLE, SAVE,  DIMENSION (:,:)  :: drainage_per_soil !! Drainage per soil type exported for soil C calculations
!$OMP THREADPRIVATE(drainage_per_soil)
  REAL(r_std),ALLOCATABLE, SAVE,  DIMENSION (:,:)  :: runoff_per_soil !! Runoff per soil type exported for soil C calculations
!$OMP THREADPRIVATE(runoff_per_soil)
  REAL, ALLOCATABLE, DIMENSION(:,:,:) :: soil_mc
!$OMP THREADPRIVATE(soil_mc)
  REAL, ALLOCATABLE, DIMENSION(:,:) :: litter_mc
!$OMP THREADPRIVATE(litter_mc)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: floodout       !! Flow out of floodplains from hydrol
!$OMP THREADPRIVATE(floodout)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: runoff         !! Surface runoff calculated by hydrol or hydrolc 
                                                                     !! @tex $(kg m^{-2})$ @endtex
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: DOC_EXP_agg    !! DOC exports, diffrenet paths (nexp), in  
                                                                     !! @tex $(gC m^{-2} dt_slow^{-1})$ @endtex 
!$OMP THREADPRIVATE(runoff)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: drainage       !! Deep drainage calculatedd by hydrol or hydrolc 
                                                                     !! @tex $(kg m^{-2})$ @endtex
!$OMP THREADPRIVATE(drainage)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: returnflow     !! Water flow from lakes and swamps which returns to 
                                                                     !! the grid box @tex $(kg m^{-2})$ @endtex
!$OMP THREADPRIVATE(returnflow)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: reinfiltration !! Routed water which returns into the soil
!$OMP THREADPRIVATE(reinfiltration)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: irrigation     !! Irrigation flux taken from the routing reservoirs and 
                                                                     !! being put into the upper layers of the soil 
                                                                     !! @tex $(kg m^{-2})$ @endtex                                           
!$OMP THREADPRIVATE(irrigation)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: DOC_to_topsoil !! Sum of DOC fluxes from irrirgation and reinfiltration 
                                                                     !! @tex $(g m^{-2})$ @endtex                                           
!$OMP THREADPRIVATE(DOC_to_topsoil)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: DOC_to_subsoil !! DOC fluxes from returnflow in lakes and swamps 
                                                                     !! @tex $(g m^{-2})$ @endtex                                           
!$OMP THREADPRIVATE(DOC_to_subsoil)
  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)      :: precip2canopy  !! Precipitation onto the canopy
!$OMP THREADPRIVATE(precip2canopy)
  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)      :: precip2ground  !! Precipitation not intercepted by canopy
!$OMP THREADPRIVATE(precip2ground)
  REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:)      :: canopy2ground  !! Water flux from canopy to the ground
!$OMP THREADPRIVATE(canopy2ground)  
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: emis           !! Surface emissivity (unitless)
!$OMP THREADPRIVATE(emis)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: z0             !! Surface roughness (m)
!$OMP THREADPRIVATE(z0)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: roughheight    !! Effective height for roughness (m)
!$OMP THREADPRIVATE(roughheight)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: reinf_slope    !! slope coefficient (reinfiltration) 
!$OMP THREADPRIVATE(reinf_slope)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag       !! Mean relative soil moisture in the different levels used 
                                                                     !! by thermosoil.f90 (unitless, 0-1)
!$OMP THREADPRIVATE(shumdiag)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag_perma !! Saturation degree of the soil 
!$OMP THREADPRIVATE(shumdiag_perma)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: k_litt         !! litter cond.
!$OMP THREADPRIVATE(k_litt)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: litterhumdiag  !! Litter dryness factor (unitless, 0-1)
!$OMP THREADPRIVATE(litterhumdiag)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: stempdiag      !! Temperature which controls canopy evolution (K)
!$OMP THREADPRIVATE(stempdiag)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintveg       !! Water on vegetation due to interception 
                                                                     !! @tex $(kg m^{-2})$ @endtex
!$OMP THREADPRIVATE(qsintveg)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta2         !! Interception resistance (unitless,0-1)
!$OMP THREADPRIVATE(vbeta2)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta3         !! Vegetation resistance (unitless,0-1)
!$OMP THREADPRIVATE(vbeta3)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vbeta3pot      !! Potential vegetation resistance
!$OMP THREADPRIVATE(vbeta3pot)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: gsmean         !! Mean stomatal conductance for CO2 (umol m-2 s-1) 
!$OMP THREADPRIVATE(gsmean) 
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: cimean         !! STOMATE: mean intercellular CO2 concentration (ppm)
!$OMP THREADPRIVATE(cimean)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vevapwet       !! Interception loss over each PFT 
                                                                     !! @tex $(kg m^{-2} days^{-1})$ @endtex
!$OMP THREADPRIVATE(vevapwet)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: transpir       !! Transpiration @tex $(kg m^{-2} days^{-1})$ @endtex
!$OMP THREADPRIVATE(transpir)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: transpot       !! Potential Transpiration (needed for irrigation)
!$OMP THREADPRIVATE(transpot)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: qsintmax       !! Maximum amount of water in the canopy interception 
                                                                     !! reservoir @tex $(kg m^{-2})$ @endtex
!$OMP THREADPRIVATE(qsintmax)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rveget         !! Surface resistance for the vegetation 
                                                                     !! @tex $(s m^{-1})$ @endtex
!$OMP THREADPRIVATE(rveget)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rstruct        !! Vegetation structural resistance
!$OMP THREADPRIVATE(rstruct)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: snow_nobio     !! Snow mass of non-vegetative surfaces 
                                                                     !! @tex $(kg m^{-2})$ @endtex
!$OMP THREADPRIVATE(snow_nobio)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: snow_nobio_age !! Snow age on non-vegetative surfaces (days)
!$OMP THREADPRIVATE(snow_nobio_age)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: frac_nobio     !! Fraction of non-vegetative surfaces (continental ice, 
                                                                     !! lakes, ...) (unitless, 0-1)
!$OMP THREADPRIVATE(frac_nobio)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: albedo         !! Surface albedo for visible and near-infrared 
                                                                     !! (unitless, 0-1)
!$OMP THREADPRIVATE(albedo)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: assim_param    !! min+max+opt temps, vcmax, vjmax for photosynthesis
!$OMP THREADPRIVATE(assim_param)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lai            !! Surface foliaire
!$OMP THREADPRIVATE(lai)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: gpp            !! STOMATE: GPP. gC/m**2 of total area
!$OMP THREADPRIVATE(gpp)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)     :: temp_growth       !! Growth temperature (°C) - Is equal to t2m_month 
!$OMP THREADPRIVATE(temp_growth) 
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: humrel         !! Relative humidity
!$OMP THREADPRIVATE(humrel)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: vegstress      !! Vegetation moisture stress (only for vegetation growth)
!$OMP THREADPRIVATE(vegstress)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: frac_age       !! Age efficacity from STOMATE for isoprene 
!$OMP THREADPRIVATE(frac_age)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: soiltile       !! Fraction of each soil tile (0-1, unitless)
!$OMP THREADPRIVATE(soiltile)
  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: njsc           !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
!$OMP THREADPRIVATE(njsc)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta1         !! Snow resistance 
!$OMP THREADPRIVATE(vbeta1)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta4         !! Bare soil resistance
!$OMP THREADPRIVATE(vbeta4)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta5         !! Floodplains resistance
!$OMP THREADPRIVATE(vbeta5)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilcap        !!
!$OMP THREADPRIVATE(soilcap)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: soilflx        !!
!$OMP THREADPRIVATE(soilflx)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: temp_sol       !! Soil temperature
!$OMP THREADPRIVATE(temp_sol)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: qsurf          !! near soil air moisture
!$OMP THREADPRIVATE(qsurf)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: flood_res      !! flood reservoir estimate
!$OMP THREADPRIVATE(flood_res)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: fastr          !! fast reservoir estimate
!$OMP THREADPRIVATE(fastr)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: flood_frac     !! flooded fraction
!$OMP THREADPRIVATE(flood_frac)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: stream_frac     !! stream fraction
!$OMP THREADPRIVATE(stream_frac)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: streamfl_frac   !! flooded by swollen stream fraction
!$OMP THREADPRIVATE(streamfl_frac)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: snow           !! Snow mass [Kg/m^2]
!$OMP THREADPRIVATE(snow)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: snow_age       !! Snow age
!$OMP THREADPRIVATE(snow_age)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: drysoil_frac   !! Fraction of visibly (albedo) Dry soil (Between 0 and 1)
!$OMP THREADPRIVATE(drysoil_frac)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: rsol           !! resistance to bare soil evaporation
!$OMP THREADPRIVATE(rsol)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evap_bare_lim  !! Bare soil stress
!$OMP THREADPRIVATE(evap_bare_lim)

  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: co2_flux       !! CO2 flux (gC/m**2 of average ground/s)
!$OMP THREADPRIVATE(co2_flux)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evapot         !! Soil Potential Evaporation
!$OMP THREADPRIVATE(evapot)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evapot_corr    !! Soil Potential Evaporation Correction (Milly 1992)
!$OMP THREADPRIVATE(evapot_corr)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapflo       !! Floodplains evaporation
!$OMP THREADPRIVATE(vevapflo)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapsno       !! Snow evaporation
!$OMP THREADPRIVATE(vevapsno)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vevapnu        !! Bare soil evaporation
!$OMP THREADPRIVATE(vevapnu)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: t2mdiag        !! 2 meter temperature
!$OMP THREADPRIVATE(t2mdiag)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: tot_melt       !! Total melt
!$OMP THREADPRIVATE(tot_melt)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: vbeta          !! Resistance coefficient
!$OMP THREADPRIVATE(vbeta)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: valpha         !! Resistance coefficient
!$OMP THREADPRIVATE(valpha)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: fusion         !!
!$OMP THREADPRIVATE(fusion)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: rau            !! Density
!$OMP THREADPRIVATE(rau)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: deadleaf_cover !! Fraction of soil covered by dead leaves
!$OMP THREADPRIVATE(deadleaf_cover)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: ptnlev1        !! 1st level Different levels soil temperature
!$OMP THREADPRIVATE(ptnlev1)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: mc_layh        !! Volumetric soil moisture for each layer in hydrol(liquid + ice) (m3/m3)
!$OMP THREADPRIVATE(mc_layh)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: mcl_layh       !! Volumetric soil moisture for each layer in hydrol(liquid) (m3/m3)
!$OMP THREADPRIVATE(mcl_layh)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: tmc_layh       !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm)
!$OMP THREADPRIVATE(tmc_layh)

  LOGICAL, SAVE                                    :: l_first_sechiba = .TRUE. !! Flag controlling the intialisation (true/false)
!$OMP THREADPRIVATE(l_first_sechiba)

  ! Variables related to snow processes calculations  

  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowrho      !! snow density for each layer
!$OMP THREADPRIVATE(snowrho)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowheat     !! snow heat content for each layer (J/m2)
!$OMP THREADPRIVATE(snowheat)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowgrain    !! snow grain size (m)
!$OMP THREADPRIVATE(snowgrain)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowtemp     !! snow temperature profile (K)
!$OMP THREADPRIVATE(snowtemp)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: snowdz       !! snow layer thickness (m)
!$OMP THREADPRIVATE(snowdz)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: gtemp        !! soil surface temperature
!$OMP THREADPRIVATE(gtemp)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: pgflux       !! net energy into snow pack
!$OMP THREADPRIVATE(pgflux)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: cgrnd_snow   !! Integration coefficient for snow numerical scheme
!$OMP THREADPRIVATE(cgrnd_snow)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:)  :: dgrnd_snow   !! Integration coefficient for snow numerical scheme
!$OMP THREADPRIVATE(dgrnd_snow)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: lambda_snow  !! Coefficient of the linear extrapolation of surface temperature 
                                                                  !! from the first and second snow layers
!$OMP THREADPRIVATE(lambda_snow)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:)    :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
!$OMP THREADPRIVATE(temp_sol_add)

  LOGICAL, SAVE                                   :: ok_doc= .TRUE.
!$OMP THREADPRIVATE(ok_doc)   
  CHARACTER(LEN=80)    :: var_name       !! To store variables names for I/O (unitless)
  CHARACTER(LEN=10)    :: flow_suff      !! To assign a suffix to the variables name indicating the transported specie
  CHARACTER(LEN=80)    :: var_long_name  !! To store variables long names for I/O (unitless) 
CONTAINS

!!  =============================================================================================================================
!! SUBROUTINE:    sechiba_initialize
!!
!>\BRIEF	  Initialize all prinicipal modules by calling their "_initialize" subroutines
!!
!! DESCRIPTION:	  Initialize all prinicipal modules by calling their "_initialize" subroutines
!!
!! \n
!_ ==============================================================================================================================

  SUBROUTINE sechiba_initialize( &
       kjit,         kjpij,        kjpindex,     index,      date0,       &
       lalo,         contfrac,     neighbours,   resolution, zlev,        &
       u,            v,            qair,         t2m,        temp_air,    &
       petAcoef,     peqAcoef,     petBcoef,     peqBcoef,                &
       precip_rain,  precip_snow,  lwdown,       swnet,      swdown,      &
       pb,           rest_id,      hist_id,      hist2_id,                &
       rest_id_stom, hist_id_stom, hist_id_stom_IPCC,                     &
       coastalflow,  riverflow,    tsol_rad,     vevapp,       qsurf_out, &
       z0_out,       albedo_out,   fluxsens,     fluxlat,      emis_out,  &
       netco2flux,   fco2_lu,      temp_sol_new, tq_cdrag)

!! 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 
                                                                                  !! (unitless)
    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only 
                                                                                  !! (unitless)
    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
    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)                               :: rest_id_stom      !! STOMATE's _Restart_ file identifier 
                                                                                  !! (unitless)
    INTEGER(i_std),INTENT (in)                               :: hist_id_stom      !! STOMATE's _History_ file identifier 
                                                                                  !! (unitless)
    INTEGER(i_std),INTENT(in)                                :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file 
                                                                                  !! identifier (unitless)
    REAL(r_std), INTENT (in)                                 :: date0             !! Initial date (??unit??)
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)          :: lalo              !! Geographic coordinates (latitude,longitude)
                                                                                  !! for grid cells (degrees)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: contfrac          !! Fraction of continent in the grid 
                                                                                  !! (unitless, 0-1)
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map. 
                                                                                  !! Sechiba uses a reduced grid excluding oceans
                                                                                  !! ::index contains the indices of the 
                                                                                  !! terrestrial pixels only! (unitless)
    INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in)       :: neighbours        !! Neighboring grid points if land!(unitless)
    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)          :: resolution        !! Size in x and y of the grid (m)
    
    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)            :: zlev              !! Height of first layer (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair              !! Lowest level specific humidity 
                                                                                  !! @tex $(kg kg^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: t2m               !! 2m air temperature (K)
    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)


!! 0.2 Output variables
    REAL(r_std),DIMENSION (kjpindex,nflow), INTENT (out)     :: coastalflow       !! Outflow on coastal points by small basins.
                                                                                  !! This is the water which flows in a disperse 
                                                                                  !! way into the ocean
                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex,nflow), INTENT (out)     :: riverflow         !! Outflow of the major rivers.
                                                                                  !! The flux will be located on the continental 
                                                                                  !! grid but this should be a coastal point  
                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad          !! Radiative surface temperature 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp            !! Total of evaporation 
                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
    
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf_out         !! Surface specific humidity 
                                                                                  !! @tex $(kg kg^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0_out            !! Surface roughness (output diagnostic, m)
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out)         :: albedo_out        !! VIS and NIR albedo (output diagnostic, 
                                                                                  !! unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens          !! Sensible heat flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat           !! Latent heat flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: emis_out          !! Emissivity (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: netco2flux        !! Sum CO2 flux over PFTs 
                                                                                  !! ??(gC m^{-2} s^{-1})??
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fco2_lu           !! Land Cover Change CO2 flux 
                                                                                  !! ??(gC m^{-2} s^{-1})??
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: temp_sol_new      !! New ground temperature (K)

!! 0.3 Modified
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: tq_cdrag          !! Surface drag coefficient 
                                                                                  !! @tex $(m.s^{-1})$ @endtex 

!! 0.4 Local variables
    INTEGER(i_std)                                           :: ji, jv		  !! Index (unitless)
    REAL(r_std), DIMENSION(kjpindex)                         :: histvar           !! Computations for history files (unitless)
    CHARACTER(LEN=80)                                        :: var_name          !! To store variables names for I/O (unitless)
    REAL(r_std), DIMENSION(kjpindex)                         :: frac_snow_veg     !! Snow cover fraction on vegetation, 
                                                                                  !! only for diagnostics (unitless)
    REAL(r_std), DIMENSION(kjpindex,nnobio)                  :: frac_snow_nobio   !! Snow cover fraction on continental ice, lakes, etc
                                                                                  !! only for diagnostics (unitless)

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

    IF (printlev>=3) WRITE(numout,*) ' sechiba kjpindex =',kjpindex

    !! 1. Initialize variables on first call

    !! 1.1 Initialize most of sechiba's variables
    CALL sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)
    
    !! 1.3 Initialize stomate's variables

    CALL slowproc_initialize (kjit,          kjpij,        kjpindex,      date0,              &
                              rest_id,       rest_id_stom, hist_id_stom,   hist_id_stom_IPCC, &
                              index,         indexveg,     lalo,           neighbours,        &
                              resolution,    contfrac,     t2m,                               &
                              soiltile,      reinf_slope,  deadleaf_cover, assim_param,       &
                              lai,           frac_age,     height,         veget,             &
                              frac_nobio,    njsc,         veget_max,      tot_bare_soil,     &
                              totfrac_nobio, qsintmax,     co2_flux,       fco2_lu, temp_growth)


    netco2flux(:) = zero
    DO jv = 2,nvm
       netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
    ENDDO
    
    !! 1.4 Initialize diffusion coefficients
    CALL diffuco_initialize (kjit,    kjpindex, index,                  &
                             rest_id, lalo,     neighbours, resolution, &
                             rstruct, tq_cdrag)
    
    !! 1.5 Initialize remaining variables of energy budget
    CALL enerbil_initialize (kjit,     kjpindex,     index,    rest_id,  &
                             temp_air, qair,                             &
                             temp_sol, temp_sol_new, tsol_rad, t2mdiag,  &
                             evapot,   evapot_corr,  qsurf,    fluxsens, &
                             fluxlat,  vevapp )
    
    
    !! 1.7 Initialize remaining hydrological variables
    IF ( .NOT. hydrol_cwrr ) THEN
       ! 1.7.1 Initialize remaining hydrological variables from Choisnel module (2 soil layers)
       CALL hydrolc_initialize( kjit,      kjpindex,     index,          rest_id, &
                                veget,     veget_max,    tot_bare_soil,           &
                                rsol,      drysoil_frac, snow,                    &
                                snow_age,  snow_nobio,   snow_nobio_age, humrel,  &
                                vegstress, qsintveg,     shumdiag,       snowrho, &
                                snowtemp,  snowgrain,    snowdz,         snowheat)
     
       evap_bare_lim(:) = -un
       k_litt(:) = huit
       
       ! No specific calculation for shumdiag_perma. We assume it to shumdiag.
       shumdiag_perma(:,:)=shumdiag(:,:)
    ELSE
       !! 1.7.2 Initialize remaining hydrological variables from CWRR module (11 soil layers)
       CALL hydrol_initialize ( kjit,           kjpindex,  index,         rest_id,          &
            njsc,           soiltile,  veget,         veget_max,        &
            humrel,         vegstress, drysoil_frac,                    &
            shumdiag_perma, k_litt,    qsintveg,                        &
            evap_bare_lim,  snow,      snow_age,      snow_nobio,       &
            snow_nobio_age, snowrho,   snowtemp,      snowgrain,        &
            snowdz,         snowheat,  &
            mc_layh,    mcl_layh,  tmc_layh)

    ENDIF
    
    !! 1.9 Initialize surface parameters (emissivity, albedo and roughness)
    CALL condveg_initialize (kjit, kjpindex, index, rest_id, &
         lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
         zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
         drysoil_frac, height, snowdz,snowrho, tot_bare_soil, &
         emis, albedo, z0, roughheight, &
         frac_snow_veg,frac_snow_nobio)
    
    !! 1.10 Initialization of soil thermodynamics
    
    IF (hydrol_cwrr) THEN
       CALL thermosoil_initialize (kjit, kjpindex, rest_id,  &
            temp_sol_new, snow,       shumdiag_perma,        &
            soilcap,      soilflx,    stempdiag,             &
            gtemp,               &
            mc_layh,  mcl_layh,   tmc_layh,       njsc ,     &
            frac_snow_veg,frac_snow_nobio,totfrac_nobio,     &
            snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb)
    ELSE
       CALL thermosoilc_initialize (kjit, kjpindex, rest_id,  &
            temp_sol_new, snow,       shumdiag_perma,        &
            soilcap,      soilflx,    stempdiag,             &
            gtemp, &
            frac_snow_veg,frac_snow_nobio,totfrac_nobio,     &
            snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb)

    END IF

    !! 1.12 Initialize river routing
    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
       !! 1.12.1 Initialize river routing
       CALL routing_initialize( kjit,        kjpindex,       index,                 &
            rest_id,     hist_id,        hist2_id,   lalo,      &
            neighbours,  resolution,     contfrac,   stempdiag, &
            returnflow,  reinfiltration, irrigation, riverflow, &
            coastalflow, flood_frac,     flood_res , streamfl_frac, stream_frac ,fastr )
       DO ji=1,kjpindex
    !      flood_frac(ji) = MIN((flood_frac(ji) + streamfl_frac(ji) + stream_frac(ji)), un)
       ENDDO !ig = 1, nbpt
    ELSE
       !! 1.12.2 No routing, set variables to zero
       riverflow(:,:) = zero
       coastalflow(:,:) = zero
       returnflow(:,:) = zero
       reinfiltration(:,:) = zero
       irrigation(:,:) = zero
       flood_frac(:) = zero
       streamfl_frac(:) = zero
       stream_frac(:) = zero
       flood_res(:) = zero
       fastr(:) = zero
    ENDIF
    
    !! 1.13 Write internal variables to output fields
    z0_out(:) = z0(:)
    emis_out(:) = emis(:)
    albedo_out(:,:) = albedo(:,:) 
    qsurf_out(:) = qsurf(:)
    
  END SUBROUTINE sechiba_initialize

!! ==============================================================================================================================\n
!! SUBROUTINE 	: sechiba_main
!!
!>\BRIEF        Main routine for the sechiba module performing three functions:
!! calculating temporal evolution of all variables and preparation of output and 
!! restart files (during the last call only)
!!
!!\n DESCRIPTION : Main routine for the sechiba module. 
!! One time step evolution consists of:
!! - call sechiba_var_init to do some initialization,
!! - call slowproc_main to do some daily calculations
!! - call diffuco_main for diffusion coefficient calculation,
!! - call enerbil_main for energy budget calculation,
!! - call hydrolc_main (or hydrol_main) for hydrologic processes calculation,
!! - call enerbil_fusion : last part with fusion,
!! - call condveg_main for surface conditions such as roughness, albedo, and emmisivity,
!! - call thermosoil_main(for cwrr) or thermosoilc_main(for choisnel) for soil thermodynamic calculation,
!! - call sechiba_end to swap previous to new fields.
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): Hydrological variables (:: coastalflow and :: riverflow),
!! components of the energy budget (:: tsol_rad, :: vevapp, :: fluxsens, 
!! :: temp_sol_new and :: fluxlat), surface characteristics (:: z0_out, :: emis_out, 
!! :: tq_cdrag and :: albedo_out) and land use related CO2 fluxes (:: netco2flux and 
!! :: fco2_lu)            
!!
!! REFERENCE(S)	: 
!!
!! FLOWCHART    : 
!! \latexonly 
!! \includegraphics[scale = 0.5]{sechibamainflow.png}
!! \endlatexonly
!! \n
!_ ================================================================================================================================

  SUBROUTINE sechiba_main (kjit, kjpij, kjpindex, index, date0, &
       & ldrestart_read, ldrestart_write, &
       & lalo, contfrac, neighbours, resolution,&
       & zlev, u, v, qair, q2m, t2m, temp_air, epot_air, ccanopy, &
       & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
       & precip_rain, precip_snow, lwdown, swnet, swdown, coszang, pb, &
       & vevapp, fluxsens, fluxlat, coastalflow, riverflow, netco2flux, fco2_lu, &
       & tsol_rad, temp_sol_new, qsurf_out, albedo_out, emis_out, z0_out, &
       & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
!       & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, &
!       hist_id_stom_IPCC,soil_mc,litter_mc)
!! 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 
                                                                                  !! (unitless)
    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only 
                                                                                  !! (unitless)
    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
    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)                               :: rest_id_stom      !! STOMATE's _Restart_ file identifier 
                                                                                  !! (unitless)
    INTEGER(i_std),INTENT (in)                               :: hist_id_stom      !! STOMATE's _History_ file identifier 
                                                                                  !! (unitless)
    INTEGER(i_std),INTENT(in)                                :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file 
                                                                                  !! identifier (unitless)
    REAL(r_std), INTENT (in)                                 :: date0             !! Initial date (??unit??)
    LOGICAL, INTENT(in)                                      :: ldrestart_read    !! Logical for _restart_ file to read 
                                                                                  !! (true/false)
    LOGICAL, INTENT(in)                                      :: ldrestart_write   !! Logical for _restart_ file to write 
                                                                                  !! (true/false)
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)          :: lalo              !! Geographic coordinates (latitude,longitude)
                                                                                  !! for grid cells (degrees)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: contfrac          !! Fraction of continent in the grid 
                                                                                  !! (unitless, 0-1)
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map. 
                                                                                  !! Sechiba uses a reduced grid excluding oceans
                                                                                  !! ::index contains the indices of the 
                                                                                  !! terrestrial pixels only! (unitless)
    INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in)       :: neighbours        !! Neighboring grid points if land!(unitless)
    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)          :: resolution        !! Size in x and y of the grid (m)
    
    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)            :: zlev              !! Height of first layer (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair              !! Lowest level specific humidity 
                                                                                  !! @tex $(kg kg^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q2m               !! 2m specific humidity 
                                                                                  !! @tex $(kg kg^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: t2m               !! 2m air temperature (K)
    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)            :: coszang           !! Cosine of the solar zenith angle (unitless)
    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)            :: epot_air          !! Air potential energy (??J)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ccanopy           !! CO2 concentration in the canopy (ppm)
    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)


!! 0.2 Output variables

    REAL(r_std),DIMENSION (kjpindex,nflow), INTENT (out)     :: coastalflow       !! Outflow on coastal points by small basins.
                                                                                  !! This is the water which flows in a disperse 
                                                                                  !! way into the ocean
                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex,nflow), INTENT (out)     :: riverflow         !! Outflow of the major rivers.
                                                                                  !! The flux will be located on the continental 
                                                                                  !! grid but this should be a coastal point  
                                                                                  !! @tex $(kg dt_routing^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: tsol_rad          !! Radiative surface temperature 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: vevapp            !! Total of evaporation 
                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
    
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: qsurf_out         !! Surface specific humidity 
                                                                                  !! @tex $(kg kg^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: z0_out            !! Surface roughness (output diagnostic, m)
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out)         :: albedo_out        !! VIS and NIR albedo (output diagnostic, 
                                                                                  !! unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxsens          !! Sensible heat flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fluxlat           !! Latent heat flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: emis_out          !! Emissivity (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: netco2flux        !! Sum CO2 flux over PFTs 
                                                                                  !! ??(gC m^{-2} s^{-1})??
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: fco2_lu           !! Land Cover Change CO2 flux 
                                                                                  !! ??(gC m^{-2} s^{-1})??

!! 0.3 Modified

    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: tq_cdrag          !! Surface drag coefficient 
                                                                                  !! @tex $(m.s^{-1})$ @endtex 
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)         :: temp_sol_new      !! New ground temperature (K)
!    REAL(r_std),DIMENSION (kjpindex,nbdl,nstm), INTENT(inout):: soil_mc           !! soil moisture content \f($m^3 \times m^3$)\f
!    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT(inout)     :: litter_mc         !! litter moisture content \f($m^3 \times m^3$)\f

!! 0.4 local variables

    INTEGER(i_std)                                           :: ji, jv, iflow     !! Index (unitless)
    REAL(r_std), DIMENSION(kjpindex)                         :: histvar           !! Computations for history files (unitless)
    CHARACTER(LEN=80)                                        :: var_name          !! To store variables names for I/O (unitless)
    REAL(r_std), DIMENSION(kjpindex)                         :: sum_treefrac      !! Total fraction occupied by trees (0-1, uniless) 
    REAL(r_std), DIMENSION(kjpindex)                         :: sum_grassfrac     !! Total fraction occupied by grasses (0-1, unitless)
    REAL(r_std), DIMENSION(kjpindex)                         :: sum_cropfrac      !! Total fraction occcupied by crops (0-1, unitess)
    REAL(r_std), DIMENSION(kjpindex)                         :: grndflux          !! Net energy into soil (W/m2)
    REAL(r_std), DIMENSION(kjpindex,nsnow)                   :: snowliq           !! Liquid water content (m)
    REAL(r_std), DIMENSION(kjpindex)                         :: frac_snow_veg     !! Snow cover fraction on vegetation, 
                                                                                  !! only for diagnostics (unitless)
    REAL(r_std), DIMENSION(kjpindex,nnobio)                  :: frac_snow_nobio   !! Snow cover fraction on continental ice, lakes, etc
                                                                                  !! only for diagnostics (unitless)


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

    IF (printlev>=3) WRITE(numout,*) ' sechiba kjpindex =',kjpindex

    !! 1. Initialize variables at each time step
    CALL sechiba_var_init (kjpindex, rau, pb, temp_air) 

    !! 2. Compute diffusion coefficients
    CALL diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
         & zlev, z0, roughheight, temp_sol, temp_air, temp_growth, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb ,  &
         & rsol, evap_bare_lim, evapot, evapot_corr, snow, flood_frac, flood_res, frac_nobio, snow_nobio, totfrac_nobio, &
         & swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
         & vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, &
         & lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, &
         & hist_id, hist2_id)
   
    !! 3. Compute energy balance
    CALL enerbil_main (kjit, kjpindex, &
         & index, indexveg, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
         & qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, &
         & emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
         & vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, t2mdiag, temp_sol, tsol_rad, &
         & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id, &
         & precip_rain,  pgflux, snowdz, snowrho, temp_sol_add)

 
    !! 4. Compute hydrology
    IF ( .NOT. hydrol_cwrr ) THEN
       ! 4.1 Water balance from Choisnel module (2 soil layers)
       CALL hydrolc_main (kjit, kjpindex, index, indexveg, &
            & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
            & qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age, tot_melt, transpir, &
            & precip_rain, precip_snow, returnflow, reinfiltration, irrigation, humrel, vegstress, rsol, drysoil_frac, &
            & evapot, evapot_corr, flood_frac, flood_res, shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id, &
            & temp_air, pb, u, v, swnet, pgflux, &
            & snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
            & grndflux,gtemp,tot_bare_soil, &
            & lambda_snow,cgrnd_snow,dgrnd_snow,temp_sol_add)

       evap_bare_lim(:) = -un
       k_litt(:) = huit
       
       ! No specific calculation for shumdiag_perma. We assume it to shumdiag.
       shumdiag_perma(:,:)=shumdiag(:,:)
    ELSE
       !! 4.1 Water balance from CWRR module (11 soil layers)
       CALL hydrol_main (kjit, kjpindex, &
            & index, indexveg, indexsoil, indexlayer, indexnbdl, &
            & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
            & qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age,  &
            & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
            & humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, flood_frac, flood_res, &
            & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, reinf_slope,&
            & rest_id, hist_id, hist2_id,&
            & stempdiag, &
            & temp_air, pb, u, v, swnet, pgflux, &
            & snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
            & grndflux,gtemp,tot_bare_soil, &
            & lambda_snow,cgrnd_snow,dgrnd_snow,temp_sol_add, &
            & mc_layh, mcl_layh, tmc_layh, runoff_per_soil, drainage_per_soil, soil_mc, &
            & litter_mc, wat_flux0, wat_flux, precip2canopy, precip2ground, canopy2ground) 

       rsol(:) = -un

    ENDIF
         
    !! 5. Compute remaining components of the energy balance
    IF ( .NOT. ok_explicitsnow ) THEN
       CALL enerbil_fusion (kjpindex, tot_melt, soilcap, snowdz, &
            temp_sol_new, fusion)
    END IF

    !! 6. Compute surface variables (emissivity, albedo and roughness)
    CALL condveg_main (kjit, kjpindex, index, rest_id, hist_id, hist2_id, &
         lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
         zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
         drysoil_frac, height, snowdz, snowrho, tot_bare_soil, &
         emis, albedo, z0, roughheight, &
         frac_snow_veg, frac_snow_nobio)

    !! 7. Compute soil thermodynamics
    IF (hydrol_cwrr) THEN
       CALL thermosoil_main (kjit, kjpindex, &
            index, indexgrnd, &
            temp_sol_new, snow, soilcap, soilflx, &
            shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
            snowdz,snowrho,snowtemp,gtemp,pb,&
            mc_layh, mcl_layh, tmc_layh, njsc,frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
	    lambda_snow, cgrnd_snow, dgrnd_snow)
    ELSE
       CALL thermosoilc_main (kjit, kjpindex, &
            index, indexgrnd, indexnbdl, &
            temp_sol_new, snow, soilcap, soilflx, &
            shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, &
            snowdz,snowrho,snowtemp,gtemp,pb,&
            frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, &
	    lambda_snow, cgrnd_snow, dgrnd_snow)
    END IF

    !! 8. Compute slow processes (i.e. 'daily' and annual time step)
    ! ::ok_co2 and ::ok_stomate are flags that determine whether the
    ! forcing files are written.
    CALL slowproc_main (kjit, kjpij, kjpindex, date0, &
         index, indexveg, lalo, neighbours, resolution, contfrac, soiltile, &
         t2mdiag, t2mdiag, temp_sol, stempdiag, &
         vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
         deadleaf_cover, &
         assim_param, &
         lai, frac_age, height, veget, frac_nobio, njsc, veget_max, totfrac_nobio, qsintmax, &
         rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
         co2_flux, fco2_lu, temp_growth, tot_bare_soil, &
         soil_mc, litter_mc, floodout, runoff, drainage, wat_flux0, wat_flux, &
         drainage_per_soil, runoff_per_soil, DOC_EXP_agg, DOC_to_topsoil, DOC_to_subsoil, &
         (flood_frac+streamfl_frac), stream_frac, precip2canopy, precip2ground, canopy2ground, fastr)

    !! 9. Compute river routing 
    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
       !! 8.1 River routing
       CALL routing_main (kjit, kjpindex, index, &
            & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
            & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, stream_frac, streamfl_frac, &
            & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, &
            & rest_id, hist_id, hist2_id, DOC_EXP_agg, temp_sol, fastr)
    ELSE
       !! 8.2 No routing, set variables to zero
       riverflow(:,:) = zero
       coastalflow(:,:) = zero
       returnflow(:,:) = zero
       reinfiltration(:,:) = zero
       irrigation(:,:) = zero
       flood_frac(:) = zero
       stream_frac(:) = zero
       streamfl_frac(:) = zero
       flood_res(:) = zero
       fastr(:) = zero
    ENDIF
    
    ! Calculate DOC inputs to soil column from stomate_soilcarbon.f90, to be exported via slowprocesses
    DOC_to_topsoil(:,:) = (reinfiltration(:,:)+irrigation(:,:))*mille*one_day/dt_sechiba
    DOC_to_subsoil(:,:) = returnflow(:,:)*mille*one_day/dt_sechiba
        
    !! 9.2 Compute global CO2 flux
    !! Here,  far, netco2flux doesn't do anything (not an output any more)
    !! It could become an output again, and it is here where I could add the CO2 evasion from waters
    netco2flux(:) = zero
    DO jv = 2,nvm
       netco2flux(:) = netco2flux(:) + co2_flux(:,jv)*veget_max(:,jv)
    ENDDO
 
    !! 10. Update the temperature (temp_sol) with newly computed values
    CALL sechiba_end (kjpindex, temp_sol_new, temp_sol)

   
    !! 11. Write internal variables to output fields
    z0_out(:) = z0(:)
    emis_out(:) = emis(:)
    albedo_out(:,:) = albedo(:,:) 
    qsurf_out(:) = qsurf(:)
 
    !! 12. Write global variables to history files
    sum_treefrac(:) = zero
    sum_grassfrac(:) = zero
    sum_cropfrac(:) = zero
    DO jv = 2, nvm 
       IF (is_tree(jv) .AND. natural(jv)) THEN
          sum_treefrac(:) = sum_treefrac(:) + veget_max(:,jv)
       ELSE IF ((.NOT. is_tree(jv))  .AND. natural(jv)) THEN
          sum_grassfrac(:) = sum_grassfrac(:) + veget_max(:,jv)
       ELSE 
          sum_cropfrac = sum_cropfrac(:) + veget_max(:,jv)
       ENDIF
    ENDDO          

    CALL xios_orchidee_send_field("evapnu",vevapnu*one_day/dt_sechiba)
    CALL xios_orchidee_send_field("snow",snow)
    CALL xios_orchidee_send_field("flood_frac",flood_frac)
    CALL xios_orchidee_send_field("stream_frac",stream_frac)
    CALL xios_orchidee_send_field("streamfl_frac",streamfl_frac)
    CALL xios_orchidee_send_field("snowage",snow_age)
    CALL xios_orchidee_send_field("snownobio",snow_nobio)
    CALL xios_orchidee_send_field("snownobioage",snow_nobio_age)
    CALL xios_orchidee_send_field("reinf_slope",reinf_slope)
    CALL xios_orchidee_send_field("soilindex",REAL(njsc, r_std))
    CALL xios_orchidee_send_field("vegetfrac",veget)
    CALL xios_orchidee_send_field("maxvegetfrac",veget_max)
    CALL xios_orchidee_send_field("nobiofrac",frac_nobio)
    CALL xios_orchidee_send_field("soiltile",soiltile)
    CALL xios_orchidee_send_field("rstruct",rstruct)
    IF (ok_co2) CALL xios_orchidee_send_field("gpp",gpp/dt_sechiba)
    CALL xios_orchidee_send_field("nee",co2_flux/dt_sechiba)
    CALL xios_orchidee_send_field("drysoil_frac",drysoil_frac)
    CALL xios_orchidee_send_field("evapflo",vevapflo*one_day/dt_sechiba)
    CALL xios_orchidee_send_field("evapflo_alma",vevapflo/dt_sechiba)
    CALL xios_orchidee_send_field("k_litt",k_litt)
    CALL xios_orchidee_send_field("beta",vbeta)
    CALL xios_orchidee_send_field("vbeta1",vbeta1)
    CALL xios_orchidee_send_field("vbeta2",vbeta2)
    CALL xios_orchidee_send_field("vbeta3",vbeta3)
    CALL xios_orchidee_send_field("vbeta4",vbeta4)
    CALL xios_orchidee_send_field("vbeta5",vbeta5)
    CALL xios_orchidee_send_field("gsmean",gsmean)
    CALL xios_orchidee_send_field("cimean",cimean)
    CALL xios_orchidee_send_field("rveget",rveget)
    CALL xios_orchidee_send_field("rsol",rsol)
 
    histvar(:)=SUM(vevapwet(:,:),dim=2)
    CALL xios_orchidee_send_field("evspsblveg",histvar/dt_sechiba)
    histvar(:)= vevapnu(:)+vevapsno(:)
    CALL xios_orchidee_send_field("evspsblsoi",histvar/dt_sechiba)
    histvar(:)=SUM(transpir(:,:),dim=2)
    CALL xios_orchidee_send_field("tran",histvar/dt_sechiba)
    histvar(:)= sum_treefrac(:)*100*contfrac(:)
    CALL xios_orchidee_send_field("treeFrac",histvar)
    histvar(:)= sum_grassfrac(:)*100*contfrac(:)
    CALL xios_orchidee_send_field("grassFrac",histvar)
    histvar(:)= sum_cropfrac(:)*100*contfrac(:)
    CALL xios_orchidee_send_field("cropFrac",histvar)
    histvar(:)=veget_max(:,1)*100*contfrac(:)
    CALL xios_orchidee_send_field("baresoilFrac",histvar)
    histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)*100*contfrac(:)
    CALL xios_orchidee_send_field("residualFrac",histvar)

    CALL xios_orchidee_send_field("tsol_rad",tsol_rad-273.15)
    CALL xios_orchidee_send_field("qsurf",qsurf)
    CALL xios_orchidee_send_field("albedo",albedo)
    CALL xios_orchidee_send_field("emis",emis)
    CALL xios_orchidee_send_field("z0",z0)
    CALL xios_orchidee_send_field("roughheight",roughheight)
    CALL xios_orchidee_send_field("lai",lai)
    histvar(:)=zero   
    DO ji = 1, kjpindex
       IF (SUM(veget_max(ji,:)) > zero) THEN
         DO jv=2,nvm
            histvar(ji) = histvar(ji) + veget_max(ji,jv)*lai(ji,jv)/SUM(veget_max(ji,:))
         END DO
       END IF
    END DO

    CALL xios_orchidee_send_field("LAImean",histvar)
    CALL xios_orchidee_send_field("subli",vevapsno*one_day/dt_sechiba)
    CALL xios_orchidee_send_field("vevapnu",vevapnu*one_day/dt_sechiba)
    CALL xios_orchidee_send_field("vevapnu_alma",vevapnu/dt_sechiba)
    CALL xios_orchidee_send_field("transpir",transpir*one_day/dt_sechiba)
    CALL xios_orchidee_send_field("inter",vevapwet*one_day/dt_sechiba)
    CALL xios_orchidee_send_field("Qf",fusion)
    histvar(:)=zero
    DO jv=1,nvm
      histvar(:) = histvar(:) + vevapwet(:,jv)
    ENDDO
    CALL xios_orchidee_send_field("ECanop",histvar/dt_sechiba)
    histvar(:)=zero
    DO jv=1,nvm
      histvar(:) = histvar(:) + transpir(:,jv)
    ENDDO
    CALL xios_orchidee_send_field("TVeg",histvar/dt_sechiba)
    CALL xios_orchidee_send_field("ACond",tq_cdrag)

    IF ( .NOT. almaoutput ) THEN
       ! Write history file in IPSL-format
       CALL histwrite_p(hist_id, 'beta', kjit, vbeta, kjpindex, index)
       CALL histwrite_p(hist_id, 'z0', kjit, z0, kjpindex, index)
       CALL histwrite_p(hist_id, 'roughheight', kjit, roughheight, kjpindex, index)
       CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
       CALL histwrite_p(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'subli', kjit, vevapsno, kjpindex, index)
       CALL histwrite_p(hist_id, 'evapnu', kjit, vevapnu, kjpindex, index)
       CALL histwrite_p(hist_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
       CALL histwrite_p(hist_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'vbeta4', kjit, vbeta4, kjpindex, index)    
       CALL histwrite_p(hist_id, 'vbeta5', kjit, vbeta5, kjpindex, index)    
       CALL histwrite_p(hist_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
       CALL histwrite_p(hist_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)

       IF ( .NOT. hydrol_cwrr ) THEN
          CALL histwrite_p(hist_id, 'rsol', kjit, rsol, kjpindex, index)
       ENDIF
       CALL histwrite_p(hist_id, 'snow', kjit, snow, kjpindex, index)
       CALL histwrite_p(hist_id, 'snowage', kjit, snow_age, kjpindex, index)
       CALL histwrite_p(hist_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
       CALL histwrite_p(hist_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)

       IF (ok_explicitsnow) THEN
          CALL histwrite_p(hist_id, 'grndflux', kjit, grndflux, kjpindex,index)
          CALL histwrite_p(hist_id, 'snowtemp',kjit,snowtemp,kjpindex*nsnow,indexsnow)
          CALL histwrite_p(hist_id, 'snowliq', kjit,snowliq,kjpindex*nsnow,indexsnow)
          CALL histwrite_p(hist_id, 'snowdz', kjit,snowdz,kjpindex*nsnow,indexsnow)
          CALL histwrite_p(hist_id, 'snowrho', kjit,snowrho,kjpindex*nsnow,indexsnow)
          CALL histwrite_p(hist_id, 'snowgrain',kjit,snowgrain,kjpindex*nsnow,indexsnow)
          CALL histwrite_p(hist_id, 'snowheat',kjit,snowheat,kjpindex*nsnow,indexsnow)
       END IF

       CALL histwrite_p(hist_id,'frac_snow_veg',kjit,frac_snow_veg,kjpindex,index)
       CALL histwrite_p(hist_id, 'frac_snow_nobio', kjit,frac_snow_nobio,kjpindex*nnobio, indexnobio)
       CALL histwrite_p(hist_id, 'pgflux',kjit,pgflux,kjpindex,index)
       CALL histwrite_p(hist_id, 'soiltile',  kjit, soiltile, kjpindex*nstm, indexsoil)
       !
       IF ( hydrol_cwrr ) THEN
          CALL histwrite_p(hist_id, 'soilindex',  kjit, REAL(njsc, r_std), kjpindex, index)
          CALL histwrite_p(hist_id, 'reinf_slope',  kjit, reinf_slope, kjpindex, index)
          CALL histwrite_p(hist_id, 'k_litt', kjit, k_litt, kjpindex, index)
       ENDIF
       IF ( do_floodplains ) THEN
          CALL histwrite_p(hist_id, 'evapflo', kjit, vevapflo, kjpindex, index)
          CALL histwrite_p(hist_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
       ENDIF
       IF ( ok_co2 ) THEN
          CALL histwrite_p(hist_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)    
          CALL histwrite_p(hist_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)    
       ENDIF
       IF ( ok_stomate ) THEN
          CALL histwrite_p(hist_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)    
       ENDIF

       histvar(:)=SUM(vevapwet(:,:),dim=2)
       CALL histwrite_p(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)

       histvar(:)= vevapnu(:)+vevapsno(:)
       CALL histwrite_p(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)

       histvar(:)=SUM(transpir(:,:),dim=2)
       CALL histwrite_p(hist_id, 'tran', kjit, histvar, kjpindex, index)

       histvar(:)= sum_treefrac(:)*100*contfrac(:)
       CALL histwrite_p(hist_id, 'treeFrac', kjit, histvar, kjpindex, index) 

       histvar(:)= sum_grassfrac(:)*100*contfrac(:)
       CALL histwrite_p(hist_id, 'grassFrac', kjit, histvar, kjpindex, index) 

       histvar(:)= sum_cropfrac(:)*100*contfrac(:)
       CALL histwrite_p(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)

       histvar(:)=veget_max(:,1)*100*contfrac(:)
       CALL histwrite_p(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)

       histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)*100*contfrac(:)
       CALL histwrite_p(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
    ELSE
       ! Write history file in ALMA format 
       CALL histwrite_p(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
       CALL histwrite_p(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
       CALL histwrite_p(hist_id, 'Qf', kjit, fusion, kjpindex, index)
       CALL histwrite_p(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
       CALL histwrite_p(hist_id, 'EWater', kjit, vevapflo, kjpindex, index)
       CALL histwrite_p(hist_id, 'SWE', kjit, snow, kjpindex, index)
       histvar(:)=zero
       DO jv=1,nvm
          histvar(:) = histvar(:) + transpir(:,jv)
       ENDDO
       CALL histwrite_p(hist_id, 'TVeg', kjit, histvar, kjpindex, index)
       histvar(:)=zero
       DO jv=1,nvm
          histvar(:) = histvar(:) + vevapwet(:,jv)
       ENDDO
       CALL histwrite_p(hist_id, 'ECanop', kjit, histvar, kjpindex, index)
       CALL histwrite_p(hist_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
       CALL histwrite_p(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
       !
       CALL histwrite_p(hist_id, 'Z0', kjit, z0, kjpindex, index)
       CALL histwrite_p(hist_id, 'EffectHeight', kjit, roughheight, kjpindex, index)
       !
       IF ( do_floodplains ) THEN
          CALL histwrite_p(hist_id, 'Qflood', kjit, vevapflo, kjpindex, index)
          CALL histwrite_p(hist_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
       ENDIF
       !
       IF ( ok_co2 ) THEN
          CALL histwrite_p(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
       ENDIF
       IF ( ok_stomate ) THEN
             CALL histwrite_p(hist_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)    
       ENDIF
    ENDIF ! almaoutput
    
    !! 13. Write additional output file with higher frequency
    IF ( hist2_id > 0 ) THEN
       IF ( .NOT. almaoutput ) THEN
          ! Write history file in IPSL-format
          CALL histwrite_p(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
          CALL histwrite_p(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
          CALL histwrite_p(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
          CALL histwrite_p(hist2_id, 'emis', kjit, emis, kjpindex, index)
          CALL histwrite_p(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
          CALL histwrite_p(hist2_id, 'z0', kjit, z0, kjpindex, index)
          CALL histwrite_p(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
          CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
          CALL histwrite_p(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
          IF ( do_floodplains ) THEN
             CALL histwrite_p(hist2_id, 'vevapflo', kjit, vevapflo, kjpindex, index)
             CALL histwrite_p(hist2_id, 'flood_frac', kjit, flood_frac, kjpindex, index)
          ENDIF
          CALL histwrite_p(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
          CALL histwrite_p(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
          CALL histwrite_p(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)    
          CALL histwrite_p(hist2_id, 'vbeta5', kjit, vbeta5, kjpindex, index)    
          CALL histwrite_p(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
          CALL histwrite_p(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
          IF ( .NOT. hydrol_cwrr ) THEN
             CALL histwrite_p(hist2_id, 'rsol', kjit, rsol, kjpindex, index)
          ENDIF
          CALL histwrite_p(hist2_id, 'snow', kjit, snow, kjpindex, index)
          CALL histwrite_p(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
          CALL histwrite_p(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
          CALL histwrite_p(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
          !
          IF (  hydrol_cwrr ) THEN
             CALL histwrite_p(hist2_id, 'soilindex',  kjit, REAL(njsc, r_std), kjpindex, index)
             CALL histwrite_p(hist2_id, 'reinf_slope',  kjit, reinf_slope, kjpindex, index)
          ENDIF
          !
          IF ( ok_co2 ) THEN
             CALL histwrite_p(hist2_id, 'gsmean', kjit, gsmean, kjpindex*nvm, indexveg)    
             CALL histwrite_p(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
             CALL histwrite_p(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)    
          ENDIF
          IF ( ok_stomate ) THEN
             CALL histwrite_p(hist2_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)    
          ENDIF
       ELSE
          ! Write history file in ALMA format
          CALL histwrite_p(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
          CALL histwrite_p(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
          CALL histwrite_p(hist2_id, 'Qf', kjit, fusion, kjpindex, index)
          CALL histwrite_p(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
          IF ( do_floodplains ) THEN
             CALL histwrite_p(hist2_id, 'EWater', kjit, vevapflo, kjpindex, index)
             CALL histwrite_p(hist2_id, 'FloodFrac', kjit, flood_frac, kjpindex, index)
          ENDIF
          CALL histwrite_p(hist2_id, 'SWE', kjit, snow, kjpindex, index)
          histvar(:)=zero
          DO jv=1,nvm
             histvar(:) = histvar(:) + transpir(:,jv)
          ENDDO
          CALL histwrite_p(hist2_id, 'TVeg', kjit, histvar, kjpindex, index)
          histvar(:)=zero
          DO jv=1,nvm
             histvar(:) = histvar(:) + vevapwet(:,jv)
          ENDDO
          CALL histwrite_p(hist2_id, 'ECanop', kjit, histvar, kjpindex, index)
          CALL histwrite_p(hist2_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
          CALL histwrite_p(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
          IF ( ok_co2 ) THEN
             CALL histwrite_p(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
          ENDIF
          IF ( ok_stomate ) THEN
             CALL histwrite_p(hist2_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)    
          ENDIF
       ENDIF ! almaoutput
    ENDIF ! hist2_id


    !! Change the vegetation fractions if a new map was read in slowproc. This is done 
    !! after lcchange has been done in stomatelpj
    IF (done_stomate_lcchange) THEN
       CALL slowproc_change_frac(kjpindex, lai, &
                                 veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile)
       done_stomate_lcchange=.FALSE.
    END IF

    !! 14. If it is the last time step, write restart files
    IF (ldrestart_write) THEN
       CALL sechiba_finalize( &
            kjit,     kjpij,  kjpindex, index,   rest_id, &
            tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad)
    END IF

  END SUBROUTINE sechiba_main


!!  =============================================================================================================================
!! SUBROUTINE:    sechiba_finalize
!!
!>\BRIEF	  Finalize all modules by calling their "_finalize" subroutines.
!!
!! DESCRIPTION:	  Finalize all modules by calling their "_finalize" subroutines. These subroutines will write variables to 
!!                restart file. 
!!
!! \n
!_ ==============================================================================================================================

  SUBROUTINE sechiba_finalize( &
       kjit,     kjpij,  kjpindex, index,   rest_id, &
       tq_cdrag, vevapp, fluxsens, fluxlat, tsol_rad)

!! 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 
                                                                                  !! (unitless)
    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size - terrestrial pixels only 
                                                                                  !! (unitless)
    INTEGER(i_std),INTENT (in)                               :: rest_id           !! _Restart_ file identifier (unitless)
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: index             !! Indices of the pixels on the map. 
                                                                                  !! Sechiba uses a reduced grid excluding oceans
                                                                                  !! ::index contains the indices of the 
                                                                                  !! terrestrial pixels only! (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: tsol_rad           !! Radiative surface temperature 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: vevapp             !! Total of evaporation 
                                                                                  !! @tex $(kg m^{-2} days^{-1})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: fluxsens           !! Sensible heat flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: fluxlat            !! Latent heat flux 
                                                                                  !! @tex $(W m^{-2})$ @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)           :: tq_cdrag           !! Surface drag coefficient 
                                                                                  !! @tex $(m.s^{-1})$ @endtex 

!! 0.2 Local variables
    INTEGER(i_std)                                          :: ji, jv		  !! Index (unitless)
    REAL(r_std), DIMENSION(kjpindex)                        :: histvar            !! Computations for history files (unitless)
    CHARACTER(LEN=80)                                       :: var_name           !! To store variables names for I/O (unitless)


    !! Write restart file for the next simulation from SECHIBA and other modules

    IF (printlev>=3) WRITE (numout,*) 'Write restart file'

    !! 1. Call diffuco_finalize to write restart files
    CALL diffuco_finalize (kjit, kjpindex, rest_id, rstruct )
    
    !! 2. Call energy budget to write restart files
    CALL enerbil_finalize (kjit,   kjpindex,    rest_id,            &
                           evapot, evapot_corr, temp_sol, tsol_rad, &
                           qsurf,  fluxsens,    fluxlat,  vevapp )
    
    !! 3. Call hydrology to write restart files
    IF ( .NOT. hydrol_cwrr ) THEN
       !! 3.1 Call water balance from Choisnel module (2 soil layers) to write restart file
       CALL hydrolc_finalize( kjit,      kjpindex,   rest_id,        snow,       &
                              snow_age,  snow_nobio, snow_nobio_age, humrel,     &
                              vegstress, qsintveg,   snowrho,        snowtemp,   &
                              snowdz,    snowheat,   snowgrain,       &
                              drysoil_frac, rsol, shumdiag)

       evap_bare_lim(:) = -un
       k_litt(:) = huit
       shumdiag_perma(:,:)=shumdiag(:,:)
    ELSE
       !! 3.2 Call water balance from CWRR module (11 soil layers) to write restart file
       CALL hydrol_finalize( kjit,           kjpindex, rest_id,  vegstress,  &
                             qsintveg,       humrel,   snow,     snow_age, snow_nobio, &
                             snow_nobio_age, snowrho,  snowtemp, snowdz,     &
                             snowheat,       snowgrain,    &
                             drysoil_frac, evap_bare_lim)
    ENDIF
    
    !! 4. Call condveg to write surface variables to restart files
    CALL condveg_finalize (kjit, kjpindex, rest_id, z0, roughheight)
    
    !! 5. Call soil thermodynamic to write restart files
    IF (hydrol_cwrr) THEN
       CALL thermosoil_finalize (kjit,    kjpindex, rest_id,   gtemp, &
            soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
    ELSE
       CALL thermosoilc_finalize (kjit,    kjpindex, rest_id,   gtemp, &
            soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)
    
    END IF ! ldrestart_write 

    !! 6. Add river routing to restart files  
    IF ( river_routing .AND. nbp_glo .GT. 1) THEN
       !! 6.1 Call river routing to write restart files 
       CALL routing_finalize( kjit, kjpindex, rest_id, flood_frac, flood_res, streamfl_frac, fastr )
    ELSE
       !! 6.2 No routing, set variables to zero
       reinfiltration(:,:) = zero
       returnflow(:,:) = zero
       irrigation(:,:) = zero
       flood_frac(:) = zero
       streamfl_frac(:) = zero
       flood_res(:) = zero
       fastr(:) = zero
    ENDIF
    
    !! 7. Call slowproc_main to add 'daily' and annual variables to restart file
    CALL slowproc_finalize (kjit,       kjpindex,  rest_id,  index,      &
                            njsc,       lai,       height,   veget,      &
                            frac_nobio, veget_max, reinf_slope,          &
                            assim_param, frac_age, soiltile)
    
    IF (printlev>=3) WRITE (numout,*) 'sechiba_finalize done'
    
  END SUBROUTINE sechiba_finalize

  
!! ==============================================================================================================================\n
!! SUBROUTINE 	: sechiba_init
!!
!>\BRIEF        Dynamic allocation of the variables, the dimensions of the 
!! variables are determined by user-specified settings. 
!! 
!! DESCRIPTION  : The domain size (:: kjpindex) is used to allocate the correct
!! dimensions to all variables in sechiba. Depending on the variable, its 
!! dimensions are also determined by the number of PFT's (::nvm), number of 
!! soil types (::nstm), number of non-vegetative surface types (::nnobio),
!! number of soil levels (::ngrnd), number of soil layers in the hydrological 
!! model (i.e. cwrr) (::nslm). Values for these variables are set in
!! constantes_soil.f90 and constantes_veg.f90.\n
!!
!! Memory is allocated for all Sechiba variables and new indexing tables
!! are build making use of both (::kjpij) and (::kjpindex). New indexing tables 
!! are needed because a single pixel can contain several PFTs, soil types, etc.
!! The new indexing tables have separate indices for the different
!! PFTs, soil types, etc.\n
!!
!! RECENT CHANGE(S): None
!! 
!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output 
!! variables. However, the routine allocates memory and builds new indexing 
!! variables for later use.
!!
!! REFERENCE(S)	: None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================ 

  SUBROUTINE sechiba_init (kjit, kjpij, kjpindex, index, rest_id, lalo)

!! 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 (unitless)
    INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size - terrestrial pixels only (unitless)
    INTEGER(i_std), INTENT (in)                         :: rest_id            !! _Restart_ file identifier (unitless)
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: index              !! Indeces of the points on the map (unitless)
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)     :: lalo               !! Geographical coordinates (latitude,longitude) 
                                                                              !! for pixels (degrees)
!! 0.2 Output variables

!! 0.3 Modified variables

!! 0.4 Local variables

    INTEGER(i_std)                                      :: ier                !! Check errors in memory allocation (unitless)
    INTEGER(i_std)                                      :: ji, jv             !! Indeces (unitless)
!_ ==============================================================================================================================

!! 1. Initialize variables 
    
    ! Dynamic allocation with user-specified dimensions on first call
    IF (l_first_sechiba) THEN 
       l_first_sechiba=.FALSE.
    ELSE 
       CALL ipslerr_p(3,'sechiba_init',' l_first_sechiba false . we stop ','','')
    ENDIF

    !! Initialize local printlev
    printlev_loc=get_printlev('sechiba')
    

    !! 1.1 Initialize 3D vegetation indexation table
    ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexveg','','')

    ALLOCATE (indexlai(kjpindex*(nlai+1)),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlai','','')

    ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsoil','','')

    ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnobio','','')

    ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexgrnd','','')

    ALLOCATE (indexsnow(kjpindex*nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexsnow','','')

    ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexlayer','','')

    ALLOCATE (indexnbdl(kjpindex*nbdl),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexnbdl','','')

    ALLOCATE (indexalb(kjpindex*2),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for indexalb','','')

    !! 1.2  Initialize 1D array allocation with restartable value
    ALLOCATE (flood_res(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for flood_res','','')
    flood_res(:) = undef_sechiba

    ALLOCATE (fastr(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fastr','','')
    fastr(:) = undef_sechiba

    ALLOCATE (flood_frac(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for flood_frac','','')
    flood_frac(:) = undef_sechiba

    ALLOCATE (stream_frac(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stream_frac','','')
    stream_frac(:) = undef_sechiba

    ALLOCATE (streamfl_frac(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for streamfl_frac','','')
    streamfl_frac(:) = undef_sechiba

    ALLOCATE (snow(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow','','')
    snow(:) = undef_sechiba

    ALLOCATE (snow_age(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_age','','')
    snow_age(:) = undef_sechiba

    ALLOCATE (drysoil_frac(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drysoil_frac','','')

    ALLOCATE (rsol(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rsol','','')

    ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim','','')

    ALLOCATE (evapot(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot','','')
    evapot(:) = undef_sechiba

    ALLOCATE (evapot_corr(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evapot_corr','','')

    ALLOCATE (humrel(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for humrel','','')
    humrel(:,:) = undef_sechiba

    ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vegstress','','')
    vegstress(:,:) = undef_sechiba

    ALLOCATE (njsc(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for njsc','','')
    njsc(:)= undef_int

    ALLOCATE (soiltile(kjpindex,nstm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soiltile','','')

    ALLOCATE (reinf_slope(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','')

    ALLOCATE (vbeta1(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta1','','')

    ALLOCATE (vbeta4(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta4','','')

    ALLOCATE (vbeta5(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta5','','')

    ALLOCATE (soilcap(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilcap','','')

    ALLOCATE (soilflx(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for soilflx','','')

    ALLOCATE (temp_sol(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol','','')
    temp_sol(:) = undef_sechiba

    ALLOCATE (qsurf(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsurf','','')
    qsurf(:) = undef_sechiba

    !! 1.3 Initialize 2D array allocation with restartable value
    ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintveg','','')
    qsintveg(:,:) = undef_sechiba

    ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta2','','')

    ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3','','')

    ALLOCATE (vbeta3pot(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta3pot','','')

    ALLOCATE (gsmean(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gsmean','','')

    ALLOCATE (cimean(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cimean','','')

    ALLOCATE (gpp(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gpp','','')
    gpp(:,:) = undef_sechiba

 
    ALLOCATE (temp_growth(kjpindex),stat=ier) 
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_growth','','')
    temp_growth(:) = undef_sechiba 

    ALLOCATE (veget(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget','','')
    veget(:,:)=undef_sechiba

    ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for veget_max','','')

    ALLOCATE (tot_bare_soil(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_bare_soil','','')

    ALLOCATE (lai(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lai','','')
    lai(:,:)=undef_sechiba

    ALLOCATE (frac_age(kjpindex,nvm,nleafages),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_age','','')
    frac_age(:,:,:)=undef_sechiba

    ALLOCATE (height(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for height','','')
    height(:,:)=undef_sechiba

    ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for frac_nobio','','')
    frac_nobio(:,:) = undef_sechiba

    ALLOCATE (albedo(kjpindex,2),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for albedo','','')

    ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio','','')
    snow_nobio(:,:) = undef_sechiba

    ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snow_nobio_age','','')
    snow_nobio_age(:,:) = undef_sechiba

    ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for assim_param','','')

    !! 1.4 Initialize 1D array allocation 
    ALLOCATE (vevapflo(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapflo','','')
    vevapflo(:)=zero

    ALLOCATE (vevapsno(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapsno','','')

    ALLOCATE (vevapnu(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapnu','','')

    ALLOCATE (t2mdiag(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for t2mdiag','','')

    ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for totfrac_nobio','','')

    ALLOCATE (floodout(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for floodout','','')

    ALLOCATE (runoff(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for runoff','','')

   ALLOCATE (wat_flux0(kjpindex,nstm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in wat_flux0 allocation. We stop. We need kjpindex*nstm words = ',kjpindex*nstm
       STOP 'sechiba_init'
    END IF
 
    ALLOCATE (wat_flux(kjpindex,nslm, nstm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in wat_flux allocation. We stop. We need kjpindex*nslm*nstm words = ',kjpindex*nslm*nstm
       STOP 'sechiba_init'
    END IF

    ALLOCATE (drainage_per_soil(kjpindex,nstm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in drainage_per_soil allocation. We stop. We need kjpindex*nstm words = ',kjpindex*nstm
       STOP 'sechiba_init'
    END IF

    ALLOCATE (runoff_per_soil(kjpindex,nstm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in runoff_per_soil allocation. We stop. We need kjpindex*nstm words = ',kjpindex*nstm
       STOP 'sechiba_init'
    END IF

    ALLOCATE (litter_mc(kjpindex,nstm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in litter_mc allocation. We stop. We need kjpindex*nstm words = ',kjpindex*nstm
       STOP 'sechiba_init'
    END IF

    ALLOCATE (soil_mc(kjpindex,nbdl,nstm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in soil_mc allocation. We stop. We need kjpindex*nbdl*nstm words = ',kjpindex*nbdl*nstm
       STOP 'sechiba_init'
    END IF

    ALLOCATE (drainage(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for drainage','','')

    ALLOCATE (returnflow(kjpindex,nflow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for returnflow','','')
    returnflow(:,:) = zero
    
    ALLOCATE (DOC_EXP_agg(kjpindex,nexp,nflow),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in DOC_EXP_agg allocation. We stop. We need kjpindex*nexp*nflow words = ',kjpindex*nexp*nflow
       STOP 'sechiba_init'
    END IF
    DOC_EXP_agg(:,:,:) = zero    

    ALLOCATE (reinfiltration(kjpindex,nflow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinfiltration','','')
    reinfiltration(:,:) = zero

    ALLOCATE (irrigation(kjpindex,nflow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for irrigation','','')
    irrigation(:,:) = zero
    
    ALLOCATE (DOC_to_topsoil(kjpindex,nflow),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in DOC_to_topsoil allocation. We stop. We need kjpindex*nflow words = ',kjpindex*nflow
       STOP 'sechiba_init'
    END IF
    DOC_to_topsoil(:,:) = zero    
    
    ALLOCATE (DOC_to_subsoil(kjpindex,nflow),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in DOC_to_subsoil allocation. We stop. We need kjpindex*nflow words = ',kjpindex*nflow
       STOP 'sechiba_init'
    END IF
    DOC_to_subsoil(:,:) = zero

    ALLOCATE (precip2canopy(kjpindex,nvm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in precip2canopy allocation. We stop. We need kjpindex*nvm words = ',kjpindex*nvm
       STOP 'sechiba_init'
    END IF

    ALLOCATE (precip2ground(kjpindex,nvm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in precip2ground allocation. We stop. We need kjpindex*nvm words = ',kjpindex*nvm
       STOP 'sechiba_init'
    END IF

    ALLOCATE (canopy2ground(kjpindex,nvm),stat=ier)
    IF (ier.NE.0) THEN
       WRITE (numout,*) ' error in canopy2ground allocation. We stop. We need kjpindex*nvm words = ',kjpindex*nvm
       STOP 'sechiba_init'
    END IF
    
    ALLOCATE (z0(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for z0','','')

    ALLOCATE (roughheight(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for roughheight','','')

    ALLOCATE (emis(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for emis','','')

    ALLOCATE (tot_melt(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tot_melt','','')

    ALLOCATE (valpha(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for valpha','','')

    ALLOCATE (vbeta(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vbeta','','')

    ALLOCATE (fusion(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for fusion','','')

    ALLOCATE (rau(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rau','','')

    ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for deadleaf_cover','','')

    ALLOCATE (stempdiag(kjpindex, nbdl),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for stempdiag','','')


    ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for co2_flux','','')
    co2_flux(:,:)=zero

    ALLOCATE (shumdiag(kjpindex,nbdl),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag','','')
    
    ALLOCATE (shumdiag_perma(kjpindex,nbdl),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag_perma','','')

    ALLOCATE (litterhumdiag(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for litterhumdiag','','')

    ALLOCATE (ptnlev1(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ptnlev1','','')

    ALLOCATE (k_litt(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for k_litt','','')

    !! 1.5 Initialize 2D array allocation
    ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for vevapwet','','')
    vevapwet(:,:)=undef_sechiba

    ALLOCATE (transpir(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpir','','')

    ALLOCATE (transpot(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for transpot','','')

    ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for qsintmax','','')

    ALLOCATE (rveget(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rveget','','')

    ALLOCATE (rstruct(kjpindex,nvm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for rstruct','','')

    ALLOCATE (pgflux(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for pgflux','','')
    pgflux(:)= 0.0

    ALLOCATE (cgrnd_snow(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for cgrnd_snow','','')
    cgrnd_snow(:,:) = 0

    ALLOCATE (dgrnd_snow(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for dgrnd_snow','','')
    dgrnd_snow(:,:) = 0

    ALLOCATE (lambda_snow(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for lambda_snow','','')
    lambda_snow(:) = 0

    ALLOCATE (temp_sol_add(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for temp_sol_add','','')

    ALLOCATE (gtemp(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for gtemp','','')

    ALLOCATE (snowrho(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowrho','','')

    ALLOCATE (snowheat(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowheat','','')

    ALLOCATE (snowgrain(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowgrain','','')

    ALLOCATE (snowtemp(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowtemp','','')

    ALLOCATE (snowdz(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for snowdz','','')

    ALLOCATE (mc_layh(kjpindex, nslm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mc_layh','','')

    ALLOCATE (mcl_layh(kjpindex, nslm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcl_layh','','')

    ALLOCATE (tmc_layh(kjpindex, nslm),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for tmc_layh','','')

    !! 1.6 Initialize indexing table for the vegetation fields. 
    ! In SECHIBA we work on reduced grids but to store in the full 3D filed vegetation variable 
    ! we need another index table : indexveg, indexsoil, indexnobio and indexgrnd
    DO ji = 1, kjpindex
       !
       DO jv = 1, nlai+1
          indexlai((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
       ENDDO
       !
       DO jv = 1, nvm
          indexveg((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
       ENDDO
       !      
       DO jv = 1, nstm
          indexsoil((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
       ENDDO
       !      
       DO jv = 1, nnobio
          indexnobio((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
       ENDDO
       !
       DO jv = 1, ngrnd
          indexgrnd((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
       ENDDO
       !
       DO jv = 1, nsnow
          indexsnow((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
       ENDDO

       DO jv = 1, nbdl
          indexnbdl((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij
       ENDDO

       DO jv = 1, nslm
          indexlayer((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
       ENDDO
       !
       DO jv = 1, 2
          indexalb((jv-1)*kjpindex + ji) = INDEX(ji) + (jv-1)*kjpij + offset_omp - offset_mpi
       ENDDO
       !
    ENDDO

!! 2. Read the default value that will be put into variable which are not in the restart file
    CALL ioget_expval(val_exp)
    
    IF (printlev>=3) WRITE (numout,*) ' sechiba_init done '

  END SUBROUTINE sechiba_init
  

!! ==============================================================================================================================\n
!! SUBROUTINE 	: sechiba_clear
!!
!>\BRIEF        Deallocate memory of sechiba's variables
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): None 
!!
!! REFERENCE(S)	: None
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================ 

  SUBROUTINE sechiba_clear (forcing_name,cforcing_name)

    CHARACTER(LEN=100), INTENT(in)           :: forcing_name       !! Name of forcing file (unitless)
    CHARACTER(LEN=100), INTENT(in)           :: cforcing_name      !! Name of forcing file with carbon related variables (unitless)
!_ ================================================================================================================================
    
!! 1. Initialize first run

    l_first_sechiba=.TRUE.

!! 2. Deallocate dynamic variables of sechiba

    IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
    IF ( ALLOCATED (indexlai)) DEALLOCATE (indexlai)
    IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
    IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
    IF ( ALLOCATED (indexsnow)) DEALLOCATE (indexsnow)
    IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
    IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
    IF ( ALLOCATED (indexnbdl)) DEALLOCATE (indexnbdl)
    IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
    IF ( ALLOCATED (flood_res)) DEALLOCATE (flood_res)
    IF ( ALLOCATED (fastr)) DEALLOCATE (fastr)
    IF ( ALLOCATED (flood_frac)) DEALLOCATE (flood_frac)
    IF ( ALLOCATED (streamfl_frac)) DEALLOCATE (streamfl_frac)
    IF ( ALLOCATED (snow)) DEALLOCATE (snow)
    IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
    IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
    IF ( ALLOCATED (rsol)) DEALLOCATE (rsol)
    IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
    IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
    IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
    IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
    IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
    IF ( ALLOCATED (soiltile)) DEALLOCATE (soiltile)
    IF ( ALLOCATED (njsc)) DEALLOCATE (njsc)
    IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope)
    IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
    IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
    IF ( ALLOCATED (vbeta5)) DEALLOCATE (vbeta5)
    IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
    IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
    IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
    IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
    IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
    IF ( ALLOCATED (vbeta2))  DEALLOCATE (vbeta2)
    IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
    IF ( ALLOCATED (vbeta3pot)) DEALLOCATE (vbeta3pot)
    IF ( ALLOCATED (gsmean)) DEALLOCATE (gsmean)
    IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
    IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
    IF ( ALLOCATED (temp_growth)) DEALLOCATE (temp_growth) 
    IF ( ALLOCATED (veget)) DEALLOCATE (veget)
    IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
    IF ( ALLOCATED (tot_bare_soil)) DEALLOCATE (tot_bare_soil)
    IF ( ALLOCATED (lai)) DEALLOCATE (lai)
    IF ( ALLOCATED (frac_age)) DEALLOCATE (frac_age)
    IF ( ALLOCATED (height)) DEALLOCATE (height)
    IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
    IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
    IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
    IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
    IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
    IF ( ALLOCATED (vevapflo)) DEALLOCATE (vevapflo)
    IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
    IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
    IF ( ALLOCATED (t2mdiag)) DEALLOCATE (t2mdiag)
    IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
    IF ( ALLOCATED (floodout)) DEALLOCATE (floodout)
    IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
    IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
    IF ( ALLOCATED (reinfiltration)) DEALLOCATE (reinfiltration)
    IF ( ALLOCATED (DOC_EXP_agg)) DEALLOCATE (DOC_EXP_agg) 
    IF ( ALLOCATED (DOC_to_topsoil)) DEALLOCATE (DOC_to_topsoil) 
    IF ( ALLOCATED (DOC_to_subsoil)) DEALLOCATE (DOC_to_subsoil) 
    IF ( ALLOCATED (precip2canopy)) DEALLOCATE (precip2canopy)
    IF ( ALLOCATED (precip2ground)) DEALLOCATE (precip2ground)
    IF ( ALLOCATED (canopy2ground)) DEALLOCATE (canopy2ground) 
    IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
    IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
    IF ( ALLOCATED (valpha)) DEALLOCATE (valpha)
    IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
    IF ( ALLOCATED (fusion)) DEALLOCATE (fusion)
    IF ( ALLOCATED (rau)) DEALLOCATE (rau)
    IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
    IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
    IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
    IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
    IF ( ALLOCATED (shumdiag_perma)) DEALLOCATE (shumdiag_perma)
    IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
    IF ( ALLOCATED (ptnlev1)) DEALLOCATE (ptnlev1)
    IF ( ALLOCATED (k_litt)) DEALLOCATE (k_litt)
    IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
    IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
    IF ( ALLOCATED (transpot)) DEALLOCATE (transpot)
    IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
    IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
    IF ( ALLOCATED (rstruct)) DEALLOCATE (rstruct)
    IF ( ALLOCATED (snowrho)) DEALLOCATE (snowrho)
    IF ( ALLOCATED (snowgrain)) DEALLOCATE (snowgrain)
    IF ( ALLOCATED (snowtemp)) DEALLOCATE (snowtemp)
    IF ( ALLOCATED (snowdz)) DEALLOCATE (snowdz)
    IF ( ALLOCATED (snowheat)) DEALLOCATE (snowheat)
    IF ( ALLOCATED (cgrnd_snow)) DEALLOCATE (cgrnd_snow)
    IF ( ALLOCATED (dgrnd_snow)) DEALLOCATE (dgrnd_snow)
    IF ( ALLOCATED (lambda_snow)) DEALLOCATE(lambda_snow) 
    IF ( ALLOCATED (gtemp)) DEALLOCATE (gtemp)
    IF ( ALLOCATED (pgflux)) DEALLOCATE (pgflux)
    IF ( ALLOCATED (mc_layh)) DEALLOCATE (mc_layh)
    IF ( ALLOCATED (mcl_layh)) DEALLOCATE (mcl_layh)
    IF ( ALLOCATED (tmc_layh)) DEALLOCATE (tmc_layh)

!! 3. Clear all allocated memory

    CALL pft_parameters_clear
    CALL slowproc_clear 
    CALL diffuco_clear 
    CALL enerbil_clear  
    IF ( hydrol_cwrr ) THEN
       CALL hydrol_clear 
       CALL thermosoil_clear
    ELSE
       CALL hydrolc_clear  
       CALL thermosoilc_clear
    ENDIF
    CALL condveg_clear 
    CALL routing_clear

  END SUBROUTINE sechiba_clear


!! ==============================================================================================================================\n
!! SUBROUTINE 	: sechiba_var_init
!!
!>\BRIEF        Calculate air density as a function of air temperature and 
!! pressure for each terrestrial pixel.
!! 
!! RECENT CHANGE(S): None
!! 
!! MAIN OUTPUT VARIABLE(S): air density (::rau, kg m^{-3}).
!! 
!! REFERENCE(S)	: None
!! 
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air) 

!! 0.1 Input variables

    INTEGER(i_std), INTENT (in)                    :: kjpindex        !! Domain size - terrestrial pixels only (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: pb              !! Surface pressure (hPa)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: temp_air        !! Air temperature (K)
    
!! 0.2 Output variables

    REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau             !! Air density @tex $(kg m^{-3})$ @endtex

!! 0.3 Modified variables

!! 0.4 Local variables

    INTEGER(i_std)                                 :: ji              !! Indices (unitless)
!_ ================================================================================================================================
    
!! 1. Calculate intial air density (::rau)
   
    DO ji = 1,kjpindex
       rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
    END DO

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

  END SUBROUTINE sechiba_var_init


!! ==============================================================================================================================\n
!! SUBROUTINE 	: sechiba_end
!!
!>\BRIEF        Swap old for newly calculated soil temperature.
!! 
!! RECENT CHANGE(S): None
!! 
!! MAIN OUTPUT VARIABLE(S): soil temperature (::temp_sol; K)
!! 
!! REFERENCE(S)	: None
!! 
!! FLOWCHART    : None
!! \n
!! ================================================================================================================================ 

  SUBROUTINE sechiba_end (kjpindex, temp_sol_new, temp_sol)
                         

!! 0.1 Input variables

    INTEGER(i_std), INTENT (in)                       :: kjpindex           !! Domain size - terrestrial pixels only (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: temp_sol_new       !! New soil temperature (K)
    
    !! 0.2 Output variables
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)    :: temp_sol           !! Soil temperature (K)

!_ ================================================================================================================================
    
!! 1. Swap temperature

    temp_sol(:) = temp_sol_new(:)
    
    IF (printlev>=3) WRITE (numout,*) ' sechiba_end done '

  END SUBROUTINE sechiba_end

END MODULE sechiba