Index: /trunk/LMDZ.MARS/changelog.txt
===================================================================
--- /trunk/LMDZ.MARS/changelog.txt	(revision 3229)
+++ /trunk/LMDZ.MARS/changelog.txt	(revision 3230)
@@ -4511,2 +4511,10 @@
 Remove interactive checking with XIOS whether a field should be sent to it;
 some yet unresolved issues arise when using this in mixed MPI-OpenMP mode...
+
+== 20/02/2024 == LL
+* Move soil_tifeedback into a module waterice_tifeedback_mod.F90
+* The flag to call tifeedback has been changed from "tifeedback" to surfaceice_tifeedback for clarity
+* Add the possibility to change the thermal inertia while pore ices is forming in the soil: to do so, use the flag poreice_tifeedback. The computation is done i
+n waterice_tifeedback_mod.F90.
+For now, surfaceice_tifeedback and poreice_tifeedback can not be called together
+, we might think about how to merge later.
Index: /trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 3229)
+++ /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 3230)
@@ -12,5 +12,5 @@
      &   ,callg2d,linear,rayleigh                                       &
      &   ,scavenging,sedimentation                                      &
-     &   ,activice,water,tifeedback,microphys,supersat,caps,photochem   &
+     &   ,activice,water,microphys,supersat,caps,photochem              &
      &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds &
      &   ,co2useh2o,meteo_flux,activeco2ice,CLFvaryingCO2,spantCO2      &
@@ -18,5 +18,6 @@
      &   ,latentheat_surfwater,gwd_convective_source,startphy_file      &
      &   ,hdo,hdofrac,cst_cap_albedo,temp_dependent_m,refill_watercap   &
-     &   ,cloud_adapt_ts, callatke
+     &   ,cloud_adapt_ts,callatke,surfaceice_tifeedback                 &
+     &   , poreice_tifeedback
 !$OMP THREADPRIVATE(/callkeys_l/)
 
@@ -79,5 +80,6 @@
       logical cloud_adapt_ts ! adaptative timestep for cloud microphysics
       logical sedimentation
-      logical activice,tifeedback,supersat,caps
+      logical activice,surfaceice_tifeedback,supersat,caps
+      logical poreice_tifeedback
       logical co2clouds,co2useh2o,meteo_flux,CLFvaryingCO2,satindexco2
       logical activeco2ice
Index: /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 3229)
+++ /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 3230)
@@ -30,12 +30,12 @@
 
 ! Subsurface tracers:
-logical, save                                  :: adsorption_soil ! boolean to call adosrption (or not)
-real,    save                                  :: choice_ads      ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
-real,    save, allocatable, dimension(:,:,:,:) :: qsoil           ! subsurface tracers (kg/m^3 of regol)
-integer, parameter :: nqsoil = 3                                  ! number of subsurface tracers, only three when working with water
+logical, save                                  :: adsorption_soil      ! boolean to call adosrption (or not)
+real,    save                                  :: choice_ads           ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
+real,    save, allocatable, dimension(:,:,:,:) :: qsoil                ! subsurface tracers (kg/m^3 of regol)
+integer, parameter :: nqsoil = 3                                       ! number of subsurface tracers, only three when working with water
 integer, parameter :: igcm_h2o_vap_soil = 1
 integer, parameter :: igcm_h2o_ice_soil = 2
 integer, parameter :: igcm_h2o_vap_ads  = 3
-
+REAL, parameter :: porosity_reg  = 0.45
 !$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads)
 
Index: /trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 3229)
+++ /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 3230)
@@ -683,22 +683,22 @@
 ! thermal inertia feedback
          write(*,*) "Activate the thermal inertia feedback ?"
-         tifeedback=.false. ! default value
-         call getin_p("tifeedback",tifeedback)
-         write(*,*) " tifeedback = ",tifeedback
+         surfaceice_tifeedback=.false. ! default value
+         call getin_p("surfaceice_tifeedback",surfaceice_tifeedback)
+         write(*,*) " surfaceice_tifeedback = ",surfaceice_tifeedback
 
 ! Test of incompatibility:
 
