Index: /trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90	(revision 3526)
+++ /trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90	(revision 3527)
@@ -1,2 +1,8 @@
+MODULE dyn_ss_ice_m_mod
+
+implicit none
+
+CONTAINS
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -15,6 +21,11 @@
 ! orbital elements remain constant
 !***********************************************************************
-  use miscparameters, only : pi, d2r, NMAX, marsDay, solsperyear 
-  !use allinterfaces
+  use constants_marspem_mod, only: sec_per_sol
+  use fast_subs_mars, only: psv, icelayer_mars, NMAX
+#ifndef CPP_STD
+    use comcstfi_h,   only: pi
+#else
+    use comcstfi_mod, only: pi
+#endif
   implicit none
   integer, parameter :: NP=1   ! # of sites
@@ -105,5 +116,5 @@
  !    print *,'  latitude (deg)',latitude(k),' rho*c (J/m^3/K)',rhoc(k),' thIn=',thIn(k)
  !    print *,'  total pressure=',p0(k),'partial pressure=',pfrost(k)
-     delta = thIn(k)/rhoc(k)*sqrt(marsDay/pi)
+     delta = thIn(k)/rhoc(k)*sqrt(sec_per_sol/pi)
  !    print *,'  skin depths (m)',delta,delta*sqrt(solsperyear)
      call soilthprop(porosity,1.d0,rhoc(k),thIn(k),1,newrhoc,newti,icefrac)
@@ -170,2 +181,113 @@
   
 end subroutine dyn_ss_ice_m
