! =================================================================================================================================
! MODULE       : thermosoilc
!
! CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
! LICENCE      : IPSL (2006)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF        Calculates the soil temperatures by solving the heat
!! diffusion equation within the soil. This module is only used with Choisnel hydrology.
!!
!!\n DESCRIPTION : General important informations about the numerical scheme and
!!                 the soil vertical discretization:\n
!!               - the soil is divided into "ngrnd" (=7 by default) layers, reaching to as
!!                 deep as 5.5m down within the soil, with thiscknesses
!!                 following a geometric series of ration 2.\n
!!               - "jg" is usually used as the index going from 1 to ngrnd to describe the
!!                  layers, from top (jg=1) to bottom (jg=ngrnd)\n
!!               - the thermal numerical scheme is implicit finite differences.\n
!!                 -- When it is resolved in thermosoil_profile at the present timestep t, the
!!                 dependancy from the previous timestep (t-1) is hidden in the
!!                 integration coefficients cgrnd and dgrnd, which are therefore
!!                 calculated at the very end of thermosoil_main (call to
!!                 thermosoil_coef) for use in the next timestep.\n
!!                 -- At timestep t, the system becomes :\n 
!!
!!                              T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
!!                                      -- EQ1 -- \n
!!
!!                 (the bottom boundary condition has been used to obtained this equation).\n
!!                 To solve it, the uppermost soil temperature T(1) is required.
!!                 It is obtained from the surface temperature Ts, which is
!!                 considered a linear extrapolation of T(1) and T(2)\n
!!
!!                           Ts=(1-lambda)*T(1) -lambda*T(2) \n 
!!                                      -- EQ2--\n
!!
!!                 -- caveat 1 : Ts is called 'temp_soil_new' in this routine,
!!                 don' t act.\n
!!                 -- caveat 2 : actually, the surface temperature at time t Ts
!!                 depends on the soil temperature at time t through the
!!                 ground heat flux. This is again implicitly solved, with Ts(t)
!!                 expressed as :\n
!!
!!                 soilcap*(Ts(t)-Ts(t-1))/dt=soilflux+otherfluxes(Ts(t))\n 
!!                                      -- EQ3 --\n
!!
!!                 and the dependency from the previous timestep is hidden in
!!                 soilcap and soilflux (apparent surface heat capacity and heat
!!                 flux respectively). Soilcap and soilflux are therefore
!!                 calculated at the previsou timestep, at the very end of thermosoil
!!                 (final call to thermosoil_coef) and stored to be used at the next time step.
!!                 At timestep t, EQ3 is solved for Ts in enerbil, and Ts
!!                 is used in thermosoil to get T(1) and solve EQ1.\n
!!
!! - lambda is the @tex $\mu$ @endtex of F. Hourdin' s PhD thesis, equation (A28); ie the
!! coefficient of the linear extrapolation of Ts (surface temperature) from T1 and T2 (ptn(jg=1) and ptn(jg=2)), so that:\n
!! Ts= (1+lambda)*T(1)-lambda*T(2) --EQ2-- \n
!! lambda = (zz_coeff(1))/((zz_coef(2)-zz_coef(1))) \n
!!
!! - cstgrnd is the attenuation depth of the diurnal temperature signal
!! (period : one_day) as a result of the heat conduction equation
!! with no coefficients :
!!\latexonly
!!\input{thermosoil_var_init0.tex}
!!\endlatexonly
!!  -- EQ4 --\n
!! This equation results from the change of variables :
!! z' =z*sqrt(Cp/K) where z' is the new depth (homogeneous
!! to sqrt(time) ), z the real depth (in m), Cp and K the soil heat
!! capacity and conductivity respectively.\n
!!
!! the attenuation depth of a diurnal thermal signal for EQ4 is therefore homogeneous to sqrt(time) and
!! equals : \n
!! cstgrnd = sqrt(oneday/Pi)
!!
!! - lskin is the attenuation depth of the diurnal temperature signal
!! (period : one_day) within the soil for the complete heat conduction equation
!! (ie : with coefficients)
!!\latexonly
!!\input{thermosoil_var_init00.tex}
!!\endlatexonly
!! -- EQ5 --  \n
!! it can be retrieved from cstgrnd using the change of variable  z' =z*sqrt(Cp/K):\n
!! lskin = sqrt(K/Cp)*cstgrnd =  sqrt(K/Cp)*sqrt(oneday//Pi)\n
!! 
!! In thermosoil, the ratio lskin/cstgrnd is frequently used as the
!! multiplicative factor to go from
!!'adimensional' depths (like z' ) to real depths (z). z' is not really
!! adimensional but is reffered to like this in the code.
!!
!!
!! RECENT CHANGE(S) : thermosoilc module is a copy of thermosoil before the vertical discretization was changed.
!!                    This module is used only for the Choisnel hydrology. 
!!
!! REFERENCE(S) : None
!!
!! SVN          :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/tags/ORCHIDEE_2_0/ORCHIDEE/src_sechiba/thermosoilc.f90 $
!! $Date: 2018-03-14 11:33:02 +0100 (mer. 14 mars 2018) $
!! $Revision: 5096 $
!! \n
!_ ================================================================================================================================

MODULE thermosoilc

  ! modules used :
  USE ioipsl
  USE ioipsl_para
  USE xios_orchidee
  USE constantes
  USE time, ONLY : one_day, dt_sechiba
  USE constantes_soil
  USE sechiba_io_p
  USE grid


  IMPLICIT NONE

  !private and public routines :
  PRIVATE
  PUBLIC :: thermosoilc_main, thermosoilc_clear, thermosoilc_levels, thermosoilc_initialize, thermosoilc_finalize

  REAL(r_std), SAVE                               :: lambda, cstgrnd, lskin   !! See Module description
!$OMP THREADPRIVATE(lambda, cstgrnd, lskin)

  REAL(r_std), SAVE                               :: fz1, zalph               !! usefull constants for diverse use
!$OMP THREADPRIVATE(fz1, zalph)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn                      !! vertically discretized 
                                                                              !! soil temperatures @tex ($K$) @endtex. 
!$OMP THREADPRIVATE(ptn)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: zz                       !! depths of the soil thermal numerical nodes. 
                                                                              !! Caveats: they are not exactly the centers of the
                                                                              !! thermal layers, see the calculation in 
                                                                              !! ::thermosoilc_var_init  @tex ($m$) @endtex.
!$OMP THREADPRIVATE(zz)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: zz_coef		      !! depths of the boundaries of the thermal layers,
                                                                              !! see the calculation in 
                                                                              !! thermosoilc_var_init  @tex ($m$) @endtex.
!$OMP THREADPRIVATE(zz_coef)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: dz1                      !! numerical constant used in the thermal numerical
                                                                              !! scheme  @tex ($m^{-1}$) @endtex. ; it corresponds
                                                                              !! to the coefficient  @tex $d_k$ @endtex of equation
                                                                              !! (A.12) in F. Hourdin PhD thesis.
!$OMP THREADPRIVATE(dz1)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: dz2                      !! thicknesses of the thermal layers  @tex ($m$)
                                                                              !! @endtex; typically: 
                                                                              !! dz2(jg)=zz_coef(jg+1)-zz_coef(jg); calculated once 
                                                                              !! and for all in thermosoilc_var_init
!$OMP THREADPRIVATE(dz2)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd                    !! integration coefficient for the numerical scheme,
                                                                              !! see eq.1
!$OMP THREADPRIVATE(cgrnd)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd                    !! integration coefficient for the numerical scheme,
                                                                              !! see eq.1
!$OMP THREADPRIVATE(dgrnd)

  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa                    !! volumetric vertically discretized soil heat 
                                                                              !! capacity  @tex ($J K^{-1} m^{-3}$) @endtex. 
                                                                              !! It depends on the soil
                                                                              !! moisture content (shum_ngrnd_perma) and is calculated at 
                                                                              !! each time step in thermosoilc_coef.
!$OMP THREADPRIVATE(pcapa)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa                   !! vertically discretized soil thermal conductivity 
                                                                              !!  @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa.
!$OMP THREADPRIVATE(pkappa)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_snow               !! volumetric vertically discretized snow heat 
                                                                              !! capacity @tex ($J K^{-1} m^{-3}$) @endtex.
!$OMP THREADPRIVATE(pcapa_snow)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa_snow              !! vertically discretized snow thermal conductivity 
                                                                              !! @tex ($W K^{-1} m^{-1}$) @endtex.
!$OMP THREADPRIVATE(pkappa_snow)

  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en                 !! heat capacity used for surfheat_incr and 
                                                                              !! coldcont_incr 
!$OMP THREADPRIVATE(pcapa_en)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg                  !! vertically discretized temperature at the 
                                                                              !! beginning of the time step  @tex ($K$) @endtex; 
                                                                              !! is used in 
                                                                              !! thermosoilc_energy for energy-related diagnostic of
                                                                              !! the routine.
!$OMP THREADPRIVATE(ptn_beg)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: temp_sol_beg             !! Surface temperature at the beginning of the 
                                                                              !! timestep  @tex ($K$) @endtex
!$OMP THREADPRIVATE(temp_sol_beg)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: surfheat_incr            !! Change in soil heat content during the timestep 
                                                                              !!  @tex ($J$) @endtex.
!$OMP THREADPRIVATE(surfheat_incr)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: coldcont_incr            !! Change in snow heat content  @tex ($J$) @endtex.
!$OMP THREADPRIVATE(coldcont_incr)
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shum_ngrnd_perma         !! Saturation degree on the thermal axes (0-1, dimensionless)
!$OMP THREADPRIVATE(shum_ngrnd_perma)

  !  Variables related to soil freezing
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: profil_froz              !! Frozen fraction of the soil on hydrological levels (-)
!$OMP THREADPRIVATE(profil_froz)
  REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat                  !! Accumulated latent heat for the whole soil (J)
!$OMP THREADPRIVATE(e_soil_lat)
  REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp               !! Additional heat capacity due to soil freezing for each soil layer (J/K)
!$OMP THREADPRIVATE(pcappa_supp)    


CONTAINS
  !!  =============================================================================================================================
  !! SUBROUTINE		 		    : thermosoilc_initialize
  !!
  !>\BRIEF			            Allocate module variables, read from restart file or initialize with default values
  !!
  !! DESCRIPTION			    : Allocate module variables, read from restart file or initialize with default values.
  !!                                          Call thermosoilc_var_init to calculate physical constants. 
  !!                                          Call thermosoilc_coef to calculate thermal soil properties.
  !!
  !! RECENT CHANGE(S)			    : None
  !!
  !! REFERENCE(S)			    : None
  !! 
  !! FLOWCHART                              : None
  !! \n
  !_ ==============================================================================================================================
  SUBROUTINE 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)

    !! 0. Variable and parameter declaration
    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                            :: kjit             !! Time step number (unitless) 
    INTEGER(i_std), INTENT(in)                            :: kjpindex         !! Domain size (unitless)
    INTEGER(i_std),INTENT (in)                            :: rest_id          !! Restart file identifier (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: temp_sol_new     !! Surface temperature at the present time-step,
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: snow             !! Snow mass (kg)
    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in)    :: shumdiag_perma   !! Soil saturation degree on the diagnostic axis (0-1, unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: frac_snow_veg    !! Snow cover fraction on vegeted area
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in)   :: frac_snow_nobio  !! Snow cover fraction on non-vegeted area
    REAL(r_std),DIMENSION (kjpindex),INTENT(in)           :: totfrac_nobio    !! Total fraction of continental ice+lakes+cities+...
                                                                              !! (unitless,0-1)
    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in)    :: snowdz           !! Snow depth
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)    :: snowrho         !! Snow density
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp        !! Snow temperature
    REAL(r_std), DIMENSION (kjpindex), INTENT (in)         :: pb              !! Surface presure (hPa)


    !! 0.2 Output variables
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)        :: soilcap          !! apparent surface heat capacity (J m-2 K-1)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)        :: soilflx          !! apparent soil heat flux (W m-2)
    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)   :: stempdiag        !! diagnostic temperature profile (K)
    REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: gtemp            !! First soil layer temperature

    !! 0.3 Modified variables
    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)       :: lambda_snow     !! Coefficient of the linear extrapolation of surface temperature 
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow      !! Integration coefficient for snow numerical scheme
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow      !! Integration coefficient for snow numerical scheme

    !! 0.4 Local variables
    INTEGER(i_std)                                        :: ier, i
    LOGICAL                                               :: calculate_coef   !! Local flag to initialize variables by call to thermosoilc_coef   