-         if (tifeedback.and..not.water) then
-           print*,'if tifeedback is used,'
+         if (surfaceice_tifeedback.and..not.water) then
+           print*,'if surfaceice_tifeedback is used,'
            print*,'water should be used too'
            call abort_physic(modname,
-     &          "tifeedback requires water",1)
-         endif
-
-         if (tifeedback.and..not.callsoil) then
-           print*,'if tifeedback is used,'
+     &          "surfaceice_tifeedback requires water",1)
+         endif
+
+         if (surfaceice_tifeedback.and..not.callsoil) then
+           print*,'if surfaceice_tifeedback is used,'
            print*,'callsoil should be used too'
            call abort_physic(modname,
-     &          "tifeedback requires callsoil",1)
+     &          "surfaceice_tifeedback requires callsoil",1)
          endif
 
@@ -1041,5 +1041,15 @@
          if (adsorption_soil) choice_ads = 1
          call getin_p("choice_ads",choice_ads) 
-           
+         
+         poreice_tifeedback = .false.
+         call getin_p("poreice_tifeedback",poreice_tifeedback)
+         write(*,*) 'poreice_tifeedback=',poreice_tifeedback 
+         if (poreice_tifeedback .and. (.not. adsorption_soil)) then
+              write(*,*)"Pore ice TI feedback can be run only if 
+     &                   adsorption_soil = True"
+              call abort_physic(modname,
+     &        "Pore ice TI feedback be used with adsorption_soil
+     &         = true",1)
+         endif 
 c        ----------------------------------------------------------
 
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 3229)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 3230)
@@ -123,4 +123,5 @@
       use pbl_parameters_mod, only: pbl_parameters
       use lmdz_atke_turbulence_ini, only : atke_ini
+      use waterice_tifeedback_mod, only : waterice_tifeedback
       IMPLICIT NONE
 c=======================================================================
@@ -282,5 +283,5 @@
       REAL inertiesoil_tifeedback(ngrid,nsoilmx,nslope) ! Time varying subsurface
                                       ! thermal inertia (J.s-1/2.m-2.K-1)
-                                      ! (used only when tifeedback=.true.)
+                                      ! (used only when tifeedback surface or pore =.true.)
 c     Variables used by the CO2 clouds microphysical scheme:
       DOUBLE PRECISION riceco2(ngrid,nlayer)   ! co2 ice geometric mean radius (m)
@@ -540,14 +541,14 @@
 
 ! Variable for ice table
-      REAL :: rhowater_surf(ngrid,nslope)         ! Water density at the surface [kg/m^3]
-      REAL :: rhowater_surf_sat(ngrid,nslope)     ! Water density at the surface at saturation [kg/m^3]
-      REAL :: rhowater_soil(ngrid,nsoilmx,nslope) ! Water density in soil layers [kg/m^3]
-      REAL,PARAMETER :: alpha_clap_h2o = 28.9074  ! Coeff for Clapeyron law [/]
-      REAL,PARAMETER :: beta_clap_h2o = -6143.7   ! Coeff for Clapeyron law [K]
-      REAL :: pvap_surf(ngrid)                    ! Water vapor partial pressure in first layer [Pa]
-      REAL,PARAMETER :: m_co2 = 44.01E-3          ! CO2 molecular mass [kg/mol]
-      REAL,PARAMETER :: m_noco2 = 33.37E-3        ! Non condensible mol mass [kg/mol]
-      REAL :: ztmp1,ztmp2                         ! intermediate variables to compute the mean molar mass of the layer
-
+      REAL :: rhowater_surf(ngrid,nslope)            ! Water density at the surface [kg/m^3]
+      REAL :: rhowater_surf_sat(ngrid,nslope)        ! Water density at the surface at saturation [kg/m^3]
+      REAL :: rhowater_soil(ngrid,nsoilmx,nslope)    ! Water density in soil layers [kg/m^3]
+      REAL,PARAMETER :: alpha_clap_h2o = 28.9074     ! Coeff for Clapeyron law [/]
+      REAL,PARAMETER :: beta_clap_h2o = -6143.7      ! Coeff for Clapeyron law [K]
+      REAL :: pvap_surf(ngrid)                       ! Water vapor partial pressure in first layer [Pa]
+      REAL,PARAMETER :: m_co2 = 44.01E-3             ! CO2 molecular mass [kg/mol]
+      REAL,PARAMETER :: m_noco2 = 33.37E-3           ! Non condensible mol mass [kg/mol]
+      REAL :: ztmp1,ztmp2                            ! intermediate variables to compute the mean molar mass of the layer
+      REAL :: pore_icefraction(ngrid,nsoilmx,nslope) ! ice filling fraction in the pores
 ! Variable for the computation of the TKE with parameterization from ATKE working group
       REAL :: viscom                              ! kinematic molecular viscosity for momentum