+
+!=======================================================================
+
+subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, &
+     &     newrhoc,newti,icefrac)
+!***********************************************************************
+! soilthprop: assign thermal properties of icy soil or dirty ice
+!
+!     porositiy = void space / total volume
+!     rhof = density of free ice in space not occupied by regolith [kg/m^3]
+!     fill = rhof/icedensity <=1 (only relevant for layertype 1)
+!     rhocobs = heat capacity per volume of dry regolith [J/m^3]
+!     tiobs = thermal inertia of dry regolith [SI-units]
+!     layertype: 1=interstitial ice, 2=pure ice or ice with dirt
+!                3=pure ice, 4=ice-cemented soil, 5=custom values
+!     icefrac: fraction of ice in icelayer (only relevant for layertype 2)
+!     output are newti and newrhoc
+!***********************************************************************
+  implicit none
+  integer, intent(IN) :: layertype
+  real(8), intent(IN) :: porosity, fill, rhocobs, tiobs
+  real(8), intent(OUT) :: newti, newrhoc
+  real(8), intent(IN) :: icefrac
+  real(8) kobs, cice, icedensity, kice
+  !parameter (cice=2000.d0, icedensity=926.d0, kice=2.4d0) ! unaffected by scaling
+  parameter (cice=1540.d0, icedensity=927.d0, kice=3.2d0) ! at 198 Kelvin
+  real(8) fA, ki0, ki, k
+  real(8), parameter :: kw=3. ! Mellon et al., JGR 102, 19357 (1997)
+
+  kobs = tiobs**2/rhocobs
+  ! k, rhoc, and ti are defined in between grid points
+  ! rhof and T are defined on grid points
+
+  newrhoc = -9999.
+  newti  = -9999.
+
+  select case (layertype)
+  case (1) ! interstitial ice
+     newrhoc = rhocobs + porosity*fill*icedensity*cice
+     if (fill>0.) then
+        !--linear addition (option A)
+        k = porosity*fill*kice + kobs
+        !--Mellon et al. 1997 (option B)
+        ki0 = porosity/(1/kobs-(1-porosity)/kw)
+        fA = sqrt(fill)
+        ki = (1-fA)*ki0 + fA*kice
+        !k = kw*ki/((1-porosity)*ki+porosity*kw)
+     else
+        k = kobs
+     endif
+     newti = sqrt(newrhoc*k)
+     
+  case (2)  ! massive ice (pure or dirty ice)
+     newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice
+     k = icefrac*kice + (1.-icefrac)*kw
+     newti = sqrt(newrhoc*k)
+  
+  case (3)  ! all ice, special case of layertype 2, which doesn't use porosity
+     newrhoc = icedensity*cice
+     k = kice 
+     newti = sqrt(newrhoc*k)
+  
+  case (4)  ! pores completely filled with ice, special case of layertype 1
+     newrhoc = rhocobs + porosity*icedensity*cice
+     k = porosity*kice + kobs ! option A, end-member case of type 1, option A 
+     !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average
+     newti = sqrt(newrhoc*k)
+
+  case (5)  ! custom values
+     ! values from Mellon et al. (2004) for ice-cemented soil
+     newrhoc = 2018.*1040.
+     k = 2.5
+     newti = sqrt(newrhoc*k)
+
+  case default
+     error stop 'invalid layer type'
+     
+  end select
+  
+end subroutine soilthprop
+
+
+!=======================================================================
+      
+      real*8 function frostpoint(p)
+!     inverse of psv
+!     input is partial pressure [Pascal]
+!     output is temperature [Kelvin]
+      implicit none
+      real*8 p
+      
+!-----inverse of parametrization 1
+!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
+!      parameter (DHmelt=6008.,DHvap=45050.)
+!      parameter (DHsub=DHmelt+DHvap)
+!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
+!      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
+      
+!-----inverse of parametrization 2
+!     inverse of eq. (2) in Murphy & Koop (2005)
+      real*8 A,B
+      parameter (A=-6143.7, B=28.9074)
+      frostpoint = A / (log(p) - B)
+
+!-----approximate inverse of parametrization 3
+!     eq. (8) in Murphy & Koop (2005)
+!      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
+      
+      end function frostpoint
+
+END MODULE dyn_ss_ice_m_mod
Index: unk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90	(revision 3526)
+++ 	(revision )
@@ -1,125 +1,0 @@
-module allinterfaces
-  ! interfaces from Fortran 90 subroutines and functions
-  ! comments have been stripped
-
-  !begin fast_subs_mars.f90
-
-  interface
-     subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, &
-          & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, &
-          & albedo,p0,icefrac,zdepthT,avrho1, &
-          & avrho1prescribed)
-       use miscparameters, only : d2r, NMAX, icedensity
-       implicit none
-       integer, intent(IN) :: nz, NP
-       real(8), intent(IN) :: bigstep
-       real(8), intent(IN), dimension(NP) :: thIn, rhoc
-       real(8), intent(IN) :: z(NMAX)
-       real(8), intent(IN) :: porosity, pfrost(NP)
-       real(8), intent(INOUT) :: Tb(NP), porefill(nz,NP), zdepthF(NP), zdepthT(NP)
-       real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG
-       real(8), intent(IN), dimension(NP) ::  albedo, p0
-       real(8), intent(IN) :: icefrac
-       real(8), intent(OUT) :: avrho1(NP)
-       real(8), intent(IN), optional :: avrho1prescribed(NP)
-     end subroutine icelayer_mars
-  end interface
-
-  interface 
-     subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, &
-          &  fracIR, fracDust, p0,  avdrho, avdrhoP, avrho1, &
-          &  Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, &
-          &  B, typeG, zdepthG, avrho1prescribed)
-       use miscparameters
-       implicit none
-       integer, intent(IN) :: nz, typeT 
-       real(8), intent(IN) :: albedo0, pfrost, z(NMAX)
-       real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, p0
-       real(8), intent(IN) ::  porefill(nz)
-       real(8), intent(OUT) :: avdrho, avdrhoP 
-       real(8), intent(OUT) :: avrho1
-       real(8), intent(INOUT) :: Tb, Tmean1
-       integer, intent(OUT) :: typeF 
-       real(8), intent(OUT) :: zdepthE, zdepthF 
-       real(8), intent(OUT) :: ypp(nz) 
-       real(8), intent(OUT) :: Tmean3, zdepthG
-       real(8), intent(IN) :: B
-       integer, intent(OUT) :: typeG
-       real(8), intent(IN), optional :: avrho1prescribed
-     end subroutine ajsub_mars
-  end interface
-  
-  !end of fast_subs_mars.f90
-  !begin fast_subs_univ.f90
-
-  interface
-     pure function zint(y1,y2,z1,z2)
-       implicit none
-       real(8), intent(IN) :: y1,y2,z1,z2
-       real(8) zint
-     end function zint
-  end interface
-
-  interface
-     pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1)
-       implicit none
-       integer, intent(IN) :: nz
-       real(8), intent(IN) :: z(nz), rhosatav(nz)
-       real(8), intent(IN) :: rhosatav0, avrho1
-       real(8) zint, equildepth
-       external zint
-     end function equildepth
-  end interface
-
-  interface
-     subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1,  &
-          & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG)
-       use miscparameters, only : icedensity
-       implicit none
-       integer, intent(IN) :: nz, typeT
-       real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill
-       real(8), intent(IN) :: rhosatav0, rlow, avrho1
-       integer, intent(OUT) :: typeF
-       real(8), intent(INOUT) :: zdepthF
-       real(8), intent(IN) :: B 
-       real(8), intent(OUT) :: ypp(nz), zdepthG
-       integer, intent(INOUT) :: typeG
-       real(8), external :: zint
-     end subroutine depths_avmeth
-  end interface
-
-  interface
-     pure function constriction(porefill)
-       implicit none
-       real(8), intent(IN) :: porefill
-       real(8) constriction
-     end function constriction
-  end interface
-
-  interface
-     pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, &
-          & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG)
-       implicit none
-       integer, intent(IN) :: nz, typeF, typeG
-       real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP
-       real(8), intent(IN) :: Diff, porosity, icefrac, bigstep
-       real(8), intent(INOUT) :: zdepthT, porefill(nz)
-     end subroutine icechanges
-  end interface
-
-  !end of fast_subs_univ
-  ! Other
-  interface
-   subroutine dyn_ss_ice_m(ssi_depth_in,T1,T_in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)
-           implicit none
-           integer, intent(IN) :: nsoil
-           real(8),  intent(IN) :: thIn,ssi_depth_in,T1 
-           real(8),  intent(IN) :: p0(1), pfrost(1)
-           real(8),  intent(IN) :: T_in(nsoil)
-           real(8), intent(INOUT) :: porefill(nsoil,1)
-           real(8), intent(INOUT) :: porefill_in(nsoil,1)
-           real(8), intent(INOUT) :: ssi_depth
-   end subroutine dyn_ss_ice_m
-  end interface
-
-end module allinterfaces
Index: /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90	(revision 3526)
+++ /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90	(revision 3527)
@@ -1,3 +1,6 @@
-module thermalmodelparam_mars
+MODULE fast_subs_mars
+
+implicit none
+
   ! parameters for thermal model
   ! they are only used in the subroutines below
