!
! $Header$
!
SUBROUTINE soil_hetero(ptimestep, indice, knon, snow, ptsrf, qsol, &
     lon, lat, ptsoil, pcapcal, pfluxgrd, ztherm_i, conv_ratio)

  USE yomcst_mod_h, ONLY: RTT, RPI
  USE dimsoil_mod_h, ONLY: nsoilmx
  USE comsoil_mod_h
  USE compbl_mod_h
  USE dimpft_mod_h
  USE dimphy
  USE mod_phys_lmdz_para
  USE indice_sol_mod
  USE print_control_mod, ONLY: lunout
  USE phys_state_var_mod, ONLY: alpha_soil_tersrf, period_tersrf 
  USE surf_param_mod, ONLY: eff_surf_param

  IMPLICIT NONE

!=======================================================================
!
!   Auteur:  Frederic Hourdin     30/01/92
!   -------
!
!   Object:  Computation of : the soil temperature evolution
!   -------                   the surfacic heat capacity "Capcal"
!                            the surface conduction flux pcapcal
!
!   Update: 2021/07 : soil thermal inertia, formerly a constant value,
!   ------   can also be now a function of soil moisture (F Cheruy's idea)
!            depending on iflag_inertie, read from physiq.def via conf_phys_m.F90
!            ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy)
!           2025/04 : A. Maison, adapting the routine for heterogeneous continental sub-surfaces
!
!   Method: Implicit time integration
!   -------
!   Consecutive ground temperatures are related by:
!           T(k+1) = C(k) + D(k)*T(k)  (*)
!   The coefficients C and D are computed at the t-dt time-step.
!   Routine structure:
!   1) C and D coefficients are computed from the old temperature
!   2) new temperatures are computed using (*)
!   3) C and D coefficients are computed from the new temperature
!      profile for the t+dt time-step
!   4) the coefficients A and B are computed where the diffusive
!      fluxes at the t+dt time-step is given by
!             Fdiff = A + B Ts(t+dt)
!      or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
!             with F0 = A + B (Ts(t))
!                 Capcal = B*dt
!           
!   Interface:
!   ----------
!
!   Arguments:
!   ----------
!   ptimestep            physical timestep (s)
!   indice               sub-surface index
!   snow(klon)           snow
!   ptsrf(klon)          surface temperature at time-step t (K)
!   qsol(klon)           soil moisture (kg/m2 or mm)
!   lon(klon)            longitude in radian
!   lat(klon)            latitude in radian
!   ptsoil(klon,nsoilmx) temperature inside the ground (K)
!   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
!   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
!   ztherm_i(klon)       soil thermal inertia (J.m-2.K.s-1/2)
!   conv_ratio(klon)     ratio to convert soil depths in meters (-)
!   
!=======================================================================
!-----------------------------------------------------------------------
! Arguments
! ---------
  REAL, INTENT(IN)                     :: ptimestep
  INTEGER, INTENT(IN)                  :: indice, knon !, knindex
  REAL, DIMENSION(klon), INTENT(IN)    :: snow
  REAL, DIMENSION(klon), INTENT(IN)    :: ptsrf
  REAL, DIMENSION(klon), INTENT(IN)    :: qsol
  REAL, DIMENSION(klon), INTENT(IN)    :: lon
  REAL, DIMENSION(klon), INTENT(IN)    :: lat
  REAL, DIMENSION(klon), INTENT(IN)    :: ztherm_i
  REAL, DIMENSION(klon), INTENT(IN)    :: conv_ratio

  REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil
  REAL, DIMENSION(klon), INTENT(OUT)           :: pcapcal
  REAL, DIMENSION(klon), INTENT(OUT)           :: pfluxgrd

!-----------------------------------------------------------------------
! Local variables
! ---------------
  INTEGER                             :: ig, jk, ierr
  REAL, DIMENSION(nsoilmx)            :: zdz2
  REAL                                :: z1s
  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: C_coef, D_coef

! Local saved variables
! ---------------------
  REAL, SAVE                     :: lambda
!$OMP THREADPRIVATE(lambda)
  REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2
!$OMP THREADPRIVATE(dz1,dz2)
  LOGICAL, SAVE                  :: firstcall=.TRUE.
!$OMP THREADPRIVATE(firstcall)
    
!-----------------------------------------------------------------------
!   Depthts:
!   --------
  REAL fz,rk,fz1,rk1,rk2
  fz(rk)=fz1*(alpha_soil_tersrf**rk-1.)/(alpha_soil_tersrf-1.)