@@ -694,9 +695,9 @@
          IF (callsoil) THEN
 c           Thermal inertia feedback:
-            IF (tifeedback) THEN
+            IF (surfaceice_tifeedback.or.poreice_tifeedback) THEN
              DO islope = 1,nslope
-                CALL soil_tifeedback(ngrid,nsoilmx,
-     s               qsurf(:,:,islope),
-     s               inertiesoil_tifeedback(:,:,islope))
+                CALL waterice_tifeedback(ngrid,nsoilmx,nslope,
+     s               qsurf(:,igcm_h2o_ice,:),pore_icefraction,
+     s               inertiesoil_tifeedback)
              ENDDO
                 CALL soil(ngrid,nsoilmx,firstcall,
@@ -1498,9 +1499,9 @@
 c        Calling vdif (Martian version WITH CO2 condensation)
          dwatercap_dif(:,:) = 0.
-
          CALL vdifc(ngrid,nlayer,nsoilmx,nq,nqsoil,zpopsk,
      $        ptimestep,capcal,lwrite,
      $        zplay,zplev,zzlay,zzlev,z0,
-     $        pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,qsoil,
+     $        pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,
+     $        qsoil,pore_icefraction,
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
@@ -2385,10 +2386,10 @@
       IF (callsoil) THEN
 c       Thermal inertia feedback
-        IF (tifeedback) THEN
-         DO islope = 1,nslope
-           CALL soil_tifeedback(ngrid,nsoilmx,
-     s               qsurf(:,:,islope),
-     s               inertiesoil_tifeedback(:,:,islope))
-         ENDDO
+        IF (surfaceice_tifeedback.or.poreice_tifeedback) THEN
+        
+           CALL waterice_tifeedback(ngrid,nsoilmx,nslope,
+     s               qsurf(:,igcm_h2o_ice,:),pore_icefraction,
+     s               inertiesoil_tifeedback(:,:,:))
+     
          CALL soil(ngrid,nsoilmx,.false.,inertiesoil_tifeedback,
      s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
@@ -3551,5 +3552,5 @@
      &                         '',albedo(:,1,islope))
          ENDDO
-              if (tifeedback) then
+              if (surfaceice_tifeedback.or.poreice_tifeedback) then
                  call write_output("soiltemp",
      &                              "Soil temperature","K",
@@ -3558,4 +3559,5 @@
      &                       'Soil Thermal Inertia',
      &         'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat))
+
          do islope=1,nslope
           write(str2(1:2),'(i2.2)') islope
Index: /trunk/LMDZ.MARS/libf/phymars/soil.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 3229)
+++ /trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 3230)
@@ -41,5 +41,5 @@
 
 ! 0. Initialisations and preprocessing step
-      if (firstcall.or.tifeedback) then
+      if(firstcall.or.surfaceice_tifeedback.or.poreice_tifeedback) then
       ! note: firstcall is set to .true. or .false. by the caller
       !       and not changed by soil.F 
Index: unk/LMDZ.MARS/libf/phymars/soil_tifeedback.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/soil_tifeedback.F	(revision 3229)
+++ 	(revision )
@@ -1,97 +1,0 @@
-      SUBROUTINE soil_tifeedback(ngrid,nsoil,icecover,newtherm_i)
-
-      use tracer_mod, only: nqmx, igcm_h2o_ice, rho_ice
-      use comsoil_h, only: layer, inertiedat
-      use surfdat_h, only: watercaptag, inert_h2o_ice
-      IMPLICIT NONE
-
-c=======================================================================
-c   Description :
-c       Surface water ice / Thermal inertia feedback.
-c
-c   When surface water-ice is thick enough, this routine creates a new
-c   soil thermal inertia with three different layers :
-c   - One layer of surface water ice (the thickness is given
-c     by the variable icecover (in kg of ice per m2) and the thermal
-c     inertia is prescribed by inert_h2o_ice (see surfdat_h)); 
-c   - A transitional layer of mixed thermal inertia;
-c   - A last layer of regolith below the ice cover whose thermal inertia
-c     is equal to inertiedat.
-c
-c   To use the model :
-c       SET THE tifeedback LOGICAL TO ".true." in callphys.def.
-c
-c   Author: J.-B. Madeleine Mars 2008 - Updated November 2012
-c=======================================================================
-
-c Local variables
-c ---------------
-
-      INTEGER :: ig                     ! Grid point (ngrid)
-      INTEGER :: ik                     ! Grid point (nsoil)
-      INTEGER :: iref                   ! Ice/Regolith boundary index
-      INTEGER :: ngrid                  ! Number of horizontal grid points 
-      INTEGER :: nsoil                  ! Number of soil layers 
-      REAL :: icedepth                  ! Ice cover thickness (m)
-
-c Inputs
-c ------
-
-      REAL icecover(ngrid,nqmx)         ! tracer on the surface (kg.m-2)
-                                        ! (iq=igcm_h2o_ice) is surface
-                                        ! water ice
-c Outputs
-c -------
-
-      REAL newtherm_i(ngrid,nsoil)    ! New soil thermal inertia
-
-c Initialization
-c --------------
-
-      newtherm_i(1:ngrid,1:nsoil) = 0
-
-c Creating the new soil thermal inertia table
-c -------------------------------------------
-      DO ig=1,ngrid
-c       Calculating the ice cover thickness
-        icedepth=icecover(ig,igcm_h2o_ice)/rho_ice
-c       If the ice cover is too thick or watercaptag=true,
-c         the entire column is changed :
-        IF ((icedepth.ge.layer(nsoil)).or.(watercaptag(ig))) THEN
-          DO ik=1,nsoil
-               newtherm_i(ig,ik)=inert_h2o_ice
-          ENDDO 
-c       We neglect the effect of a very thin ice cover :
-        ELSE IF (icedepth.lt.layer(1)) THEN
-          DO ik=1,nsoil
-               newtherm_i(ig,ik)=inertiedat(ig,ik)
-          ENDDO 
-        ELSE
-c         Ice/regolith boundary index :
-          iref=1
-c         Otherwise, we find the ice/regolith boundary:
-          DO ik=1,nsoil-1
-              IF ((icedepth.ge.layer(ik)).and.
-     &        (icedepth.lt.layer(ik+1))) THEN
-                  iref=ik+1
-                  EXIT
-              ENDIF
-          ENDDO
-c         And we change the thermal inertia:
-          DO ik=1,iref-1
-            newtherm_i(ig,ik)=inert_h2o_ice
-          ENDDO
-c         Transition (based on the equations of thermal conduction):
-          newtherm_i(ig,iref)=sqrt( (layer(iref)-layer(iref-1)) /
-     &        ( ((icedepth-layer(iref-1))/inert_h2o_ice**2) +
-     &        ((layer(iref)-icedepth)/inertiedat(ig,ik)**2) ) )
-c         Underlying regolith:
-          DO ik=iref+1,nsoil
-              newtherm_i(ig,ik)=inertiedat(ig,ik)
-          ENDDO
-        ENDIF ! icedepth
-      ENDDO ! ig
-
-c=======================================================================
-      RETURN
-      END
Index: /trunk/LMDZ.MARS/libf/phymars/soilwater.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/soilwater.F90	(revision 3229)
+++ /trunk/LMDZ.MARS/libf/phymars/soilwater.F90	(revision 3230)
@@ -1,8 +1,8 @@
 subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep,  & 
       exchange, qsat_surf, pq, pa, pb, pc, pd, pdqsdifpot, pqsurf,  & 
-      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2) 
-
-
-      use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads
+      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2, saturation_water_ice) 
+
+
+      use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg
       use comcstfi_h
       use tracer_mod
