Index: /trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90
===================================================================
--- /trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90	(revision 3482)
+++ /trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90	(revision 3483)
@@ -8,5 +8,5 @@
                      ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
                      day_ini,time,tsurf,tsoil, &
-                     emis,q2,qsurf)
+                     emis,q2,qsurf,therm_inertia)
                     !  ,cloudfrac,totcloudfrac,hice, &
                     !  rnat,pctsrf_sic)
@@ -18,4 +18,5 @@
   use write_field_phy, only: Writefield_phy
 !!
+  use comsoil_h, only: nsoilmx
   use tabfi_mod, only: tabfi
   USE tracer_h, ONLY: noms
@@ -46,6 +47,7 @@
   real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
   real,intent(out) :: emis(ngrid) ! surface emissivity
-  real,intent(out) :: q2(ngrid,nlayer+1) ! 
+  real,intent(out) :: q2(ngrid,nlayer+1) !
   real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
+  real,intent(out) :: therm_inertia(ngrid,nsoilmx) ! thermal inertia
 ! real n2ice(ngrid) ! n2 ice cover
 
@@ -62,5 +64,5 @@
       INTEGER nid, nvarid
       INTEGER ierr, i, nsrf
-!      integer isoil 
+!      integer isoil
 !      INTEGER length
 !      PARAMETER (length=100)
@@ -69,5 +71,5 @@
       CHARACTER*1 yes
 !
-      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
+      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,TI
       INTEGER nqold
 
@@ -76,8 +78,8 @@
       integer :: count
       character(len=30) :: txt ! to store some text
-      
+
       INTEGER :: indextime=1 ! index of selected time, default value=1
       logical :: found
-      
+
       character(len=8) :: modname="phyetat0"
 
@@ -265,4 +267,8 @@
     call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
   endif
+else
+  call getin_p("therm_inertia",TI)
+  therm_inertia(:,:) = TI
+    !AF24
 endif ! of if (startphy_file)
 
Index: /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
===================================================================
--- /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3482)
+++ /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3483)
@@ -28,5 +28,5 @@
       use comdiurn_h, only: coslat, sinlat, coslon, sinlon
       use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
-      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
+      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
       use geometry_mod, only: latitude, longitude, cell_area
       USE comgeomfi_h, only: totarea, totarea_planet
@@ -263,5 +263,4 @@
       REAL,save  :: tcond1p4Pa
       DATA tcond1p4Pa/38/
-      real :: tidat(ngrid,nsoilmx)    ! thermal inertia soil
 
 
@@ -523,4 +522,10 @@
 #endif
 
+      ! local variables for skin depth check
+      real :: inertia_min,inertia_max
+      real :: diurnal_skin ! diurnal skin depth (m)
+      real :: annual_skin ! anuual skin depth (m)
+
+
       ! flags to trigger extra sanity checks
       logical, save ::  check_physics_inputs=.false.
@@ -585,5 +590,5 @@
          call phyetat0(startphy_file,                                 &
                        ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
-                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
+                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
 
 #else
@@ -999,10 +1004,47 @@
       end if ! if ngrid ne 1
 
-      call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,   &
-      capcal,adjust,dist_star,albedo,emis,flusurfold,ptimestep,   &
-      zls)
+      call surfprop(ngrid,nq,fract,qsurf,tsurf,        &
+      capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
+      albedo,emis,inertiedat)
       ! do i=2,ngrid
       !    albedo(i,:) = albedo(1,:)
       ! enddo
+
+      if (callsoil) then
+         ! AF24 Originally in soil.F, but inertiedat is modified by surfprop
+         ! 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=",longitude(ig)*(180./pi)
+            write(*,*) "             latitude=",latitude(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)"
+            call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
+            endif
+            if (2.*annual_skin>layer(nsoilmx)) then
+            ! one should have the full soil be at least twice the diurnal skin depth
+            write(*,*) "soil Error: grid point ig=",ig
+            write(*,*) "            longitude=",longitude(ig)*(180./pi)
+            write(*,*) "             latitude=",latitude(ig)*(180./pi)
+            write(*,*) "  total soil layer depth ",layer(nsoilmx)
+            write(*,*) "  not large enough for an annual skin depth of ", &
+                        annual_skin
+            write(*,*) " change soil layer distribution (comsoil_h.F90)"
+            call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1)
+            endif
+         enddo ! of do ig=1,ngrid
+
+      end if ! callsoil
 
       if (callrad) then
