Index: /trunk/LMDZ.MARS/changelog.txt
===================================================================
--- /trunk/LMDZ.MARS/changelog.txt	(revision 3108)
+++ /trunk/LMDZ.MARS/changelog.txt	(revision 3109)
@@ -4297,2 +4297,8 @@
 == 26/10/2023 == JBC
 Few small fixes following r3098.
+
+== 02/11/2023 == LL
+Correction of a small bug in soil_settings: the interpolation on the new grid
+for the thermal inertia was not used with the old thermal inertia vector but
+only with the last value of the oldthermal inertia
+Cleaning of the code
Index: /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 3108)
+++ /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 3109)
@@ -3,5 +3,4 @@
 implicit none
 ! nsoilmx : number of subterranean layers
-!EM: old soil routine:      integer, parameter :: nsoilmx = 10
   integer, parameter :: nsoilmx = 18 
 
Index: /trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3108)
+++ /trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3109)
@@ -476,5 +476,5 @@
     ice_depth = -1 ! default value: no ice
     call getin("subsurface_ice_depth",ice_depth)
-
+    write(*,*) " subsurface_ice_depth = ",ice_depth
 
     z0(1) = z0_default ! default value for roughness
@@ -654,5 +654,4 @@
 
     inertiesoil(1,:,1) = inertiedat(1,:)
-
     tsoil(:,:,1) = tsurf(1,1) ! soil temperature
 endif !(.not. startfiles_1D)
Index: /trunk/LMDZ.MARS/libf/phymars/soil_settings.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/soil_settings.F	(revision 3108)
+++ /trunk/LMDZ.MARS/libf/phymars/soil_settings.F	(revision 3109)
@@ -89,18 +89,4 @@
 ! 1.1 Start by reading how many layers of soil there are
 
-!	ierr=nf90_inq_dimid(nid,"subsurface_layers",dimid)
-!	if (ierr.ne.nf90_noerr) then
-!	 write(*,*)'soil_settings: Problem reading <subsurface_layers>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!	 call abort
-!	endif
-
-!	ierr=nf90_inquire_dimension(nid,dimid,len=dimlen)
-!	if (ierr.ne.nf90_noerr) then
-!	 write(*,*)'soil_settings: Problem getting <subsurface_layers>',
-!     &             'length'
-!	 write(*,*)trim(nf90_strerror(ierr))
-!         call abort
-!	endif
         dimlen=inquire_dimension_length("subsurface_layers")
 
@@ -124,15 +110,5 @@
 ! 1.2 Find out the # of dimensions <inertiedat> was defined as using
 !     (in ye old days, thermal inertia was only given at the "surface")
-      ! Look for thermal inertia data
-!      ierr=nf90_inq_varid(nid,"inertiedat",nvarid)
-!      if (ierr.NE.nf90_noerr) then
-!         write(*,*)'soil_settings: Field <inertiedat> not found!'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         call abort
-!      endif
-!
-!      ! Read the # of dimensions <inertidat> was defined as using
-!      ierr=nf90_inquire_variable(nid,nvarid,ndims=ndims)
-!      ! if (ndims.eq.1) then we have the "old 2D-surface" format
+
       ndims=inquire_field_ndims("inertiedat")
 
@@ -148,18 +124,6 @@
 	enddo
       else ! Look for depth 
-!        ierr=nf90_inq_varid(nid,"soildepth",nvarid)
-!        if (ierr.ne.nf90_noerr) then
-!          write(*,*)'soil_settings: Field <soildepth> not found!'
-!          write(*,*)trim(nf90_strerror(ierr))
-!	  call abort
-!        endif
         ! read <depth> coordinate
         if (interpol) then !put values in oldmlayer
-!          ierr=nf90_get_var(nid,nvarid,oldmlayer)
-!          if (ierr.ne.nf90_noerr) then
-!           write(*,*)'soil_settings: Problem while reading <soildepth>'
-!           write(*,*)trim(nf90_strerror(ierr))
-!           call abort
-!          endif
           call get_var("soildepth",oldmlayer,found)
           if (.not.found) then
@@ -167,10 +131,4 @@
           endif
         else ! put values in mlayer
-!          ierr=nf90_get_var(nid,nvarid,mlayer)
-!          if (ierr.ne.nf90_noerr) then
-!           write(*,*)'soil_settings: Problem while reading <soildepth>'
-!           write(*,*)trim(nf90_strerror(ierr))
-!           call abort
-!          endif
           call get_var("soildepth",mlayer,found)
           if (.not.found) then
