Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 1515)
+++ /trunk/LMDZ.GENERIC/README	(revision 1516)
@@ -1117,2 +1117,6 @@
 - Added a generic routine (in utilities directory) to extrapolate ice fields.
 - Added a generic version of streamfunction (in utilities directory).
+
+== 17/03/2016 == EM
+- Added tests to ensure that soil model settings are adequate to resolve
+  sub-surface diurnal and annual thermal waves.
Index: /trunk/LMDZ.GENERIC/libf/phystd/comsoil_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/comsoil_h.F90	(revision 1515)
+++ /trunk/LMDZ.GENERIC/libf/phystd/comsoil_h.F90	(revision 1516)
@@ -6,4 +6,8 @@
 !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
 
   real,save,allocatable,dimension(:) :: layer      ! soil layer depths
Index: /trunk/LMDZ.GENERIC/libf/phystd/soil.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/soil.F	(revision 1515)
+++ /trunk/LMDZ.GENERIC/libf/phystd/soil.F	(revision 1516)
@@ -4,5 +4,8 @@
      &          capcal,fluxgrd)
 
-      use comsoil_h, only: layer, mlayer, volcapa
+      use comsoil_h, only: layer, mlayer, volcapa, inertiedat
+      use comcstfi_mod, only: pi, daysec
+      use planete_mod, only: year_day
+      use comgeomfi_h, only: long, lati
 
       implicit none
@@ -49,5 +52,7 @@
 ! local variables:
       integer ig,ik
-
+      real :: inertia_min,inertia_max
+      real :: diurnal_skin ! diurnal skin depth (m)
+      real :: annual_skin ! anuual skin depth (m)
 
 ! 0. Initialisations and preprocessing step
@@ -118,5 +123,38 @@
       !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
       enddo ! of do ig=1,ngrid
-            
+      
+      ! Additional checks: is the vertical discretization sufficient
+      ! to resolve diurnal and annual waves?
+      do ig=1,ngrid
+        ! extreme inertia for this column
+        inertia_min=minval(inertiedat(ig,:))
+        inertia_max=maxval(inertiedat(ig,:))
+        ! diurnal and annual skin depth
+        diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
+        annual_skin=(inertia_max/volcapa)*sqrt(year_day*daysec/pi)
+        if (0.5*diurnal_skin<layer(1)) then
+        ! one should have the fist layer be at least half of diurnal skin depth
+          write(*,*) "soil Error: grid point ig=",ig
+          write(*,*) "            longitude=",long(ig)*(180./pi)
+          write(*,*) "             latitude=",lati(ig)*(180./pi)
+          write(*,*) "  first soil layer depth ",layer(1)
+          write(*,*) "  not small enough for a diurnal skin depth of ",
+     &                diurnal_skin
+          write(*,*) " change soil layer distribution (comsoil_h.F90)"
+          stop
+        endif
+        if (2.*annual_skin>layer(nsoil)) then
+        ! one should have the full soil be at least twice the diurnal skin depth
+          write(*,*) "soil Error: grid point ig=",ig
+          write(*,*) "            longitude=",long(ig)*(180./pi)
+          write(*,*) "             latitude=",lati(ig)*(180./pi)
+          write(*,*) "  total soil layer depth ",layer(nsoil)
+          write(*,*) "  not large enough for an annual skin depth of ",
+     &                annual_skin
+          write(*,*) " change soil layer distribution (comsoil_h.F90)"
+          stop
+        endif
+      enddo ! of do ig=1,ngrid
+      
       else ! of if (firstcall)
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F	(revision 1515)
+++ /trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F	(revision 1516)
@@ -2,5 +2,6 @@
 
 !      use netcdf
-      use comsoil_h, only: layer, mlayer, inertiedat, volcapa
+      use comsoil_h, only: layer, mlayer, inertiedat, volcapa,
+     &                     lay1, alpha
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
@@ -70,5 +71,5 @@
       real,dimension(:),allocatable :: newval
 
-      real alpha,lay1 ! coefficients for building layers
+      real malpha,mlay1 ! coefficients for building layers
       real xmin,xmax ! to display min and max of a field
 
@@ -133,6 +134,4 @@
       ! default mlayer distribution, following a power law:
       !  mlayer(k)=lay1*alpha**(k-1/2)
-        lay1=2.e-4
-	alpha=2
         do iloop=0,nsoil-1
 	  mlayer(iloop)=lay1*(alpha**(iloop-0.5))
@@ -142,9 +141,21 @@
       ! Assuming layer distribution follows mid-layer law:
       ! layer(k)=lay1*alpha**(k-1)
-      lay1=sqrt(mlayer(0)*mlayer(1))
-      alpha=mlayer(1)/mlayer(0)
+      mlay1=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
+         stop
+      endif
+      if ((abs(malpha-alpha)/alpha).gt.1.e-6) then
+         write(*,*) "soil_settings error: malpha=",malpha
+         write(*,*) " does not match comsoil_h alpha=",alpha
+         stop
+      endif
       do iloop=1,nsoil
-        layer(iloop)=lay1*(alpha**(iloop-1))
+        layer(iloop)=mlay1*(malpha**(iloop-1))
       enddo
+
 ! 2. Volumetric heat capacity (note: it is declared in comsoil_h)
 ! ---------------------------
