Index: trunk/LMDZ.COMMON/libf/evolution/NS_conductionQ.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_conductionQ.F90	(revision 3525)
+++ 	(revision )
@@ -1,113 +1,0 @@
-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/NS_conductionT.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_conductionT.F90	(revision 3525)
+++ 	(revision )
@@ -1,70 +1,0 @@
-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/NS_derivs.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_derivs.F90	(revision 3525)
+++ 	(revision )
@@ -1,115 +1,0 @@
-
-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/NS_dyn_ss_ice_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90	(revision 3525)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90	(revision 3526)
@@ -1,3 +1,2 @@
-
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -171,5 +170,2 @@
   
 end subroutine dyn_ss_ice_m
-
-
-
Index: trunk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90	(revision 3525)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90	(revision 3526)
@@ -1,3 +1,2 @@
-
 module allinterfaces
   ! interfaces from Fortran 90 subroutines and functions
@@ -86,5 +85,5 @@
        real(8), intent(OUT) :: ypp(nz), zdepthG
        integer, intent(INOUT) :: typeG
-       real(8), external :: zint, deriv1_onesided, colint
+       real(8), external :: zint
      end subroutine depths_avmeth
   end interface
@@ -96,26 +95,4 @@
        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
 
@@ -131,110 +108,6 @@
   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)
@@ -250,16 +123,3 @@
   end interface
 
-    interface
-      subroutine dyn_ss_ice_m_wrapper(ngrid,nsoil,tHIn,p0,pfrost,T_in,ssi_depth_in,porefill_in,porefill,ssi_depth)
-           implicit none
-           integer, intent(IN) :: nsoil,ngrid
-           real(8),  intent(IN) :: thIn(ngrid),ssi_depth_in(ngrid)
-           real(8),  intent(IN) :: p0(ngrid), pfrost(ngrid)
-           real(8),  intent(IN) :: T_in(nsoil,ngrid)
-           real(8), intent(OUT) :: porefill(nsoil,ngrid)
-           real(8), intent(IN) :: porefill_in(nsoil,ngrid)
-           real(8), intent(OUT) :: ssi_depth(ngrid)
-   end subroutine dyn_ss_ice_m_wrapper
-  end interface
-
 end module allinterfaces
Index: trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90	(revision 3525)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90	(revision 3526)
@@ -294,17 +294,9 @@
 
 
-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
-
-
-
+real*8 function tfrostco2(p)
+!     the inverse of function psvco2
+!     input is pressure in Pascal, output is temperature in Kelvin
+implicit none
+real*8 p
+tfrostco2 = 3182.48/(23.3494+log(100./p))
+end function
Index: trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90	(revision 3525)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90	(revision 3526)
@@ -55,4 +55,5 @@
 !***********************************************************************
   use allinterfaces, except_this_one => depths_avmeth
+  use math_mod, only: deriv2_simple, deriv1_onesided, deriv1, colint
   implicit none
   integer, intent(IN) :: nz, typeT
@@ -118,5 +119,5 @@
   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(:))
+     call deriv1_onesided(typeP,z,nz,eta(:),ap_one)
      ! print *,typeP,ap(typeP),ap_one
      ap(typeP)=ap_one
@@ -151,9 +152,9 @@
   endif
   if (typeG>0 .and. typeT<0) then
-     cumfillabove = colint(porefill(:)/eta(:),z,nz,typeG-1,nz) 
+     call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove) 
      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)
+        call colint(porefill(:)/eta(:),z,nz,i,nz,cumfill)
         if (cumfill<yp(i)*18./8314.*B) then  ! usually executes on i=typeG
            if (i>typeG) then
@@ -194,55 +195,4 @@
 
 
-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)
@@ -254,4 +204,5 @@
   use miscparameters, only : icedensity
   use allinterfaces, except_this_one => icechanges
+  use math_mod, only: colint
   implicit none
   integer, intent(IN) :: nz, typeF, typeG
@@ -293,5 +244,5 @@
         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)
+        call colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j,integ)
         erase = j
         if (integ>B*avdrhoP*18./8314.) exit