@@ -184,9 +142,9 @@
       ! 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))
-	enddo
+         lay1=2.e-4
+         alpha=2
+         do iloop=0,nsoil-1
+           mlayer(iloop)=lay1*(alpha**(iloop-0.5))
+         enddo
       endif
 
@@ -214,36 +172,11 @@
 	volcapa=default_volcapa
       endif
-! Look for it
-!      ierr=NF_INQ_VARID(nid,"volcapa",nvarid)
-!      if (ierr.NE.nf90_noerr) then
-!         write(*,*)'soil_settings: Field <volcapa> not found!'
-!         write(*,*)'Initializing Volumetric heat capacity to ',
-!     &             default_volcapa
-!         volcapa=default_volcapa
-!      else
-!#ifdef NC_DOUBLE
-!       ierr = NF_GET_VAR_DOUBLE(nid,nvarid,volcapa)
-!#else
-!       ierr = NF_GET_VAR_REAL(nid,nvarid,volcapa)
-!#endif
-!        if (ierr.ne.nf90_noerr) then
-!	  write(*,*)'soil_settings: Problem while reading <volcapa>'
-!         call abort
-!	endif
-!      endif
-
-! 3. Thermal inertia for presetend climate -inertiedat- (note: it is declared in comsoil_h)
+
+! 3. Thermal inertia for present day climate -inertiedat- & the one given by the pem -ineritesoil-(note: it is declared in comsoil_h)
 ! ------------------
-
-! 3.1 Look (again) for thermal inertia data (to reset nvarid)
-!      ierr=nf90_inq_varid(nid,"inertiedat",nvarid)
-!      if (ierr.NE.nf90_noerr) then
-!         write(*,*)'soil_settings: Field <inertiedat> not found!'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         call abort
-!      endif
-
-! 3.2 Knowing the # of dimensions <inertidat> was defined as using,
+!     Knowing the # of dimensions <inertidat> was defined as using,
 !     read/build thermal inertia
+
+! 3.1 Present day - inertiedat
 
       if (ndims.eq.1) then ! "old 2D-surface" format
@@ -283,10 +216,5 @@
            endif
          endif ! of if (.not.allocated(oldinertiedat))
-!         ierr=nf90_get_var(nid,nvarid,oldinertiedat)
-!        if (ierr.NE.nf90_noerr) then
-!         write(*,*)'soil_settings: Problem while reading <inertiedat>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         call abort
-!        endif
+
         call get_field("inertiedat",oldinertiedat,found)
         if (.not.found) then
@@ -296,9 +224,4 @@
         endif
        else ! put values in therm_i
-!        ierr=nf90_get_var(nid,nvarid,inertiedat)
-!        if (ierr.NE.nf90_noerr) then
-!         write(*,*)'soil_settings: Problem while reading <inertiedat>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!         call abort
          call get_field("inertiedat",inertiedat,found)
          if (.not.found) then
@@ -307,15 +230,11 @@
      &        "failed loading <inertiedat>",1)
          endif
-!        endif
        endif ! of if (interpol)
       endif ! of if (ndims.eq.1)
 
-
-! 4. Read soil temperatures
-! -------------------------
-
-!      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
+! 3.2 Inertie from the PEM
+
       ok=inquire_field("inertiesoil")
-!      if (ierr.ne.nf90_noerr) then
+
       if (.not.ok) then
         write(*,*)'soil_settings: Field <inertiesoil> not found!'
@@ -326,5 +245,5 @@
         enddo
       else ! <inertiesoil> found
-       if (interpol) then ! put values in oldtsoil
+       if (interpol) then ! put values in oldinertiesoil
          if (.not.allocated(oldinertiesoil)) then
            allocate(oldinertiesoil(ngrid,dimlen,nslope),stat=ierr)
@@ -336,10 +255,5 @@
            endif
          endif
-!        ierr=nf90_get_var(nid,nvarid,oldtsoil)
-!        if (ierr.ne.nf90_noerr) then
-!	 write(*,*)'soil_settings: Problem while reading <tsoil>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!	 call abort
-!	endif
+
          call get_field("inertiesoil",oldinertiesoil,found)
          if (.not.found) then
@@ -348,18 +262,5 @@
      &          "failed loading <inertiesoil>",1)
          endif
