SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,   newtherm_i)

      use comsoil_h_PEM, only: layer_PEM, inertiedat_PEM

      IMPLICIT NONE

!=======================================================================
!   Description :
!       Surface water ice / Thermal inertia feedback.
!
!   When surface water-ice is thick enough, this routine creates a new
!   soil thermal inertia with three different layers :
!   - One layer of surface water ice (the thickness is given
!     by the variable icecover (in kg of ice per m2) and the thermal
!     inertia is prescribed by inert_h2o_ice (see surfdat_h)); 
!   - A transitional layer of mixed thermal inertia;
!   - A last layer of regolith below the ice cover whose thermal inertia
!     is equal to inertiedat.
!
!  To use the model :
!      SET THE tifeedback LOGICAL TO ".true." in callphys.def.
!
!  Author: Adapted from J.-B. Madeleine Mars 2008 ( Updated November 2012) by LL, 2022
!=======================================================================

!Local variables
!---------------

      INTEGER :: ig                     ! Grid point (ngrid)
      INTEGER :: ik                     ! Grid point (nsoil)
      INTEGER :: iref                   ! Ice/Regolith boundary index
      INTEGER, INTENT(IN) :: ngrid                  ! Number of horizontal grid points 
      INTEGER, INTENT(IN) :: nsoil                  ! Number of soil layers 
      REAL :: icedepth                  ! Ice cover thickness (m)
      REAL :: inert_h2o_ice = 800.      ! surface water ice thermal inertia [SI]
      REAL :: rho_ice = 920.            ! density of water ice [kg/m^3]
      REAL :: prev_thermi(ngrid,nsoil)  ! previous thermal inertia [SI]
!Inputs
!------

      REAL ,INTENT(IN):: icecover(ngrid)         ! tracer on the surface (kg.m-2)

!Outputs
!-------

      REAL,INTENT(INOUT) :: newtherm_i(ngrid,nsoil)    ! New soil thermal inertia

      prev_thermi(:,:) = newtherm_i(:,:)

!Creating the new soil thermal inertia table
!-------------------------------------------
      DO ig=1,ngrid
!      Calculating the ice cover thickness
       
        icedepth=icecover(ig)/rho_ice
       
!      If the ice cover is too thick or watercaptag=true,
!        the entire column is changed :
        IF (icedepth.ge.layer_PEM(nsoil)) THEN
          DO ik=1,nsoil
               newtherm_i(ig,ik)=max(inert_h2o_ice,prev_thermi(ig,ik))
          ENDDO 
!      We neglect the effect of a very thin ice cover :
        ELSE IF (icedepth.lt.layer_PEM(1)) THEN
          DO ik=1,nsoil
               newtherm_i(ig,ik)=inertiedat_PEM(ig,ik)
          ENDDO 
        ELSE
!        Ice/regolith boundary index :
          iref=1
!        Otherwise, we find the ice/regolith boundary:
          DO ik=1,nsoil-1
              IF ((icedepth.ge.layer_PEM(ik)).and. (icedepth.lt.layer_PEM(ik+1))) THEN
                  iref=ik+1
                  EXIT
              ENDIF
          ENDDO
!        And we change the thermal inertia:
          DO ik=1,iref-1
            newtherm_i(ig,ik)=max(inert_h2o_ice,prev_thermi(ig,ik))
          ENDDO
!        Transition (based on the equations of thermal conduction):
          newtherm_i(ig,iref)=sqrt( (layer_PEM(iref)-layer_PEM(iref-1)) / &
           ( ((icedepth-layer_PEM(iref-1))/newtherm_i(ig,iref-1)**2) + &
             ((layer_PEM(iref)-icedepth)/inertiedat_PEM(ig,ik)**2) ) )
!        Underlying regolith:
          DO ik=iref+1,nsoil
              newtherm_i(ig,ik)=inertiedat_PEM(ig,ik)
          ENDDO
        ENDIF ! icedepth
      ENDDO ! ig

!=======================================================================
      RETURN
      END