@@ -332,52 +283,2 @@
   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/NS_flux_mars.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_flux_mars.F90	(revision 3525)
+++ 	(revision )
@@ -1,125 +1,0 @@
-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/NS_generalorbit.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_generalorbit.F90	(revision 3525)
+++ 	(revision )
@@ -1,54 +1,0 @@
-!=============================================================
-! Subroutine to return distance, longitude, and declination of 
-! the sun in planetocentric coordinates from orbital elements
-!
-! INPUTS: 
-! edays = time since perihelion (earth days)
-! a = semimajor axis (AU)
-! ecc = eccentricity
-! omega = Ls of perihelion, relative to equinox (radians)
-! eps = obliquity (radians)
-!
-! OUTPUTS:
-! Ls = areocentric longitude (radians)
-! dec = planetocentric solar declination (radians)
-! r = heliocentric distance (AU)
-!=============================================================
-
-      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
-
-!     T = orbital period (days)
-      T = sqrt(4*pi**2/(6.674e-11*1.989e30)*(a*149.598e9)**3)/86400.
-
-!     M = mean anomaly (radians)
-      M = 2.*pi*edays/T  ! M=0 at perihelion
-
-!     E = eccentric anomaly 
-!     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
-
-!     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/NS_grids.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_grids.F90	(revision 3525)
+++ 	(revision )
@@ -1,96 +1,0 @@
-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/NS_misc_params.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_misc_params.F90	(revision 3525)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_misc_params.F90	(revision 3526)
@@ -12,3 +12,2 @@
   !   length of Earth year in days = 365.24
 end module miscparameters
-
Index: trunk/LMDZ.COMMON/libf/evolution/NS_psvco2.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_psvco2.F90	(revision 3525)
+++ 	(revision )
@@ -1,50 +1,0 @@
-      real*8 function psvco2(T)
-!     solid-vapor transition for CO2
-!     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)
-!     the inverse of function psvco2
-!     input is pressure in Pascal, output is temperature in Kelvin
-      implicit none
-      real*8 p
-      tfrostco2 = 3182.48/(23.3494+log(100./p))
-      end
-
-
-
-
-!------------------------------------------------------------------------------
-!     Antoine equation parameters from NIST, 154K-196K 
-!     based on Giauque and Egan (1937)
-!     A=6.81228, B=1301.679, C=-3.494
-!     p = 10**(A-(B/(T+C)))*1.e5
-
-!     Expressions from Int. Crit. Tabl. 3.207, based on many references
-!     mm2Pa = 133.32
-!     -135C to -56.7C (138K-216K)
-!     p = 10**(-0.05223*26179.3/T + 9.9082)*mm2Pa
-!     -183C to -135C (90K-138K)
-!     p = 10**(-1275.62/T + 0.006833*T + 8.3071)*mm2Pa
-!     Expressions from Int. Crit. Tabl. 3.208, based on Henning
-!     -110 to -80C (163K-193K)
-!     p = 10**(- 1279.11/T + 1.75*log10(T) - 0.0020757*T + 5.85242)*mm2Pa
-!     p = 10**(- 1352.6/T + 9.8318)*mm2Pa
-      
-!     Mars book (1992), p959, fit by chapter authors
-!     p = exp(23.3494 - 3182.48/T)*100.   ! 120-160K
-!     p = exp(25.2194 - 3311.57/T - 6.71e-3*T)*100  ! 100-170K
-!     Mars book (1992), p960, based on Miller & Smythe (1970)
-!     p = exp(26.1228 - 3385.26/T - 9.4461e-3*T)*100 ! 123-173K
-
-!     Fray & Schmitt, PSS 57, 2053 (2009)
-!     A0=1.476e1, A1=-2.571e3, A2=-7.781e4, A3=4.325e6, A4=-1.207e8, A5=1.350e9
-!     p = exp(A0+A1/T+A2/T**2+A3/T**3+A4/T**4+A5/T**5)*1e5 ! 40K-194.7K
-!     A0=1.861e1, A1=-4.154e3, A2=1.041e5
-!     p = exp(A0 + A1/T + A2/T**2)*1e5  ! 194.7K-216.58K
-!------------------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90	(revision 3525)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90	(revision 3526)
@@ -76,77 +76,2 @@
   
 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/NS_tridag.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_tridag.F90	(revision 3525)