!_ ================================================================================================================================
    
    IF (printlev >= 3) WRITE (numout,*) 'Start thermosoilc_initialize '

    !! 1. Allocate soil temperatures variables
    ALLOCATE (ptn(kjpindex,ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of ptn','','')

    ALLOCATE (zz(ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of zz','','')

    ALLOCATE (zz_coef(ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of zz_coef','','')

    ALLOCATE (dz1(ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of dz1','','')

    ALLOCATE (dz2(ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of dz2','','')

    ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of cgrnd','','')

    ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of dgrnd','','')

    ALLOCATE (pkappa_snow(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pkappa_snow','','')

    ALLOCATE (pcapa_snow(kjpindex,nsnow),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pcapa_snow','','')

    ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pcapa','','')

    ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pkappa','','')

    ! Temporary fix: Initialize following variable because they are output to xios before the first calculation
    pcapa  = 0
    pkappa = 0
    pcapa_snow  = 0
    pkappa_snow = 0

    ALLOCATE (surfheat_incr(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of surfheat_incr','','')

    ALLOCATE (coldcont_incr(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of coldcont_incr','','')

    ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of pcapa_en','','')

    ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of ptn_beg','','')

    ALLOCATE (temp_sol_beg(kjpindex),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of temp_sol_beg','','')

    ALLOCATE (shum_ngrnd_perma(kjpindex,ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of shum_ngrnd_perma','','')

    ALLOCATE (profil_froz(kjpindex,ngrnd),stat=ier)
    IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of profil_froz','','')

    IF (ok_freeze_thermix) THEN
       ALLOCATE (pcappa_supp(kjpindex,ngrnd),stat=ier)
       IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of ok_freeze_termix','','')
    END IF
    IF (ok_Ecorr) THEN
       ALLOCATE (e_soil_lat(kjpindex),stat=ier)
       IF (ier /= 0) CALL ipslerr_p(3,'thermosoilc_initialize', 'Error in allocation of e_soil_lat','','')
    END IF

    !! 2. Initialize variable from restart file or with default values 
    
    !! Reads restart files for soil temperatures only. If no restart file is
    !! found,  the initial soil temperature is by default set to 280K at all depths. The user
    !! can decide to initialize soil temperatures at an other value, in which case he should set the flag THERMOSOIL_TPRO
    !! to this specific value in the run.def.
    IF (printlev>=3) WRITE (numout,*) 'Read restart file for THERMOSOIL variables'

    CALL ioconf_setatt_p('UNITS', 'K')
    CALL ioconf_setatt_p('LONG_NAME','Soil Temperature profile')
    CALL restget_p (rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g)

    ! Initialize ptn if it was not found in restart file
    IF (ALL(ptn(:,:)==val_exp)) THEN 
       ! ptn was not found in restart file

       IF (read_reftemp) THEN
          ! Read variable ptn from file
          CALL thermosoilc_read_reftempfile(kjpindex,lalo,ptn)
       ELSE
          ! Initialize ptn with a constant value which can be set in run.def
          
          !Config Key   = THERMOSOIL_TPRO
          !Config Desc  = Initial soil temperature profile if not found in restart
          !Config Def   = 280.
          !Config If    = OK_SECHIBA
          !Config Help  = The initial value of the temperature profile in the soil if 
          !Config         its value is not found in the restart file. Here
          !Config         we only require one value as we will assume a constant 
          !Config         throughout the column.
          !Config Units = Kelvin [K]
          CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std)
       END IF
    END IF
    
    ! Initialize ptn_beg (variable needed in thermosoilc_coef before calucation in thermosoilc_energy)
    ptn_beg(:,:) = ptn(:,:)
    
    ! Initialize temp_sol_beg with values from previous time-step
    temp_sol_beg(:) = temp_sol_new(:) 
    
    ! Read e_soil_lat from restart file or initialize
    IF (ok_Ecorr) THEN
       CALL restget_p (rest_id, 'e_soil_lat', nbp_glo, 1, 1, kjit, .TRUE., &
            e_soil_lat, "gather", nbp_glo, index_g)
       CALL setvar_p (e_soil_lat, val_exp,'NO_KEYWORD',zero)
    END IF

    ! Read gtemp from restart file
    CALL restget_p (rest_id, 'gtemp', nbp_glo, 1, 1, kjit, .TRUE., &
         gtemp, "gather", nbp_glo, index_g)
    CALL setvar_p (gtemp, val_exp,'NO_KEYWORD',zero)
    
    ! Read variables calculated in thermosoilc_coef from restart file
    ! If the variables were not found in the restart file, the logical 
    ! calculate_coef will be true and thermosoilc_coef will be called further below.
    ! These variables need to be in the restart file to avoid a time shift that
    ! would be done using thermosoilc_coef at this stage.
    calculate_coef=.FALSE.
    CALL ioconf_setatt_p('UNITS', 'J m-2 K-1')
    CALL ioconf_setatt_p('LONG_NAME','Apparent surface heat capacity')
    CALL restget_p (rest_id, 'soilcap', nbp_glo, 1, 1, kjit, .TRUE., &
         soilcap, "gather", nbp_glo, index_g)
    IF (ALL(soilcap(:)==val_exp)) calculate_coef=.TRUE.

    CALL ioconf_setatt_p('UNITS', 'W m-2')
    CALL ioconf_setatt_p('LONG_NAME','Apparent soil heat flux')
    CALL restget_p (rest_id, 'soilflx', nbp_glo, 1, 1, kjit, .TRUE., &
         soilflx, "gather", nbp_glo, index_g)
    IF (ALL(soilflx(:)==val_exp)) calculate_coef=.TRUE.

    CALL ioconf_setatt_p('UNITS', '')
    CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
    CALL restget_p (rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
         cgrnd, "gather", nbp_glo, index_g)
    IF (ALL(cgrnd(:,:)==val_exp)) calculate_coef=.TRUE.

    CALL ioconf_setatt_p('UNITS', '')
    CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
    CALL restget_p (rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., &
         dgrnd, "gather", nbp_glo, index_g)
    IF (ALL(dgrnd(:,:)==val_exp)) calculate_coef=.TRUE.

    CALL ioconf_setatt_p('UNITS', '')
    CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
    CALL restget_p (rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
         cgrnd_snow, "gather", nbp_glo, index_g)
    IF (ALL(cgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.

    CALL ioconf_setatt_p('UNITS', '')
    CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme')
    CALL restget_p (rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., &
         dgrnd_snow, "gather", nbp_glo, index_g)
    IF (ALL(dgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE.

    CALL ioconf_setatt_p('UNITS', '')
    CALL ioconf_setatt_p('LONG_NAME','Coefficient of the linear extrapolation of surface temperature')
    CALL restget_p (rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, .TRUE., &
         lambda_snow, "gather", nbp_glo, index_g)
    IF (ALL(lambda_snow(:)==val_exp)) calculate_coef=.TRUE.


    !! 2.2.Computes physical constants and arrays; initializes soil thermal properties; produces the first stempdiag
    !!  Computes some physical constants and arrays depending on the soil vertical discretization 
    !! (lskin, cstgrnd, zz, zz_coef, dz1, dz2); get the vertical humidity onto the thermal levels
    CALL thermosoilc_var_init (kjpindex, shumdiag_perma, snow, &
         zz, zz_coef, dz1, dz2, stempdiag)
    
    !! 2.3. Computes cgrd, dgrd, soilflx and soilcap coefficients from restart values or initialisation values.
    !!      This is done only if they were not found in restart file.
    IF (calculate_coef) THEN
       IF (printlev>=3) WRITE (numout,*) 'thermosoilc_coef will be called in the intialization phase'
       CALL thermosoilc_coef (&
            kjpindex, temp_sol_new, snow, &
            ptn,      soilcap, soilflx,      zz, &
            dz1,      dz2,           &
            cgrnd,    dgrnd, &
            frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
            snowdz,snowrho,  snowtemp,     pb,         &
	    lambda_snow,    cgrnd_snow,    dgrnd_snow)
    END IF

  END SUBROUTINE thermosoilc_initialize


!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_main
!!
!>\BRIEF        Thermosoilc_main computes the soil thermal properties and dynamics, ie solves
!! the heat diffusion equation within the soil. The soil temperature profile is
!! then interpolated onto the diagnostic axis.
!!
!! DESCRIPTION : The resolution of the soil heat diffusion equation 
!! relies on a numerical finite-difference implicit scheme
!! fully described in the reference and in the header of the thermosoilc module.
!! - The dependency of the previous timestep hidden in the 
!! integration coefficients cgrnd and dgrnd (EQ1), calculated in thermosoilc_coef, and 
!! called at the end of the routine to prepare for the next timestep.
!! - The effective computation of the new soil temperatures is performed in thermosoilc_profile. 
!!
!! - thermosoilc_coef calculates the coefficients for the numerical scheme for the very first iteration of thermosoilc;
!! after that, thermosoilc_coef is called only at the end of the module to calculate the coefficients for the next timestep.
!! - thermosoilc_profile solves the numerical scheme.\n
!!
!! - Flags : one unique flag : THERMOSOIL_TPRO (to be set to the desired initial soil in-depth temperature in K; by default 280K)
!!
!! RECENT CHANGE(S) : None
!!
!! MAIN OUTPUT VARIABLE(S): vertically discretized soil temperatures ptn, soil
!! thermal properties (pcapa, pkappa), apparent surface heat capacity (soilcap)
!! and heat flux (soilflux) to be used in enerbil at the next timestep to solve
!! the surface energy balance.
!!
!! REFERENCE(S) : 
!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
!!  Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin' s PhD thesis relative to the thermal
!!  integration scheme has been scanned and is provided along with the documentation, with name : 
!!  Hourdin_1992_PhD_thermal_scheme.pdf
!!
!! FLOWCHART    : 
!! \latexonly
!! \includegraphics[scale = 1]{thermosoil_flowchart.png}
!! \endlatexonly
!! 
!! \n
!_ ================================================================================================================================

  SUBROUTINE thermosoilc_main (kjit, kjpindex, &
       index, indexgrnd, &
       indexnslm, &
       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)

    !! 0. Variable and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                            :: kjit             !! Time step number (unitless) 
    INTEGER(i_std), INTENT(in)                            :: kjpindex         !! Domain size (unitless)
    INTEGER(i_std),INTENT (in)                            :: rest_id,hist_id  !! Restart_ file and history file identifier 
                                                                              !! (unitless)
    INTEGER(i_std),INTENT (in)                            :: hist2_id         !! history file 2 identifier (unitless)
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)      :: index            !! Indeces of the points on the map (unitless)
    INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd        !! Indeces of the points on the 3D map (vertical 
                                                                              !! dimension towards the ground) (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)      :: temp_sol_new     !! Surface temperature at the present time-step,
                                                                              !! temp_sol_new is only modified for the case ok_explicitsnow
                                                                              !! Ts @tex ($K$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: snow             !! Snow mass @tex ($kg$) @endtex.
                                                                              !! Caveat: when there is snow on the
                                                                              !! ground, the snow is integrated into the soil for
                                                                              !! the calculation of the thermal dynamics. It means
                                                                              !! that the uppermost soil layers can completely or 
                                                                              !! partially consist in snow. In the second case, zx1
                                                                              !! and zx2 are the fraction of the soil layer 
                                                                              !! consisting in snow and 'normal' soil, respectively
                                                                              !! This is calculated in thermosoilc_coef.
    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in)    :: shumdiag_perma   !! Soil saturation degree on the diagnostic axis (0-1, unitless)
    INTEGER(i_std),DIMENSION (kjpindex*nslm), INTENT (in) :: indexnslm        !! Indeces of the points on the 3D map
    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in)    :: snowdz           !! Snow depth
    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in)    :: snowrho          !! Snow density
    REAL(r_std), DIMENSION (kjpindex), INTENT (in)        :: pb               !! Surface presure (hPa)

    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (inout) :: snowtemp         !! Snow temperature
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)          :: frac_snow_veg    !! Snow cover fraction on vegeted area
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in)   :: frac_snow_nobio  !! Snow cover fraction on non-vegeted area
    REAL(r_std),DIMENSION (kjpindex),INTENT(in)           :: totfrac_nobio    !! Total fraction of continental ice+lakes+cities+...
                                                                              !!(unitless,0-1)
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)      :: temp_sol_add     !! additional surface temperature due to the melt of first layer
                                                                              !!at the present time-step @tex ($K$) @endtex

    !! 0.2 Output variables

    REAL(r_std),DIMENSION (kjpindex), INTENT (out)        :: ptnlev1          !! 1st level soil temperature   
    REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: gtemp            !! First soil layer temperature


    !! 0.3 Modified variables

    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)      :: soilcap          !! apparent surface heat capacity
                                                                              !! @tex ($J m^{-2} K^{-1}$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)      :: soilflx          !! apparent soil heat flux @tex ($W m^{-2}$) @endtex
                                                                              !! , positive 
                                                                              !! towards the soil, writen as Qg (ground heat flux) 
                                                                              !! in the history files, and computed at the end of 
                                                                              !! thermosoilc for the calculation of Ts in enerbil, 
                                                                              !! see EQ3.
    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (inout) :: stempdiag        !! diagnostic temperature profile @tex ($K$) @endtex
                                                                              !! , eg on the 
                                                                              !! diagnostic axis (levels:1:nslm). The soil 
                                                                              !! temperature is put on this diagnostic axis to be
                                                                              !! used by other modules (slowproc.f90; routing.f90;
                                                                              !! hydrol or hydrolc when a frozen soil 
                                                                              !! parametrization is used..)
    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)       :: lambda_snow     !! Coefficient of the linear extrapolation of surface temperature 
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow      !! Integration coefficient for snow numerical scheme
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow      !! Integration coefficient for snow numerical scheme

    !! 0.4 Local variables

    INTEGER(i_std)                                        :: jv,ji,ii
    REAL(r_std),DIMENSION (kjpindex)                      :: snowtemp_weighted!! Snow temperature weighted by snow density (K)

!_ ================================================================================================================================
    
  !! 3. Put the soil wetness diagnostic on the levels of the soil temperature

    !!?? this could logically be put just before the last call to
    !!thermosoil_coef, as the results are used there...
    CALL thermosoilc_humlev(kjpindex, shumdiag_perma, snow)

    
  !! 4. Effective computation of the soil temperatures profile, using the cgrd and dgrd coefficients from previsou tstep.
    
    CALL thermosoilc_profile (kjpindex, temp_sol_new, ptn,  &
         stempdiag,snowtemp, frac_snow_veg, frac_snow_nobio,totfrac_nobio, &
         cgrnd_snow, dgrnd_snow)

  !! 5. Call to thermosoilc_energy, still to be clarified..

    CALL thermosoilc_energy (kjpindex, temp_sol_new, soilcap)

  !! 6. Writing the history files according to the ALMA standards (or not..)

    !in only one file (hist2_id <=0) or in 2 different files (hist2_id >0).
    CALL xios_orchidee_send_field("ptn",ptn)
    CALL xios_orchidee_send_field("soilflx",soilflx)
    CALL xios_orchidee_send_field("surfheat_incr",surfheat_incr)
    CALL xios_orchidee_send_field("coldcont_incr",coldcont_incr)
    CALL xios_orchidee_send_field("pkappa",pkappa)
    CALL xios_orchidee_send_field("pkappa_snow",pkappa_snow)
    CALL xios_orchidee_send_field("pcapa",pcapa)
    CALL xios_orchidee_send_field("pcapa_snow",pcapa_snow)
    CALL xios_orchidee_send_field("snowtemp",snowtemp)

    IF (ok_explicitsnow) THEN
       DO ji=1,kjpindex 
          ! Use min_sechiba instead of zero to avoid problem with division by zero
          IF (snow(ji) .GT. min_sechiba) THEN
             snowtemp_weighted(ji) = SUM(snowtemp(ji,:)*snowrho(ji,:))/SUM(snowrho(ji,:))
          ELSE
             snowtemp_weighted(ji) = xios_default_val
          END IF
       END DO
       CALL xios_orchidee_send_field("snowtemp_weighted",snowtemp_weighted)
    END IF

    IF ( .NOT. almaoutput ) THEN
      CALL histwrite_p(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
      CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
    ELSE
      CALL histwrite_p(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
      CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index)
      CALL histwrite_p(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
      CALL histwrite_p(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
    ENDIF
    IF ( hist2_id > 0 ) THEN
       IF ( .NOT. almaoutput ) THEN
          CALL histwrite_p(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd)
       ELSE
          CALL histwrite_p(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd)
          CALL histwrite_p(hist2_id, 'Qg', kjit, soilflx, kjpindex, index)
          CALL histwrite_p(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index)
          CALL histwrite_p(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index)
       ENDIF
    ENDIF
    
  !! 7. A last final call to thermosoilc_coef
 
    !! A last final call to thermosoilc_coef, which calculates the different
    !!coefficients (cgrnd, dgrnd, dz1, soilcap, soilflx) from this time step to be
    !!used at the next time step, either in the surface temperature calculation
    !!(soilcap, soilflx) or in the soil thermal numerical scheme.
    CALL thermosoilc_coef (&
         kjpindex, temp_sol_new, snow, &
         ptn,      soilcap, soilflx,      zz,   &
         dz1,      dz2,                     &
         cgrnd,   dgrnd, &
         frac_snow_veg,frac_snow_nobio,totfrac_nobio,&
         snowdz,snowrho,  snowtemp,     pb, &
	 lambda_snow,    cgrnd_snow,    dgrnd_snow)

    ! Save variables for explicit snow model
    gtemp(:) = ptn(:,1)

    !! Initialize output arguments to be used in sechiba
    ptnlev1(:) = ptn(:,1)


    !! Surface temperature is forced to zero celcius if its value is larger than melting point, only for explicit snow scheme
    IF  (ok_explicitsnow) THEN
       DO ji=1,kjpindex
          IF  (SUM(snowdz(ji,:)) .GT. 0.0) THEN
             IF (temp_sol_new(ji) .GE. tp_00) THEN
                temp_sol_new(ji) = tp_00
             ENDIF
          END IF
       END DO
    ENDIF


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

  END SUBROUTINE thermosoilc_main

  !!  =============================================================================================================================
  !! SUBROUTINE		 		    : thermosoilc_finalize
  !!
  !>\BRIEF                                    Write to restart file
  !!
  !! DESCRIPTION			    : This subroutine writes the module variables and variables calculated in thermosoilc
  !!                                          to restart file
  !! \n
  !_ ==============================================================================================================================
  SUBROUTINE thermosoilc_finalize (kjit,    kjpindex, rest_id,   gtemp,     &
                                  soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow)

    !! 0. Variable and parameter declaration
    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                            :: kjit             !! Time step number (unitless) 
    INTEGER(i_std), INTENT(in)                            :: kjpindex         !! Domain size (unitless)
    INTEGER(i_std),INTENT (in)                            :: rest_id          !! Restart file identifier(unitless)
    REAL(r_std),DIMENSION (kjpindex),INTENT(in)           :: gtemp            !! First soil layer temperature
    REAL(r_std),DIMENSION (kjpindex),INTENT(in)           :: soilcap
    REAL(r_std),DIMENSION (kjpindex),INTENT(in)           :: soilflx
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)         :: lambda_snow      !! Coefficient of the linear extrapolation of surface temperature 
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in)  :: cgrnd_snow       !! Integration coefficient for snow numerical scheme
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in)  :: dgrnd_snow       !! Integration coefficient for snow numerical scheme

!_ ================================================================================================================================
    
    !! 1. Write variables to restart file to be used for the next simulation
    IF (printlev>=3) WRITE (numout,*) 'Write restart file with THERMOSOIL variables'
    
    CALL restput_p(rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g)
    
    IF (ok_Ecorr) THEN
       CALL restput_p(rest_id, 'e_soil_lat', nbp_glo, 1 , 1, kjit, e_soil_lat, 'scatter', nbp_glo, index_g)
    END IF
    
    CALL restput_p(rest_id, 'gtemp', nbp_glo, 1, 1, kjit, gtemp, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'soilcap', nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'soilflx', nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, cgrnd_snow, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, dgrnd_snow, 'scatter', nbp_glo, index_g)
    CALL restput_p(rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, lambda_snow, 'scatter', nbp_glo, index_g)

  END SUBROUTINE thermosoilc_finalize


!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_clear
!!
!>\BRIEF        Desallocates the allocated arrays.
!! The call of thermosoilc_clear originates from sechiba_clear but the calling sequence and 
!! its purpose require further investigation.
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S) : None 
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) : None 
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

 SUBROUTINE thermosoilc_clear()

        IF ( ALLOCATED (ptn)) DEALLOCATE (ptn)
        IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd) 
        IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd) 
        IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa)
        IF ( ALLOCATED (pkappa))  DEALLOCATE (pkappa)
        IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en)
        IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg)
        IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg)
        IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr)
        IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr)
        IF ( ALLOCATED (shum_ngrnd_perma)) DEALLOCATE (shum_ngrnd_perma)
        IF ( ALLOCATED (profil_froz)) DEALLOCATE (profil_froz)

  END SUBROUTINE thermosoilc_clear


!! ================================================================================================================================
!! FUNCTION     : fz
!!
!>\BRIEF        fz(rk), the function's result, is the "rk"th element of a geometric series 
!! with first element fz1 and ration zalph.
!!
!! DESCRIPTION  : This function is used to calculate the depths of the boudaries of the thermal layers (zz_coef) and 
!! of the numerical nodes (zz) of the thermal scheme. Formulae to get the adimensional depths are followings :
!!      zz(jg)      = fz(REAL(jg,r_std) - undemi); \n
!!      zz_coef(jg) = fz(REAL(jg,r_std))
!!
!! RECENT CHANGE(S) : None
!!
!! RETURN VALUE : fz(rk)
!!
!! REFERENCE(S) : None 
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  FUNCTION fz(rk) RESULT (fz_result)

  !! 0. Variables and parameter declaration

    !! 0.1 Input variables

    REAL(r_std), INTENT(in)                        :: rk
    
    !! 0.2 Output variables

    REAL(r_std)                                    :: fz_result
    
    !! 0.3 Modified variables

    !! 0.4 Local variables

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

    fz_result = fz1 * (zalph ** rk - un) / (zalph - un)

  END FUNCTION fz


!! ================================================================================================================================
!! FUNCTION     : thermosoilc_levels
!!
!>\BRIEF          Depth of nodes for the thermal layers in meters. 
!!
!! DESCRIPTION  : Calculate and return the depth in meters of the nodes of the soil layers. This calculation is the same 
!!                as done in thermosoilc_var_init for zz. See thermosoilc_var_init for more details. 
!!
!! RECENT CHANGE(S) : None
!!
!! RETURN VALUE : Vector of soil depth for the nodes in meters
!!
!! REFERENCE(S) : None 
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  FUNCTION thermosoilc_levels() RESULT (zz_out)
    
    !! 0.1 Return variable
    
    REAL(r_std), DIMENSION (ngrnd)  :: zz_out      !! Depth of soil layers in meters
    
    !! 0.2 Local variables
    INTEGER(i_std)                  :: jg
    REAL(r_std)                     :: so_capa
    REAL(r_std)                     :: so_cond
    
!_ ================================================================================================================================

    !! 1. Define some parameters
    so_capa = (so_capa_dry + so_capa_wet)/deux
    so_cond = (so_cond_dry + so_cond_wet)/deux
    
    cstgrnd=SQRT(one_day / pi)
    lskin = SQRT(so_cond / so_capa * one_day / pi)

    !! Parameters needed by fz function
    fz1 = 0.3_r_std * cstgrnd
    zalph = deux

    !!  2. Get adimentional depth of the numerical nodes
    DO jg=1,ngrnd
       zz_out(jg) = fz(REAL(jg,r_std) - undemi)
    ENDDO
    
    !! 3. Convert to meters
    DO jg=1,ngrnd
       zz_out(jg) = zz_out(jg) /  cstgrnd * lskin
    END DO

  END FUNCTION thermosoilc_levels


!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_var_init
!!
!>\BRIEF        Define and initializes the soil thermal parameters
!!		  
!! DESCRIPTION	: This routine\n
!! 1. Defines the parameters ruling the vertical grid of the thermal scheme (fz1, zalpha).\n
!! 2. Defines the scaling coefficients for adimensional depths (lskin, cstgrnd, see explanation in the 
!!    variables description of thermosoilc_main). \n
!! 3. Calculates the vertical discretization of the soil (zz, zz_coef, dz2) and the constants used
!!    in the numerical scheme and which depend only on the discretization (dz1, lambda).\n
!! 4. Initializes the soil thermal parameters (capacity, conductivity) based on initial soil moisture content.\n
!! 5. Produces a first temperature diagnostic based on temperature initialization.\n
!!
!! The scheme comprizes ngrnd=7 layers by default.
!! The layer' s boundaries depths (zz_coef) follow a geometric series of ratio zalph=2 and first term fz1.\n
!! zz_coef(jg)=fz1.(1-zalph^jg)/(1-zalph) \n
!! The layers' boudaries depths are first calculated 'adimensionally', ie with a
!! discretization adapted to EQ5. This discretization is chosen for its ability at
!! reproducing a thermal signal with periods ranging from days to centuries. (see
!! Hourdin, 1992). Typically, fz1 is chosen as : fz1=0.3*cstgrnd (with cstgrnd the
!! adimensional attenuation depth). \n
!! The factor lskin/cstgrnd is then used to go from adimensional depths to
!! depths in m.\n
!! zz(real)=lskin/cstgrnd*zz(adimensional)\n
!! Similarly, the depths of the numerical nodes are first calculated
!! adimensionally, then the conversion factor is applied.\n
!! the numerical nodes (zz) are not exactly the layers' centers : their depths are calculated as follows:\n
!! zz(jg)=fz1.(1-zalph^(jg-1/2))/(1-zalph)\n
!! The values of zz and zz_coef used in the default thermal discretization are in the following table.
!! \latexonly
!! \includegraphics{thermosoilc_var_init1.jpg}
!! \endlatexonly\n
!!
!! RECENT CHANGE(S) : None
!!
!! MAIN OUTPUT VARIABLE(S) : None
!!
!! REFERENCE(S)	:
!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of 
!! planetary atmospheres, Ph.D. thesis, Paris VII University.
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE thermosoilc_var_init(kjpindex, shumdiag_perma, snow, &
                                 zz, zz_coef, dz1, dz2, stempdiag)

  !! 0. Variables and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                               :: kjpindex          !! Domain size (unitless)
    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (in)      :: shumdiag_perma    !! Relative soil humidity on the diagnostic axis 
                                                                                  !! (unitless), [0,1]. (see description of the 
                                                                                  !! variables of thermosoilc_main for more 
                                                                                  !! explanations) 
    REAL(r_std), DIMENSION (kjpindex), INTENT (in)     	     :: snow              !! Snow quantity   
   
    !! 0.2 Output variables

    REAL(r_std), DIMENSION (ngrnd), INTENT(out)              :: zz                !! depths of the layers'numerical nodes 
                                                                                  !! @tex ($m$)@endtex 
    REAL(r_std), DIMENSION (ngrnd), INTENT(out)              :: zz_coef		  !! depths of the layers'boundaries 
                                                                                  !! @tex ($m$)@endtex 
    REAL(r_std), DIMENSION (ngrnd), INTENT(out)              :: dz1               !! numerical constant depending on the vertical
                                                                                  !! thermal grid only @tex  ($m^{-1}$) @endtex. 
                                                                                  !! (see description
                                                                                  !! of the variables of thermosoilc_main for more
                                                                                  !! explanations)
    REAL(r_std), DIMENSION (ngrnd), INTENT(out)              :: dz2               !! thicknesses of the soil thermal layers 
                                                                                  !! @tex ($m$) @endtex
    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: stempdiag         !! Diagnostic temperature profile @tex ($K$)
                                                                                  !! @endtex

    !! 0.3 Modified variables

    !! 0.4 Local variables

    INTEGER(i_std)                                           :: ier, ji, jg       !! Index (unitless)
    REAL(r_std)                                              :: sum               !! Temporary variable
    REAL(r_std)                                              :: so_capa           !! Average Thermal Conductivity of soils 
                                                                                  !! @tex $(W.m^{-2}.K^{-1})$ @endtex
    REAL(r_std)                                              :: so_cond           !! Average Thermal Conductivity of soils 
                                                                                  !! @tex $(W.m^{-2}.K^{-1})$ @endtex

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

  !! 1. Initialization of the parameters of the vertical discretization and of the attenuation depths
    
    !! so_capa and so_cond are temporary variables which contain the average values of soil conductivity 
    !! and soil conductivity and are only needed in thermosoilc_var_init to set the vertical layering. 
    so_capa = (so_capa_dry + so_capa_wet)/deux
    so_cond = (so_cond_dry + so_cond_wet)/deux

    cstgrnd=SQRT(one_day / pi)
    lskin = SQRT(so_cond / so_capa * one_day / pi)
    fz1 = 0.3_r_std * cstgrnd
    zalph = deux

    !! Calculate so_capa_ice
    so_capa_ice = so_capa_dry + poros*capa_ice*rho_ice
    WRITE(numout,*) 'Calculation of so_capa_ice=', so_capa_ice,' using poros=',poros,' and capa_ice=',capa_ice
    
  !! 2.  Computing the depth of the thermal levels (numerical nodes) and the layers boundaries
   
    !! Computing the depth of the thermal levels (numerical nodes) and 
    !! the layers boundariesusing the so-called
    !! adimentional variable z' = z/lskin*cstgrnd (with z in m)
    
    !! 2.1 adimensional thicknesses of the layers
    DO jg=1,ngrnd

    !!?? code simplification hopefully possible here with up-to-date compilers !
    !!! This needs to be solved soon. Either we allow CPP options in SECHIBA or the VPP
    !!! fixes its compiler 
    !!!#ifdef VPP5000
      dz2(jg) = fz(REAL(jg,r_std)-undemi+undemi) - fz(REAL(jg-1,r_std)-undemi+undemi)
    !!!#else
    !!!      dz2(jg) = fz(REAL(jg,r_std)) - fz(REAL(jg-1,r_std))
    !!!#endif
    ENDDO
    
    !! 2.2 adimentional depth of the numerical nodes and layers' boudaries
    DO jg=1,ngrnd
      zz(jg)      = fz(REAL(jg,r_std) - undemi)
      zz_coef(jg) = fz(REAL(jg,r_std)-undemi+undemi)
    ENDDO

    !! 2.3 Converting to meters
    DO jg=1,ngrnd
      zz(jg)      = zz(jg) /  cstgrnd * lskin
      zz_coef(jg) = zz_coef(jg) / cstgrnd * lskin 
      dz2(jg)     = dz2(jg) /  cstgrnd * lskin
    ENDDO

    !! 2.4 Computing some usefull constants for the numerical scheme
    DO jg=1,ngrnd-1
      dz1(jg)  = un / (zz(jg+1) - zz(jg))
    ENDDO
    lambda = zz(1) * dz1(1)

    !! 2.5 Get the wetness profile on the thermal vertical grid from the diagnostic axis
    CALL thermosoilc_humlev(kjpindex, shumdiag_perma, snow)

  !! 3. Diagnostics : consistency checks on the vertical grid.

    WRITE (numout,*) 'diagnostic des niveaux dans le sol' !!?? to be changed,
    WRITE (numout,*) 'niveaux intermediaires et pleins'
    sum = zero
    DO jg=1,ngrnd
      sum = sum + dz2(jg)
      WRITE (numout,*) zz(jg),sum
    ENDDO

    
  !! 4. Compute a first diagnostic temperature profile
    
    CALL thermosoilc_diaglev(kjpindex, stempdiag)

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

  END SUBROUTINE thermosoilc_var_init
  

!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_coef
!!
!>\BRIEF        Calculate soil thermal properties, integration coefficients, apparent heat flux,
!! surface heat capacity,  
!!
!! DESCRIPTION	: This routine computes : \n
!!		1. the soil thermal properties. \n 
!!		2. the integration coefficients of the thermal numerical scheme, cgrnd and dgrnd,
!!              which depend on the vertical grid and on soil properties, and are used at the next 
!!              timestep.\n
!!              3. the soil apparent heat flux and surface heat capacity soilflux
!!              and soilcap, used by enerbil to compute the surface temperature at the next
!!              timestep.\n
!!             -  The soil thermal properties depend on water content (shum_ngrnd_perma, shumdiag_perma) and on the presence 
!!              of snow : snow is integrated into the soil for the thermal calculations, ie if there 
!!              is snow on the ground, the first thermal layer(s) consist in snow, depending on the 
!!              snow-depth. If a layer consists out of snow and soil, wheighed soil properties are 
!!              calculated\n
!!             - The coefficients cgrnd and dgrnd are the integration
!!              coefficients for the thermal scheme \n
!!                              T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n
!!                                      -- EQ1 -- \n
!!              They correspond respectively to $\beta$ and $\alpha$ from F. Hourdin\'s thesis and 
!!              their expression can be found in this document (eq A19 and A20)
!!             - soilcap and soilflux are the apparent surface heat capacity and flux
!!               used in enerbil at the next timestep to solve the surface
!!               balance for Ts (EQ3); they correspond to $C_s$ and $F_s$ in F.
!!               Hourdin\'s PhD thesis and are expressed in eq. A30 and A31. \n
!!                 soilcap*(Ts(t)-Ts(t-1))/dt=soilflux+otherfluxes(Ts(t)) \n
!!                                      -- EQ3 --\n
!!
!! RECENT CHANGE(S) : None
!!
!! MAIN OUTPUT VARIABLE(S): cgrnd, dgrnd, pcapa, pkappa, soilcap, soilflx
!!
!! REFERENCE(S) :
!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
!! integration scheme has been scanned and is provided along with the documentation, with name : 
!! Hourdin_1992_PhD_thermal_scheme.pdf
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

  SUBROUTINE thermosoilc_coef (kjpindex, temp_sol_new, snow,            &
                              ptn,      soilcap, soilflx,      zz,     &
                              dz1,      dz2,                           &
                              cgrnd,    dgrnd,                         &
                              frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
                              snowdz,snowrho,  snowtemp,     pb,         &
			      lambda_snow,    cgrnd_snow,    dgrnd_snow)

    !! 0. Variables and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                             :: kjpindex     !! Domain size (unitless)
    REAL(r_std), DIMENSION (kjpindex), INTENT (in)         :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex
    REAL(r_std), DIMENSION (kjpindex), INTENT (in)         :: snow         !! snow mass @tex ($Kg$) @endtex
    REAL(r_std), DIMENSION (ngrnd), INTENT(in)             :: zz           !! depths of the soil thermal numerical nodes 
                                                                           !! @tex ($m$) @endtex 
    REAL(r_std), DIMENSION (ngrnd), INTENT(in)             :: dz1          !! numerical constant depending on the vertical 
                                                                           !! thermal grid only @tex ($m^{-1}$) @endtex 
    REAL(r_std), DIMENSION (ngrnd), INTENT(in)             :: dz2          !! thicknesses of the soil thermal layers
                                                                           !! @tex ($m$) @endtex 
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)           :: frac_snow_veg   !! Snow cover fraction on vegeted area
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in)    :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
    REAL(r_std),DIMENSION (kjpindex),INTENT(in)            :: totfrac_nobio   !! Total fraction of continental ice+lakes+cities+...
                                                                              !!(unitless,0-1)
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)    :: snowdz          !! Snow depth
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)    :: snowrho         !! Snow density
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)    :: snowtemp        !! Snow temperature
    REAL(r_std), DIMENSION (kjpindex), INTENT (in)         :: pb              !! Surface presure (hPa)
    
    !! 0.2 Output variables

    REAL(r_std), DIMENSION (kjpindex), INTENT (out)        :: soilcap      !! surface heat capacity
                                                                           !! @tex ($J m^{-2} K^{-1}$) @endtex
    REAL(r_std), DIMENSION (kjpindex), INTENT (out)        :: soilflx      !! surface heat flux @tex ($W m^{-2}$) @endtex,
                                                                           !! positive towards the 
                                                                           !! soil, writen as Qg (ground heat flux) in the history 
                                                                           !! files.
    REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd        !! matrix coefficient for the computation of soil 
                                                                           !! temperatures (beta in F. Hourdin thesis)
    REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd        !! matrix coefficient for the computation of soil 
                                                                           !! temperatures (alpha in F. Hourdin thesis)


    !! 0.3 Modified variable

    REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (inout):: ptn          !! vertically discretized soil temperatures
                                                                           !! @tex ($K$) @endtex
    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)       :: lambda_snow  !! Coefficient of the linear extrapolation of surface temperature 
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow   !! Integration coefficient for snow numerical scheme
    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow   !! Integration coefficient for snow numerical scheme

    !! 0.4 Local variables

    INTEGER(i_std)                                         :: ji, jg
    REAL(r_std), DIMENSION (kjpindex,ngrnd-1)              :: zdz1         !! numerical (buffer) constant 
                                                                           !! @tex ($W m^{-1} K^{-1}$) @endtex
    REAL(r_std), DIMENSION (kjpindex,ngrnd)                :: zdz2         !! numerical (buffer) constant  
                                                                           !! @tex ($W m^{-1} K^{-1}$) @endtex
    REAL(r_std), DIMENSION (kjpindex)                      :: z1           !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex
    REAL(r_std), DIMENSION (kjpindex)                      :: soilcap_nosnow !! surface heat capacity
                                                                             !! @tex ($J m^{-2} K^{-1}$) @endtex
    REAL(r_std), DIMENSION (kjpindex)                      :: soilflx_nosnow !! surface heat flux @tex ($W m^{-2}$) @endtex,
                                                                             !! positive towards the soil, written as Qg (ground heat flux in the history files).
    REAL(r_std), DIMENSION (kjpindex)                      :: snowcap        !! apparent snow heat capacity @tex ($J m^{-2} K^{-1}$)
    REAL(r_std), DIMENSION (kjpindex)                      :: snowflx        !! apparent snow-atmosphere heat flux @tex ($W m^{-2}$) @endtex
    REAL(r_std), DIMENSION (kjpindex,nsnow)                :: dz1_snow
    REAL(r_std), DIMENSION (kjpindex,nsnow)                :: ZSNOWDZM
    REAL(r_std), DIMENSION (kjpindex,nsnow)                :: dz2_snow
    REAL(r_std), DIMENSION (kjpindex,nsnow)                :: zdz1_snow
    REAL(r_std), DIMENSION (kjpindex,nsnow)                :: zdz2_snow
    REAL(r_std), DIMENSION (kjpindex)                      :: z1_snow

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

  !! 1. Computation of the soil thermal properties
   
    ! Computation of the soil thermal properties; snow properties are also accounted for
    IF (ok_freeze_thermix) THEN
       CALL thermosoilc_getdiff( kjpindex, snow, ptn, snowrho, snowtemp, pb )
    ELSE IF (ok_explicitsnow) THEN
       ! Special case with explicit snow model without soil freezing
       CALL thermosoilc_getdiff_old_thermix_without_snow( kjpindex, snowrho, snowtemp, pb )
    ELSE
       ! Special case with old snow without soil freezing
       CALL thermosoilc_getdiff_old_thermix_with_snow( kjpindex, snow )
    ENDIF

    ! Energy conservation : Correction to make sure that the same latent heat is released and 
    ! consumed during freezing and thawing
    IF (ok_Ecorr) THEN
       CALL thermosoilc_readjust(kjpindex, ptn)
    ENDIF
    

    !! 2. Computation of the coefficients of the numerical integration scheme for the soil layers

    !! 2.1 Calculate numerical coefficients zdz1 and zdz2
    DO jg=1,ngrnd
      DO ji=1,kjpindex
        zdz2(ji,jg)=pcapa(ji,jg) * dz2(jg)/dt_sechiba
      ENDDO
    ENDDO
    
    DO jg=1,ngrnd-1
      DO ji=1,kjpindex
        zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg)
      ENDDO
    ENDDO
    
    !! 2.2 Calculate coefficients cgrnd and dgrnd for soil
    DO ji = 1,kjpindex
      z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1)
      cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji)
      dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji)
    ENDDO

    DO jg = ngrnd-1,2,-1
      DO ji = 1,kjpindex
        z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg)))
        cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji)
        dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji)
      ENDDO
    ENDDO


    !! 3. Computation of the coefficients of the numerical integration scheme for the snow layers

    !! 3.1 Calculate numerical coefficients zdz1_snow, zdz2_snow and lambda_snow
    DO ji = 1, kjpindex

       dz2_snow(ji,:) = 0

       IF ( ok_explicitsnow ) THEN

          ! Calculate internal values
          DO jg = 1, nsnow
             ZSNOWDZM(ji,jg) = MAX(snowdz(ji,jg),psnowdzmin)
          ENDDO
          dz2_snow(ji,:)=ZSNOWDZM(ji,:)

          DO jg = 1, nsnow-1
             dz1_snow(ji,jg)  = 2.0 / (dz2_snow(ji,jg+1)+dz2_snow(ji,jg))
          ENDDO

          lambda_snow(ji) = dz2_snow(ji,1)/2.0 * dz1_snow(ji,1)

          DO jg=1,nsnow
             zdz2_snow(ji,jg)=pcapa_snow(ji,jg) * dz2_snow(ji,jg)/dt_sechiba
          ENDDO

          DO jg=1,nsnow-1
             zdz1_snow(ji,jg) = dz1_snow(ji,jg) * pkappa_snow(ji,jg)
          ENDDO

          ! the bottom snow
          zdz1_snow(ji,nsnow) = pkappa_snow(ji,nsnow) / ( zz_coef(1) + dz2_snow(ji,nsnow)/2 )

       ELSE
          ! There is no snow, these values will never be used
          lambda_snow(ji) = lambda
       ENDIF

    ENDDO

    !! 3.2 Calculate coefficients cgrnd_snow and dgrnd_snow for snow
    DO ji = 1,kjpindex
       IF ( ok_explicitsnow ) THEN
          ! bottom level
          z1_snow(ji) = zdz2(ji,1)+(un-dgrnd(ji,1))*zdz1(ji,1)+zdz1_snow(ji,nsnow)
          cgrnd_snow(ji,nsnow) = (zdz2(ji,1) * ptn(ji,1) + zdz1(ji,1) * cgrnd(ji,1) ) / z1_snow(ji)
          dgrnd_snow(ji,nsnow) = zdz1_snow(ji,nsnow) / z1_snow(ji)

          ! next-to-bottom level
          z1_snow(ji) = zdz2_snow(ji,nsnow)+(un-dgrnd_snow(ji,nsnow))*zdz1_snow(ji,nsnow)+zdz1_snow(ji,nsnow-1)
          cgrnd_snow(ji,nsnow-1) = (zdz2_snow(ji,nsnow)*snowtemp(ji,nsnow)+&
               zdz1_snow(ji,nsnow)*cgrnd_snow(ji,nsnow))/z1_snow(ji)
          dgrnd_snow(ji,nsnow-1) = zdz1_snow(ji,nsnow-1) / z1_snow(ji)

          DO jg = nsnow-1,2,-1
             z1_snow(ji) = un / (zdz2_snow(ji,jg) + zdz1_snow(ji,jg-1) + zdz1_snow(ji,jg) * (un - dgrnd_snow(ji,jg)))
             cgrnd_snow(ji,jg-1) = (snowtemp(ji,jg) * zdz2_snow(ji,jg) + zdz1_snow(ji,jg) * cgrnd_snow(ji,jg)) * z1_snow(ji)
             dgrnd_snow(ji,jg-1) = zdz1_snow(ji,jg-1) * z1_snow(ji)
          ENDDO
       ELSE
          ! There is no snow, these values will never be used
          cgrnd_snow(ji,:) = cgrnd(ji,1)
          dgrnd_snow(ji,:) = dgrnd(ji,1)
       ENDIF
    ENDDO


    !! 4. Computation of the apparent ground heat flux 
    !! Computation of apparent snow-atmosphere flux  
    DO ji = 1,kjpindex
       IF ( ok_explicitsnow ) THEN
          snowflx(ji) = zdz1_snow(ji,1) * (cgrnd_snow(ji,1) + (dgrnd_snow(ji,1)-1.) * snowtemp(ji,1))
          snowcap(ji) = (zdz2_snow(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd_snow(ji,1)) * zdz1_snow(ji,1))
          z1_snow(ji) = lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un 
          snowcap(ji) = snowcap(ji) / z1_snow(ji)
          snowflx(ji) = snowflx(ji) + &
               & snowcap(ji) * (snowtemp(ji,1) * z1_snow(ji) - lambda_snow(ji) * cgrnd_snow(ji,1) - temp_sol_new(ji)) / dt_sechiba
       ELSE
          snowflx(ji) = zero
          snowcap(ji) = zero
       ENDIF
    ENDDO


    !! Computation of the apparent ground heat flux (> towards the soil) and
    !! apparent surface heat capacity, used at the next timestep by enerbil to
    !! compute the surface temperature.
    DO ji = 1,kjpindex
      soilflx_nosnow(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1))
      soilcap_nosnow(ji) = (zdz2(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd(ji,1)) * zdz1(ji,1))
      z1(ji) = lambda * (un - dgrnd(ji,1)) + un
      soilcap_nosnow(ji) = soilcap_nosnow(ji) / z1(ji)
      soilflx_nosnow(ji) = soilflx_nosnow(ji) + &
         & soilcap_nosnow(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dt_sechiba
    ENDDO

    !! Add snow fraction
    IF ( ok_explicitsnow ) THEN
       ! Using an effective heat capacity and heat flux by a simple pondering of snow and soil fraction
       DO ji = 1, kjpindex
          soilcap(ji) = snowcap(ji)*frac_snow_veg(ji)*(1-totfrac_nobio(ji))+    & ! weights related to snow cover fraction on vegetation  
               soilcap_nosnow(ji)*SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji)+ & ! weights related to SCF on nobio
               soilcap_nosnow(ji)*(1-(frac_snow_veg(ji)*(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
          soilflx(ji) = snowflx(ji)*frac_snow_veg(ji)*(1-totfrac_nobio(ji))+    & ! weights related to snow cover fraction on vegetation  
               soilflx_nosnow(ji)*SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji)+ & ! weights related to SCF on nobio
               soilflx_nosnow(ji)*(1-(frac_snow_veg(ji)*(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
       ENDDO
    ELSE
       ! Do not consider snow fraction
       soilcap(:)=soilcap_nosnow(:)
       soilflx(:)=soilflx_nosnow(:)
    END IF

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

  END SUBROUTINE thermosoilc_coef
 
 
!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_profile
!!
!>\BRIEF        In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile; 
!! This profile is then exported onto the diagnostic axis (call thermosoilc_diaglev)
!!
!! DESCRIPTION	: The calculation of the new soil temperature profile is based on
!! the cgrnd and dgrnd values from the previous timestep and the surface temperature Ts aka temp_sol_new. (see detailed
!! explanation in the header of the thermosoilc module or in the reference).\n
!!                              T(k+1)=cgrnd(k)+dgrnd(k)*T(k)\n
!!                                      -- EQ1 --\n
!!                           Ts=(1-lambda)*T(1) -lambda*T(2)\n 
!!                                      -- EQ2--\n
!!
!! RECENT CHANGE(S) : None
!! 
!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis), 
!!                          stempdiag (soil temperature profile on the diagnostic axis)
!!
!! REFERENCE(S) :
!! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres,
!! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal
!! integration scheme has been scanned and is provided along with the documentation, with name : 
!! Hourdin_1992_PhD_thermal_scheme.pdf
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================

 SUBROUTINE thermosoilc_profile (kjpindex, temp_sol_new, ptn, stempdiag, &
                                 snowtemp, frac_snow_veg,frac_snow_nobio,totfrac_nobio, &
                                 cgrnd_snow,dgrnd_snow)

  !! 0. Variables and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                               :: kjpindex       !! Domain size (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_sol_new   !! Surface temperature at the present time-step 
                                                                               !! @tex ($K$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)             :: frac_snow_veg  !! Snow cover fraction on vegeted area
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in)      :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
    REAL(r_std),DIMENSION (kjpindex),INTENT(in)              :: totfrac_nobio  !! Total fraction of continental ice+lakes+cities+...(unitless,0-1)
    REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in)       :: snowtemp       !! Snow temperature
    REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in)       :: cgrnd_snow     !! Integration coefficient for snow numerical scheme
    REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in)       :: dgrnd_snow     !! Integration coefficient for snow numerical scheme

   
    !! 0.2 Output variables
    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)      :: stempdiag      !! diagnostic temperature profile 
                                                                               !! @tex ($K$) @endtex
    REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (out)     :: ptn            !! vertically discretized soil temperatures 
                                                                               !! @tex ($K$) @endtex

    !! 0.3 Modified variables


    !! 0.4 Local variables

    INTEGER(i_std)                                           :: ji, jg
    REAL(r_std)                                              :: temp_sol_eff
     