!
!-----------------------------------------------------------------------
! Calculation of some constants
! NB! These constants do not depend on the sub-surfaces
!-----------------------------------------------------------------------

  IF (firstcall) THEN
!-----------------------------------------------------------------------
!   ground levels 
!   grnd=z/l where l is the skin depth of the diurnal cycle:
!-----------------------------------------------------------------------
!   la premiere couche represente un dixieme de cycle diurne
     fz1=SQRT(period_tersrf/RPI)
     
     DO jk=1,nsoilmx
        rk1=jk
        rk2=jk-1
        dz2(jk)=fz(rk1)-fz(rk2)
     ENDDO
     DO jk=1,nsoilmx-1
        rk1=jk+.5
        rk2=jk-.5
        dz1(jk)=1./(fz(rk1)-fz(rk2))
     ENDDO
     lambda=fz(.5)*dz1(1)
     WRITE(lunout,*) 'surface index:', indice
     WRITE(lunout,*)'full layers, intermediate layers (seconds)'
     DO jk=1,nsoilmx
        rk=jk
        rk1=jk+.5
        rk2=jk-.5
        WRITE(lunout,*)'fz=', &
             fz(rk1)*fz(rk2)*RPI,fz(rk)*fz(rk)*RPI
     ENDDO
     WRITE(lunout,*)'full layers, intermediate layers (meters)'
     DO jk=1,nsoilmx
        rk=jk
        rk2=jk-.5
        WRITE(lunout,*)'fz=', &
             fz(rk2)*conv_ratio, fz(rk)*conv_ratio
     ENDDO
     firstcall =.FALSE.
  END IF

!-----------------------------------------------------------------------
! 1)
! Calculation of Cgrf and Dgrd coefficients using soil temperature from 
! previous time step.
!
! These variables are recalculated on the local compressed grid instead 
! of saved in restart file.
!-----------------------------------------------------------------------
  DO jk=1,nsoilmx
     zdz2(jk)=dz2(jk)/ptimestep
  ENDDO
  
  DO ig=1,knon
     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
     C_coef(ig,nsoilmx-1,indice)= &
          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
     D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
  ENDDO
  
  DO jk=nsoilmx-1,2,-1
     DO ig=1,knon
        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
             *(1.-D_coef(ig,jk,indice)))
        C_coef(ig,jk-1,indice)= &
             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
        D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s
     ENDDO
  ENDDO

!-----------------------------------------------------------------------
! 2)
! Computation of the soil temperatures using the Cgrd and Dgrd
! coefficient computed above
!
!-----------------------------------------------------------------------

!    Surface temperature
  DO ig=1,knon
     ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/  &
          (lambda*(1.-D_coef(ig,1,indice))+1.)
  ENDDO
  
!   Other temperatures
  DO jk=1,nsoilmx-1
     DO ig=1,knon
        ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
             *ptsoil(ig,jk)
     ENDDO
  ENDDO

  IF (indice == is_sic) THEN
     DO ig = 1 , knon
        ptsoil(ig,nsoilmx) = RTT - 1.8
     END DO
  ENDIF

!-----------------------------------------------------------------------
! 3)
! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil 
! temperature
!-----------------------------------------------------------------------
  DO ig=1,knon
     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
     C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
     D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s
  ENDDO
  
  DO jk=nsoilmx-1,2,-1
     DO ig=1,knon
        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
             *(1.-D_coef(ig,jk,indice)))
        C_coef(ig,jk-1,indice) = &
             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
        D_coef(ig,jk-1,indice) = dz1(jk-1)*z1s
     ENDDO
  ENDDO

!-----------------------------------------------------------------------
! 4)
! Computation of the surface diffusive flux from ground and
! calorific capacity of the ground
!-----------------------------------------------------------------------
  DO ig=1,knon
     pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
          (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
     pcapcal(ig)  = ztherm_i(ig)* &
          (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
     z1s = lambda*(1.-D_coef(ig,1,indice))+1.
     pcapcal(ig)  = pcapcal(ig)/z1s
     pfluxgrd(ig) = pfluxgrd(ig) &
          + pcapcal(ig) * (ptsoil(ig,1) * z1s &
          - lambda * C_coef(ig,1,indice) &
          - ptsrf(ig)) &
          /ptimestep
  ENDDO
    
END SUBROUTINE soil_hetero