@@ -8,7 +11,8 @@
   real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
   integer, parameter :: EQUILTIME =15 ! [Mars years]
-end module thermalmodelparam_mars
-
-
+  
+  integer, parameter :: NMAX = 1000
+
+CONTAINS
 !*****************************************************
 ! Subroutines for fast method
@@ -25,6 +29,6 @@
 ! latitude  [degree]
 !*************************************************************************
-  use miscparameters, only : d2r, NMAX, icedensity
-  use allinterfaces, except_this_one => icelayer_mars
+  use ice_table_mod, only: rho_ice
+  use fast_subs_univ, only: icechanges
   !use omp_lib
   implicit none
@@ -53,5 +57,6 @@
      Diff = 4e-4*600./p0(k)
      fracIR = 0.04*p0(k)/600.; fracDust = 0.02*p0(k)/600.
-     B = Diff*bigstep*86400.*365.24/(porosity*icedensity)
+     B = Diff*bigstep*86400.*365.24/(porosity*927.)
+     !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'h2o'))
      
      typeT = -9
@@ -140,7 +145,6 @@
 !  Tb>0 initializes temperatures with Tb 
 !***********************************************************************
-  use miscparameters
-  use thermalmodelparam_mars
-  use allinterfaces, except_this_one => ajsub_mars
+  use fast_subs_univ, only: depths_avmeth, equildepth
+  use constants_marspem_mod, only: sols_per_my, sec_per_sol
   implicit none
   integer, intent(IN) :: nz, typeT
@@ -166,8 +170,7 @@
   real(8) Tsurf, Tco2frost, albedo, Fsurf, m, dE, emiss, T(NMAX)
   real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2
-  real(8) rhosatav0, rhosatav(nz), rlow 
-  real(8), external :: psv, tfrostco2
+  real(8) rhosatav0, rhosatav(nz), rlow
   
-  tmax = EQUILTIME*solsperyear
+  tmax = EQUILTIME*sols_per_my
   nsteps = int(tmax/dt)     ! calculate total number of timesteps
 
@@ -207,5 +210,5 @@
   do n=0,nsteps-1
      time = (n+1)*dt         !   time at n+1 
-     tdays = time*(marsDay/earthDay) ! parenthesis may improve roundoff
+   !  tdays = time*(sec_per_sol/earthDay) ! parenthesis may improve roundoff
    !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
    !  HA = 2.*pi*mod(time,1.d0)  ! hour angle
@@ -217,15 +220,15 @@
      Told(1:nz) = T(1:nz)
      if (m<=0. .or. Tsurf>Tco2frost) then