@@ -1795,5 +1837,5 @@
       ! tidat_out(:,:)=0.
       ! DO l=1,min(nlayermx,nsoilmx)
-      !    tidat_out(:,l)=tidat(:,l)
+      !    tidat_out(:,l)=inertiedat(:,l)
       ! ENDDO
 
@@ -2217,5 +2259,5 @@
             !      ptimestep,pdaypal, &
             !      ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, &
-            !      cell_area,albedodat,tidat,zmea,zstd,zsig, &
+            !      cell_area,albedodat,inertiedat,zmea,zstd,zsig, &
             !      zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal,  &
             !      peri_daypal)
Index: /trunk/LMDZ.PLUTO/libf/phypluto/soil.F90
===================================================================
--- /trunk/LMDZ.PLUTO/libf/phypluto/soil.F90	(revision 3482)
+++ /trunk/LMDZ.PLUTO/libf/phypluto/soil.F90	(revision 3483)
@@ -1,6 +1,6 @@
-      subroutine soil(ngrid,nsoil,firstcall,lastcall,
-     &          therm_i,
-     &          timestep,tsurf,tsoil,
-     &          capcal,fluxgrd)
+      subroutine soil(ngrid,nsoil,firstcall,lastcall,   &
+               therm_i, &
+               timestep,tsurf,tsoil,    &
+               capcal,fluxgrd)
 
       use comsoil_h, only: layer, mlayer, volcapa, inertiedat
@@ -21,5 +21,5 @@
 !-----------------------------------------------------------------------
 
-c-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
 !  arguments
 !  ---------
@@ -49,7 +49,4 @@
 ! 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
@@ -76,7 +73,7 @@
       do ig=1,ngrid
         do ik=1,nsoil-1
-      thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)
-     &                +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))
-     &                    /(mlayer(ik)-mlayer(ik-1))
+      thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)  &
+                     +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))   &
+                         /(mlayer(ik)-mlayer(ik-1))
 !	write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
 	enddo
@@ -91,6 +88,6 @@
 	! q_{k+1/2}
         do ik=1,nsoil-1
-          coefq(ik)=volcapa*(layer(ik+1)-layer(ik))
-     &                 /timestep
+          coefq(ik)=volcapa*(layer(ik+1)-layer(ik)) &
+                      /timestep
 	enddo
 
@@ -102,54 +99,21 @@
 
 	! alph_{N-1}
-	alph(ig,nsoil-1)=coefd(ig,nsoil-1)/
-     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1))
+	alph(ig,nsoil-1)=coefd(ig,nsoil-1)/ &
+                       (coefq(nsoil-1)+coefd(ig,nsoil-1))
         ! alph_k
         do ik=nsoil-2,1,-1
-	  alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*
-     &                              (1.-alph(ig,ik+1))+coefd(ig,ik))
+	  alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*   &
+                                   (1.-alph(ig,ik+1))+coefd(ig,ik))
 	enddo
 
         ! capcal
 ! Cstar
-        capcal(ig)=volcapa*layer(1)+
-     &              (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
-     &              (timestep*(1.-alph(ig,1)))
+        capcal(ig)=volcapa*layer(1)+    &
+                   (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* &
+                   (timestep*(1.-alph(ig,1)))
 ! Cs
-        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))
+        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*  &
+                              thermdiff(ig,1)/mthermdiff(ig,0))
       !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=",longitude(ig)*(180./pi)
-          write(*,*) "             latitude=",latitude(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=",longitude(ig)*(180./pi)
-          write(*,*) "             latitude=",latitude(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
 
@@ -160,8 +124,8 @@
 ! First layer:
       do ig=1,ngrid
-        tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))/
-     &              (1.+mu*(1.0-alph(ig,1))*
-     &               thermdiff(ig,1)/mthermdiff(ig,0))
+        tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*   &
+                              thermdiff(ig,1)/mthermdiff(ig,0))/    &
+                   (1.+mu*(1.0-alph(ig,1))* &
+                    thermdiff(ig,1)/mthermdiff(ig,0))
       enddo
 ! Other layers:
@@ -177,14 +141,14 @@
 ! Bottom layer, beta_{N-1}
       do ig=1,ngrid
-        beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
-     &                   /(coefq(nsoil-1)+coefd(ig,nsoil-1))
+        beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil) &
+                        /(coefq(nsoil-1)+coefd(ig,nsoil-1))
       enddo
 ! Other layers
       do ik=nsoil-2,1,-1
         do ig=1,ngrid