!_ ================================================================================================================================
    
  !! 1. Computes the soil temperatures ptn.

    !! 1.1. ptn(jg=1) using EQ1 and EQ2
    DO ji = 1,kjpindex

       IF ( ok_explicitsnow ) THEN
          ! Soil temperature calculation with explicit snow if there is snow on the ground
          ! using an effective surface temperature by a simple pondering 
          temp_sol_eff=snowtemp(ji,nsnow)*frac_snow_veg(ji)*(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation  
               temp_sol_new(ji)*SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji)+ & ! weights related to SCF on nobio
               temp_sol_new(ji)*(1-(frac_snow_veg(ji)*(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
          ! Soil temperature calculation with explicit snow if there is snow on the ground
          ptn(ji,1) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * temp_sol_eff
       ELSE
          ! Standard soil temperature calculation
          ptn(ji,1) = (lambda * cgrnd(ji,1) + temp_sol_new(ji)) / (lambda *(un - dgrnd(ji,1)) + un)
       ENDIF
    ENDDO

    !! 1.2. ptn(jg=2:ngrnd) using EQ1.
    DO jg = 1,ngrnd-1
      DO ji = 1,kjpindex
        ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg)
      ENDDO
    ENDDO

  !! 2. Put the soil temperatures onto the diagnostic axis 
  
    !! Put the soil temperatures onto the diagnostic axis for convenient
    !! use in other routines (stomate..)
    CALL thermosoilc_diaglev(kjpindex, stempdiag)

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

  END SUBROUTINE thermosoilc_profile


!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_diaglev
!!
!>\BRIEF        Interpolation of the soil in-depth temperatures onto the diagnostic profile.
!!
!! DESCRIPTION  : This is a very easy linear interpolation, with intfact(jsl, jg) the fraction
!! the thermal layer jg comprised within the diagnostic layer jsl. The depths of
!! the diagnostic levels are diaglev(1:nslm), computed in slowproc.f90.
!!
!! RECENT CHANGE(S) : None
!!
!! MAIN OUTPUT VARIABLE(S): stempdiag (soil temperature profile on the diagnostic axis)
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================

  SUBROUTINE thermosoilc_diaglev(kjpindex, stempdiag)
    
  !! 0. Variables and parameter declaration

    !! 0.1 Input variables
 
    INTEGER(i_std), INTENT(in)                          :: kjpindex       !! Domain size (unitless)
   
    !! 0.2 Output variables

    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag      !! Diagnostoc soil temperature profile @tex ($K$) @endtex
    
    !! 0.3 Modified variables

    !! 0.4 Local variables

    INTEGER(i_std)                                      :: ji, jsl, jg
    REAL(r_std)                                         :: lev_diag, prev_diag, lev_prog, prev_prog
    REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)      :: intfact
!$OMP THREADPRIVATE(intfact)
    LOGICAL, PARAMETER                                  :: check=.FALSE.
!_ ================================================================================================================================
    
  !! 1. Computes intfact(jsl, jg)

    !! Computes intfact(jsl, jg), the fraction
    !! the thermal layer jg comprised within the diagnostic layer jsl.
    IF ( .NOT. ALLOCATED(intfact)) THEN
        
        ALLOCATE(intfact(nslm, ngrnd))
        
        prev_diag = zero
        DO jsl = 1, nslm
          lev_diag = diaglev(jsl)
          prev_prog = zero
          DO jg = 1, ngrnd
             IF ( jg == ngrnd .AND. (prev_prog + dz2(jg)) < lev_diag ) THEN
                lev_prog = lev_diag
             ELSE
                lev_prog = prev_prog + dz2(jg)
             ENDIF
            intfact(jsl,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
            prev_prog = lev_prog
          ENDDO
           prev_diag = lev_diag
        ENDDO
        
        IF ( check ) THEN
           WRITE(numout,*) 'thermosoilc_diagev -- thermosoilc_diaglev -- thermosoilc_diaglev --' 
           DO jsl = 1, nslm
              WRITE(numout,*) jsl, '-', intfact(jsl,1:ngrnd)
           ENDDO
           WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM"
           DO jsl = 1, nslm
              WRITE(numout,*) jsl, '-', SUM(intfact(jsl,1:ngrnd))
           ENDDO
           WRITE(numout,*) 'thermosoilc_diaglev -- thermosoilc_diaglev -- thermosoilc_diaglev --' 
        ENDIF
        
    ENDIF

 !! 2. does the interpolation

    stempdiag(:,:) = zero
    DO jg = 1, ngrnd
      DO jsl = 1, nslm
        DO ji = 1, kjpindex
          stempdiag(ji,jsl) = stempdiag(ji,jsl) + ptn(ji,jg)*intfact(jsl,jg)
        ENDDO
      ENDDO
    ENDDO

  END SUBROUTINE thermosoilc_diaglev
  

!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_humlev
!!
!>\BRIEF        Interpolates the diagnostic soil humidity profile shumdiag_perma(nslm, diagnostic axis) onto 
!!              the thermal axis, which gives shum_ngrnd_perma(ngrnd, thermal axis).
!!
!! DESCRIPTION  : Same as in thermosoilc_diaglev : This is a very easy linear interpolation, with intfactw(jsl, jg) the fraction
!! the thermal layer jsl comprised within the diagnostic layer jg. 
!! 
!! The depths of the diagnostic levels are diaglev(1:nslm), computed in slowproc.f90.
!!
!! RECENT CHANGE(S) : None
!!
!! MAIN OUTPUT VARIABLE(S): shum_ngrnd_perma (soil saturation degree on the thermal axis)
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================
  SUBROUTINE thermosoilc_humlev(kjpindex, shumdiag_perma, snow)
  
  !! 0. Variables and parameter declaration

    !! 0.1 Input variables
 
    INTEGER(i_std), INTENT(in)                            :: kjpindex    !! Domain size (unitless)
    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in)    :: shumdiag_perma  
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: snow 
    
    !! 0.2 Output variables

    !! 0.3 Modified variables

    !! 0.4 Local variables
    INTEGER(i_std)                                       :: ji, jsl, jg
    REAL(r_std)                                          :: lev_diag, prev_diag, lev_prog, prev_prog
    REAL(r_std), DIMENSION(ngrnd,nslm)                   :: intfactw     !! fraction of each diagnostic layer (jsl) comprized within
                                                                         !! a given thermal layer (jg)(0-1, unitless) 
    REAL(r_std), DIMENSION(kjpindex)               :: snow_h
    LOGICAL, PARAMETER :: check=.FALSE.

!_ ================================================================================================================================
    
  !! 1. computes intfactw(jsl,jg), the fraction of each diagnostic layer (jg) comprized within a given thermal layer (jsl)
    IF ( check ) &
         WRITE(numout,*) 'thermosoilc_humlev --' 

    ! Snow height
    snow_h(:)=snow(:)/sn_dens
    !
    shum_ngrnd_perma(:,:) = zero
    DO ji=1,kjpindex
       prev_diag = zero
       DO jsl = 1, ngrnd
          lev_diag = prev_diag + dz2(jsl)
          prev_prog = snow_h(ji)
          DO jg = 1, nslm
             IF ( jg == nslm .AND. diaglev(jg)+snow_h(ji) < lev_diag ) THEN
                lev_prog = lev_diag+snow_h(ji)
             ELSE
                lev_prog = diaglev(jg)+snow_h(ji)
             ENDIF
             intfactw(jsl,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag)
             prev_prog = lev_prog
          ENDDO
          prev_diag = lev_diag
       ENDDO

       DO jsl = 1, nslm
          DO jg = 1, ngrnd
             shum_ngrnd_perma(ji,jg) = shum_ngrnd_perma(ji,jg) + shumdiag_perma(ji,jsl)*intfactw(jg,jsl)
          ENDDO
       ENDDO
       !
       IF ( check ) THEN
          DO jg = 1, ngrnd
             WRITE(numout,*) ji,jg, '-', intfactw(jg,1:nslm),'-sum-', SUM(intfactw(jg,1:nslm))
          ENDDO
       ENDIF
    ENDDO

    IF ( check ) &
         WRITE(numout,*) 'thermosoilc_humlev --' 

  END SUBROUTINE thermosoilc_humlev


!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_energy
!!
!>\BRIEF         Energy check-up.
!!
!! DESCRIPTION  : I didn\'t comment this routine since at do not understand its use, please
!! ask initial designers (Jan ? Nathalie ?).
!!
!! RECENT CHANGE(S) : None
!!
!! MAIN OUTPUT VARIABLE(S) : ??
!!
!! REFERENCE(S) : None
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================

  SUBROUTINE thermosoilc_energy(kjpindex, temp_sol_new, soilcap)
  
   !! 0. Variables and parameter declaration

    !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                     :: kjpindex     !! Domain size (unitless)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: temp_sol_new !! Surface temperature at the present time-step, Ts 
                                                                   !! @tex ($K$) @endtex
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)  :: soilcap      !! Apparent surface heat capacity 
                                                                   !! @tex ($J m^{-2} K^{-1}$) @endtex, 
                                                                   !! see eq. A29 of F. Hourdin\'s PhD thesis.
    
    !! 0.2 Output variables

    !! 0.3 Modified variables
    
    !! 0.4 Local variables
    
    INTEGER(i_std)                                 :: ji, jg