-       ! call conductionQ(nz,z,dt*marsDay,Qn,Qnp1,T,ti,rhocv,emiss, &
+       ! call conductionQ(nz,z,dt*sec_per_sol,Qn,Qnp1,T,ti,rhocv,emiss, &
        !      &           Tsurf,Fgeotherm,Fsurf)
      endif
      if (Tsurf<Tco2frost .or. m>0.) then ! CO2 condensation
         T(1:nz) = Told(1:nz)
-        !call conductionT(nz,z,dt*marsDay,T,Tsurfold,Tco2frost,ti, &
+        !call conductionT(nz,z,dt*sec_per_sol,T,Tsurfold,Tco2frost,ti, &
              !&              rhocv,Fgeotherm,Fsurf) 
         Tsurf = Tco2frost
     !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
     !         &      emiss*sigSB*(Tsurfold**4+Tsurf**4)/2.
-        m = m + dt*marsDay*dE/Lco2frost
+        m = m + dt*sec_per_sol*dE/Lco2frost
      endif
      if (Tsurf>Tco2frost .or. m<=0.) then
@@ -238,5 +241,5 @@
      !Qn=Qnp1
      
-     if (time>=tmax-solsperyear) then
+     if (time>=tmax-sols_per_my) then
         Tmean1 = Tmean1 + Tsurf
         Tmean3 = Tmean3 + T(nz)
@@ -301,2 +304,33 @@
 tfrostco2 = 3182.48/(23.3494+log(100./p))
 end function
+
+!=======================================================================
+
+      real*8 function psv(T)
+!     saturation vapor pressure of H2O ice [Pascal]
+!     input is temperature [Kelvin]
+      implicit none
+      real*8 T
+
+!-----parametrization 1
+!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
+!      parameter (DHmelt=6008.,DHvap=45050.)
+!      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
+!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
+!      C = (DHsub/R)*(1./T - 1./Tt)
+!      psv = pt*exp(-C)
+
+!-----parametrization 2
+!     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
+!     differs from parametrization 1 by only 0.1%
+      real*8 A,B
+      parameter (A=-6143.7, B=28.9074)
+      psv = exp(A/T+B)  ! Clapeyron
+
+!-----parametrization 3      
+!     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
+!     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
+      
+      end function psv
+
+END MODULE fast_subs_mars
Index: /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90	(revision 3526)
+++ /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90	(revision 3527)
@@ -1,2 +1,8 @@
+MODULE fast_subs_univ
+
+implicit none
+
+CONTAINS
+
 !*****************************************************
 ! Commonly used subroutines for fast method
@@ -19,5 +25,4 @@
 !  this is not the true (final) equilibrium depth
 !***********************************************************************
-  use allinterfaces, except_this_one => equildepth
   implicit none
   integer, intent(IN) :: nz
@@ -54,6 +59,6 @@
 !  B = Diff*bigstep/(porosity*icedensity)  [SI units]
 !***********************************************************************
-  use allinterfaces, except_this_one => depths_avmeth
   use math_mod, only: deriv2_simple, deriv1_onesided, deriv1, colint
+  use ice_table_mod, only: constriction
   implicit none
   integer, intent(IN) :: nz, typeT
@@ -178,21 +183,4 @@
 
 
-pure function constriction(porefill)
-! specify constriction function here, 0<=eta<=1
-  implicit none
-  real(8), intent(IN) :: porefill
-  real(8) eta, constriction
-  if (porefill<=0.) eta = 1.
-  if (porefill>0. .and. porefill<1.) then
-     ! eta = 1.
-     ! eta = 1-porefill
-     eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
-  endif
-  if (porefill>=1.) eta = 0.
-  constriction = eta
-end function constriction
-
-
-
 pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, &
      & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG)
@@ -202,7 +190,6 @@
 ! a crucial subroutine
 !***********************************************************
-  use miscparameters, only : icedensity
-  use allinterfaces, except_this_one => icechanges
   use math_mod, only: colint
+  use ice_table_mod, only: rho_ice
   implicit none
   integer, intent(IN) :: nz, typeF, typeG
@@ -213,5 +200,6 @@
   real(8) B, beta, integ
 
-  B = Diff*bigstep*86400.*365.24/(porosity*icedensity)
+  B = Diff*bigstep*86400.*365.24/(porosity*927.)
+  !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(T,'h2o'))
 
   ! advance ice table, avdrho>0 is retreat
@@ -232,5 +220,6 @@
      if (typeP==typeT) then   ! new 2011-09-01
         beta = (1-icefrac)/(1-porosity)/icefrac
-        beta = Diff*bigstep*beta*86400*365.24/icedensity
+        beta = Diff*bigstep*beta*86400*365.24/927.
+        !beta = Diff*bigstep*beta*86400*365.24/rho_ice(T,'h2o')
         zdepthT = sqrt(2*beta*avdrho*18./8314. + zdepthT**2)
      endif
@@ -283,2 +272,4 @@
   end if
 end subroutine icechanges
