Index: trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/newstart.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/newstart.F	(revision 1537)
+++ trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/newstart.F	(revision 1538)
@@ -36,4 +36,5 @@
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
       use inifis_mod, only: inifis
+      use control_mod, only: nday 
       implicit none
 
@@ -95,5 +96,5 @@
 c------------------
       REAL tsurf(ngridmx)	! surface temperature
-      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
+      REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
 !      REAL co2ice(ngridmx)	! CO2 ice layer
       REAL emis(ngridmx)	! surface emissivity
@@ -103,5 +104,5 @@
 !      REAL rnaturfi(ngridmx)
       real alb(iip1,jjp1),albfi(ngridmx) ! albedos
-      real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D)
+      real,allocatable :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
       real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
       REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
@@ -113,5 +114,5 @@
 
       REAL rnat(ngridmx)
-      REAL tslab(ngridmx,nsoilmx) ! slab ocean temperature
+      REAL,ALLOCATABLE :: tslab(:,:) ! slab ocean temperature
       REAL pctsrf_sic(ngridmx) ! sea ice cover
       REAL tsea_ice(ngridmx) ! temperature sea_ice
@@ -274,4 +275,13 @@
 ! ... this has to be done before phyetat0
       call initracer(ngridmx,nqtot,tname)
+
+      ! default nsoilmx is set in comsoil_h
+      call getin_p("nsoilmx",nsoilmx)
+
+      ! allocate some local arrays now that we know all dimensions
+      allocate(tsoil(ngridmx,nsoilmx))
+      allocate(ith(iip1,jjp1,nsoilmx))
+      allocate(ithfi(ngridmx,nsoilmx))
+      allocate(tslab(ngridmx,nsoilmx))
 
 ! Take care of arrays in common modules
@@ -475,5 +485,5 @@
 !    (for instance initracer needs to know about some flags, and/or
 !      'datafile' path may be changed by user)
-      call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
+      call inifis(ngridmx,llm,nqtot,day_ini,daysec,nday,dtphys,
      &                latfi,lonfi,airefi,rad,g,r,cpp)
 
Index: trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/start2archive.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/start2archive.F	(revision 1537)
+++ trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/start2archive.F	(revision 1538)
@@ -19,4 +19,5 @@
 c=======================================================================
 
+      use ioipsl_getin_p_mod, only: getin_p
       use infotrac, only: infotrac_init, nqtot, tname
       USE comsoil_h
@@ -69,5 +70,5 @@
 c ------------------------------------
       REAL tsurf(ngridmx)	! Surface temperature
-      REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
+      REAL,ALLOCATABLE :: tsoil(:,:) ! Soil temperature
       REAL co2ice(ngridmx)	! CO2 ice layer
       REAL q2(ngridmx,llm+1)
@@ -93,6 +94,6 @@
       REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
       REAL tsurfS(ip1jmp1)
-      REAL tsoilS(ip1jmp1,nsoilmx)
-      REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
+      REAL,ALLOCATABLE :: tsoilS(:,:)
+      REAL,ALLOCATABLE :: ithS(:,:) ! Soil Thermal Inertia
       REAL co2iceS(ip1jmp1)
       REAL q2S(ip1jmp1,llm+1)
@@ -147,4 +148,12 @@
       CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
       call initcomgeomphy
+
+! get value of nsoilmx (default is set in comsoil_h)
+      call getin_p("nsoilmx",nsoilmx)
+      
+      ! allocate some local arrays now that we know all dimensions
+      allocate(tsoil(ngridmx,nsoilmx))
+      allocate(tsoilS(ip1jmp1,nsoilmx))
+      allocate(ithS(ip1jmp1,nsoilmx))
 
       ! ALLOCATE ARRAYS IN comgeomfi_h (usually done in inifis)
Index: trunk/LMDZ.GENERIC/libf/phystd/comsoil_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/comsoil_h.F90	(revision 1537)
+++ trunk/LMDZ.GENERIC/libf/phystd/comsoil_h.F90	(revision 1538)
@@ -5,9 +5,10 @@
 !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case 
 !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
-  integer, parameter :: nsoilmx = 18
-! Full soil layer depths are set as: layer(k)=lay1*alpha**(k-1) , k=1,nsoil
-! Mid soil layer depths are set as: mlayer(k)=lay1*alpha**(k-1/2) , k=0,nsoil-1 
-  real,parameter :: lay1=2.e-4 ! depth of first full soil layer
-  real,parameter :: alpha=2
+  integer,save :: nsoilmx = 18 ! default, but may be set in callphys.def
+! Full soil layer depths are set as: layer(k)=lay1_soil*alpha_soil**(k-1) , k=1,nsoil
+! Mid soil layer depths are set as: mlayer(k)=lay1_soil*alpha_soil**(k-1/2) , k=0,nsoil-1 
+  real,save :: lay1_soil=2.e-4 ! depth (m) of first full soil layer (may be set in callphys.def)
+  real,save :: alpha_soil=2 ! coefficient for soil layer thickness (may be set in callphys.def)
+!$OMP THREADPRIVATE(nsoilmx,lay1_soil,alpha_soil)
 
   real,save,allocatable,dimension(:) :: layer      ! soil layer depths
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 1537)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 1538)
@@ -14,5 +14,5 @@
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
   use comgeomfi_h, only: long, lati, area, totarea, totarea_planet