!_ ================================================================================================================================
   

     DO ji = 1, kjpindex
      surfheat_incr(ji) = zero
      coldcont_incr(ji) = zero
     ENDDO
    
    !  Sum up the energy content of all layers in the soil.
    DO ji = 1, kjpindex
   
       IF (pcapa_en(ji,1) .LE. sn_capa) THEN
          
          ! Verify the energy conservation in the surface layer
          coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
          surfheat_incr(ji) = zero
       ELSE
          
          ! Verify the energy conservation in the surface layer
          surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji))
          coldcont_incr(ji) = zero
       ENDIF
    ENDDO
    
    ptn_beg(:,:)      = ptn(:,:)
    temp_sol_beg(:)   = temp_sol_new(:)

  END SUBROUTINE thermosoilc_energy



!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_readjust
!!
!>\BRIEF        
!!
!! DESCRIPTION	: Energy conservation : Correction to make sure that the same latent heat is released and 
!!                consumed during freezing and thawing  
!!
!! RECENT CHANGE(S) : None
!! 
!! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis), 
!!                          
!! REFERENCE(S) :
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================

  SUBROUTINE thermosoilc_readjust(kjpindex, ptn)

   !! 0. Variables and parameter declaration

    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in)                             :: kjpindex

    !! 0.2 Modified variables
    REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(inout)    :: ptn

    !! 0.3 Local variables
    INTEGER(i_std)  :: ji, jg
    REAL(r_std) :: ptn_tmp

    DO jg=1, ngrnd
       DO ji=1, kjpindex
          ! All soil latent energy is put into e_soil_lat(ji)
          ! because the variable soil layers make it difficult to keep track of all
          ! layers in this version
          ! NOTE : pcapa has unit J/K/m3 and pcappa_supp has J/K
          e_soil_lat(ji)=e_soil_lat(ji)+pcappa_supp(ji,jg)*(ptn(ji,jg)-ptn_beg(ji,jg))
       END DO
   END DO

   DO ji=1, kjpindex
      IF (e_soil_lat(ji).GT.min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
         ! The soil is thawed: we spread the excess of energy over the uppermost 6 levels e.g. 2.7m
         ! Here we increase the temperatures
         DO jg=1,6
            ptn_tmp=ptn(ji,jg)
            
            ptn(ji,jg)=ptn(ji,jg)+MIN(e_soil_lat(ji)/pcapa(ji,jg)/zz_coef(6), 0.5)
            e_soil_lat(ji)=e_soil_lat(ji)-(ptn(ji,jg)-ptn_tmp)*pcapa(ji,jg)*dz2(jg)
         ENDDO
      ELSE IF (e_soil_lat(ji).LT.-min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN
         ! The soil is thawed
         ! Here we decrease the temperatures
         DO jg=1,6
            ptn_tmp=ptn(ji,jg)
            ptn(ji,jg)=MAX(ZeroCelsius+fr_dT/2., ptn_tmp+e_soil_lat(ji)/pcapa(ji,jg)/zz_coef(6))
            e_soil_lat(ji)=e_soil_lat(ji)+(ptn_tmp-ptn(ji,jg))*pcapa(ji,jg)*dz2(jg)
         END DO
      END IF
   END DO

  END SUBROUTINE thermosoilc_readjust
   
!-------------------------------------------------------------------



!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_getdiff
!!
!>\BRIEF          Computes soil and snow heat capacity and conductivity    
!!
!! DESCRIPTION	: Computation of the soil and snow thermal properties; snow properties are also accounted for
!!
!! RECENT CHANGE(S) : None
!! 
!! MAIN OUTPUT VARIABLE(S): 
!!                          
!! REFERENCE(S) :
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================

  SUBROUTINE thermosoilc_getdiff( kjpindex, snow, ptn, snowrho, snowtemp, pb )

   !! 0. Variables and parameter declaration

    !! 0.1 Input variables
    INTEGER(i_std),INTENT(in)				:: kjpindex
    REAL(r_std),DIMENSION(kjpindex),INTENT (in)	        :: snow       !! Snow mass
    REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in)    :: snowrho    !! Snow density
    REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in)    :: snowtemp   !! Snow temperature
    REAL(r_std),DIMENSION(kjpindex),INTENT(in)          :: pb         !! Surface pressure
    REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in)	:: ptn        !! Soil temperature profile

    !! 0.3 Local variables
    REAL						:: xx         !! Unfrozen fraction of the soil
    REAL(r_std), DIMENSION(kjpindex)             	:: snow_h
    REAL(r_std), DIMENSION(kjpindex,ngrnd) 		:: zx1, zx2    
    REAL			   			:: cap_iw     !! Heat capacity of ice/water mixture 
    REAL						:: csat       !! Thermal conductivity for saturated soil
    INTEGER						:: ji,jg


    DO ji = 1,kjpindex
       IF (.NOT. ok_explicitsnow) THEN      
          ! 1. Determine the fractions of snow and soil
          snow_h(ji) = snow(ji) / sn_dens
      
          !
          !  1.1. The first level
          !
          IF ( snow_h(ji) .GT. zz_coef(1) ) THEN
             ! the 1st level is in the snow => the 1st layer is entirely snow
             zx1(ji,1) = 1.
             zx2(ji,1) = 0.
          ELSE IF ( snow_h(ji) .GT. zero ) THEN      
             ! the 1st level is beyond the snow and the snow is present
             zx1(ji,1) = snow_h(ji) / zz_coef(1)
             zx2(ji,1) = ( zz_coef(1) - snow_h(ji)) / zz_coef(1)	
          ELSE
             ! there is no snow at all, quoi ;-)
             zx1(ji,1) = 0.
             zx2(ji,1) = 1.       
          ENDIF
      
          !
          !  1.2. The other levels
          !
          DO jg = 2, ngrnd
             IF ( snow_h(ji) .GT. zz_coef(jg) ) THEN
                ! the current level is in the snow => the current layer is entirely snow
                zx1(ji,jg) = 1.
                zx2(ji,jg) = 0.
             ELSE IF ( snow_h(ji) .GT. zz_coef(jg-1) ) THEN
                ! the current layer is partially snow
                zx1(ji,jg) = (snow_h(ji) - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1))
                zx2(ji,jg) = ( zz_coef(jg) - snow_h(ji)) / (zz_coef(jg) - zz_coef(jg-1))
             ELSE
                ! both levels are out of snow => the current layer is entirely soil	  
                zx1(ji,jg) = 0.
                zx2(ji,jg) = 1.       
             ENDIF
          ENDDO
       ELSE
          zx1(ji,:) = 0.
          zx2(ji,:) = 1.
       END IF

      DO jg = 1, ngrnd
         !
         ! 2. Calculate heat capacity with allowance for permafrost
         ! 2.1. soil heat capacity depending on temperature and humidity
	 
         IF (ptn(ji,jg) .LT. ZeroCelsius-fr_dT/2.) THEN
	    ! frozen soil
            pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_ice - so_capa_dry)!Isa : old version, proved to be correct
	    profil_froz(ji,jg) = 1.
 	    pcappa_supp(ji,jg)= 0.

     	 ELSEIF (ptn(ji,jg) .GT. ZeroCelsius+fr_dT/2.) THEN
	    ! unfrozen soil	 
     	    pcapa(ji,jg) =  so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)!Isa : old version, proved to be correct
	    profil_froz(ji,jg) = 0.
 	    pcappa_supp(ji,jg)= 0.
     	 ELSE
     	   ! xx is the unfrozen fraction of soil water	   	   
     	   xx = (ptn(ji,jg)-(ZeroCelsius-fr_dT/2.)) / fr_dT
           profil_froz(ji,jg) = (1. - xx)

	   ! net heat capacity of the ice/water mixture
     	   cap_iw = xx * so_capa_wet + (1.-xx) * so_capa_ice
	   pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(cap_iw-so_capa_dry) + shum_ngrnd_perma(ji,jg)*poros*lhf*rho_water/fr_dT
 	   pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*poros*lhf*rho_water/fr_dT*zx2(ji,jg)*dz2(jg)

         ENDIF

         !
	 ! 2.2. Take into account the snow and soil fractions in the layer
	 !
         pcapa(ji,jg) = zx1(ji,jg) * sn_capa + zx2(ji,jg) * pcapa(ji,jg)

	 !
	 ! 2.3. Calculate the heat capacity for energy conservation check 
	 IF ( zx1(ji,jg).GT.0. ) THEN
            pcapa_en(ji,jg) = sn_capa
	 ELSE
            pcapa_en(ji,jg) = pcapa(ji,jg)
	 ENDIF
 
         !
         ! 3. Calculate the heat conductivity with allowance for permafrost (Farouki, 1981, Cold Reg. Sci. Technol.)
         !
	 ! 3.1. unfrozen fraction
         xx = (ptn(ji,jg)-(ZeroCelsius-fr_dT/2.)) / fr_dT * poros
         xx = MIN( poros, MAX( 0., xx ) )

         ! 3.2. saturated conductivity
      	 csat = cond_solid**(1.-poros) * cond_ice**(poros-xx) * cond_water**xx
     	
       	 ! 3.3. unsaturated conductivity
     	 pkappa(ji,jg) =(csat - so_cond_dry)*shum_ngrnd_perma(ji,jg) + so_cond_dry

         !
	 ! 3.4. Take into account the snow and soil fractions in the layer
         pkappa(ji,jg) = un / ( zx1(ji,jg) / sn_cond + zx2(ji,jg) / pkappa(ji,jg) )
      END DO            
    ENDDO   

    ! 4. Calculate snow heat capacity and conductivity
     DO ji = 1,kjpindex
         pcapa_snow(ji,:) = snowrho(ji,:) * xci
         pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2*snowrho(ji,:)*snowrho(ji,:)) +      &
               MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
                ZSNOWTHRMCOND_CVAP)))*(XP00/(pb(ji)*100.)))
    END DO


   
   END SUBROUTINE thermosoilc_getdiff

