Index: trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90	(revision 3347)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90	(revision 3525)
@@ -9,4 +9,13 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! Constants:
+real, parameter :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
+real, parameter :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
+real, parameter :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
+real, parameter :: P610 = 610.0                ! current average pressure on Mars [Pa]
+real, parameter :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless]
+real, parameter :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa]
+real, parameter :: Pa2Torr = 1./133.3          ! Conversion from Pa to tor [Pa/Torr]
 
 !=======================================================================
@@ -50,5 +59,5 @@
     ice_thermalinertia = inertie_purewaterice
 else
-    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2)  ! Siegler et al., 2012
+    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) ! Siegler et al., 2012
 endif
 
@@ -56,5 +65,5 @@
 !=======================================================================
 
-SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,ice_thickness,TI_PEM)
+SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil,tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
 
 use comsoil_h,             only: volcapa
@@ -65,23 +74,15 @@
 
 ! Input:
-integer,                       intent(in) :: ngrid, nslope, nsoil_PEM ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
-real,                          intent(in) :: p_avg_new                ! Global average surface pressure [Pa]
-real, dimension(ngrid,nslope), intent(in) :: tendencies_waterice      ! Tendencies on the water ice [kg/m^2/year]
-real, dimension(ngrid,nslope), intent(in) :: waterice                 ! Surface Water ice [kg/m^2]
-real, dimension(ngrid,nslope), intent(in) :: ice_depth                ! Ice table depth [m]
-real, dimension(ngrid,nslope), intent(in) :: ice_thickness            ! Ice table thickness [m]
+integer,                             intent(in) :: ngrid, nslope, nsoil ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
+real,                                intent(in) :: p_avg_new            ! Global average surface pressure [Pa]
+real, dimension(ngrid,nslope),       intent(in) :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
+real, dimension(ngrid,nslope),       intent(in) :: waterice             ! Surface Water ice [kg/m^2]
+real, dimension(ngrid,nslope),       intent(in) :: icetable_depth       ! Ice table depth [m]
+real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness   ! Ice table thickness [m]
+real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling      ! Ice porefilling [m^3/m^3]
+logical,                             intent(in) :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
 
 ! Outputs:
-real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
-
-! Constants:
-real :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
-real :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
-real :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
-real :: ice_inertia                 ! Inertia of water ice [SI]
-real :: P610 = 610.0                ! current average pressure on Mars [Pa]
-real :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless]
-real :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa]
-real :: Pa2Tor = 1./133.3           ! Conversion from Pa to tor [Pa/tor]
+real, dimension(ngrid,nsoil,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
 
 ! Local variables:
@@ -91,4 +92,5 @@
 real    :: ice_bottom_depth                       ! Bottom depth of the subsurface ice [m]
 real    :: d_part                                 ! Regolith particle size [micrometer]
+real    :: ice_inertia                            ! Inertia of water ice [SI]
 
 write(*,*) "Update soil properties"
@@ -101,6 +103,6 @@
         if (reg_thprop_dependp) then
             if (TI_PEM(ig,1,islope) < reg_inertie_thresold) then
-                d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Tor)**(0.6)))**(-1./(0.11*log10(P610*Pa2Tor/K))) ! compute particle size, in micrometer
-                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K)))
+                d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Torr)**(0.6)))**(-1./(0.11*log10(P610*Pa2Torr/K))) ! compute particle size, in micrometer
+                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Torr/K)))
                 if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
                 if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
@@ -135,5 +137,5 @@
                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
-        do isoil = index_bedrock + 2,nsoil_PEM
+        do isoil = index_bedrock + 2,nsoil
             TI_PEM(ig,isoil,islope) = TI_bedrock
         enddo
@@ -141,76 +143,79 @@
 enddo ! islope
 
-!  3. Build new TI in case of ice table
+! 3. Build new TI in case of ice table
 do ig = 1,ngrid
     do islope = 1,nslope
-        if (ice_depth(ig,islope) > -1.e-6) then
+        if (icetable_depth(ig,islope) > -1.e-6) then
         ! 3.0 Case where it is perennial ice
-            if (ice_depth(ig,islope) < 1.e-10) then
+            if (icetable_depth(ig,islope) < 1.e-10) then
                 call ice_thermal_properties(.true.,1.,0.,ice_inertia)
-                do isoil = 1,nsoil_PEM
+                do isoil = 1,nsoil
                     TI_PEM(ig,isoil,islope) = ice_inertia
                 enddo
             else
