Index: /trunk/LMDZ.COMMON/libf/evolution/conductionQ.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/conductionQ.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/conductionQ.f90	(revision 3470)
@@ -0,0 +1,113 @@
+subroutine conductionQ(nz,z,dt,Qn,Qnp1,T,ti,rhoc,emiss,Tsurf, &
+     &     Fgeotherm,Fsurf)
+!***********************************************************************
+!   conductionQ:  program to calculate the diffusion of temperature 
+!                 into the ground and thermal emission at the surface 
+!                 with variable thermal properties on irregular grid
+!   Crank-Nicolson scheme, flux conservative
+!                          uses Samar's radiation formula
+!   Eqn: rhoc*T_t = (k*T_z)_z 
+!   BC (z=0): Q(t) + kT_z = em*sig*T^4
+!   BC (z=L): heat flux = Fgeotherm
+!
+!   nz = number of grid points
+!   dt = time step
+!   Qn,Qnp1 = net solar insolation at time steps n and n+1 [W/m^2]
+!   T = vertical temperature profile [K]  (in- and output)  
+!   ti = thermal inertia [J m^-2 K^-1 s^-1/2]  VECTOR
+!   rhoc = rho*c  VECTOR where rho=density [kg/m^3] and 
+!                              c=specific heat [J K^-1 kg^-1]
+!   ti and rhoc are not allowed to vary in the layers immediately 
+!               adjacent to the surface or the bottom
+!   emiss = emissivity
+!   Tsurf = surface temperature [K]  (in- and output)
+!   Fgeotherm = geothermal heat flux at bottom boundary [W/m^2]
+!   Fsurf = heat flux at surface [W/m^2]  (output)
+!
+!   Grid: surface is at z=0
+!         z(0)=0, z(2)=3*z(1), i.e., the top layer has half the width
+!         T(1) is at z(1); ...; T(i) is at z(i)
+!         k(i), rhoc(i), ti(i) are midway between z(i-1) and z(i)
+!     
+!   originally written by Samar Khatiwala, 2001
+!   extended to variable thermal properties
+!         and irregular grid by Norbert Schorghofer
+!   added predictor-corrector 9/2019  
+!***********************************************************************
+
+  implicit none
+  real*8, parameter :: sigSB=5.6704d-8
+  
+  integer, intent(IN) :: nz
+  real*8, intent(IN) :: z(nz), dt, Qn, Qnp1, ti(nz),rhoc(nz)
+  real*8, intent(IN) :: emiss, Fgeotherm
+  real*8, intent(INOUT) :: T(nz), Tsurf
+  real*8, intent(OUT) :: Fsurf
+  integer i, iter
+  real*8 k(nz), k1, alpha(nz), gamma(nz), Tr
+  real*8 a(nz), b(nz), c(nz), r(nz), Told(nz)
+  real*8 arad, brad, ann, annp1, bn, buf, dz, beta
+  
+  ! set some constants
+  k(:) = ti(:)**2/rhoc(:) ! thermal conductivity
+  dz = 2.*z(1)
+  beta = dt/rhoc(1)/(2.*dz**2)   ! assumes rhoc(0)=rhoc(1)
+  alpha(1) = beta*k(2)
+  gamma(1) = beta*k(1)
+  do i=2,nz-1
+     buf = dt/(z(i+1)-z(i-1))
+     alpha(i) = 2.*k(i+1)*buf/(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i))
+     gamma(i) = 2.*k(i)*buf/(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1))
+  enddo
+  buf = dt/(z(nz)-z(nz-1))**2
+  gamma(nz) = k(nz)*buf/(2*rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz)
+  
+  k1 = k(1)/dz
+  
+  ! elements of tridiagonal matrix
+  a(:) = -gamma(:)   !  a(1) is not used
+  b(:) = 1. + alpha(:) + gamma(:) !  b(1) has to be reset at every timestep
+  c(:) = -alpha(:)   !  c(nz) is not used
+  b(nz) = 1. + gamma(nz)
+
+  Tr = Tsurf            ! 'reference' temperature
+  iter = 0
+  Told(:) = T(:)
+30 continue  ! in rare case of iterative correction
+
+  ! Emission
+  arad = -3.*emiss*sigSB*Tr**4
+  brad = 2.*emiss*sigSB*Tr**3
+  ann = (Qn-arad)/(k1+brad)
+  annp1 = (Qnp1-arad)/(k1+brad)
+  bn = (k1-brad)/(k1+brad)
+  b(1) = 1. + alpha(1) + gamma(1) - gamma(1)*bn
+  
+  ! Set RHS         
+  r(1) = gamma(1)*(annp1+ann) + &
+       &     (1.-alpha(1)-gamma(1)+gamma(1)*bn)*T(1) + alpha(1)*T(2)
+  do concurrent (i=2:nz-1)
+     r(i) = gamma(i)*T(i-1) + (1.-alpha(i)-gamma(i))*T(i)+ alpha(i)*T(i+1)
+  enddo
+  r(nz) = gamma(nz)*T(nz-1) + (1.-gamma(nz))*T(nz) &
+       &     + dt/rhoc(nz)*Fgeotherm/(z(nz)-z(nz-1)) ! assumes rhoc(nz+1)=rhoc(nz)
+
+  ! Solve for T at n+1
+  call tridag(a,b,c,r,T,nz) ! update by tridiagonal inversion
+  
+  Tsurf = 0.5*(annp1 + bn*T(1) + T(1)) ! (T0+T1)/2
+
+  ! iterative predictor-corrector
+  if ((Tsurf > 1.2*Tr .or. Tsurf < 0.8*Tr) .and. iter<10) then  ! linearization error expected
+     ! redo until Tr is within 20% of new surface temperature
+     ! (under most circumstances, the 20% threshold is never exceeded)
+     iter = iter+1
+     Tr = sqrt(Tr*Tsurf)  ! linearize around an intermediate temperature
+     T(:) = Told(:)
+     goto 30
+  endif
+  !if (iter>=5) print *,'consider taking shorter time steps',iter
+  !if (iter>=10) print *,'warning: too many iterations'
+  
+  Fsurf = - k(1)*(T(1)-Tsurf)/z(1) ! heat flux into surface
+end subroutine conductionQ
Index: /trunk/LMDZ.COMMON/libf/evolution/conductionT.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/conductionT.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/conductionT.f90	(revision 3470)
@@ -0,0 +1,70 @@
+subroutine conductionT(nz,z,dt,T,Tsurf,Tsurfp1,ti,rhoc,Fgeotherm,Fsurf)
+!***********************************************************************
+!   conductionT:  program to calculate the diffusion of temperature 
+!                 into the ground with prescribed surface temperature 
+!                 and variable thermal properties on irregular grid
+!   Crank-Nicholson scheme, flux conservative
+!
+!   Eqn: rhoc*T_t = (k*T_z)_z 
+!   BC (z=0): T=T(t)
+!   BC (z=L): heat flux = Fgeotherm
+!
+!   nz = number of grid points
+!   dt = time step
+!   T = vertical temperature profile [K]  (in- and output)
+!   Tsurf, Tsurfp1 = surface temperatures at times n and n+1  
+!   ti = thermal inertia [J m^-2 K^-1 s^-1/2]  VECTOR
+!   rhoc = rho*c  heat capacity per volume [J m^-3 K^-1]  VECTOR
+!   ti and rhoc are not allowed to vary in the layers immediately
+!               adjacent to the surface or the bottom
+!   Fgeotherm = geothermal heat flux at bottom boundary [W/m^2]
+!   Fsurf = heat flux at surface [W/m^2]  (output)
+!
+!   Grid: surface is at z=0
+!         T(1) is at z(1); ...; T(i) is at z(i)
+!         k(i) is midway between z(i-1) and z(i)
+!         rhoc(i) is midway between z(i-1) and z(i)
+!***********************************************************************
+  implicit none
+
+  integer, intent(IN) :: nz
+  real*8, intent(IN) :: z(nz), dt, Tsurf, Tsurfp1, ti(nz), rhoc(nz)
+  real*8, intent(IN) :: Fgeotherm
+  real*8, intent(INOUT) :: T(nz)
+  real*8, intent(OUT) :: Fsurf
+  integer i
+  real*8 alpha(nz), k(nz), gamma(nz), buf
+  real*8 a(nz), b(nz), c(nz), r(nz)
+  
+  ! set some constants
+  k(:) = ti(:)**2/rhoc(:) ! thermal conductivity
+  alpha(1) = k(2)*dt/rhoc(1)/(z(2)-z(1))/z(2) 
+  gamma(1) = k(1)*dt/rhoc(1)/z(1)/z(2) 
+  do i=2,nz-1
+     buf = dt/(z(i+1)-z(i-1))
+     alpha(i) = k(i+1)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i))
+     gamma(i) = k(i)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1))
+  enddo
+  buf = dt/(z(nz)-z(nz-1))**2
+  gamma(nz) = k(nz)*buf/(rhoc(nz)+rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz)
+  
+  ! elements of tridiagonal matrix
+  a(:) = -gamma(:)   !  a(1) is not used
+  b(:) = 1. + alpha(:) + gamma(:)
+  c(:) = -alpha(:)   !  c(nz) is not used
+  b(nz) = 1. + gamma(nz)
+  
+  ! Set RHS         
+  r(1)= alpha(1)*T(2) + (1.-alpha(1)-gamma(1))*T(1) + gamma(1)*(Tsurf+Tsurfp1)
+  do concurrent (i=2:nz-1)
+     r(i) = gamma(i)*T(i-1) + (1.-alpha(i)-gamma(i))*T(i) + alpha(i)*T(i+1)
+  enddo
+  r(nz) = gamma(nz)*T(nz-1) + (1.-gamma(nz))*T(nz) + &
+       &     dt/rhoc(nz)*Fgeotherm/(z(nz)-z(nz-1)) ! assumes rhoc(nz+1)=rhoc(nz)
+
+  ! Solve for T at n+1
+  call tridag(a,b,c,r,T,nz) ! update by tridiagonal inversion
+  
+  Fsurf = -k(1)*(T(1)-Tsurfp1)/z(1) ! heat flux into surface
+  
+end subroutine conductionT
Index: /trunk/LMDZ.COMMON/libf/evolution/derivs.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/derivs.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/derivs.f90	(revision 3470)
@@ -0,0 +1,115 @@
+
+subroutine deriv1(z,nz,y,y0,yNp1,yp)
+  ! first derivative of a function y(z) on irregular grid
+  ! upper b.c.: y(0)=y0
+  ! lower b.c.: yp =0
+  implicit none
+  integer, intent(IN) :: nz
+  real(8), intent(IN) :: z(nz),y(nz),y0,yNp1
+  real(8), intent(OUT) :: yp(nz)
+  integer j
+  real(8) hm,hp,c1,c2,c3
+
+  hp = z(2)-z(1)
+  !hm = z(1)
+  c1 = z(1)/(hp*z(2))
+  c2 = 1/z(1) - 1/(z(2)-z(1))
+  c3 = -hp/(z(1)*z(2))
+  yp(1) = c1*y(2) + c2*y(1) + c3*y0
+  do j=2,nz-1
+     hp = z(j+1)-z(j)
+     hm = z(j)-z(j-1)
+     c1 = +hm/(hp*(z(j+1)-z(j-1)))
+     c2 = 1/hm - 1/hp
+     c3 = -hp/(hm*(z(j+1)-z(j-1)))
+     yp(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
+  enddo
+  yp(nz) = (yNp1 - y(nz-1))/(2.*(z(nz)-z(nz-1)))
+end subroutine deriv1
+
+
+
+subroutine deriv2_full(z,nz,a,b,a0,b0,bNp1,yp2)
+  ! second derivative (a*b_z)_z on irregular grid
+  ! upper b.c.: a(0)=a0, b(0)=b0
+  ! 2nd order, without cross-terms
+  implicit none
+  integer, intent(IN) :: nz
+  real(8), intent(IN) :: z(nz),a(nz),b(nz),a0,b0,bNp1
+  real(8), intent(OUT) :: yp2(nz)
+  integer j
+  real(8) hm,hp,c1,c2,c3
+  
+  c1 = 1./(z(1)*z(2))
+  c2 = 1./((z(2)-z(1))*z(1))
+  c3 = 1./((z(2)-z(1))*z(2))
+  yp2(1) = -c2*a(1)*b(1) &   
+       &  -c1*a0*b(1) + c1*(a(1)+a0)*b0 &
+       &  -c3*a(2)*b(1) + c3*(a(1)+a(2))*b(2)
+  do j=2,nz-1
+     hp = z(j+1)-z(j)
+     hm = z(j)-z(j-1)
+     c1 = 1./(hm*(z(j+1)-z(j-1)))
+     c2 = 1./(hp*hm)
+     c3 = 1./(hp*(z(j+1)-z(j-1)))
+     yp2(j) = -c2*a(j)*b(j) &   
+          &  -c1*a(j-1)*b(j) + c1*(a(j)+a(j-1))*b(j-1) &
+          &  -c3*a(j+1)*b(j) + c3*(a(j)+a(j+1))*b(j+1)
+  enddo
+
+  ! lower bc: b_z = 0
+  !yp(nz)= 2.*a(nz)*(b(nz-1)-b(nz))/(z(nz)-z(nz-1))**2
+  !more general lower bc
+  yp2(nz) = a(nz)*(bNp1 - 2*b(nz) + b(nz-1))/(z(nz)-z(nz-1))**2
+end subroutine deriv2_full
+
+
+
+subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
+  ! second derivative y_zz on irregular grid
+  ! boundary conditions: y(0)=y0, y(nz+1)=yNp1
+  implicit none
+  integer, intent(IN) :: nz
+  real(8), intent(IN) :: z(nz),y(nz),y0,yNp1
+  real(8), intent(OUT) :: yp2(nz)
+  integer j
+  real(8) hm,hp,c1,c2,c3
+
+  c1 = +2./((z(2)-z(1))*z(2))
+  c2 = -2./((z(2)-z(1))*z(1))
+  c3 = +2./(z(1)*z(2))
+  yp2(1) = c1*y(2) + c2*y(1) + c3*y0
+
+  do j=2,nz-1
+     hp = z(j+1)-z(j)
+     hm = z(j)-z(j-1)
+     c1 = +2./(hp*(z(j+1)-z(j-1)))
+     c2 = -2./(hp*hm)
+     c3 = +2./(hm*(z(j+1)-z(j-1)))
+     yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
+  enddo
+
+  yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2
+end subroutine deriv2_simple
+
+
+
+real(8) function deriv1_onesided(j,z,nz,y)
+  ! first derivative of function y(z) at z(j)
+  ! one-sided derivative on irregular grid
+  implicit none
+  integer, intent(IN) :: nz,j
+  real(8), intent(IN) :: z(nz),y(nz)
+  real(8) h1,h2,c1,c2,c3
+  if (j<1 .or. j>nz-2) then
+     deriv1_onesided = -9999.
+  else
+     h1 = z(j+1)-z(j)
+     h2 = z(j+2)-z(j+1)
+     c1 = -(2*h1+h2)/(h1*(h1+h2))
+     c2 =  (h1+h2)/(h1*h2)
+     c3 = -h1/(h2*(h1+h2))
+     deriv1_onesided = c1*y(j) + c2*y(j+1) + c3*y(j+2)
+  endif
+end function deriv1_onesided
+
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 3470)
@@ -0,0 +1,177 @@
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose:  Retreat and growth of subsurface ice on Mars
+!!!           orbital elements remain constant
+!!!
+!!!
+!!! Author: EV, updated NS MSIM dynamical program for the PEM 
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+
+SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,T_in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)
+
+!***********************************************************************
+! Retreat and growth of subsurface ice on Mars
+! orbital elements remain constant
+!***********************************************************************
+  use miscparameters, only : pi, d2r, NMAX, marsDay, solsperyear 
+  !use allinterfaces
+  implicit none
+  integer, parameter :: NP=1   ! # of sites
+  integer nz, i, k, nsoil, iloop
+  real(8)  zmax, delta, z(NMAX), icetime, porosity, icefrac
+  real(8), dimension(NP) ::  albedo, thIn, rhoc
+  real(8), dimension(NP) :: pfrost, p0
+  real(8) newti, stretch, newrhoc, ecc, omega, eps, timestep
+  real(8)  ssi_depth_in, ssi_depth, T1
+  real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG
+  real(8), dimension(NMAX,NP) :: porefill, porefill_in 
+  real(8), dimension(NP) :: Tb, T_in, Tmean1, Tmean3, avrho1
+  real(8) tmax, tlast, avrho1prescribed(NP), l1
+  real(8), external :: smartzfac
+
+  !if (iargc() /= 1) then
+  !   stop 'USAGE: icages ext'
+  !endif
+  !call getarg( 1, ext )
+
+  if (NP>100) stop 'subroutine icelayer_mars cannot handle this many sites'
+
+  ! parameters that never ever change
+  nz=nsoil
+  porosity = 0.4d0  ! porosity of till
+  !rhoc(:) = 1500.*800.  ! will be overwritten
+  icefrac = 0.98
+  tmax = 1
+  tlast = 0.
+  avrho1prescribed(:) = pfrost/T1  ! <0 means absent
+  albedo=0.23
+  !avrho1prescribed(:) = 0.16/200.  ! units are Pa/K
+
+  !open(unit=21,file='lats.'//ext,action='read',status='old',iostat=ierr)
+  !if (ierr /= 0) then
+  !   print *,'File lats.'//ext,'not found'
+  !   stop
+  !endif
+  do k=1,NP
+     !read(21,*) latitude(k),albedo(k),thIn(k),htopo(k)
+     ! empirical relation from Mellon & Jakosky
+     rhoc(k) = 800.*(150.+100.*sqrt(34.2+0.714*thIn(k))) 
+  enddo
+  !close(21)
+
+  ! set eternal grid
+  zmax = 25.
+  !zfac = smartzfac(nz,zmax,6,0.032d0)
+  !call setgrid(nz,z,zmax,zfac)
+  l1=2.e-4
+  do iloop=0,nsoil
+    z(iloop) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
+  enddo
+  
+
+  !open(unit=30,file='z.'//ext,action='write',status='unknown')
+  !write(30,'(999(f8.5,1x))') z(1:nz)
+  !close(30)
+
+  !ecc = ecc_in;  eps = obl_in*d2r;  omega = Lp_in*d2r   ! today
+  ! total atmospheric pressure
+  !p0(:) = 600.
+  ! presently 520 Pa at zero elevation (Smith & Zuber, 1998)
+ ! do k=1,NP
+ !    p0(k)=520*exp(-htopo(k)/10800.)
+ ! enddo
+  timestep = 1 ! must be integer fraction of 1 ka
+  icetime = -tmax-timestep  ! earth years
+  
+  ! initializations 
+  Tb=T_in
+  zdepthF(:) = -9999.
+
+  !zdepthT(1:NP) = -9999.   ! reset again below
+ ! zdepthT(1:NP) = 0.
+
+ ! print *,'RUNNING MARS_FAST'
+ ! print *,'Global model parameters:'
+ ! print *,'nz=',nz,' zfac=',zfac,'zmax=',zmax
+ ! print *,'porosity=',porosity
+ ! print *,'starting at time',icetime,'years'
+ ! print *,'time step=',timestep,'years'
+ ! print *,'eps=',eps/d2r,'ecc=',ecc,'omega=',omega/d2r
+ ! print *,'number of sites=',NP
+ ! print *,'Site specific parameters:'
+  do k=1,NP
+     if (NP>1) print *,'  Site ',k
+ !    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)
+ !    print *,'  skin depths (m)',delta,delta*sqrt(solsperyear)
+     call soilthprop(porosity,1.d0,rhoc(k),thIn(k),1,newrhoc,newti,icefrac)
+     stretch = (newti/thIn(k))*(rhoc(k)/newrhoc)
+     do i=1,nz
+        if (z(i)<delta) cycle
+ !       print *,'  ',i-1,' grid points within diurnal skin depth'
+        exit
+     enddo
+ !    print *,'  ',zmax/(sqrt(solsperyear)*delta),'times seasonal dry skin depth'
+ !    print *,'  ',zmax/(sqrt(solsperyear)*delta*stretch),'times seasonal filled skin depth'
+ !    print *,'  Initial ice depth=',zdepthT(k)
+ !    print *
+  enddo
+ ! call outputmoduleparameters
+ ! print *
+
+  ! open and name all output files
+!  open(unit=34,file='subout.'//ext,action='write',status='unknown')
+!  open(unit=36,file='depthF.'//ext,action='write',status='unknown')
+!  open(unit=37,file='depths.'//ext,action='write',status='unknown')
+
+ ! print *,'Equilibrating initial temperature'
+ ! do i=1,4
+ !    call icelayer_mars(0d0,nz,NP,thIn,rhoc,z,porosity,pfrost,Tb,zdepthF, &
+ !      &  zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, &
+ !      &  latitude,albedo,p0,ecc,omega,eps,icefrac,zdepthT,avrho1, &
+ !      &  avrho1prescribed)
+ ! enddo
+
+  print *,'History begins here'
+ porefill(1:nz,1:NP) =  porefill_in(1:nz,1:NP)
+  zdepthT(1:NP) = ssi_depth_in
+  do
+  !print *,'Zt0=  ',ZdepthT
+     call icelayer_mars(timestep,nz,NP,thIn,rhoc,z,porosity,pfrost,Tb,zdepthF, &
+          & zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, & 
+          & albedo,p0,icefrac,zdepthT,avrho1, &
+          & avrho1prescribed)
+     icetime = icetime+timestep
+   !  print *,'T_after=  ',Tb(:)
+ !    print *,'z=  ',z(:)
+     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
+         !write(36,*) icetime,latitude(k),zdepthF(k),porefill(1:nz,k)
+           ! compact output format
+ !          write(36,'(f10.0,2x,f7.3,1x,f11.5,1x)',advance='no') & 
+ !               & icetime,latitude(k),zdepthF(k)
+ !          call compactoutput(36,porefill(:,k),nz)
+ !          write(37,501) icetime,latitude(k),zdepthT(k), &
+ !               & Tmean1(k),Tmean3(k),zdepthG(k),avrho1(k)
+  !      enddo
+  !   endif
+     print *,icetime
+     if (icetime>=tlast) exit
+  enddo
+
+ ! close(34)
+ ! close(36); close(37)
+
+!501 format (f10.0,2x,f7.3,2x,f10.4,2(2x,f6.2),2x,f9.3,2x,g11.4)
+  
+end subroutine dyn_ss_ice_m
+
+
+
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 3470)
@@ -0,0 +1,251 @@
+
+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, deriv1_onesided, colint
+     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
+     subroutine assignthermalproperties(nz,thIn,rhoc,ti,rhocv, &
+          &                      typeT,icefrac,porosity,porefill)
+       implicit none
+       integer, intent(IN) :: nz
+       integer, intent(IN), optional :: typeT
+       real(8), intent(IN), optional :: icefrac
+       real(8), intent(IN) :: thIn, rhoc
+       real(8), intent(IN), optional :: porosity, porefill(nz)
+       real(8), intent(OUT) :: ti(nz), rhocv(nz)
+     end subroutine assignthermalproperties
+  end interface
+
+  interface
+     pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill)
+       implicit none
+       integer, intent(IN) :: nz, typeF, typeG
+       real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B
+       real(8), intent(INOUT) :: porefill(nz)
+     end subroutine icechanges_poreonly
+  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
+
+  interface
+     subroutine compactoutput(unit,porefill,nz)
+       implicit none
+       integer, intent(IN) :: unit,nz
+       real(8), intent(IN) :: porefill(nz)
+     end subroutine compactoutput
+  end interface
+
+  !end of fast_subs_univ
+  !begin derivs.f90 
+
+  interface
+     subroutine deriv1(z,nz,y,y0,yNp1,yp)
+       implicit none
+       integer, intent(IN) :: nz
+       real(8), intent(IN) :: z(nz),y(nz),y0,yNp1
+       real(8), intent(OUT) :: yp(nz)
+     end subroutine deriv1
+  end interface
+
+  interface
+     subroutine deriv2_full(z,nz,a,b,a0,b0,bNp1,yp2)
+       implicit none
+       integer, intent(IN) :: nz
+       real(8), intent(IN) :: z(nz),a(nz),b(nz),a0,b0,bNp1
+       real(8), intent(OUT) :: yp2(nz)
+     end subroutine deriv2_full
+  end interface
+
+  interface
+     subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
+       implicit none
+       integer, intent(IN) :: nz
+       real(8), intent(IN) :: z(nz),y(nz),y0,yNp1
+       real(8), intent(OUT) :: yp2(nz)
+     end subroutine deriv2_simple
+  end interface
+  interface
+     real(8) pure function deriv1_onesided(j,z,nz,y)
+       implicit none
+       integer, intent(IN) :: nz,j
+       real(8), intent(IN) :: z(nz),y(nz)
+     end function deriv1_onesided
+  end interface
+
+  !end of derivs.f90
+  !begin fast_subs_exper.f90
+
+  interface
+     subroutine icelayer_exper(bigstep, nz, thIn, rhoc, z, porosity, pfrost, &
+          & zdepthT, icefrac, zdepthF, zdepthE, porefill, Tmean, Tampl, Diff, zdepthG)
+       use miscparameters, only : d2r, NMAX, icedensity
+       implicit none
+       integer, intent(IN) :: nz
+       real(8), intent(IN) :: bigstep
+       real(8), intent(IN) :: thIn, rhoc, z(NMAX), porosity, pfrost
+       real(8), intent(INOUT) :: zdepthT, zdepthF, porefill(nz)
+       real(8), intent(OUT) :: zdepthE
+       real(8), intent(IN) :: icefrac, Diff, Tmean, Tampl
+       real(8), intent(OUT) :: zdepthG
+     end subroutine icelayer_exper
+  end interface
+
+  interface
+     subroutine ajsub_exper(typeT, nz, z, ti, rhocv, pfrost, Tmean, Tampl, &
+          &     avdrho, avdrhoP, zdepthE, typeF, zdepthF, ypp, porefill, & 
+          &     B, typeG, zdepthG)
+       use miscparameters, only : NMAX, solsperyear, marsDay
+       implicit none
+       integer, intent(IN) :: nz, typeT
+       real(8), intent(IN) :: z(NMAX), ti(NMAX), rhocv(NMAX), pfrost
+       real(8), intent(IN) :: Tmean, Tampl
+       real(8), intent(OUT) :: avdrho, avdrhoP 
+       real(8), intent(OUT) :: zdepthE
+       integer, intent(OUT) :: typeF 
+       real(8), intent(INOUT) :: zdepthF 
+       real(8), intent(OUT) :: ypp(nz)
+       real(8), intent(IN) :: porefill(nz)
+       real(8), intent(IN) :: B 
+       integer, intent(OUT) :: typeG
+       real(8), intent(OUT) :: zdepthG
+       real(8), external :: Tsurface, psv
+       real(8), external :: equildepth 
+     end subroutine ajsub_exper
+  end interface
+  
+  !end fast_subs_exper
+
+  ! Other
+  interface
+     pure function flux_mars77(R,decl,HA,albedo,fracir,fracscat)
+       implicit none
+       real*8 flux_mars77
+       real*8, intent(IN) :: R,decl,HA,albedo,fracIR,fracScat
+     end function flux_mars77
+  end interface
+
+  interface
+     pure function colint(y,z,nz,i1,i2)
+       implicit none
+       integer, intent(IN) :: nz, i1, i2
+       real(8), intent(IN) :: y(nz),z(nz)
+       real(8) colint
+     end function colint
+  end interface
+
+  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/fast_subs_mars.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/fast_subs_mars.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/fast_subs_mars.f90	(revision 3470)
@@ -0,0 +1,307 @@
+module thermalmodelparam_mars
+  ! parameters for thermal model
+  ! they are only used in the subroutines below
+  real(8), parameter :: dt = 0.02  ! in units of Mars solar days
+  !real(8), parameter :: Fgeotherm = 0.
+  real(8), parameter :: Fgeotherm = 0 !0.028  ! [W/m^2]
+  real(8), parameter :: Lco2frost=6.0e5, co2albedo=0.60, co2emiss=1.
+  real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
+  integer, parameter :: EQUILTIME =15 ! [Mars years]
+end module thermalmodelparam_mars
+
+
+!*****************************************************
+! Subroutines for fast method
+! written by Norbert Schorghofer 2007-2011
+!*****************************************************
+
+
+subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, &
+     & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, &
+     & albedo,p0,icefrac,zdepthT,avrho1, &
+     & avrho1prescribed)
+!*************************************************************************
+! bigstep = time step [Earth years]
+! latitude  [degree]
+!*************************************************************************
+  use miscparameters, only : d2r, NMAX, icedensity
+  use allinterfaces, except_this_one => icelayer_mars
+  !use omp_lib
+  implicit none
+  integer, intent(IN) :: nz, NP
+  real(8), intent(IN) :: bigstep
+  real(8), intent(IN) :: thIn(NP), rhoc(NP), z(NMAX), 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)  ! <0 means absent
+
+  integer k, typeF, typeG, typeT, j, jump, typeP
+  real(8) fracIR, fracDust, ti(NMAX), rhocv(NMAX)
+  real(8) Diff, ypp(nz), avdrho(NP), avdrhoP(NP), B, deltaz
+  real(8), SAVE :: avdrho_old(100), zdepth_old(100)  ! NP<=100
+  logical mode2
+
+  !$omp parallel &
+  !$omp    private(Diff,fracIR,fracDust,B,typeT,j,ti,rhocv,typeF,jump,typeG)
+  !$omp do
+  do k=1,NP   ! big loop
+
+     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)
+     
+     typeT = -9
+     if (zdepthT(k)>=0. .and. zdepthT(k)<z(nz)) then
+        do j=1,nz
+           if (z(j)>zdepthT(k)) then ! ice
+              typeT = j  ! first point in ice
+              exit
+           endif
+        enddo
+     endif
+
+   !  call assignthermalproperties(nz,thIn(k),rhoc(k),ti,rhocv, &
+   !       & typeT,icefrac,porosity,porefill(:,k))
+     
+     !----run thermal model
+     call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, & 
+          &     ti, rhocv, fracIR, fracDust, p0(k), &
+          &     avdrho(k), avdrhoP(k), avrho1(k), Tb(k), zdepthE(k), typeF, &
+          &     zdepthF(k), ypp, porefill(:,k), Tmean1(k), Tmean3(k), B, &
+          &     typeG, zdepthG(k), avrho1prescribed(k))
+
+     if (typeF*zdepthF(k)<0.) stop 'error in icelayer_mars'
+     ! diagnose
+     if (zdepthT(k)>=0.) then
+        jump = 0
+        do j=1,nz
+           if (zdepth_old(k)<z(j) .and. zdepthT(k)>z(j)) jump = jump+1
+        enddo
+     else
+        jump = -9
+     endif
+     if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then 
+        write(34,*) '# zdepth arrested'
+        if (jump>1) write(34,*) '# previous step was too large',jump
+     endif
+!     write(34,'(f8.3,1x,f6.2,1x,f11.5,2(1x,g11.4),1x,i3)') &
+!          &        bigstep,latitude(k),zdepthT(k),avdrho(k),avdrhoP(k),jump
+     zdepth_old(k) = zdepthT(k)
+     avdrho_old(k) = avdrho(k)
+
+!----mode 2 growth  
+     typeP = -9;  mode2 = .FALSE.
+     do j=1,nz
+        if (porefill(j,k)>0.) then
+           typeP = j  ! first point with ice
+           exit
+        endif
+     enddo
+     if (typeT>0 .and. typeP>2 .and. zdepthE(k)>0.) then
+        if (porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0. .and. &
+             & zdepthE(k)<z(typeP) .and. &
+             & z(typeP)-zdepthE(k)>2*(z(typeP)-z(typeP-1))) then  ! trick that avoids oscillations
+           deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B  ! conservation of mass 
+           if (deltaz>z(typeP)-z(typeP-1)) then  ! also implies avdrhoP<0.
+              mode2 = .TRUE.
+           endif
+        endif
+     endif
+
+     !call icechanges_poreonly(nz,z,typeF,typeG,avdrhoP(k),ypp,B,porefill(:,k))
+     call icechanges(nz,z(:),typeF,avdrho(k),avdrhoP(k),ypp(:), &
+          & Diff,porosity,icefrac,bigstep,zdepthT(k),porefill(:,k),typeG)
+
+     if (mode2 .and. porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0.) then  ! nothing changed
+        porefill(typeP-1,k)=1.  ! paste a layer
+   !     write(34,*) '# mode 2 growth occurred',typeP,typeF,typeT
+     endif
+
+  enddo  ! end of big loop
+  !$omp end do
+  !$omp end parallel
+end subroutine icelayer_mars
+
+
+
+subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, &
+     &     fracIR, fracDust, patm, avdrho, avdrhoP, avrho1, &
+     &     Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, &
+     &     B, typeG, zdepthG, avrho1prescribed)
+!***********************************************************************
+!  A 1D thermal model that returns various averaged quantities
+!
+!  Tb<0 initializes temperatures
+!  Tb>0 initializes temperatures with Tb 
+!***********************************************************************
+  use miscparameters
+  use thermalmodelparam_mars
+  use allinterfaces, except_this_one => ajsub_mars
+  implicit none
+  integer, intent(IN) :: nz, typeT
+!  real(8), intent(IN) :: latitude  ! in radians
+  real(8), intent(IN) :: albedo0, pfrost, z(NMAX)
+  real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, patm
+  real(8), intent(IN) ::  porefill(nz)
+  real(8), intent(OUT) :: avdrho, avdrhoP  ! difference in annual mean vapor density
+  real(8), intent(OUT) :: avrho1  ! mean vapor density on surface
+  real(8), intent(INOUT) :: Tb, Tmean1
+  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
+  real(8), intent(OUT) :: zdepthE, zdepthF 
+  real(8), intent(OUT) :: ypp(nz) ! (d rho/dt)/Diff
+  real(8), intent(OUT) :: Tmean3, zdepthG
+  real(8), intent(IN) :: B  ! just passing through
+  integer, intent(OUT) :: typeG
+  real(8), intent(IN), optional :: avrho1prescribed
+  real(8), parameter :: a = 1.52366 ! Mars semimajor axis in A.U.
+  integer nsteps, n, i, nm, typeP
+  real(8) tmax, time, Qn, Qnp1, tdays
+  real(8) marsR, marsLs, marsDec, HA
+  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
+  
+  tmax = EQUILTIME*solsperyear
+  nsteps = int(tmax/dt)     ! calculate total number of timesteps
+
+  Tco2frost = tfrostco2(patm) 
+
+  !if (Tb<=0.) then  ! initialize
+     !Tmean0 = 210.15       ! black-body temperature of planet
+     !Tmean0 = (589.*(1.-albedo0)*cos(latitude)/(pi*emiss0*sigSB))**0.25 ! better estimate
+     !Tmean0 = Tmean0-5.
+     !write(34,*) '# initialized with temperature estimate at',latitude/d2r,'of',Tmean0,'K'
+     !T(1:nz) = Tmean0 
+  !else
+     T(1:nz) = Tb
+     ! not so good when Fgeotherm is on
+  !endif
+  
+  albedo = albedo0
+  emiss = emiss0
+  do i=1,nz
+     if (T(i)<Tco2frost) T(i)=Tco2frost
+  enddo
+  Tsurf = T(1)
+  m=0.; Fsurf=0.
+
+  nm=0
+  avrho1=0.; avrho2=0.
+  Tmean1=0.; Tmean3=0.
+  rhosatav0 = 0.
+  rhosatav(:) = 0.
+
+  time=0.
+  !call generalorbit(0.d0,a,ecc,omega,eps,marsLs,marsDec,marsR)
+  !HA = 2.*pi*time            ! hour angle
+!  Qn=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
+  !Qn = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
+  !----loop over time steps 
+  do n=0,nsteps-1
+     time = (n+1)*dt         !   time at n+1 
+     tdays = time*(marsDay/earthDay) ! parenthesis may improve roundoff
+   !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
+   !  HA = 2.*pi*mod(time,1.d0)  ! hour angle
+!     Qnp1=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
+     !Qnp1 = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
+     
+     Tsurfold = Tsurf
+     Fsurfold = Fsurf
+     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, &
+       !      &           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, &
+             !&              rhocv,Fgeotherm,Fsurf) 
+        Tsurf = Tco2frost
+    !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
+    !         &      emiss*sigSB*(Tsurfold**4+Tsurf**4)/2.
+        m = m + dt*marsDay*dE/Lco2frost
+     endif
+     if (Tsurf>Tco2frost .or. m<=0.) then
+        albedo = albedo0
+        emiss = emiss0
+     else
+        albedo = co2albedo
+        emiss = co2emiss
+     endif
+     !Qn=Qnp1
+     
+     if (time>=tmax-solsperyear) then
+        Tmean1 = Tmean1 + Tsurf
+        Tmean3 = Tmean3 + T(nz)
+        avrho1 = avrho1 + min(psv(Tsurf),pfrost)/Tsurf
+        rhosatav0 = rhosatav0 + psv(Tsurf)/Tsurf
+        do i=1,nz
+           rhosatav(i) = rhosatav(i) + psv(T(i))/T(i)
+        enddo
+        nm = nm+1
+     endif
+
+  enddo  ! end of time loop
+  
+  avrho1 = avrho1/nm
+  if (present(avrho1prescribed)) then
+     if (avrho1prescribed>=0.) avrho1=avrho1prescribed
+  endif
+  Tmean1 = Tmean1/nm; Tmean3 = Tmean3/nm
+  rhosatav0 = rhosatav0/nm; rhosatav(:)=rhosatav(:)/nm
+  if (typeT>0) then
+     avrho2 = rhosatav(typeT)
+  else
+     avrho2 = rhosatav(nz) ! no ice
+  endif
+  avdrho = avrho2-avrho1
+  typeP = -9 
+  do i=1,nz
+     if (porefill(i)>0.) then
+        typeP = i  ! first point with ice
+        exit
+     endif
+  enddo
+  if (typeP>0) then
+     avdrhoP = rhosatav(typeP) - avrho1
+  else  
+     avdrhoP = -9999.
+  end if
+
+  zdepthE = equildepth(nz, z, rhosatav, rhosatav0, avrho1)
+
+  if (Fgeotherm>0.) then
+     Tb = Tmean1 
+     typeG = 1   ! will be overwritten by depths_avmeth
+     rlow = 2*rhosatav(nz)-rhosatav(nz-1)
+  else
+     Tb = T(nz)
+     typeG = -9
+     rlow = rhosatav(nz-1)
+  endif
+  call depths_avmeth(typeT, nz, z, rhosatav(:), rhosatav0, rlow, avrho1,  &
+       & porefill(:), typeF, zdepthF, B, ypp(:), typeG, zdepthG)
+
+end subroutine ajsub_mars
+
+
+
+subroutine outputmoduleparameters
+  use miscparameters
+  use thermalmodelparam_mars
+  implicit none
+!  print *,'Parameters stored in modules'
+!  print *,'  Ice bulk density',icedensity,'kg/m^3'
+!  print *,'  dt=',dt,'Mars solar days'
+!  print *,'  Fgeotherm=',Fgeotherm,'W/m^2'
+!  write(6,'(2x,a27,1x,f5.3)') 'Emissivity of dry surface=',emiss0
+!  write(6,'(2x,a12,1x,f5.3,2x,a16,1x,f5.3)') 'CO2 albedo=',co2albedo,'CO2 emissivity=',co2emiss
+!  print *,'  Thermal model equilibration time',EQUILTIME,'Mars years'
+end subroutine outputmoduleparameters
+
+
+
Index: /trunk/LMDZ.COMMON/libf/evolution/fast_subs_univ.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/fast_subs_univ.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/fast_subs_univ.f90	(revision 3470)
@@ -0,0 +1,383 @@
+!*****************************************************
+! Commonly used subroutines for fast method
+! written by Norbert Schorghofer 2007-2011
+!*****************************************************
+
+pure function zint(y1,y2,z1,z2)
+  ! interpolate between two values, y1*y2<0
+  implicit none
+  real(8), intent(IN) :: y1,y2,z1,z2
+  real(8) zint
+  zint = (y1*z2 - y2*z1)/(y1-y2)
+end function zint
+
+
+
+pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1)
+!***********************************************************************
+!  returns equilibrium depth for given ice content
+!  this is not the true (final) equilibrium depth
+!***********************************************************************
+  use allinterfaces, except_this_one => equildepth
+  implicit none
+  integer, intent(IN) :: nz
+  real(8), intent(IN) :: z(nz), rhosatav(nz)
+  real(8), intent(IN) :: rhosatav0, avrho1
+  integer i, typeE
+  real(8) equildepth  ! =zdepthE
+  !real(8), external :: zint  ! defined in allinterfaces.mod
+  
+  typeE = -9; equildepth = -9999.
+  do i=1,nz 
+     if (rhosatav(i) <= avrho1) then
+        typeE=i
+        exit
+     endif
+  enddo
+  if (typeE>1) then ! interpolate
+     equildepth=zint(avrho1-rhosatav(typeE-1),avrho1-rhosatav(typeE),z(typeE-1),z(typeE))
+  end if
+  if (typeE==1) equildepth=zint(avrho1-rhosatav0,avrho1-rhosatav(1),0.d0,z(1))
+  if (equildepth>z(nz)) equildepth=-9999.   ! happens very rarely
+end function equildepth
+
+
+
+subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1,  &
+     & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG)
+!***********************************************************************
+!  returns interface depth and ypp
+!  also returns lower geothermally limited boundary, if applicable
+!
+!  this is a crucial subroutine
+!
+!  B = Diff*bigstep/(porosity*icedensity)  [SI units]
+!***********************************************************************
+  use allinterfaces, except_this_one => depths_avmeth
+  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  ! index of depth below which filling occurs
+  real(8), intent(INOUT) :: zdepthF
+  real(8), intent(IN) :: B 
+  real(8), intent(OUT) :: ypp(nz), zdepthG
+  integer, intent(INOUT) :: typeG  ! positive on input when Fgeotherm>0
+
+  integer i, typeP, nlb, newtypeG
+  real(8) eta(nz), Jpump1, help(nz), yp(nz), zdepthFold, ap_one, ap(nz)
+  real(8) leak, cumfill, cumfillabove
+
+  if (typeT<0) then
+     nlb = nz
+     do i=1,nz
+        eta(i) = constriction(porefill(i))
+     enddo
+  else
+     !nlb = typeT-1
+     nlb = typeT ! changed 2010-09-29
+     do i=1,typeT-1
+        eta(i) = constriction(porefill(i))
+     enddo
+     do i=typeT,nz
+        eta(i)=0.
+     enddo
+  end if
+
+!-fill depth
+  zdepthFold = zdepthF
+  typeF = -9;  zdepthF = -9999.
+  call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp)  ! yp also used below
+  do i=1,nlb
+     Jpump1 = (rhosatav(i)-avrho1)/z(i)  ! <0 when stable
+     ! yp is always <0
+     help(i) = Jpump1 - eta(i)*yp(i)
+     leak = porefill(i)/B*(z(i)-zdepthFold)/(18./8314.)
+     !help(i) = help(i)-leak   ! optional
+     if (help(i) <= 0.) then
+        typeF=i
+        !print *,'#',typeF,Jpump1,eta(typeF)*yp(typeF),leak
+        exit
+     endif
+  enddo
+  if (typeF>1) zdepthF = zint(help(typeF-1),help(typeF),z(typeF-1),z(typeF))
+  if (typeF==1) zdepthF=z(1)
+
+
+!-depth to shallowest perennial ice
+  typeP = -9 
+  do i=1,nz
+     if (porefill(i)>0.) then
+        typeP = i  ! first point with ice
+        exit
+     endif
+  enddo
+
+!-calculate ypp
+  !call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp)
+  call deriv1(z,nz,eta(:),1.d0,eta(nz-1),ap)
+  if (typeP>0 .and. typeP<nz-2) then
+     ap_one=deriv1_onesided(typeP,z,nz,eta(:))
+     ! print *,typeP,ap(typeP),ap_one
+     ap(typeP)=ap_one
+  endif
+  call deriv2_simple(z,nz,rhosatav(1:nz),rhosatav0,rlow,ypp(:))
+  !call deriv2_full(z,nz,eta(:),rhosatav(:),1.d0,rhosatav0,rhosatav(nz-1),ypp(:))
+  !write(40,*) rhosatav
+  !write(41,*) yp
+  !write(42,*) ypp
+  ypp(:) = ap(:)*yp(1:)+eta(:)*ypp(:)
+  !write(43,*) ypp
+  !write(44,*) eta(1:nz)
+  !write(45,*) (rhosatav(:)-avrho1)/z(:)
+  ypp(:) = ypp(:)*18./8314.
+  ! ypp values beyond nlb should not be used
+
+!-geothermal stuff
+  if (typeT>0) typeG=-9
+  if (typeG<0) zdepthG=-9999.
+  if (typeG>0 .and. typeT<0) then
+     typeG=-9
+     do i=2,nz
+        if (yp(i)>0.) then  ! first point with reversed flux
+           typeG=i
+           zdepthG=zint(yp(i-1),yp(i),z(i-1),z(i))
+           !zdepthG=z(typeG)
+           exit
+        endif
+     enddo
+  else
+     typeG = -9
+  endif
+  if (typeG>0 .and. typeT<0) then
+     cumfillabove = colint(porefill(:)/eta(:),z,nz,typeG-1,nz) 
+     newtypeG = -9
+     do i=typeG,nz
+        if (minval(eta(i:nz))<=0.) cycle  ! eta=0 means completely full
+        cumfill=colint(porefill(:)/eta(:),z,nz,i,nz)
+        if (cumfill<yp(i)*18./8314.*B) then  ! usually executes on i=typeG
+           if (i>typeG) then
+              write(34,*) '# adjustment to geotherm depth by',i-typeG
+              zdepthG = zint(yp(i-1)*18./8314.*B-cumfillabove, &  ! no good
+                   &        yp(i)*18./8314.*B-cumfill,z(i-1),z(i))
+              if (zdepthG>z(i) .or. zdepthG<z(i-1)) write(34,*) &
+                   & '# WARNING: zdepthG interpolation failed',i,z(i-1),zdepthG,z(i)
+              newtypeG=i
+           endif
+           ! otherwise leave zdepthG untouched
+           exit
+        endif
+        cumfillabove = cumfill
+     enddo
+     if (newtypeG>0) typeG=newtypeG
+  end if
+  ! if typeG>0, then all ice at and below typeG should be erased 
+end subroutine depths_avmeth
+
+
+
+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_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill)
+  use allinterfaces, except_this_one => icechanges_poreonly
+  implicit none
+  integer, intent(IN) :: nz, typeF, typeG
+  real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B
+  real(8), intent(INOUT) :: porefill(nz)
+  integer j, erase, newtypeP, ub
+  real(8) integ
+  
+  !----retreat
+  ! avdrhoP>0 is outward loss from zdepthP
+  ! avdrhoP<0 means gain at zdepthP or no ice anywhere
+  if (avdrhoP>0.) then
+     erase=0
+     do j=1,nz
+        if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF
+        integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j)
+        erase = j
+        if (integ>B*avdrhoP*18./8314.) exit
+     end do
+     if (erase>0) porefill(1:erase)=0.
+  endif
+
+  ! new depth
+  newtypeP = -9 
+  do j=1,nz
+     if (porefill(j)>0.) then
+        newtypeP = j  ! first point with ice
+        exit
+     endif
+  enddo
+
+  !----diffusive filling
+  ub = typeF
+  if (newtypeP>0 .and. typeF>0 .and. newtypeP<ub) ub=newtypeP
+  if (ub>0) then  
+     do j=ub,nz
+        ! B=Diff/(porosity*icedensity)*86400*365.24*bigstep
+        porefill(j) = porefill(j) + B*ypp(j)
+        if (porefill(j)<0.) porefill(j)=0.
+        if (porefill(j)>1.) porefill(j)=1.
+     enddo
+  end if
+  
+  !----enact bottom boundary
+  if (typeG>0) porefill(typeG:nz)=0.
+  
+end subroutine icechanges_poreonly
+
+
+
+pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, &
+     & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG)
+!***********************************************************
+! advances ice table, advances interface, grows pore ice
+!
+! a crucial subroutine
+!***********************************************************
+  use miscparameters, only : icedensity
+  use allinterfaces, except_this_one => icechanges
+  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)
+  integer j, erase, newtypeP, ub, typeP, typeT
+  real(8) B, beta, integ
+
+  B = Diff*bigstep*86400.*365.24/(porosity*icedensity)
+
+  ! advance ice table, avdrho>0 is retreat
+  if (zdepthT>=0. .and. avdrho>0.) then 
+     typeP=-9999; typeT=-9999
+     do j=1,nz
+        if (z(j)>zdepthT) then
+           typeT=j
+           exit
+        endif
+     enddo
+     do j=1,nz
+        if (porefill(j)>0.) then
+           typeP=j
+           exit
+        endif
+     enddo
+     if (typeP==typeT) then   ! new 2011-09-01
+        beta = (1-icefrac)/(1-porosity)/icefrac
+        beta = Diff*bigstep*beta*86400*365.24/icedensity
+        zdepthT = sqrt(2*beta*avdrho*18./8314. + zdepthT**2)
+     endif
+  endif
+  if (zdepthT>z(nz)) zdepthT=-9999.
+  
+  ! advance interface, avdrhoP>0 is loss from zdepthP
+  if (avdrhoP>0.) then
+     erase=0
+     do j=1,nz
+        if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF
+        if (zdepthT>=0. .and. z(j)>zdepthT) exit 
+        integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j)
+        erase = j
+        if (integ>B*avdrhoP*18./8314.) exit
+     end do
+     if (erase>0) porefill(1:erase)=0.
+  endif
+
+  ! new depth
+  newtypeP = -9 
+  do j=1,nz
+     if (zdepthT>=0. .and. z(j)>zdepthT) exit
+     if (porefill(j)>0.) then
+        newtypeP = j  ! first point with pore ice
+        exit
+     endif
+  enddo
+
+  ! diffusive filling
+  ub = typeF
+  if (newtypeP>0 .and. typeF>0 .and. newtypeP<ub) ub=newtypeP
+  if (ub>0) then  
+     do j=ub,nz
+        porefill(j) = porefill(j) + B*ypp(j)
+        if (porefill(j)<0.) porefill(j)=0.
+        if (porefill(j)>1.) porefill(j)=1.
+        if (zdepthT>=0. .and. z(j)>zdepthT) exit
+     enddo
+  end if
+
+  ! below icetable
+  if (zdepthT>=0.) then
+     do j=1,nz
+        if (z(j)>zdepthT) porefill(j) = icefrac/porosity
+     enddo
+  else
+     ! geothermal lower boundary
+     if (typeG>0) porefill(typeG:nz)=0.
+  end if
+end subroutine icechanges
+
+
+subroutine assignthermalproperties(nz,thIn,rhoc, &
+     &    ti,rhocv,typeT,icefrac,porosity,porefill)
+!*********************************************************
+! assign thermal properties of soil
+!*********************************************************
+  implicit none
+  integer, intent(IN) :: nz
+  integer, intent(IN), optional :: typeT
+  real(8), intent(IN), optional :: icefrac
+  real(8), intent(IN) :: thIn, rhoc
+  real(8), intent(IN), optional :: porosity, porefill(nz)
+  real(8), intent(OUT) :: ti(nz), rhocv(nz)
+  integer j
+  real(8) newrhoc, newti, fill
+  real(8), parameter :: NULL=0.
+
+  ti(1:nz) = thIn
+  rhocv(1:nz) = rhoc
+  if (typeT>0) then 
+     call soilthprop(porosity,NULL,rhoc,thIn,2,newrhoc,newti,icefrac)
+     rhocv(typeT:nz) = newrhoc
+     ti(typeT:nz) = newti
+  endif
+  do j=1,nz 
+     fill = porefill(j)   ! off by half point
+     if (fill>0. .and. (typeT<0 .or. (typeT>0 .and. j<typeT))) then
+        call soilthprop(porosity,fill,rhoc,thIn,1,rhocv(j),ti(j),NULL)
+     endif
+  enddo
+end subroutine assignthermalproperties
+
+
+
+subroutine compactoutput(unit,porefill,nz)
+  implicit none
+  integer, intent(IN) :: unit,nz
+  real(8), intent(IN) :: porefill(nz)
+  integer j
+  do j=1,nz
+     if (porefill(j)==0.) then
+        write(unit,'(1x,f2.0)',advance='no') porefill(j)
+     else
+        write(unit,'(1x,f7.5)',advance='no') porefill(j)
+     endif
+  enddo
+  write(unit,"('')")
+end subroutine compactoutput
+
Index: /trunk/LMDZ.COMMON/libf/evolution/flux_mars.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/flux_mars.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/flux_mars.f90	(revision 3470)
@@ -0,0 +1,125 @@
+pure function flux_mars77(R,decl,latitude,HA,albedo,fracir,fracscat)
+!***********************************************************************
+! flux_mars77: calculates insolation (incoming solar radiation) on Mars
+!     flat surface only; also works in polar regions
+!
+!     R: distance from sun [AU]
+!     decl: planetocentric solar declination [radians]
+!     latitude: [radians]
+!     HA: hour angle [radians from noon, clockwise]
+!     fracir: fraction of absorption
+!     fracscat: fraction of scattering
+!***********************************************************************
+  implicit none
+  real(8) flux_mars77
+  real(8), intent(IN) :: R,decl,latitude,HA,albedo,fracIR,fracScat
+  real(8), parameter :: So=1365. ! [W/m^2]
+  real(8), parameter :: sigSB=5.6704d-8
+  real(8) c1,s1,Qo,solarAttenuation,Q
+  real(8) sinbeta,sinbetaNoon,Qdir,Qscat,Qlw,Fout
+
+  c1 = cos(latitude)*cos(decl)
+  s1 = sin(latitude)*sin(decl)
+  ! beta = 90 minus incidence angle for horizontal surface
+  ! beta = elevation of sun above (horizontal) horizon 
+  sinbeta = c1*cos(HA) + s1
+  sinbetaNoon = c1 + s1
+  sinbetaNoon = max(sinbetaNoon,0.d0)  ! horizon
+
+  Qo = So/(R**2)  ! solar constant at Mars
+  
+  ! short-wavelength irradiance
+  if (sinbeta>0.d0) then
+     solarAttenuation = (1.- fracIR - fracScat)**(1./max(sinbeta,0.04))
+     Qdir = Qo*sinbeta*solarAttenuation
+     Qscat = 0.5*Qo*fracScat
+     Q = (1.-albedo)*(Qdir + Qscat)
+  else
+     Q = 0.
+  endif
+
+  ! atmospheric IR contribution based on Kieffer et al. (1977), JGR 82, 4249
+  Fout = 1.*sigSB*150**4  ! matters only in polar region
+  Qlw = fracIR*max(Qo*sinbetaNoon,Fout)
+
+  ! net flux = short-wavelength + long-wavelength irradiance
+  flux_mars77 = Q + Qlw
+end function flux_mars77
+
+
+
+pure subroutine flux_mars2(R,decl,latitude,HA,fracIR,fracScat, &
+     &   SlopeAngle,azFac,emax,Qdir,Qscat,Qlw)
+!*****************************************************************************
+! flux_mars2: Insolation (incoming solar radiation) at Mars on sloped surface;
+!             returns several irradiances
+!
+! INPUTS:
+!     R: distance from sun [AU]
+!     decl: planetocentric solar declination [radians]
+!     latitude: [radians]
+!     HA: hour angle (radians from noon, clockwise)
+!     fracIR: fraction of absorption
+!     fracScat: fraction of scattering
+!     SlopeAngle: >0, [radians]
+!     azFac: azimuth of topographic gradient (radians east of north)
+!            azFac=0 is south-facing  
+!     emax: maximum horizon elevation in direction of azimuth [radians]
+! OUTPUTS:
+!     Qdir: direct incoming short-wavelength irradiance [W/m^2]
+!     Qscat: diffuse short-wavelength irradiance from atmosphere [W/m^2]
+!     Qlw: diffuse long-wavelength irradiance from atmosphere [W/m^2]
+!*****************************************************************************
+  implicit none
+  real(8), parameter :: pi=3.1415926535897932, So=1365.
+  real(8), parameter :: sigSB=5.6704d-8
+  real(8), intent(IN) :: R,decl,latitude,HA,SlopeAngle,azFac,emax
+  real(8), intent(IN) :: fracIR,fracScat
+  real(8), intent(OUT) :: Qdir,Qscat,Qlw
+  real(8) c1,s1,solarAttenuation,Qo
+  real(8) sinbeta,sinbetaNoon,cosbeta,sintheta,azSun,buf,Fout
+  
+  c1 = cos(latitude)*cos(decl)
+  s1 = sin(latitude)*sin(decl)
+  ! beta = 90 minus incidence angle for horizontal surface
+  ! beta = elevation of sun above (horizontal) horizon 
+  sinbeta = c1*cos(HA) + s1
+  sinbetaNoon = c1 + s1
+  sinbetaNoon = max(sinbetaNoon,0.d0)  ! horizon
+      
+  cosbeta = sqrt(1-sinbeta**2)
+  buf = (sin(decl)-sin(latitude)*sinbeta)/(cos(latitude)*cosbeta)
+  if (buf>+1.) buf=+1.d0  ! roundoff
+  if (buf<-1.) buf=-1.d0  ! roundoff
+  azSun = acos(buf)
+  if (sin(HA)>=0) azSun = 2*pi-azSun
+  
+! theta = 90 minus incidence angle for sloped surface
+  sintheta = cos(SlopeAngle)*sinbeta - &  
+       &     sin(SlopeAngle)*cosbeta*cos(azSun-azFac)
+  if (cosbeta==0) sintheta=cos(SlopeAngle)*sinbeta  ! zenith, where azSun=NaN
+  
+  if (sintheta<0.) sintheta=0.  ! local horizon (self-shadowing)
+  if (sinbeta<0.) sintheta=0.  ! horizontal horizon at infinity
+  if (sinbeta<sin(emax)) sintheta=0. ! distant horizon
+
+! fluxes and contributions from atmosphere
+  Qo = So/R**2  ! solar constant at Mars
+  solarAttenuation = (1.- fracIR - fracScat)**(1./max(sinbeta,0.04d0))
+  Qdir = Qo*sintheta*solarAttenuation
+  Fout = 1.*sigSB*150**4  ! matters only in polar region
+  Qlw = fracIR*max(Qo*sinbetaNoon,Fout) 
+  if (sinbeta>0.d0) then
+     Qscat = 0.5*fracScat*Qo
+  else
+     Qscat = 0.
+  endif
+  
+! For a horizontal and unobstructed surface
+!   absorbed flux = (1-albedo)*(Qdir+Qscat) + emiss*Qlw
+!
+! For a tilted surface
+!   absorbed flux = (1-albedo)*(Qdir+Qscat*viewfactor) + emiss*Qlw*viewfactor
+!   then add irradiance from land in field of view  
+!   in the case of a planar slope, viewfactor = cos(SlopeAngle/2.)**2
+end subroutine flux_mars2
Index: /trunk/LMDZ.COMMON/libf/evolution/generalorbit.f
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/generalorbit.f	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/generalorbit.f	(revision 3470)
@@ -0,0 +1,54 @@
+C=============================================================
+C Subroutine to return distance, longitude, and declination of 
+C the sun in planetocentric coordinates from orbital elements
+C
+C INPUTS: 
+C edays = time since perihelion (earth days)
+C a = semimajor axis (AU)
+C ecc = eccentricity
+C omega = Ls of perihelion, relative to equinox (radians)
+C eps = obliquity (radians)
+C
+C OUTPUTS:
+C Ls = areocentric longitude (radians)
+C dec = planetocentric solar declination (radians)
+C r = heliocentric distance (AU)
+C=============================================================
+
+      SUBROUTINE generalorbit(edays,a,ecc,omega,eps,Ls,dec,r)
+      implicit none
+      real*8 edays,a,ecc,omega,eps  ! input
+      real*8 Ls,dec,r  ! output
+      real*8 pi,d2r
+      parameter (pi=3.1415926535897932,d2r=pi/180.d0)
+      integer j
+      real*8 M,E,nu,T,Eold
+
+c     T = orbital period (days)
+      T = sqrt(4*pi**2/(6.674e-11*1.989e30)*(a*149.598e9)**3)/86400.
+
+c     M = mean anomaly (radians)
+      M = 2.*pi*edays/T  ! M=0 at perihelion
+
+c     E = eccentric anomaly 
+c     solve M = E - ecc*sin(E) by newton method
+      E = M
+      do j=1,10
+         Eold = E
+         E = E - (E - ecc*sin(E) - M)/(1.-ecc*cos(E))
+         if (abs(E-Eold)<1.e-8) exit
+      enddo
+
+c     nu = true anomaly
+      !nu = acos(cos(E)-ecc/(1.-ecc*cos(E)))
+      !nu = sqrt(1-ecc^2)*sin(E)/(1.-ecc*cos(E))
+      !nu = atan(sqrt(1-ecc^2)*sin(E)/(1-cos(E)))
+      nu = 2.*atan(sqrt((1.+ecc)/(1.-ecc))*tan(E/2.))
+
+      !r = a*(1.-ecc**2)/(1.+ecc*cos(nu))
+      r = a*(1-ecc*cos(E))
+      Ls = mod(nu + omega,2.*pi)   
+      dec = asin(sin(eps)*sin(Ls))
+
+      END
+
Index: /trunk/LMDZ.COMMON/libf/evolution/grids.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/grids.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/grids.f90	(revision 3470)
@@ -0,0 +1,96 @@
+subroutine setgrid(nz,z,zmax,zfac)
+  ! construct regularly or geometrically spaced 1D grid
+  ! z(n)-z(1) = 2*z(1)*(zfac**(n-1)-1)/(zfac-1)
+  ! choice of z(1) and z(2) is compatible with conductionQ
+  implicit none
+  integer, intent(IN) :: nz
+  real(8), intent(IN) :: zfac, zmax
+  real(8), intent(OUT) :: z(nz)
+  integer i
+  real(8) dz
+  
+  dz = zmax/nz
+  do i=1,nz
+     z(i) = (i-0.5)*dz
+  enddo
+  
+  if (zfac>1.) then
+     dz = zmax/(3.+2.*zfac*(zfac**(nz-2)-1.)/(zfac-1.))
+     z(1) = dz
+     z(2) = 3*z(1)
+     do i=3,nz
+        z(i) = (1.+zfac)*z(i-1) - zfac*z(i-2)
+        ! z(i) = z(i-1) + zfac*(z(i-1)-z(i-2)) ! equivalent
+     enddo
+  endif
+  
+end subroutine setgrid
+
+
+
+real(8) function smartzfac(nz_max,zmax,nz_delta,delta)
+  ! output can be used as input to setgrid
+  ! produces zfac with desired number of points within depth delta
+  implicit none
+  integer, intent(IN) :: nz_max, nz_delta
+  real(8), intent(IN) :: zmax, delta
+  integer j
+  real(8) f, g, gprime
+  
+  if (nz_max<nz_delta .or. nz_delta<=1) then
+     stop 'inappropriate input to smartzfac'
+  endif
+  f = 1.05
+  do j=1,7  ! Newton iteration
+     !print *,j,f
+     g = (f-3+2*f**(nz_max-1))/(f-3+2*f**(nz_delta-1))-zmax/delta
+     gprime=(1+2*(nz_max-1)*f**(nz_max-2))/(f-3+2*f**(nz_delta-1)) - &
+          &        (f-3+2*f**(nz_max-1))*(1+2*(nz_delta-1)*f**(nz_delta-2))/ & 
+          &        (f-3+2*f**(nz_delta-1))**2
+     f = f-g/gprime
+  enddo
+  smartzfac = f
+  if (smartzfac<1. .or. smartzfac>2.) then
+     print *,'zfac=',smartzfac
+     stop 'unwanted result in smartzfac'
+  endif
+end function smartzfac
+
+
+      
+!-grid-dependent utility functions
+
+function colint(y,z,nz,i1,i2)
+  ! column integrates y
+  implicit none
+  integer, intent(IN) :: nz, i1, i2
+  real(8), intent(IN) :: y(nz), z(nz)
+  real(8) colint
+  integer i
+  real(8) dz(nz)
+  
+  dz(1) = (z(2)-0.)/2
+  do i=2,nz-1
+     dz(i) = (z(i+1)-z(i-1))/2.
+  enddo
+  dz(nz) = z(nz)-z(nz-1)
+  colint = sum(y(i1:i2)*dz(i1:i2))
+end function colint
+
+
+ 
+subroutine heatflux_from_temperature(nz,z,T,k,H)
+  ! calculates heat flux from temperature profile
+  ! like k, the heat flux H is defined mid-point
+  implicit none
+  integer, intent(IN) :: nz
+  real(8), intent(IN) :: z(nz), T(nz), k(nz)
+  real(8), intent(OUT) :: H(nz)
+  integer j
+  
+  ! H(1) = -k(1)*(T(1)-Tsurf)/z(1)
+  H(1) = 0. ! to avoid ill-defined value
+  do j=2,nz
+     H(j) = -k(j)*(T(j)-T(j-1))/(z(j)-z(j-1))
+  enddo
+end subroutine heatflux_from_temperature
Index: /trunk/LMDZ.COMMON/libf/evolution/misc_params.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/misc_params.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/misc_params.f90	(revision 3470)
@@ -0,0 +1,14 @@
+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: /trunk/LMDZ.COMMON/libf/evolution/psv.f
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/psv.f	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/psv.f	(revision 3470)
@@ -0,0 +1,54 @@
+      real*8 function psv(T)
+C     saturation vapor pressure of H2O ice [Pascal]
+C     input is temperature [Kelvin]
+      implicit none
+      real*8 T
+
+C-----parametrization 1
+c      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
+c      parameter (DHmelt=6008.,DHvap=45050.)
+c      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
+c      parameter (R=8.314,pt=6.11e2,Tt=273.16)
+c      C = (DHsub/R)*(1./T - 1./Tt)
+c      psv = pt*exp(-C)
+
+C-----parametrization 2
+C     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
+C     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
+
+C-----parametrization 3      
+C     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
+c     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
+      
+      end
+
+
+      
+      real*8 function frostpoint(p)
+C     inverse of psv
+C     input is partial pressure [Pascal]
+C     output is temperature [Kelvin]
+      implicit none
+      real*8 p
+      
+C-----inverse of parametrization 1
+c      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
+c      parameter (DHmelt=6008.,DHvap=45050.)
+c      parameter (DHsub=DHmelt+DHvap)
+c      parameter (R=8.314,pt=6.11e2,Tt=273.16)
+c      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
+      
+C-----inverse of parametrization 2
+C     inverse of eq. (2) in Murphy & Koop (2005)
+      real*8 A,B
+      parameter (A=-6143.7, B=28.9074)
+      frostpoint = A / (log(p) - B)
+
+C-----approximate inverse of parametrization 3
+C     eq. (8) in Murphy & Koop (2005)
+c      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
+      
+      end
Index: /trunk/LMDZ.COMMON/libf/evolution/psvco2.f
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/psvco2.f	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/psvco2.f	(revision 3470)
@@ -0,0 +1,50 @@
+      real*8 function psvco2(T)
+C     solid-vapor transition for CO2
+C     returns saturation pressure in Pascal, input is temperature in Kelvin
+      implicit none
+      real*8 T
+
+      psvco2 = exp(23.3494 - 3182.48/T)*100.
+      end
+
+
+      real*8 function tfrostco2(p)
+C     the inverse of function psvco2
+C     input is pressure in Pascal, output is temperature in Kelvin
+      implicit none
+      real*8 p
+      tfrostco2 = 3182.48/(23.3494+log(100./p))
+      end
+
+
+
+
+C------------------------------------------------------------------------------
+C     Antoine equation parameters from NIST, 154K-196K 
+C     based on Giauque and Egan (1937)
+C     A=6.81228, B=1301.679, C=-3.494
+c     p = 10**(A-(B/(T+C)))*1.e5
+
+C     Expressions from Int. Crit. Tabl. 3.207, based on many references
+C     mm2Pa = 133.32
+C     -135C to -56.7C (138K-216K)
+c     p = 10**(-0.05223*26179.3/T + 9.9082)*mm2Pa
+C     -183C to -135C (90K-138K)
+c     p = 10**(-1275.62/T + 0.006833*T + 8.3071)*mm2Pa
+C     Expressions from Int. Crit. Tabl. 3.208, based on Henning
+C     -110 to -80C (163K-193K)
+c     p = 10**(- 1279.11/T + 1.75*log10(T) - 0.0020757*T + 5.85242)*mm2Pa
+c     p = 10**(- 1352.6/T + 9.8318)*mm2Pa
+      
+C     Mars book (1992), p959, fit by chapter authors
+c     p = exp(23.3494 - 3182.48/T)*100.   ! 120-160K
+c     p = exp(25.2194 - 3311.57/T - 6.71e-3*T)*100  ! 100-170K
+C     Mars book (1992), p960, based on Miller & Smythe (1970)
+c     p = exp(26.1228 - 3385.26/T - 9.4461e-3*T)*100 ! 123-173K
+
+C     Fray & Schmitt, PSS 57, 2053 (2009)
+C     A0=1.476e1, A1=-2.571e3, A2=-7.781e4, A3=4.325e6, A4=-1.207e8, A5=1.350e9
+c     p = exp(A0+A1/T+A2/T**2+A3/T**3+A4/T**4+A5/T**5)*1e5 ! 40K-194.7K
+c     A0=1.861e1, A1=-4.154e3, A2=1.041e5
+c     p = exp(A0 + A1/T + A2/T**2)*1e5  ! 194.7K-216.58K
+C------------------------------------------------------------------------------
Index: /trunk/LMDZ.COMMON/libf/evolution/soilthprop_mars.f90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soilthprop_mars.f90	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/soilthprop_mars.f90	(revision 3470)
@@ -0,0 +1,152 @@
+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
+
+
+      
+subroutine smartgrid(nz,z,zdepth,thIn,rhoc,porosity,ti,rhocv,layertype,icefrac)
+!***********************************************************************
+! smartgrid: returns intelligently spaced grid and appropriate 
+!            values of thermal inertia ti and rho*c in icy layer
+!                  
+!     INPUTS: 
+!             nz = number of grid points
+!             z = grid spacing for dry regolith 
+!                 (will be partially overwritten)
+!             zdepth = depth where ice table starts
+!                      negative values indicate no ice
+!             rhoc = heat capacity per volume of dry regolith [J/m^3]
+!             thIn = thermal inertia of dry regolith [SI-units]
+!             porosity = void space / total volume
+!             layertypes are explained below  
+!             icefrac = fraction of ice in icelayer
+!
+!     OUTPUTS: z = grid coordinates
+!              vectors ti and rhocv
+!***********************************************************************
+  implicit none
+  integer, intent(IN) :: nz, layertype
+  real(8), intent(IN) :: zdepth, thIn, rhoc, porosity, icefrac
+  real(8), intent(INOUT) :: z(nz)
+  real(8), intent(OUT) :: ti(nz), rhocv(nz)
+  integer j, b
+  real(8) stretch, newrhoc, newti
+  real(8), parameter :: NULL=0.
+  
+  if (zdepth>0 .and. zdepth<z(nz)) then
+
+     select case (layertype)
+     case (1)  ! interstitial ice
+        call soilthprop(porosity,1.d0,rhoc,thIn,1,newrhoc,newti,NULL)
+     case (2)  ! mostly ice (massive ice)
+        call soilthprop(porosity,NULL,rhoc,thIn,2,newrhoc,newti,icefrac)
+     case (3)  ! all ice (pure ice)
+        call soilthprop(NULL,NULL,NULL,NULL,3,newrhoc,newti,NULL)
+     case (4)  ! ice + rock + nothing else (ice-cemented soil)
+        call soilthprop(porosity,NULL,rhoc,NULL,4,newrhoc,newti,NULL)
+     case default
+        error stop 'invalid layer type'
+     end select
+
+     ! Thermal skin depth is proportional to sqrt(kappa)
+     ! thermal diffusivity kappa = k/(rho*c) = I^2/(rho*c)**2
+     stretch = (newti/thIn)*(rhoc/newrhoc) ! ratio of sqrt(thermal diffusivity)
+     
+     b = 1
+     do j=1,nz
+        if (z(j)<=zdepth) then 
+           b = j+1
+           rhocv(j) = rhoc
+           ti(j) = thIn
+        else
+           rhocv(j) = newrhoc
+           ti(j) = newti
+        endif
+        ! print *,j,z(j),ti(j),rhocv(j)
+     enddo
+     do j=b+1,nz
+        z(j) = z(b) + stretch*(z(j)-z(b))
+     enddo
+     
+     ! print *,'zmax=',z(nz),' b=',b,' stretch=',stretch
+     ! print *,'depth at b-1, b ',z(b-1),z(b)
+     ! print *,'ti1=',thIn,' ti2=',newti
+     ! print *,'rhoc1=',rhoc,' rhoc2=',newrhoc 
+  endif
+  
+end subroutine smartgrid
+
Index: /trunk/LMDZ.COMMON/libf/evolution/tridag.for
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/tridag.for	(revision 3470)
+++ /trunk/LMDZ.COMMON/libf/evolution/tridag.for	(revision 3470)
@@ -0,0 +1,30 @@
+C================================================
+C Tridiagonal solver
+C================================================
+      SUBROUTINE tridag(a,b,c,r,u,n)
+      INTEGER n,NMAX
+      REAL*8 a(n),b(n),c(n),r(n),u(n)
+      PARAMETER (NMAX=1000)
+      INTEGER j
+      REAL*8 bet,gam(NMAX)
+      if(b(1).eq.0.) then
+         stop 'tridag: rewrite equations'
+      endif 
+c      if(n.gt.NMAX) then
+c         print *, 'tridag: too many points, set NMAX>',n
+c         stop
+c      endif 
+      bet=b(1)
+      u(1)=r(1)/bet
+      do 11 j=2,n
+        gam(j)=c(j-1)/bet
+        bet=b(j)-a(j)*gam(j)
+c        if(bet.eq.0.)pause 'tridag failed'
+        u(j)=(r(j)-a(j)*u(j-1))/bet
+11    continue
+      do 12 j=n-1,1,-1
+        u(j)=u(j)-gam(j+1)*u(j+1)
+12    continue
+      return
+      END
+C  (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31&#5(+.