!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_getdiff_old_thermix_with_snow
!!
!>\BRIEF          Computes soil heat capacity and conductivity    
!!
!! DESCRIPTION	: Computes soil heat capacity and conductivity
!!                Special case with old snow without soil freezing
!!
!! RECENT CHANGE(S) : None
!! 
!! MAIN OUTPUT VARIABLE(S):
!!                          
!! REFERENCE(S) :
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================


   SUBROUTINE thermosoilc_getdiff_old_thermix_with_snow( kjpindex, snow )


   !! 0. Variables and parameter declaration

    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in) :: kjpindex
    REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow

    !! 0.2 Local variables
    INTEGER					:: ji,jg
    REAL(r_std)                                 :: snow_h       !! snow_h is the snow height @tex ($m$) @endtex 
    REAL(r_std)                                 :: zx1, zx2     !! zx1 and zx2 are the layer fraction consisting in snow and soil respectively.

     
    ! Computation of the soil thermal properties; snow properties are also accounted for
    DO ji = 1,kjpindex
      snow_h = snow(ji) / sn_dens

      ! First layer
      IF ( snow_h .GT. zz_coef(1) ) THEN
          pcapa(ji,1) = sn_capa
          pcapa_en(ji,1) = sn_capa
          pkappa(ji,1) = sn_cond
      ELSE IF ( snow_h .GT. zero ) THEN
          pcapa_en(ji,1) = sn_capa
          zx1 = snow_h / zz_coef(1)
          zx2 = ( zz_coef(1) - snow_h) / zz_coef(1)
          pcapa(ji,1) = zx1 * sn_capa + zx2 * so_capa_wet
          pkappa(ji,1) = un / ( zx1 / sn_cond + zx2 / so_cond_wet )
      ELSE
          pcapa(ji,1) = so_capa_dry + shum_ngrnd_perma(ji,1)*(so_capa_wet - so_capa_dry)
          pkappa(ji,1) = so_cond_dry + shum_ngrnd_perma(ji,1)*(so_cond_wet - so_cond_dry)
          pcapa_en(ji,1) = so_capa_dry + shum_ngrnd_perma(ji,1)*(so_capa_wet - so_capa_dry)
      ENDIF

      ! Mid layers
      DO jg = 2, ngrnd - 2
        IF ( snow_h .GT. zz_coef(jg) ) THEN
            pcapa(ji,jg) = sn_capa
            pkappa(ji,jg) = sn_cond
            pcapa_en(ji,jg) = sn_capa
        ELSE IF ( snow_h .GT. zz_coef(jg-1) ) THEN
            zx1 = (snow_h - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1))
            zx2 = ( zz_coef(jg) - snow_h) / (zz_coef(jg) - zz_coef(jg-1))
            pcapa(ji,jg) = zx1 * sn_capa + zx2 * so_capa_wet
            pkappa(ji,jg) = un / ( zx1 / sn_cond + zx2 / so_cond_wet )
            pcapa_en(ji,jg) = sn_capa
        ELSE
            pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
            pkappa(ji,jg) = so_cond_dry + shum_ngrnd_perma(ji,jg)*(so_cond_wet - so_cond_dry)
            pcapa_en(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
        ENDIF
      ENDDO

      ! Last two layers: These layers can not be filled with snow
      DO jg = ngrnd - 1, ngrnd
         pcapa(ji,jg) = so_capa_dry
         pkappa(ji,jg) = so_cond_dry
         pcapa_en(ji,jg) = so_capa_dry
      END DO
      
    ENDDO ! DO ji = 1,kjpindex


    END SUBROUTINE thermosoilc_getdiff_old_thermix_with_snow



!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_getdiff_old_thermix_without_snow
!!
!>\BRIEF          Computes soil and snow heat capacity and conductivity    
!!
!! DESCRIPTION	: Calculations of soil and snow thermal properties without effect of freezing. This subroutine is only called
!!                when explicitsnow is activated.
!!                
!!
!! RECENT CHANGE(S) : None
!! 
!! MAIN OUTPUT VARIABLE(S):
!!                          
!! REFERENCE(S) :
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================

    SUBROUTINE thermosoilc_getdiff_old_thermix_without_snow( kjpindex,snowrho, snowtemp, pb )

   !! 0. Variables and parameter declaration

    !! 0.1 Input variables
      INTEGER(i_std), INTENT(in) :: kjpindex
      REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho  !! Snow density
      REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
      REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: pb       !! Surface presure (hPa)


    !! 0.1 Local variables
      INTEGER						  :: ji,jg

     ! soil
      DO jg = 1,ngrnd
         DO ji = 1,kjpindex
            pkappa(ji,jg) = so_cond_dry + shum_ngrnd_perma(ji,jg)*(so_cond_wet - so_cond_dry)
            pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
            pcapa_en(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
         ENDDO
      ENDDO

     ! snow 
     DO ji = 1,kjpindex
         pcapa_snow(ji,:) = snowrho(ji,:) * xci
         pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2*snowrho(ji,:)*snowrho(ji,:)) +      & 
              MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ &
               ZSNOWTHRMCOND_CVAP)))*(XP00/(pb(ji)*100.)))
     END DO


   END SUBROUTINE thermosoilc_getdiff_old_thermix_without_snow