+
+END MODULE fast_subs_univ
Index: unk/LMDZ.COMMON/libf/evolution/NS_misc_params.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_misc_params.F90	(revision 3526)
+++ 	(revision )
@@ -1,13 +1,0 @@
-module miscparameters
-  ! miscellaneous parameters that are very constant
-  real(8), parameter :: pi=3.1415926535897932, d2r=pi/180.
-  integer, parameter :: NMAX=1000
-  real(8), parameter :: marsDay=88775.244, solsperyear=668.60
-  real(8), parameter :: icedensity=927.
-  real(8), parameter :: earthDay=86400.
-  real(8), parameter :: sigSB=5.6704e-8
-  ! for reference here are some parameters that are hard coded
-  !   mass of H2O molecule = 18
-  !   universal gas constant = 8314
-  !   length of Earth year in days = 365.24
-end module miscparameters
Index: unk/LMDZ.COMMON/libf/evolution/NS_psv.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_psv.F90	(revision 3526)
+++ 	(revision )
@@ -1,54 +1,0 @@
-      real*8 function psv(T)
-!     saturation vapor pressure of H2O ice [Pascal]
-!     input is temperature [Kelvin]
-      implicit none
-      real*8 T
-
-!-----parametrization 1
-!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
-!      parameter (DHmelt=6008.,DHvap=45050.)
-!      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
-!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
-!      C = (DHsub/R)*(1./T - 1./Tt)
-!      psv = pt*exp(-C)
-
-!-----parametrization 2
-!     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
-!     differs from parametrization 1 by only 0.1%
-      real*8 A,B
-      parameter (A=-6143.7, B=28.9074)
-      psv = exp(A/T+B)  ! Clapeyron
-
-!-----parametrization 3      
-!     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
-!     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
-      
-      end
-
-
-      
-      real*8 function frostpoint(p)
-!     inverse of psv
-!     input is partial pressure [Pascal]
-!     output is temperature [Kelvin]
-      implicit none
-      real*8 p
-      
-!-----inverse of parametrization 1
-!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
-!      parameter (DHmelt=6008.,DHvap=45050.)
-!      parameter (DHsub=DHmelt+DHvap)
-!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
-!      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
-      
-!-----inverse of parametrization 2
-!     inverse of eq. (2) in Murphy & Koop (2005)
-      real*8 A,B
-      parameter (A=-6143.7, B=28.9074)
-      frostpoint = A / (log(p) - B)
-
-!-----approximate inverse of parametrization 3
-!     eq. (8) in Murphy & Koop (2005)
-!      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
-      
-      end
Index: unk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90	(revision 3526)
+++ 	(revision )
@@ -1,77 +1,0 @@
-subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, &
-     &     newrhoc,newti,icefrac)
-!***********************************************************************
-! soilthprop: assign thermal properties of icy soil or dirty ice
-!
-!     porositiy = void space / total volume
-!     rhof = density of free ice in space not occupied by regolith [kg/m^3]
-!     fill = rhof/icedensity <=1 (only relevant for layertype 1)
-!     rhocobs = heat capacity per volume of dry regolith [J/m^3]
-!     tiobs = thermal inertia of dry regolith [SI-units]
-!     layertype: 1=interstitial ice, 2=pure ice or ice with dirt
-!                3=pure ice, 4=ice-cemented soil, 5=custom values
-!     icefrac: fraction of ice in icelayer (only relevant for layertype 2)
-!     output are newti and newrhoc
-!***********************************************************************
-  implicit none
-  integer, intent(IN) :: layertype
-  real(8), intent(IN) :: porosity, fill, rhocobs, tiobs
-  real(8), intent(OUT) :: newti, newrhoc
-  real(8), intent(IN) :: icefrac
-  real(8) kobs, cice, icedensity, kice
-  !parameter (cice=2000.d0, icedensity=926.d0, kice=2.4d0) ! unaffected by scaling
-  parameter (cice=1540.d0, icedensity=927.d0, kice=3.2d0) ! at 198 Kelvin
-  real(8) fA, ki0, ki, k
-  real(8), parameter :: kw=3. ! Mellon et al., JGR 102, 19357 (1997)
-
-  kobs = tiobs**2/rhocobs
-  ! k, rhoc, and ti are defined in between grid points
-  ! rhof and T are defined on grid points
-
-  newrhoc = -9999.
-  newti  = -9999.
-
-  select case (layertype)
-  case (1) ! interstitial ice
-     newrhoc = rhocobs + porosity*fill*icedensity*cice
-     if (fill>0.) then
-        !--linear addition (option A)
-        k = porosity*fill*kice + kobs
-        !--Mellon et al. 1997 (option B)
-        ki0 = porosity/(1/kobs-(1-porosity)/kw)
-        fA = sqrt(fill)
-        ki = (1-fA)*ki0 + fA*kice
-        !k = kw*ki/((1-porosity)*ki+porosity*kw)
-     else
-        k = kobs
-     endif
-     newti = sqrt(newrhoc*k)
-     
-  case (2)  ! massive ice (pure or dirty ice)
-     newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice
-     k = icefrac*kice + (1.-icefrac)*kw
-     newti = sqrt(newrhoc*k)
-  
-  case (3)  ! all ice, special case of layertype 2, which doesn't use porosity
-     newrhoc = icedensity*cice
-     k = kice 
-     newti = sqrt(newrhoc*k)
-  
-  case (4)  ! pores completely filled with ice, special case of layertype 1
-     newrhoc = rhocobs + porosity*icedensity*cice
-     k = porosity*kice + kobs ! option A, end-member case of type 1, option A 
-     !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average
-     newti = sqrt(newrhoc*k)
-
-  case (5)  ! custom values
-     ! values from Mellon et al. (2004) for ice-cemented soil
-     newrhoc = 2018.*1040.
-     k = 2.5
-     newti = sqrt(newrhoc*k)
-
-  case default
-     error stop 'invalid layer type'
-     
-  end select
-  
-end subroutine soilthprop
Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3526)
+++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3527)
@@ -495,2 +495,8 @@
 == 19/11/2024 == JBC
 Removing unused or redundant Norbert Schorghofer's subroutines (follow-up of r3493) + cleaning and some modifications of related subroutines.
