Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3250)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3251)
@@ -20,5 +20,6 @@
 use iostart,                  only: open_startphy, get_var, close_startphy
 use physics_distribution_mod, only: init_physics_distribution
-use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, qsoil, &
+use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, &
+                                    adsorption_soil, qsoil, igcm_h2o_ice_soil, porosity_reg, &
                                     ini_comsoil_h_slope_var, end_comsoil_h_slope_var
 use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
@@ -110,6 +111,9 @@
 ! LL: Possibility to add subsurface ice
 real    :: ice_depth          ! depth of the ice table, ice_depth < 0. means no ice
-real    :: inertieice = 2100. ! ice thermal inertia
+real    :: inertieice         ! ice thermal inertia
 integer :: iref
+
+! LL: Subsurface water ice 
+real :: rho_H2O_ice           ! Density (kg/m^3) of water ice
 
 !=======================================================================
@@ -154,4 +158,9 @@
 dtemisice(1) = 2.      ! time scale for snow metamorphism (north)
 dtemisice(2) = 2.      ! time scale for snow metamorphism (south)
+
+! Water ice properties
+! ---------------------------------------
+inertieice = 2100.
+rho_H2O_ice = 920. 
 
 ! Mesh surface (not a very usefull quantity in 1D)
@@ -654,7 +663,10 @@
 endif
 
-! Initialize soil properties and temperature
+! Initialize soil properties, temperature and content
 ! ------------------------------------------
 volcapa = 1.e6 ! volumetric heat capacity
+
+if (.not. therestartfi) qsoil = 0.
+
 
 if (.not. therestartfi) then
@@ -675,4 +687,9 @@
             inertiedat(1,1) = sqrt(layer(1)/((ice_depth/inertiedat(1,1)**2) + ((layer(1) - ice_depth)/inertieice**2)))
             inertiedat(1,2:) = inertieice
+            if(adsorption_soil) then
+            ! add subsurface ice in qsoil if one runs with adsorption
+                qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*porosity_reg*rho_H2O_ice
+                qsoil(1,2:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
+            endif
         else ! searching for the ice/regolith boundary:
             do isoil = 1,nsoil
@@ -688,4 +705,11 @@
             ! Finally, we compute the underlying ice:
             inertiedat(1,iref + 1:) = inertieice
+            
+            if(adsorption_soil) then
+            ! add subsurface ice in qsoil if one runs with adsorption
+                qsoil(1,:iref - 1,igcm_h2o_ice_soil,:) = 0.
+                qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
+                qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref+1))/(layer(iref) - layer(iref+1))*porosity_reg*rho_H2O_ice
+            endif
         endif ! (ice_depth < layer(1))
     else ! ice_depth < 0 all is set to surface thermal inertia
@@ -706,7 +730,4 @@
 flux_geo = flux_geo(1,1)
 
-! Initialize soil content
-! -----------------------
-if (.not. therestartfi) qsoil = 0.
 
 ! Initialize traceurs