-                call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
-        ! 3.1.1 find the index of the mixed layer
-                iref = 0 ! initialize iref
-                do isoil = 1,nsoil_PEM ! loop on layers to find the beginning of the ice table
-                    if (ice_depth(ig,islope) >= layer_PEM(isoil)) then
-                        iref = isoil ! pure regolith layer up to here
-                    else
-                        exit ! correct iref was obtained in previous cycle
-                    endif
-                enddo
-        ! 3.1.2 find the index of the end of the ice table
-                iend = 0 ! initialize iend
-                ice_bottom_depth = ice_depth(ig,islope) + ice_thickness(ig,islope)
-                do isoil = 1,nsoil_PEM ! loop on layers to find the end of the ice table
-                    if (ice_bottom_depth >= layer_PEM(isoil)) then
-                        iend = isoil ! pure regolith layer up to here
-                    else
-                        exit ! correct iref was obtained in previous cycle
-                    endif
-                enddo
-        ! 3.2 Build the new ti
-                if (iref < nsoil_PEM) then
-                    if (iref == iend) then
-                        ! Ice table begins and end in the same mixture Mixtures with three components
-                        if (iref /= 0) then ! mixed layer
-                            TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                             &
-                                                         (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
-                                                         ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) +            &
-                                                         ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
-                        else ! first layer is already a mixed layer
-                            ! (ie: take layer(iref=0)=0)
-                            TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/ &
-                                                  (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
-                                                  ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) + &
-                                                  ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
-                        endif ! of if (iref /= 0)
-                    else
-                        if (iref /= 0) then ! mixed layer
-                            TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                             &
-                                                         (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
-                                                         ((layer_PEM(iref + 1) - ice_depth(ig,islope))/(ice_inertia**2))))
-                        else ! first layer is already a mixed layer
-                            ! (ie: take layer(iref=0)=0)
-                            TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                 &
-                                                  (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
-                                                  ((layer_PEM(1) - ice_depth(ig,islope))/(ice_inertia**2))))
-                        endif ! of if (iref /= 0)
-                    endif ! iref == iend
-        ! 3.3 Build the new ti
-                    do isoil = iref + 2,iend
-                        TI_PEM(ig,isoil,islope) = ice_inertia
+                if (icetable_equilibrium) then
+                    call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
+                    ! 3.1.1 find the index of the mixed layer
+                    iref = 0 ! initialize iref
+                    do isoil = 1,nsoil ! loop on layers to find the beginning of the ice table
+                        if (icetable_depth(ig,islope) >= layer_PEM(isoil)) then
+                            iref = isoil ! pure regolith layer up to here
+                        else
+                            exit ! correct iref was obtained in previous cycle
+                        endif
                     enddo
-                    if (iend < nsoil_PEM) then
-                        TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
-                                                     (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
-                                                     ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
-                    endif
-                endif ! of if (iref < nsoilmx)
+                    ! 3.1.2 find the index of the end of the ice table
+                    iend = 0 ! initialize iend
+                    ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope)
+                    do isoil = 1,nsoil ! loop on layers to find the end of the ice table
+                        if (ice_bottom_depth >= layer_PEM(isoil)) then
+                            iend = isoil ! pure regolith layer up to here
+                        else
+                            exit ! correct iref was obtained in previous cycle
+                        endif
+                    enddo
+                    ! 3.2 Build the new ti
+                    if (iref < nsoil) then
+                        if (iref == iend) then
+                            ! Ice table begins and end in the same mixture with three components
+                            if (iref /= 0) then ! mixed layer
+                                TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
+                                                             (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
+                                                             ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) +            &
+                                                             ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
+                            else ! first layer is already a mixed layer
+                                ! (ie: take layer(iref=0)=0)
+                                TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                                &
+                                                      (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
+                                                      ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + &
+                                                      ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
+                            endif ! of if (iref /= 0)
+                        else
+                            if (iref /= 0) then ! mixed layer
+                                TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
+                                                             (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
+                                                             ((layer_PEM(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2))))
+                            else ! first layer is already a mixed layer
+                                ! (ie: take layer(iref=0)=0)
+                                TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                      &
+                                                      (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
+                                                      ((layer_PEM(1) - icetable_depth(ig,islope))/(ice_inertia**2))))
+                            endif ! of if (iref /= 0)
+                        endif ! iref == iend
+
+                        TI_PEM(ig,iref + 2:iend,islope) = ice_inertia
+                        if (iend < nsoil) then
+                            TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
+                                                         (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
+                                                         ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
+                        endif
+                    endif ! of if (iref < nsoilmx)
+                else if (icetable_dynamic) then
+                    do  isoil = 1,nsoil
+                        call ice_thermal_properties(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI_PEM(ig,isoil,islope))
+                    enddo
+                endif ! of if icetable_equilibrium
             endif ! permanent glaciers
-        endif ! ice_depth(ig,islope) > 0.
-!        write(*,*) 'TI=', TI_PEM(ig,:,islope)
+        endif ! icetable_depth(ig,islope) > 0.
     enddo !islope
 enddo !ig