+++ 	(revision )
@@ -1,30 +1,0 @@
-!================================================
-! Tridiagonal solver
-!================================================
-      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 
-!      if(n.gt.NMAX) then
-!         print *, 'tridag: too many points, set NMAX>',n
-!         stop
-!      endif 
-      bet=b(1)
-      u(1)=r(1)/bet
-      do j=2,n
-        gam(j)=c(j-1)/bet
-        bet=b(j)-a(j)*gam(j)
-!        if(bet.eq.0.)pause 'tridag failed'
-        u(j)=(r(j)-a(j)*u(j-1))/bet
-      enddo
-      do j=n-1,1,-1
-        u(j)=u(j)-gam(j+1)*u(j+1)
-      enddo
-!      return
-      END SUBROUTINE tridag
-!  (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31&#5(+.
Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3525)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3526)
@@ -492,2 +492,5 @@
 == 19/11/2024 == JBC
 Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings.
+
+== 19/11/2024 == JBC
+Removing unused or redundant Norbert Schorghofer's subroutines (follow-up of r3493) + cleaning and some modifications of related subroutines.
Index: trunk/LMDZ.COMMON/libf/evolution/math_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/math_mod.F90	(revision 3525)
+++ trunk/LMDZ.COMMON/libf/evolution/math_mod.F90	(revision 3526)
@@ -2,14 +2,17 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
-!!! Purpose: The module contains all the mathematical subroutine used in the PEM
+!!! Purpose: The module contains all the mathematical SUBROUTINE used in the PEM
 !!!           
 !!! Author: Adapted from Schorgofer MSIM (N.S, Icarus 2010), impletented here by LL
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  implicit none
+
+implicit none
+
+!=======================================================================
   contains
+!=======================================================================
 
-
-subroutine deriv1(z,nz,y,y0,ybot,dzY)
+SUBROUTINE deriv1(z,nz,y,y0,ybot,dzY)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -19,38 +22,39 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        implicit none
+! first derivative of a function y(z) on irregular grid
+! upper boundary conditions: y(0) = y0
+! lower boundary condition.: yp = ybottom
 
-! first derivative of a function y(z) on irregular grid
-! upper boundary conditions: y(0)=y0
-! lower boundary condition.: yp = ybottom
-  integer, intent(IN) :: nz ! number of layer
-  real, intent(IN) :: z(nz) ! depth layer
-  real, intent(IN) :: y(nz) ! function which needs to be derived
-  real, intent(IN) :: y0,ybot ! boundary conditions
-  real, intent(OUT) :: dzY(nz) ! derivative of y w.r.t depth
+implicit none
+
+integer, intent(in) :: nz      ! number of layer
+real,    intent(in) :: z(nz)   ! depth layer
+real,    intent(in) :: y(nz)   ! function which needs to be derived
+real,    intent(in) :: y0,ybot ! boundary conditions
+real, intent(out) :: dzY(nz) ! derivative of y w.r.t depth
 ! local
-  integer :: j
-  real :: hm,hp,c1,c2,c3
+integer :: j
+real    :: hm, hp, c1, c2, c3
 