+
+== 20/11/2024 == JBC
+- Removing redundant Norbert Schorghofer's subroutines/parameters (follow-up of r3526);
+- Making all Norbert Schorghofer's subroutines with modern explicit interface via modules;
+- Cleaning of "glaciers_mod.F90";
+- Optimization for the computation of ice density according to temperature by using a function.
Index: /trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3526)
+++ /trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3527)
@@ -2,5 +2,4 @@
 
 implicit none
-
 
 ! Flags for ice management
@@ -33,28 +32,25 @@
 implicit none
 
-! arguments
-! ---------
-
 ! Inputs:
-      integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
-      real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
-      real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
-      real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
-      real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
-      real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
-      real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
-     
+!--------
+integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
+real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
+real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
+real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
+real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
+real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
+real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
+
 ! Ouputs:
-      real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
-      real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
-      real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
-
+!--------
+real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
+real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
+real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
 ! Local 
-      real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
-      real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
-
-!-----------------------------
+!------
+real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
+real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
+
 write(*,*) "Flow of CO2 glaciers"
-    
 call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
 call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
@@ -83,17 +79,16 @@
 
 ! Inputs:
-      integer,intent(in) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
-      real,intent(in) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
-      real,intent(in) :: Tice(ngrid,nslope) ! Ice Temperature [K]
+integer,                       intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
+real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Physical points x Slopes : Distribution of the subgrid slopes
+real, dimension(ngrid),        intent(in) :: def_slope_mean ! Slopes: values of the sub grid slope angles
+real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice Temperature [K]
 ! Ouputs:
-      real,intent(inout) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
-      real,intent(inout) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
-      real,intent(inout) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
+real, dimension(ngrid,nslope), intent(inout) :: h2oice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
+real, dimension(ngrid,nslope), intent(inout) :: flag_h2oflow      ! flag to see if there is flow on the subgrid slopes
+real, dimension(ngrid),        intent(inout) :: flag_h2oflow_mesh ! same but within the mesh
 ! Local 
-      real :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
-
-!-----------------------------
+real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2  glacier before flow 
+
 write(*,*) "Flow of H2O glaciers"
-
 call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
 call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
@@ -105,4 +100,5 @@
 SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
 
+use ice_table_mod, only: rho_ice
 use abort_pem_mod, only: abort_pem
 #ifndef CPP_STD
@@ -114,6 +110,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
-!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle
-!!!          before initating flow
+!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow
 !!!          
 !!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
@@ -123,51 +118,40 @@
 implicit none
 
-! arguments
-! --------
-
 ! Inputs
-      integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes
-      integer,intent(in) :: iflat        ! index of the flat subslope
-      real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
-      real,intent(in) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
-      character(len=3), intent(in) :: name_ice ! Nature of the ice
+! ------
+integer,                       intent(in) :: ngrid, nslope  ! # of grid points and subslopes
+integer,                       intent(in) :: iflat          ! index of the flat subslope
+real, dimension(nslope),       intent(in) :: def_slope_mean ! Slope field: Values of the subgrid slope angles [deg]
+real, dimension(ngrid,nslope), intent(in) :: Tice           ! Physical field:  ice temperature [K]
+character(3),                  intent(in) :: name_ice       ! Nature of ice
 ! Outputs
-      real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
+! -------
+real, dimension(ngrid,nslope), intent(out) :: hmax ! Physical grid x Slope field: maximum  thickness before flaw [m]
 ! Local
-      DOUBLE PRECISION :: tau_d    ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
-      real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
-      integer :: ig,islope ! loop variables
-      real :: slo_angle
-
-! 1. Compute rho
-    if(name_ice.eq."co2") then
-      DO ig = 1,ngrid
-        DO islope = 1,nslope
-        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
+! -----
+real    :: tau_d      ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
+integer :: ig, islope ! loop variables
+real    :: slo_angle
+
+select case (trim(adjustl(name_ice)))
+    case('h2o')
+        tau_d = 1.e5
+    case('co2')
         tau_d = 5.e3
