Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3470)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3486)
@@ -448,2 +448,6 @@
 == 15/10/2024 == JBC
 Update and correction for the visualization of the layerings over time.
+
+== 24,29/10/2024 == EV
+we added the option to use NS dynamical subsurface ice in the model to more realisticly calculate the amount of ice in the subsurface and therfore the subsurface thermal inertia
+
Index: trunk/LMDZ.COMMON/libf/evolution/dyn_ss_ice_m.f90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/dyn_ss_ice_m.f90	(revision 3470)
+++ trunk/LMDZ.COMMON/libf/evolution/dyn_ss_ice_m.f90	(revision 3486)
@@ -138,5 +138,5 @@
  ! enddo
 
-  print *,'History begins here'
+  !print *,'History begins here'
  porefill(1:nz,1:NP) =  porefill_in(1:nz,1:NP)
   zdepthT(1:NP) = ssi_depth_in
@@ -150,6 +150,6 @@
    !  print *,'T_after=  ',Tb(:)
  !    print *,'z=  ',z(:)
-     print *,'Zt=  ',ZdepthT
-      ssi_depth=ZdepthT(1)
+ !    print *,'Zt=  ',ZdepthT
+     ssi_depth=ZdepthT(1)
     ! if (abs(mod(icetime/100.,1.d0))<1.e-3) then ! output every 1000 years
     !    do k=1,NP
@@ -163,5 +163,5 @@
   !      enddo
   !   endif
-     print *,icetime
+!     print *,icetime
      if (icetime>=tlast) exit
   enddo
Index: trunk/LMDZ.COMMON/libf/evolution/fast_modules.f90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/fast_modules.f90	(revision 3470)
+++ trunk/LMDZ.COMMON/libf/evolution/fast_modules.f90	(revision 3486)
@@ -249,3 +249,17 @@
    end subroutine dyn_ss_ice_m
   end interface
+
+    interface
+      subroutine dyn_ss_ice_m_wrapper(ngrid,nsoil,tHIn,p0,pfrost,T_in,ssi_depth_in,porefill_in,porefill,ssi_depth)
+           implicit none
+           integer, intent(IN) :: nsoil,ngrid
+           real(8),  intent(IN) :: thIn(ngrid),ssi_depth_in(ngrid)
+           real(8),  intent(IN) :: p0(ngrid), pfrost(ngrid)
+           real(8),  intent(IN) :: T_in(nsoil,ngrid)
+           real(8), intent(OUT) :: porefill(nsoil,ngrid)
+           real(8), intent(IN) :: porefill_in(nsoil,ngrid)
+           real(8), intent(OUT) :: ssi_depth(ngrid)
+   end subroutine dyn_ss_ice_m_wrapper
+  end interface
+
 end module allinterfaces
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3470)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3486)
@@ -55,5 +55,5 @@
 use recomp_orb_param_mod,       only: recomp_orb_param
 use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
-                                      ini_ice_table_porefilling, icetable_equilibrium, computeice_table_equilibrium,compute_massh2o_exchange_ssi
+                                      ini_ice_table_porefilling, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium,compute_massh2o_exchange_ssi
 use soil_thermalproperties_mod, only: update_soil_thermalproperties
 use time_phylmdz_mod,           only: daysec, dtphys
@@ -226,4 +226,8 @@
 real, dimension(:),       allocatable :: delta_h2o_icetablesublim           ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
 
+! Dynamic ice table
+real, dimension(:,:),     allocatable :: ss_ice_pf                          ! the amout of porefilling in each layer in each grid [m^3/m^3]
+real, dimension(:,:),     allocatable :: ss_ice_pf_new                      ! the amout of porefilling in each layer in each grid after the ss module[m^3/m^3]
+real, dimension(:,:),     allocatable :: porefillingice_depth_new           ! new pure ice table depth
 ! Some variables for the PEM run
 real, parameter :: year_step = 1   ! Timestep for the pem
@@ -929,4 +933,13 @@
             call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
         endif
+        ! II_d.3 Update the ice table
+        if (icetable_dynamic) then
+            write(*,*) "Compute ice table"
+            porefillingice_thickness_prev_iter = porefillingice_thickness
+            call dyn_ss_ice_m_wrapper(ngrid,nsoilmx,TI_PEM,ps,mmol(igcm_h2o_vap),tsoil_PEM,porefillingice_depth,ss_ice_pf,ss_ice_pf_new,porefillingice_depth_new)
+            
+            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+        endif
+
 ! II_d.4 Update the soil thermal properties
         call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
@@ -979,4 +992,9 @@
             call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
         endif
+        if (icetable_dynamic) then
+            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
+            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
+        endif
+
         if (soil_pem) then
             call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))