-	  beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+
-     &                 coefd(ig,ik+1)*beta(ig,ik+1))/
-     &                 (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1))
-     &                  +coefd(ig,ik))
+	  beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+    &
+                      coefd(ig,ik+1)*beta(ig,ik+1))/    &
+                      (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1)) &
+                       +coefd(ig,ik))
 	enddo
       enddo
@@ -206,6 +170,6 @@
 !         print*,'tsoil=',tsoil(ig,1)
 
-        fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
-     &              (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
+        fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*    &
+                   (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
 
 !        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
@@ -213,9 +177,9 @@
 !     &                         thermdiff(ig,1)/mthermdiff(ig,0))
 ! Fs
-        fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*
-     &              (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))
-     &               -tsurf(ig)-mu*beta(ig,1)*
-     &                          thermdiff(ig,1)/mthermdiff(ig,0))
+        fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*  &
+                   (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*    &
+                              thermdiff(ig,1)/mthermdiff(ig,0)) &
+                    -tsurf(ig)-mu*beta(ig,1)*   &
+                               thermdiff(ig,1)/mthermdiff(ig,0))
       enddo
 
Index: /trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90
===================================================================
--- /trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90	(revision 3482)
+++ /trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90	(revision 3483)
@@ -1,5 +1,5 @@
-      subroutine surfprop(ngrid,nq,pfra,qsurf,tsurface,     &
-                     tidat,capcal,adjust,dist,albedo,emis,fluold, &
-                     ptimestep,zls)
+      subroutine surfprop(ngrid,nq,pfra,qsurf,tsurface,       &
+                     capcal,adjust,dist,fluold,ptimestep,zls, &
+                     albedo,emis,therm_inertia)
 
     !   use comgeomfi_h, only:
@@ -43,4 +43,5 @@
 !     albedo(ngrid)
 !     emis(ngrid)
+!     therm_inertia(ngrid,nsoilmx)
 !
 !     Both
@@ -56,10 +57,12 @@
 !     Arguments
 
+
       INTEGER ngrid, nq
+      REAL,INTENT(IN) :: pfra(ngrid)
       REAL,INTENT(IN) :: qsurf(ngrid,nq)
+      REAL,INTENT(IN) :: tsurface(ngrid)
       REAL,INTENT(IN) :: fluold(ngrid,nq)
       REAL,INTENT(IN) :: ptimestep
       REAL,INTENT(IN) :: zls
-      REAL,INTENT(IN) :: tsurface(ngrid)
       REAL,INTENT(IN) :: capcal(ngrid)
       REAL,INTENT(IN) :: adjust
@@ -67,6 +70,5 @@
       REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
       REAL,INTENT(OUT) :: emis(ngrid)
-      REAL,INTENT(OUT) :: tidat(ngrid,nsoilmx)
-      REAL,INTENT(IN) :: pfra(ngrid)
+      REAL,INTENT(OUT) :: therm_inertia(ngrid,nsoilmx) ! therm_inertia = inertiedat
 !-----------------------------------------------------------------------
 !     Local variables
@@ -448,13 +450,13 @@
            do isoil=0,nsoilmx-1
               if (mlayer(isoil).le.emin) then ! diurnal TI
-                   tidat(ig,isoil+1)=tid
+                   therm_inertia(ig,isoil+1)=tid
               else if (isoil.gt.0.and.(mlayer(isoil).gt.emin).and.(mlayer(isoil-1).lt.emin)) then ! still diurnal TI
-                   tidat(ig,isoil+1)=tid
+                   therm_inertia(ig,isoil+1)=tid
               else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax)) then ! TI N2
-                   tidat(ig,isoil+1)=ITN2
+                   therm_inertia(ig,isoil+1)=ITN2
               else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax2)) then
-                   tidat(ig,isoil+1)=ITCH4
+                   therm_inertia(ig,isoil+1)=ITCH4
               else
-                   tidat(ig,isoil+1)=ITH2O
+                   therm_inertia(ig,isoil+1)=ITH2O
               endif
 
@@ -465,5 +467,5 @@
 
         DO ig=1,ngrid
-           tidat(ig,:)=inertiedat(ig,:)
+           therm_inertia(ig,:)=inertiedat(ig,:)
         ENDDO
 