-       else ! put values in tsoil
-!        corner(1)=1
-!        corner(2)=1
-!        corner(3)=indextime
-!        edges(1)=ngrid
-!        edges(2)=nsoil
-!        edges(3)=1
-!        !ierr=nf90_get_var(nid,nvarid,tsoil,corner,edges)
-!        ierr=nf90_get_var(nid,nvarid,tsoil)
-!        if (ierr.ne.nf90_noerr) then
-!	 write(*,*)'soil_settings: Problem while reading <tsoil>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!	 call abort
-!	endif
+       else ! put values in inertiesoil
          call get_field("inertiesoil",inertiesoil,found,
      &                 timeindex=indextime)
@@ -373,13 +274,10 @@
 
 
-
-
-      
-! 5. Read soil temperatures
+      
+! 4. Read soil temperatures
 ! -------------------------
 
-!      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
       ok=inquire_field("tsoil")
-!      if (ierr.ne.nf90_noerr) then
+
       if (.not.ok) then
         write(*,*)'soil_settings: Field <tsoil> not found!'
@@ -401,10 +299,4 @@
            endif
          endif
-!        ierr=nf90_get_var(nid,nvarid,oldtsoil)
-!        if (ierr.ne.nf90_noerr) then
-!	 write(*,*)'soil_settings: Problem while reading <tsoil>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!	 call abort
-!	endif
          call get_field("tsoil",oldtsoil,found)
          if (.not.found) then
@@ -414,17 +306,4 @@
          endif
        else ! put values in tsoil
-!        corner(1)=1
-!        corner(2)=1
-!        corner(3)=indextime
-!        edges(1)=ngrid
-!        edges(2)=nsoil
-!        edges(3)=1
-!        !ierr=nf90_get_var(nid,nvarid,tsoil,corner,edges)
-!        ierr=nf90_get_var(nid,nvarid,tsoil)
-!        if (ierr.ne.nf90_noerr) then
-!	 write(*,*)'soil_settings: Problem while reading <tsoil>'
-!         write(*,*)trim(nf90_strerror(ierr))
-!	 call abort
-!	endif
          call get_field("tsoil",tsoil,found,timeindex=indextime)
          if (.not.found) then
@@ -436,5 +315,6 @@
       endif! of if (.not.ok)
 
-! 6. If necessary, interpolate soil temperatures and thermal inertias
+
+! 5. If necessary, interpolate soil temperatures and thermal inertias
 ! -------------------------------------------------------------------
 
@@ -474,24 +354,22 @@
       endif
 
-      ! thermal inertia - present day
+      ! thermal inertia - present day (inertiedat)
       do ig=1,ngrid
 	  ! copy data
-	  oldval(1:dimlen)=oldinertiedat(ig,dimlen)
+	  oldval(1:dimlen)=oldinertiedat(ig,1:dimlen)
 	  ! interpolate
 	  call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil)
 	  !copy result in inertiedat
 	  inertiedat(ig,:)=newval(:)
-	  enddo
-      ! thermal inertia - general case
-      do ig=1,ngrid
-       do islope=1,nslope
+      ! thermal inertia - general case with PEM (inertie soil)
+          do islope=1,nslope
 	  ! copy data
-	     oldval(1:dimlen)=oldinertiesoil(ig,dimlen,islope)
+	     oldval(1:dimlen)=oldinertiesoil(ig,1:dimlen,islope)
 	  ! interpolate
 	     call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil)
 	  !copy result in inertiedat
-	  inertiesoil(ig,:,islope)=newval(:)
-          enddo
-	  enddo
+	     inertiesoil(ig,:,islope)=newval(:)
+           enddo
+	enddo ! ig
 
         
@@ -510,8 +388,8 @@
 	  tsoil(ig,:,islope)=newval(:)
           enddo
-	  enddo
+	 enddo
 	
 	!cleanup
-      deallocate(oldgrid)
+         deallocate(oldgrid)
 	  deallocate(oldval)
 	  deallocate(newval)
@@ -523,5 +401,5 @@
 
       
-! 7. Load Geothermal Flux
+! 6. Load Geothermal Flux
 ! ----------------------------------------------------------------------
 
@@ -536,5 +414,5 @@
 
       
-! 8. Report min and max values of soil temperatures and thermal inertias
+! 7. Report min and max values of soil temperatures and thermal inertias
 ! ----------------------------------------------------------------------
 