-  use comsoil_h, only: ini_comsoil_h
+  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
   use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
                               init_time, daysec, dtphys
@@ -322,4 +322,20 @@
      call getin_p("graybody",graybody)
      write(*,*)" graybody = ",graybody
+
+! Soil model
+     write(*,*)"Number of sub-surface layers for soil scheme?"
+     ! default value of nsoilmx set in comsoil_h
+     call getin_p("nsoilmx",nsoilmx)
+     write(*,*)" nsoilmx=",nsoilmx
+     
+     write(*,*)"Thickness of topmost soil layer (m)?"
+     ! default value of lay1_soil set in comsoil_h
+     call getin_p("lay1_soil",lay1_soil)
+     write(*,*)" lay1_soil = ",lay1_soil
+     
+     write(*,*)"Coefficient for soil layer thickness distribution?"
+     ! default value of alpha_soil set in comsoil_h
+     call getin_p("alpha_soil",alpha_soil)
+     write(*,*)" alpha_soil = ",alpha_soil
 
 ! Slab Ocean 
Index: trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F	(revision 1537)
+++ trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F	(revision 1538)
@@ -3,5 +3,5 @@
 !      use netcdf
       use comsoil_h, only: layer, mlayer, inertiedat, volcapa,
-     &                     lay1, alpha
+     &                     lay1_soil, alpha_soil
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
@@ -71,5 +71,5 @@
       real,dimension(:),allocatable :: newval
 
-      real malpha,mlay1 ! coefficients for building layers
+      real malpha,mlay1_soil ! coefficients for building layers
       real xmin,xmax ! to display min and max of a field
 
@@ -84,16 +84,36 @@
 
         if (dimlen.ne.nsoil) then
-          write(*,*)'soil_settings: Interpolation of soil temperature ',
-     &              'and thermal inertia will be required!'
 	! if dimlen doesn't match nsoil, then interpolation of
 	! soil temperatures and thermal inertia will be requiered
 	  interpol=.true.
-        ! allocate olmlayer  
-          allocate(oldmlayer(dimlen),stat=ierr)
-          if (ierr.ne.0) then
-            write(*,*) 'soil_settings: failed allocation of oldmlayer!'
-            stop
-          endif
 	endif
+
+        ! allocate oldmlayer  
+        allocate(oldmlayer(dimlen),stat=ierr)
+        if (ierr.ne.0) then
+          write(*,*) 'soil_settings: failed allocation of oldmlayer!'
+          stop
+        endif
+
+        ! check if olmlayer distribution matches current one
+        call get_var("soildepth",oldmlayer,found)
+        if (found) then
+          malpha=oldmlayer(2)/oldmlayer(1)
+          if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then
+            ! alpha values are too different, intepolation needed
+            interpol=.true.
+          endif
+          ! check if loaded mid-layer depth value differs from current one
+          if (abs((oldmlayer(1)-lay1_soil*alpha_soil**(-1./2.))/
+     &             (lay1_soil*alpha_soil**(-1./2.))).gt.1.e-6) then
+            interpol=.true.
+          endif
+        endif
+
+        if (interpol) then
+          write(*,*)'soil_settings: Interpolation of soil temperature ',
+     &              'and thermal inertia will be required!'
+        endif
+        
 ! 1.2 Find out the # of dimensions <inertiedat> was defined as using
       ndims=inquire_field_ndims("inertiedat")
@@ -133,27 +153,27 @@
       if (interpol) then
       ! default mlayer distribution, following a power law:
-      !  mlayer(k)=lay1*alpha**(k-1/2)
+      !  mlayer(k)=lay1_soil*alpha_soil**(k-1/2)
         do iloop=0,nsoil-1
-	  mlayer(iloop)=lay1*(alpha**(iloop-0.5))
+	  mlayer(iloop)=lay1_soil*(alpha_soil**(iloop-0.5))
         enddo
       endif
 ! 1.5 Build layer(); following the same law as mlayer()
       ! Assuming layer distribution follows mid-layer law:
-      ! layer(k)=lay1*alpha**(k-1)
-      mlay1=sqrt(mlayer(0)*mlayer(1))
+      ! layer(k)=lay1_soil*alpha_soil**(k-1)
+      mlay1_soil=sqrt(mlayer(0)*mlayer(1))
       malpha=mlayer(1)/mlayer(0)
       ! Check that these values are the same as those prescibed for mlayers
-      if ((abs(mlay1-lay1)/lay1).gt.1.e-6) then
-         write(*,*) "soil_settings error: mlay1=",mlay1
-         write(*,*) " does not match comsoil_h lay1=",lay1
+      if ((abs(mlay1_soil-lay1_soil)/lay1_soil).gt.1.e-6) then
+         write(*,*) "soil_settings error: mlay1_soil=",mlay1_soil
+         write(*,*) " does not match comsoil_h lay1_soil=",lay1_soil
          stop
       endif
-      if ((abs(malpha-alpha)/alpha).gt.1.e-6) then
+      if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then
          write(*,*) "soil_settings error: malpha=",malpha
-         write(*,*) " does not match comsoil_h alpha=",alpha
+         write(*,*) " does not match comsoil_h alpha_soil=",alpha_soil
          stop
       endif
       do iloop=1,nsoil
-        layer(iloop)=mlay1*(malpha**(iloop-1))
+        layer(iloop)=mlay1_soil*(malpha**(iloop-1))
       enddo
 