@@ -92,6 +92,7 @@
 ! Outputs :
 ! ---------------------------------------------------------------------- 
-real, intent(out) :: zdqsdifrego(ngrid)         ! Flux from subsurface (positive pointing outwards)
-real, intent(out) :: zq1temp2(ngrid)            ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
+real, intent(out) :: zdqsdifrego(ngrid)                     ! Flux from subsurface (positive pointing outwards)
+real, intent(out) :: zq1temp2(ngrid)                        ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
+real*8, intent(out) :: saturation_water_ice(ngrid, nsoil)   ! Water pore ice saturation level (formerly Sw)
 
 ! Outputs for the output files
@@ -162,5 +163,4 @@
 real*8, allocatable, save :: tortuosity(:, :)         ! Tortuosity factor (formerly tortuo)
 
-real*8 :: saturation_water_ice(ngrid, nsoil)          ! Water pore ice saturation level (formerly Sw)
 real*8 :: saturation_water_ice_inter(ngrid, nsoil)          ! Water pore ice saturation level at the interlayer
 
@@ -346,5 +346,5 @@
                         
                         ! These properties are defined here in order to enable custom profiles
-                        porosity_ice_free(ig, ik) = 0.45D0
+                        porosity_ice_free(ig, ik) = porosity_reg
                         tortuosity(ig, ik) = 1.5D0
                         rho_soil(ig, ik) = 1.3D3 ! in kg/m3 of regolith (incl. porosity)
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3229)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3230)
@@ -9,5 +9,6 @@
      $                pplay,pplev,pzlay,pzlev,pz0,
      $                pu,pv,ph,pq,ptsrf,ptsoil,pemis,pqsurf,qsoil,