-        ENDDO
-      ENDDO
-    elseif (name_ice.eq."h2o") then
-      DO ig = 1,ngrid
-        DO islope = 1,nslope
-          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
-          tau_d = 1.e5
-        ENDDO
-      ENDDO
-    else
-      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
-    endif
-
-! 3. Compute max thickness
-    DO ig = 1,ngrid
-       DO islope = 1,nslope
-          if(islope.eq.iflat) then
+    case default
+        call abort_pem("compute_hmaxglaciers","Type of ice not known!",1)
+end select
+
+do ig = 1,ngrid
+    do islope = 1,nslope
+        if (islope == iflat) then
             hmax(ig,islope) = 1.e8
-          else
+        else
             slo_angle = abs(def_slope_mean(islope)*pi/180.)
-            hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle)
-          endif
-       ENDDO
-    ENDDO
+            hmax(ig,islope) = tau_d/(rho_ice(Tice(ig,islope),name_ice)*g*slo_angle)
+        endif
+    enddo
+enddo
+
 END SUBROUTINE compute_hmaxglaciers
 
@@ -183,4 +167,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+use ice_table_mod, only: rho_ice
 use abort_pem_mod, only: abort_pem
 #ifndef CPP_STD
@@ -192,80 +177,52 @@
 implicit none
 
-! arguments
-! --------
-
 ! Inputs
-      integer, intent(in) :: ngrid,nslope               !# of physical points and subslope
-      integer, intent(in) :: iflat                      ! index of the flat subslope
-      real, intent(in) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
-      real, intent(in) :: def_slope_mean(nslope)        ! values of the subgrid slopes
-      real, intent(in) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
-      real, intent(in) :: Tice(ngrid,nslope)            ! Ice temperature[K]
-      character(len=3), intent(in) :: name_ice              ! Nature of the ice
-
+!-------
+integer,                       intent(in) :: ngrid, nslope  ! # of physical points and subslope
+integer,                       intent(in) :: iflat          ! index of the flat subslope
+real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Distribution of the subgrid slopes within the mesh
+real, dimension(nslope),       intent(in) :: def_slope_mean ! values of the subgrid slopes
+real, dimension(ngrid,nslope), intent(in) :: hmax           ! maximum height of the  glaciers before initiating flow [m]
+real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice temperature[K]
+character(3),                  intent(in) :: name_ice       ! Nature of the ice
 ! Outputs
-      real, intent(inout) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
-      real, intent(inout) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
-      real, intent(inout) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
+!--------
+real, dimension(ngrid,nslope), intent(inout) :: qice          ! CO2 in the subslope [kg/m^2]
+real, dimension(ngrid,nslope), intent(inout) :: flag_flow     ! boolean to check if there is flow on a subgrid slope
+real, dimension(ngrid),        intent(inout) :: flag_flowmesh ! boolean to check if there is flow in the mesh
 ! Local
-      integer ig,islope ! loop
-      real rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
-      integer iaval ! ice will be transfered here
-
-! 0. Compute rho
-    if(name_ice.eq."co2") then
-      DO ig = 1,ngrid
-        DO islope = 1,nslope
-        rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
-        ENDDO
-      ENDDO
-    elseif (name_ice.eq."h2o") then
-      DO ig = 1,ngrid
-        DO islope = 1,nslope
-          rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
-        ENDDO
-      ENDDO
-    else
-      call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
-    endif
-
-! 1. Compute the transfer of ice
-
-       DO ig = 1,ngrid
-        DO islope = 1,nslope
-          IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
+!------
+integer :: ig, islope ! loop
+integer :: iaval      ! ice will be transfered here
+
+do ig = 1,ngrid
+    do islope = 1,nslope
+        if (islope /= iflat) then ! ice can be infinite on flat ground
 ! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
-            IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * &
-                  cos(pi*def_slope_mean(islope)/180.)) THEN
+            if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) then
 ! Second: determine the flatest slopes possible:
-                IF(islope.gt.iflat) THEN
-                  iaval=islope-1
-                ELSE
-                 iaval=islope+1
-                ENDIF
-                do while ((iaval.ne.iflat).and.  &
-                    (subslope_dist(ig,iaval).eq.0))
-                  IF(iaval.gt.iflat) THEN
-                     iaval=iaval-1
-                  ELSE
-                     iaval=iaval+1
-                  ENDIF
+                if (islope > iflat) then
+                    iaval=islope-1
+                else
+                    iaval = islope + 1
+                endif
+                do while (iaval /= iflat .and. subslope_dist(ig,iaval) == 0)
+                    if (iaval > iflat) then
+                        iaval = iaval - 1
+                    else
+                        iaval = iaval + 1
+                    endif
                 enddo