-  hp = z(2)-z(1)
-  c1 = z(1)/(hp*z(2))
-  c2 = 1/z(1) - 1/(z(2)-z(1))
-  c3 = -hp/(z(1)*z(2))
-  dzY(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)))
-     dzY(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
-  enddo
-  dzY(nz) = (ybot - y(nz-1))/(2.*(z(nz)-z(nz-1)))
- return
-end subroutine deriv1
+hp = z(2) - z(1)
+c1 = z(1)/(hp*z(2))
+c2 = 1/z(1) - 1/(z(2) - z(1))
+c3 = -hp/(z(1)*z(2))
+dzY(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)))
+    dzY(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
+enddo
+dzY(nz) = (ybot - y(nz - 1))/(2.*(z(nz) - z(nz - 1)))
 
+END SUBROUTINE deriv1
 
+!=======================================================================
 
-subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
+SUBROUTINE deriv2_simple(z,nz,y,y0,yNp1,yp2)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -60,37 +64,37 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! second derivative y_zz on irregular grid
+! boundary conditions: y(0) = y0, y(nz + 1) = yNp1
 
-  ! second derivative y_zz on irregular grid
-  ! boundary conditions: y(0)=y0, y(nz+1)=yNp1
-  implicit none
-  integer, intent(IN) :: nz
-  real, intent(IN) :: z(nz),y(nz),y0,yNp1
-  real, intent(OUT) :: yp2(nz)
-  integer j
-  real hm,hp,c1,c2,c3
+implicit none
 
-  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
+integer, intent(in) :: nz
+real,    intent(in) :: z(nz), y(nz), y0, yNp1
+real, intent(out) :: yp2(nz)
+integer :: j
+real    :: hm, hp, c1, c2, c3
 
-  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
-  return
-end subroutine deriv2_simple
+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
 
+!=======================================================================
 
-subroutine  deriv1_onesided(j,z,nz,y,dy_zj)
+SUBROUTINE  deriv1_onesided(j,z,nz,y,dy_zj)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
-!!! Purpose:   First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
+!!! Purpose: First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
 !!!           
 !!! Author: N.S (raw copy/paste from MSIM)
@@ -99,23 +103,26 @@
 
   implicit none
-  integer, intent(IN) :: nz,j
-  real, intent(IN) :: z(nz),y(nz)
-  real, intent(out) :: dy_zj
-  real h1,h2,c1,c2,c3
-  if (j<1 .or. j>nz-2) then
-     dy_zj = -1.
-  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))
-     dy_zj = c1*y(j) + c2*y(j+1) + c3*y(j+2)
-  endif
-  return
-end subroutine deriv1_onesided
 
+integer, intent(in) :: nz, j
+real,    intent(in) :: z(nz), y(nz)
+real, intent(out) :: dy_zj
+real :: h1, h2, c1, c2, c3
 
-subroutine  colint(y,z,nz,i1,i2,integral)
+if (j < 1 .or. j > nz - 2) then
+    dy_zj = -1.
+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))
+    dy_zj = c1*y(j) + c2*y(j + 1) + c3*y(j + 2)
+endif
+
+END SUBROUTINE deriv1_onesided
+
+!=======================================================================
+
+PURE SUBROUTINE colint(y,z,nz,i1,i2,integral)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -126,39 +133,41 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
-  implicit none
-  integer, intent(IN) :: nz, i1, i2
-  real, intent(IN) :: y(nz), z(nz)
-  real,intent(out) :: integral
-  integer i
-  real dz(nz)
+implicit none
   
-  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)
-  integral = sum(y(i1:i2)*dz(i1:i2))
-end subroutine colint
+integer, intent(in) :: nz, i1, i2
+real,    intent(in) :: y(nz), z(nz)
+real, intent(out) :: integral
+integer i
+real 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)
+integral = sum(y(i1:i2)*dz(i1:i2))
 
+END SUBROUTINE colint
 
-
+!=======================================================================
 
 SUBROUTINE findroot(y1,y2,z1,z2,zr)
-        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        !!!
-        !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
-        !!! 
-        !!! Author: LL
-        !!!
-        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
+!!! 
+!!! Author: LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-        implicit none
-        real,intent(in) :: y1,y2 ! difference between surface water density and at depth  [kg/m^3]
-        real,intent(in) :: z1,z2 ! depth [m]
-        real,intent(out) :: zr ! depth at which we have zero
-        zr = (y1*z2 - y2*z1)/(y1-y2)
-        RETURN
-        end subroutine findroot
+implicit none
+
+real, intent(in) :: y1, y2 ! difference between surface water density and at depth  [kg/m^3]
+real, intent(in) :: z1, z2 ! depth [m]
+real, intent(out) :: zr    ! depth at which we have zero
+
+zr = (y1*z2 - y2*z1)/(y1 - y2)
+
+END SUBROUTINE findroot
 
 end module math_mod