-     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
+     $                pore_icefraction,pdufi,pdvfi,pdhfi,
+     $                pdqfi,pfluxsrf,
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
      $                pdqdif,pdqsdif,wstar,
@@ -239,5 +240,5 @@
       REAL zdqsdif_tot(ngrid)         ! subtimestep pdqsdif for water ice
       LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
-      
+      REAL, INTENT(out) :: pore_icefraction(ngrid,nsoil,nslope) ! ice filling fraction in the pores
 !! Water buyoncy
       LOGICAL :: virtual
@@ -307,5 +308,5 @@
       h2o_ice_depth(1:ngrid,1:nslope)=1
       virtual = .false.
-      
+      pore_icefraction(:,:,:) = 0.
 c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
 c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
@@ -1084,5 +1085,6 @@
      &                     qsoil(ig,:,:,islope), pplev(ig,1), rho(ig), 
      &                     writeoutput,zdqsdif_regolith(ig,islope), 
-     &                     zq1temp_regolith(ig))
+     &                     zq1temp_regolith(ig),
+     &                     pore_icefraction(ig,:,islope))
 
 
Index: /trunk/LMDZ.MARS/libf/phymars/waterice_tifeedback_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/waterice_tifeedback_mod.F90	(revision 3230)
+++ /trunk/LMDZ.MARS/libf/phymars/waterice_tifeedback_mod.F90	(revision 3230)
@@ -0,0 +1,140 @@
+MODULE waterice_tifeedback_mod
+ 
+!======================================================================================================================!
+! Subject:
+!---------
+!   Module used to compute the thermal inertia of an icy soil (either pore filling or massive ice)
+!----------------------------------------------------------------------------------------------------------------------!
+! Reference:
+!-----------
+!  For pure ice on the surface: J.-B. Madeleine, F. Forget, James W. Head, B. Levrard, F. Montmessin, E. Millour,
+!  Amazonian northern mid-latitude glaciation on Mars: A proposed climate scenario, Icarus (10.1016/j.icarus.2009.04.037).
+!
+!  For pore-filling ice: Siegler, M., O. Aharonson, E. Carey, M. Choukroun, T. Hudson, N. Schorghofer, and S. Xu (2012), Measurements of thermal properties of icy Mars regolith analogs, JGR (10.1029/2011JE003938). 
+!
+! Originally written by JBM for pure ice (2008-2012), moved in to a module and .F90 by LL (2024). Pore filling ice included by LL (2024)
+!
+!======================================================================================================================!
+
+
+
+IMPLICIT NONE
+
+CONTAINS
+
+SUBROUTINE waterice_tifeedback(ngrid,nsoil,nslope,icecover,poreice,newtherm_i)
+
+      use tracer_mod, only: rho_ice
+      use comsoil_h, only: layer, inertiedat, porosity_reg
+      use surfdat_h, only: watercaptag, inert_h2o_ice
+      IMPLICIT NONE
+      include "callkeys.h"
+!=======================================================================
+!   Description :
+!       Surface water ice / pore filling ice thermal inertia feedback.
+!
+!   When surface water-ice is thick enough (flag surfaceice_tifeedback), 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 surfaceice_tifeedback LOGICAL TO ".true." in callphys.def.
+!
+!   When pore filling ice is present, surface ice is computed as in Siegler et al., 2012
+!   \sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) 
+!
+!   For now, the code can not run with both options
+!
+!
+!   Author: J.-B. Madeleine Mars 2008 - Updated November 2012; Added porous ice by LL February 2024
+!=======================================================================
+
+! Local variables
+! ---------------
+
+      INTEGER :: ig                       ! Grid point (ngrid)
+      INTEGER :: islope                   ! SubGrid point (nslope)
+      INTEGER :: ik                       ! Grid point (nsoil)
+      INTEGER :: iref                     ! Ice/Regolith boundary index
+      REAL :: icedepth                    ! Ice cover thickness (m)
+      REAL :: inertie_purewaterice = 2100 ! Thermal inertia of pure water ice (J/m^2/K/s^1/2)
+
+! Inputs
+! ------
+      INTEGER :: ngrid                 ! Number of horizontal grid points 
+      INTEGER :: nsoil                 ! Number of soil layers 
+      INTEGER :: nslope                ! Number of subgrid slopes 
+      REAL icecover(ngrid,nslope)      ! water iceon the surface (kg.m-2)
+      REAL poreice(ngrid,nsoil,nslope) ! pore ice filling fraction (1)
+! Outputs
+! -------
+
+      REAL newtherm_i(ngrid,nsoil,nslope)  ! New soil thermal inertia
+
+! Initialization
+! --------------
+
+  newtherm_i(1:ngrid,1:nsoil,1:nslope) = 0
+  
+  IF (surfaceice_tifeedback) THEN
+      
+! Creating the new soil thermal inertia table because of the massive surface ice
+! ------------------------------------------------------------------------------
+     DO islope = 1,nslope
+        DO ig=1,ngrid
+!       Calculating the ice cover thickness
+           icedepth=icecover(ig,islope)/rho_ice
+!       If the ice cover is too thick or watercaptag=true,
+!       the entire column is changed :
+           IF ((icedepth.ge.layer(nsoil)).or.(watercaptag(ig))) THEN
+              DO ik=1,nsoil
+                 newtherm_i(ig,ik,islope)=inert_h2o_ice
+              ENDDO 
+!          We neglect the effect of a very thin ice cover :
+           ELSE IF (icedepth.lt.layer(1)) THEN
+              DO ik=1,nsoil
+                 newtherm_i(ig,ik,islope)=inertiedat(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(ik)).and.icedepth.lt.layer(ik+1)) THEN
+                    iref=ik+1
+                 EXIT
+                 ENDIF
+              ENDDO
+!             And we change the thermal inertia:
+              DO ik=1,iref-1
+                 newtherm_i(ig,ik,islope)=inert_h2o_ice
+              ENDDO
+!             Transition (based on the equations of thermal conduction):
+              newtherm_i(ig,iref,islope)=sqrt( (layer(iref)-layer(iref-1))/(((icedepth-layer(iref-1))/inert_h2o_ice**2) + &
+              ((layer(iref)-icedepth)/inertiedat(ig,ik)**2) ) )
+!              Underlying regolith:
+              DO ik=iref+1,nsoil
+                 newtherm_i(ig,ik,islope)=inertiedat(ig,ik)
+              ENDDO
+           ENDIF ! icedepth
+        ENDDO ! ig
+     ENDDO ! islope
+     
+  ELSE IF (poreice_tifeedback) THEN
+  
+! Updating soil thermal properties to includes the effect of porous ice
+! ------------------------------------------------------------------------------ 
+     DO islope = 1,nslope
+        newtherm_i(:,:,islope) = sqrt(inertiedat(:,:)**2 + porosity_reg*poreice(:,:,islope)*inertie_purewaterice**2)
+     ENDDO
+  ENDIF
+
+!=======================================================================
+
+END SUBROUTINE waterice_tifeedback
+END MODULE waterice_tifeedback_mod