-              qice(ig,iaval) = qice(ig,iaval) + &
-               (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) *     &
-               cos(pi*def_slope_mean(islope)/180.)) *             &
-               subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
-               cos(pi*def_slope_mean(iaval)/180.) /               &
-               cos(pi*def_slope_mean(islope)/180.)
-
-              qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) *        &
-               cos(pi*def_slope_mean(islope)/180.)
-
-              flag_flow(ig,islope) = 1.
-              flag_flowmesh(ig) = 1.
-            ENDIF ! co2ice > hmax
-          ENDIF ! iflat
-        ENDDO !islope 
-       ENDDO !ig
+                qice(ig,iaval) = qice(ig,iaval) + (qice(ig,islope) - rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) &
+                               *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180.)/cos(pi*def_slope_mean(islope)/180.)
+
+                qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)
+                flag_flow(ig,islope) = 1.
+                flag_flowmesh(ig) = 1.
+            endif ! co2ice > hmax
+        endif ! iflat
+    enddo !islope 
+enddo !ig
+
 END SUBROUTINE
 
@@ -281,5 +238,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
-use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2
+use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2
 
 implicit none 
@@ -302,5 +259,4 @@
 real    :: ave    ! Intermediate to compute average
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 do ig = 1,ngrid
     ave = 0
Index: /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3526)
+++ /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3527)
@@ -184,6 +184,5 @@
         do islope = 1,nslope
             call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
-            rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
-            delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
+            delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho_ice(Tice,'h2o')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
                               *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
         enddo
@@ -195,6 +194,5 @@
             do isoil = 1,nsoil
                 call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),mlayer_PEM(isoil - 1),Tice)
-                rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
-                delta_m_h2o(ig) = delta_m_h2o(ig) + rho*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
+                delta_m_h2o(ig) = delta_m_h2o(ig) + rho_ice(Tice,'h2o')*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
                                   *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
             enddo
@@ -258,10 +256,10 @@
     index_tmp = nsoilmx_PEM
     do ilay = 1,nsoilmx_PEM
-        call constriction(porefill(ilay),eta(ilay))
+        eta(ilay) = constriction(porefill(ilay))
     enddo
 else
     index_tmp = index_IS
     do ilay = 1,index_IS - 1
-        call constriction(porefill(ilay),eta(ilay))
+        eta(ilay) = constriction(porefill(ilay))
     enddo
     do ilay = index_IS,nsoilmx_PEM
@@ -350,5 +348,5 @@
 
 !-----------------------------------------------------------------------
-SUBROUTINE constriction(porefill,eta)
+FUNCTION constriction(porefill) RESULT(eta)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -363,5 +361,5 @@
 
 real, intent(in) :: porefill ! pore filling fraction
-real, intent(out) :: eta ! constriction
+real :: eta ! constriction
 
 if (porefill <= 0.) then
@@ -373,5 +371,5 @@
 endif
 
-END SUBROUTINE constriction
+END FUNCTION constriction
 
 !-----------------------------------------------------------------------
@@ -386,4 +384,5 @@
 
 use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
+use abort_pem_mod, only: abort_pem
 
 implicit none
@@ -422,5 +421,5 @@
     endif
     if (indexice < 0) then
-        call abort_physic("compute_Tice_pem","subsurface ice is below the last soil layer",1)
+        call abort_pem("compute_Tice_pem","Subsurface ice is below the last soil layer!",1)
     else
         if(indexice >= 1) then ! Linear inteprolation between soil temperature
@@ -434,3 +433,32 @@
 END SUBROUTINE compute_Tice_pem
 
+!-----------------------------------------------------------------------
+FUNCTION rho_ice(T,ice) RESULT(rho)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute subsurface ice density
+!!!
+!!! Author: JBC
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+use abort_pem_mod, only: abort_pem
+
+implicit none
+
+real,         intent(in) :: T
+character(3), intent(in) :: ice
+real :: rho
+
+select case (trim(adjustl(ice)))
+    case('h2o')
+        rho = -3.5353e-4*T**2 + 0.0351*T + 933.5030 ! Rottgers 2012
+    case('co2')
+        rho = (1.72391 - 2.53e-4*T-2.87*1e-7*T**2)*1.e3 ! Mangan et al. 2017
+    case default
+        call abort_pem("rho_ice","Type of ice not known!",1)
+end select
+
+END FUNCTION rho_ice
+
 END MODULE ice_table_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3526)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3527)
@@ -71,4 +71,5 @@
 use co2condens_mod,             only: CO2cond_ps
 use layering_mod,               only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
+use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
 
 #ifndef CPP_STD