!! ================================================================================================================================
!! SUBROUTINE   : thermosoilc_read_reftempfile
!!
!>\BRIEF          
!!
!! DESCRIPTION	: Read file with longterm temperature
!!                
!!
!! RECENT CHANGE(S) : None
!! 
!! MAIN OUTPUT VARIABLE(S): reftemp : Reference temerature
!!                          
!! REFERENCE(S) :
!!
!! FLOWCHART    : None 
!! \n 
!_ ================================================================================================================================

 SUBROUTINE thermosoilc_read_reftempfile(kjpindex,lalo,reftemp)
    
    USE interpweight

    IMPLICIT NONE
    !! 0. Variables and parameter declaration

    !! 0.1 Input variables
    INTEGER(i_std), INTENT(in) :: kjpindex
    REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo
    
    !! 0.2 Output variables
    REAL(r_std), DIMENSION(kjpindex, ngrnd), INTENT(out) :: reftemp
    REAL(r_std), DIMENSION(kjpindex)                     :: areftemp         !! Availability of data for  the interpolation

    !! 0.3 Local variables
    INTEGER(i_std) :: ib
    CHARACTER(LEN=80) :: filename
    REAL(r_std),DIMENSION(kjpindex)                      :: reftemp_file     !! Horizontal temperature field interpolated from file [C]
    REAL(r_std)                                          :: vmin, vmax       !! min/max values to use for the
                                                                             !!   renormalization
    CHARACTER(LEN=80)                                    :: variablename     !! Variable to interpolate
                                                                             !!   the file
    CHARACTER(LEN=80)                                    :: lonname, latname !! lon, lat names in input file
    REAL(r_std), DIMENSION(:), ALLOCATABLE               :: variabletypevals !! Values for all the types of the variable
                                                                             !!   (variabletypevals(1) = -un, not used)
    CHARACTER(LEN=50)                                    :: fractype         !! method of calculation of fraction
                                                                             !!   'XYKindTime': Input values are kinds 
                                                                             !!     of something with a temporal 
                                                                             !!     evolution on the dx*dy matrix'
    LOGICAL                                              :: nonegative       !! whether negative values should be removed
    CHARACTER(LEN=50)                                    :: maskingtype      !! Type of masking
                                                                             !!   'nomask': no-mask is applied
                                                                             !!   'mbelow': take values below maskvals(1)
                                                                             !!   'mabove': take values above maskvals(1)
                                                                             !!   'msumrange': take values within 2 ranges;
                                                                             !!      maskvals(2) <= SUM(vals(k)) <= maskvals(1)
                                                                             !!      maskvals(1) < SUM(vals(k)) <= maskvals(3)
                                                                             !!       (normalized by maskvals(3))
                                                                             !!   'var': mask values are taken from a 
                                                                             !!     variable inside the file (>0)
    REAL(r_std), DIMENSION(3)                            :: maskvals         !! values to use to mask (according to 
                                                                             !!   `maskingtype') 
    CHARACTER(LEN=250)                                   :: namemaskvar      !! name of the variable to use to mask 
    REAL(r_std)                                          :: reftemp_norefinf
    REAL(r_std)                                          :: reftemp_default  !! Default value
       
    
    !Config Key   = SOIL_REFTEMP_FILE
    !Config Desc  = File with climatological soil temperature
    !Config If    = READ_REFTEMP
    !Config Def   = reftemp.nc
    !Config Help  = 
    !Config Units = [FILE]
    filename = 'reftemp.nc'
    CALL getin_p('REFTEMP_FILE',filename)

    variablename = 'temperature'

    IF (printlev >= 1) WRITE(numout,*) "thermosoilc_read_reftempfile: Read and interpolate file " &
         // TRIM(filename) //" for variable " //TRIM(variablename)

    ! For this case there are not types/categories. We have 'only' a continuos
    ! field
    ! Assigning values to vmin, vmax

    vmin = 0.
    vmax = 9999.

!   For this file we do not need neightbours!
    neighbours = 0

    !! Variables for interpweight
    ! Type of calculation of cell fractions
    fractype = 'default'
    ! Name of the longitude and latitude in the input file
    lonname = 'nav_lon'
    latname = 'nav_lat'
    ! Default value when no value is get from input file
    reftemp_default = 1.
    ! Reference value when no value is get from input file
    reftemp_norefinf = 1.
    ! Should negative values be set to zero from input file?
    nonegative = .FALSE.
    ! Type of mask to apply to the input data (see header for more details)
    maskingtype = 'nomask'
    ! Values to use for the masking (here not used)
    maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
    ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
    namemaskvar = ''

    CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours,                            &
      contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,        &
      maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf,                         &
      reftemp_file, areftemp)
    IF (printlev >= 5) WRITE(numout,*)'  thermosoilc_read_reftempfile after interpweight2D_cont'

    ! Copy reftemp_file temperature to all ground levels and transform into Kelvin
    DO ib=1, kjpindex
      reftemp(ib, :) = reftemp_file(ib)+ZeroCelsius
    END DO

    ! Write diagnostics
    CALL xios_orchidee_send_field("areftemp",areftemp)
       
  END SUBROUTINE thermosoilc_read_reftempfile
  

END MODULE thermosoilc
