Index: trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90	(revision 3532)
@@ -11,5 +11,5 @@
 !!!
 !!!
-!!! Author: EV, updated NS MSIM dynamical program for the PEM 
+!!! Author: EV, updated NS MSIM dynamical program for the PEM
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -37,5 +37,5 @@
   real(8)  ssi_depth_in, ssi_depth, T1
   real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG
-  real(8), dimension(NMAX,NP) :: porefill, porefill_in 
+  real(8), dimension(NMAX,NP) :: porefill, porefill_in
   real(8), dimension(nz) :: Tb
   real(8), dimension(NP) :: Tmean1, Tmean3, avrho1
@@ -68,5 +68,5 @@
      !read(21,*) latitude(k),albedo(k),thIn(k),htopo(k)
      ! empirical relation from Mellon & Jakosky
-     rhoc(k) = 800.*(150.+100.*sqrt(34.2+0.714*thIn(k))) 
+     rhoc(k) = 800.*(150.+100.*sqrt(34.2+0.714*thIn(k)))
   enddo
   !close(21)
@@ -80,5 +80,5 @@
     z(iloop + 1) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
   enddo
-  
+
 
   !open(unit=30,file='z.'//ext,action='write',status='unknown')
@@ -95,6 +95,6 @@
   timestep = 1 ! must be integer fraction of 1 ka
   icetime = -tmax-timestep  ! earth years
-  
-  ! initializations 
+
+  ! initializations
   !Tb = -9999.
   zdepthF(:) = -9999.
@@ -152,5 +152,5 @@
   !print *,'Zt0=  ',ZdepthT
      call icelayer_mars(timestep,nz,NP,thIn,rhoc,z,porosity,pfrost,Tb,zdepthF, &
-          & zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, & 
+          & zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, &
           & albedo,p0,icefrac,zdepthT,avrho1, &
           & avrho1prescribed)
@@ -164,5 +164,5 @@
          !write(36,*) icetime,latitude(k),zdepthF(k),porefill(1:nz,k)
            ! compact output format
- !          write(36,'(f10.0,2x,f7.3,1x,f11.5,1x)',advance='no') & 
+ !          write(36,'(f10.0,2x,f7.3,1x,f11.5,1x)',advance='no') &
  !               & icetime,latitude(k),zdepthF(k)
  !          call compactoutput(36,porefill(:,k),nz)
@@ -179,5 +179,5 @@
 
 !501 format (f10.0,2x,f7.3,2x,f10.4,2(2x,f6.2),2x,f9.3,2x,g11.4)
-  
+
 end subroutine dyn_ss_ice_m
 
@@ -232,18 +232,18 @@
      endif
      newti = sqrt(newrhoc*k)
-     
+
   case (2)  ! massive ice (pure or dirty ice)
      newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice
      k = icefrac*kice + (1.-icefrac)*kw
      newti = sqrt(newrhoc*k)
-  
+
   case (3)  ! all ice, special case of layertype 2, which doesn't use porosity
      newrhoc = icedensity*cice
-     k = kice 
-     newti = sqrt(newrhoc*k)
-  
+     k = kice
+     newti = sqrt(newrhoc*k)
+
   case (4)  ! pores completely filled with ice, special case of layertype 1
      newrhoc = rhocobs + porosity*icedensity*cice
-     k = porosity*kice + kobs ! option A, end-member case of type 1, option A 
+     k = porosity*kice + kobs ! option A, end-member case of type 1, option A
      !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average
      newti = sqrt(newrhoc*k)
@@ -257,12 +257,12 @@
   case default
      error stop 'invalid layer type'
-     
+
   end select
-  
+
 end subroutine soilthprop
 
 
 !=======================================================================
-      
+
       real*8 function frostpoint(p)
 !     inverse of psv
@@ -271,5 +271,5 @@
       implicit none
       real*8 p
-      
+
 !-----inverse of parametrization 1
 !      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
@@ -278,5 +278,5 @@
 !      parameter (R=8.314,pt=6.11e2,Tt=273.16)
 !      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
-      
+
 !-----inverse of parametrization 2
 !     inverse of eq. (2) in Murphy & Koop (2005)
@@ -288,5 +288,5 @@
 !     eq. (8) in Murphy & Koop (2005)
 !      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
-      
+
       end function frostpoint
 
Index: trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90	(revision 3532)
@@ -11,5 +11,5 @@
   real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
   integer, parameter :: EQUILTIME =15 ! [Mars years]
-  
+
   integer, parameter :: NMAX = 1000
 
@@ -59,5 +59,5 @@
      B = Diff*bigstep*86400.*365.24/(porosity*927.)
      !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'h2o'))
-     
+
      typeT = -9
      if (zdepthT(k)>=0. .and. zdepthT(k)<z(nz)) then
@@ -72,7 +72,7 @@
    !  call assignthermalproperties(nz,thIn(k),rhoc(k),ti,rhocv, &
    !       & typeT,icefrac,porosity,porefill(:,k))
-     
+
      !----run thermal model
-     call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, & 
+     call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, &
           &     ti, rhocv, fracIR, fracDust, p0(k), &
           &     avdrho(k), avdrhoP(k), avrho1(k), Tb(k), zdepthE(k), typeF, &
@@ -90,5 +90,5 @@
         jump = -9
      endif
-     if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then 
+     if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then
         write(34,*) '# zdepth arrested'
         if (jump>1) write(34,*) '# previous step was too large',jump
@@ -99,5 +99,5 @@
      avdrho_old(k) = avdrho(k)
 
-!----mode 2 growth  
+!----mode 2 growth
      typeP = -9;  mode2 = .FALSE.
      do j=1,nz
@@ -111,5 +111,5 @@
              & zdepthE(k)<z(typeP) .and. &
              & z(typeP)-zdepthE(k)>2*(z(typeP)-z(typeP-1))) then  ! trick that avoids oscillations
-           deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B  ! conservation of mass 
+           deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B  ! conservation of mass
            if (deltaz>z(typeP)-z(typeP-1)) then  ! also implies avdrhoP<0.
               mode2 = .TRUE.
@@ -143,5 +143,5 @@
 !
 !  Tb<0 initializes temperatures
-!  Tb>0 initializes temperatures with Tb 
+!  Tb>0 initializes temperatures with Tb
 !***********************************************************************
   use fast_subs_univ, only: depths_avmeth, equildepth
@@ -158,5 +158,5 @@
   real(8), intent(INOUT) :: Tb(nz)
   integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
-  real(8), intent(OUT) :: zdepthE, zdepthF 
+  real(8), intent(OUT) :: zdepthE, zdepthF
   real(8), intent(OUT) :: ypp(nz) ! (d rho/dt)/Diff
   real(8), intent(OUT) :: Tmean3, zdepthG
@@ -171,9 +171,9 @@
   real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2
   real(8) rhosatav0, rhosatav(nz), rlow
-  
+
   tmax = EQUILTIME*sols_per_my
   nsteps = int(tmax/dt)     ! calculate total number of timesteps
 
-  Tco2frost = tfrostco2(patm) 
+  Tco2frost = tfrostco2(patm)
 
   !if (Tb<=0.) then  ! initialize
@@ -182,10 +182,10 @@
      !Tmean0 = Tmean0-5.
      !write(34,*) '# initialized with temperature estimate at',latitude/d2r,'of',Tmean0,'K'
-     !T(1:nz) = Tmean0 
+     !T(1:nz) = Tmean0
   !else
      T(1:nz) = Tb
      ! not so good when Fgeotherm is on
   !endif
-  
+
   albedo = albedo0
   emiss = emiss0
@@ -207,7 +207,7 @@
 !  Qn=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
   !Qn = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
-  !----loop over time steps 
+  !----loop over time steps
   do n=0,nsteps-1
-     time = (n+1)*dt         !   time at n+1 
+     time = (n+1)*dt         !   time at n+1
    !  tdays = time*(sec_per_sol/earthDay) ! parenthesis may improve roundoff
    !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
@@ -215,5 +215,5 @@
 !     Qnp1=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
      !Qnp1 = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
-     
+
      Tsurfold = Tsurf
      Fsurfold = Fsurf
@@ -226,5 +226,5 @@
         T(1:nz) = Told(1:nz)
         !call conductionT(nz,z,dt*sec_per_sol,T,Tsurfold,Tco2frost,ti, &
-             !&              rhocv,Fgeotherm,Fsurf) 
+             !&              rhocv,Fgeotherm,Fsurf)
         Tsurf = Tco2frost
     !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
@@ -240,5 +240,5 @@
      endif
      !Qn=Qnp1
-     
+
      if (time>=tmax-sols_per_my) then
         Tmean1 = Tmean1 + Tsurf
@@ -253,5 +253,5 @@
 
   enddo  ! end of time loop
-  
+
   avrho1 = avrho1/nm
   if (present(avrho1prescribed)) then
@@ -266,5 +266,5 @@
   endif
   avdrho = avrho2-avrho1
-  typeP = -9 
+  typeP = -9
   do i=1,nz
      if (porefill(i)>0.) then
@@ -275,5 +275,5 @@
   if (typeP>0) then
      avdrhoP = rhosatav(typeP) - avrho1
-  else  
+  else
      avdrhoP = -9999.
   end if
@@ -282,5 +282,5 @@
 
   if (Fgeotherm>0.) then
-     Tb = Tmean1 
+     Tb = Tmean1
      typeG = 1   ! will be overwritten by depths_avmeth
      rlow = 2*rhosatav(nz)-rhosatav(nz-1)
@@ -328,8 +328,8 @@
       psv = exp(A/T+B)  ! Clapeyron
 
-!-----parametrization 3      
+!-----parametrization 3
 !     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
 !     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
-      
+
       end function psv
 
Index: trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90	(revision 3532)
@@ -32,7 +32,7 @@
   real(8) equildepth  ! =zdepthE
   !real(8), external :: zint  ! defined in allinterfaces.mod
-  
+
   typeE = -9; equildepth = -9999.
-  do i=1,nz 
+  do i=1,nz
      if (rhosatav(i) <= avrho1) then
         typeE=i
@@ -67,5 +67,5 @@
   integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
   real(8), intent(INOUT) :: zdepthF
-  real(8), intent(IN) :: B 
+  real(8), intent(IN) :: B
   real(8), intent(OUT) :: ypp(nz), zdepthG
   integer, intent(INOUT) :: typeG  ! positive on input when Fgeotherm>0
@@ -112,5 +112,5 @@
 
 !-depth to shallowest perennial ice
-  typeP = -9 
+  typeP = -9
   do i=1,nz
      if (porefill(i)>0.) then
@@ -157,5 +157,5 @@
   endif
   if (typeG>0 .and. typeT<0) then
-     call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove) 
+     call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove)
      newtypeG = -9
      do i=typeG,nz
@@ -178,5 +178,5 @@
      if (newtypeG>0) typeG=newtypeG
   end if
-  ! if typeG>0, then all ice at and below typeG should be erased 
+  ! if typeG>0, then all ice at and below typeG should be erased
 end subroutine depths_avmeth
 
@@ -204,5 +204,5 @@
 
   ! advance ice table, avdrho>0 is retreat
-  if (zdepthT>=0. .and. avdrho>0.) then 
+  if (zdepthT>=0. .and. avdrho>0.) then
      typeP=-9999; typeT=-9999
      do j=1,nz
@@ -226,5 +226,5 @@
   endif
   if (zdepthT>z(nz)) zdepthT=-9999.
-  
+
   ! advance interface, avdrhoP>0 is loss from zdepthP
   if (avdrhoP>0.) then
@@ -232,5 +232,5 @@
      do j=1,nz
         if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF
-        if (zdepthT>=0. .and. z(j)>zdepthT) exit 
+        if (zdepthT>=0. .and. z(j)>zdepthT) exit
         call colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j,integ)
         erase = j
@@ -241,5 +241,5 @@
 
   ! new depth
-  newtypeP = -9 
+  newtypeP = -9
   do j=1,nz
      if (zdepthT>=0. .and. z(j)>zdepthT) exit
@@ -253,5 +253,5 @@
   ub = typeF
   if (newtypeP>0 .and. typeF>0 .and. newtypeP<ub) ub=newtypeP
-  if (ub>0) then  
+  if (ub>0) then
      do j=ub,nz
         porefill(j) = porefill(j) + B*ypp(j)
Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3532)
@@ -501,2 +501,5 @@
 - Cleaning of "glaciers_mod.F90";
 - Optimization for the computation of ice density according to temperature by using a function.
+
+== 04/12/2024 == JBC
+Removing unecessary module/subroutine "interpol_TI_PEM2PCM.F90" + Few small corrections/cleanings.
Index: trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90	(revision 3532)
@@ -4,7 +4,7 @@
 !-----------------------------------------------------------------------
 !  Author: LL
-!  Purpose: This module gathers the different routines used in the PEM to compute the soil temperature evolution and initialisation. 
-!  
-!  Note: depths of layers and mid-layers, soil thermal inertia and 
+!  Purpose: This module gathers the different routines used in the PEM to compute the soil temperature evolution and initialisation.
+!
+!  Note: depths of layers and mid-layers, soil thermal inertia and
 !        heat capacity are commons in comsoil_PEM.h
 !-----------------------------------------------------------------------
@@ -12,6 +12,4 @@
 !=======================================================================
 
-
-
 SUBROUTINE compute_tsoil_pem(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
 
@@ -24,6 +22,6 @@
 !  Author: LL
 !  Purpose: Compute soil temperature using an implict 1st order scheme
-!  
-!  Note: depths of layers and mid-layers, soil thermal inertia and 
+!
+!  Note: depths of layers and mid-layers, soil thermal inertia and
 !        heat capacity are commons in comsoil_PEM.h
 !-----------------------------------------------------------------------
@@ -31,19 +29,18 @@
 #include "dimensions.h"
 
-!-----------------------------------------------------------------------
-!  arguments
-!  ---------
-!  inputs:
-integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points 
-integer,                      intent(in) :: nsoil     ! number of soil layers  
-logical,                      intent(in) :: firstcall ! identifier for initialization call 
+! Inputs:
+! -------
+integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points
+integer,                      intent(in) :: nsoil     ! number of soil layers
+logical,                      intent(in) :: firstcall ! identifier for initialization call
 real, dimension(ngrid,nsoil), intent(in) :: therm_i   ! thermal inertia [SI]
 real,                         intent(in) :: timestep  ! time step [s]
 real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
- 
-! outputs:
+! Outputs:
+!---------
 real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
-! local variables:
-integer :: ig, ik    
+! Local:
+!-------
+integer :: ig, ik
 
 ! 0. Initialisations and preprocessing step
@@ -52,5 +49,5 @@
     do ig = 1,ngrid
         do ik = 0,nsoil - 1
-            mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa    
+            mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
         enddo
     enddo
@@ -64,5 +61,5 @@
     enddo
 
-! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
+! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEM
     ! mu_PEM
     mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
@@ -99,5 +96,5 @@
 ! Other layers:
         do ik = 1,nsoil - 1
-                tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik) 
+                tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
         enddo
     enddo
@@ -139,16 +136,15 @@
 #include "dimensions.h"
 
-!-----------------------------------------------------------------------
-!  arguments
-!  ---------
-!  inputs:
+! Inputs:
+!--------
 integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
 integer,                      intent(in) :: nsoil   ! number of soil layers
 real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
 real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]
-
-! outputs:
+! Outputs:
+!---------
 real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
-! local variables:
+! Local:
+!-------
 integer :: ig, ik, iloop
 
@@ -169,5 +165,5 @@
 enddo
 
-! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
+! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEM
 ! mu_PEM
 mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
@@ -179,5 +175,5 @@
 do ig = 1,ngrid
     ! d_k
-    do ik = 1,nsoil-1
+    do ik = 1,nsoil - 1
         coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
     enddo
@@ -207,14 +203,13 @@
 !  2. Compute soil temperatures
 do iloop = 1,10 !just convergence
-    ! First layer:
-    do ig = 1,ngrid
-        tsoil(ig,1)=(tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
-                   (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
-    ! Other layers:
-    do ik = 1,nsoil - 1
-        tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
-    enddo
-enddo
-
+    do ig = 1,ngrid
+        ! First layer:
+        tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
+                      (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
+        ! Other layers:
+        do ik = 1,nsoil - 1
+            tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
+        enddo
+    enddo
 enddo ! iloop
 
Index: trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 3532)
@@ -3,26 +3,26 @@
 implicit none
 
-integer, parameter                        :: nsoilmx_PEM = 68   ! number of layers in the PEM
-real, allocatable, dimension(:),     save :: layer_PEM          ! soil layer depths [m]
-real, allocatable, dimension(:),     save :: mlayer_PEM         ! soil mid-layer depths [m]
-real, allocatable, dimension(:,:,:), save :: TI_PEM             ! soil thermal inertia [SI]
-real, allocatable, dimension(:,:),   save :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
+integer, parameter                  :: nsoilmx_PEM = 68   ! number of layers in the PEM
+real, allocatable, dimension(:)     :: layer_PEM          ! soil layer depths [m]
+real, allocatable, dimension(:)     :: mlayer_PEM         ! soil mid-layer depths [m]
+real, allocatable, dimension(:,:,:) :: TI_PEM             ! soil thermal inertia [SI]
+real, allocatable, dimension(:,:)   :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
 ! Variables (FC: built in firstcall in soil.F)
-real, allocatable, dimension(:,:,:), save :: tsoil_PEM          ! sub-surface temperatures [K]
-real, allocatable, dimension(:,:),   save :: mthermdiff_PEM     ! (FC) mid-layer thermal diffusivity [SI]
-real, allocatable, dimension(:,:),   save :: thermdiff_PEM      ! (FC) inter-layer thermal diffusivity [SI]
-real, allocatable, dimension(:),     save :: coefq_PEM          ! (FC) q_{k+1/2} coefficients [SI]
-real, allocatable, dimension(:,:),   save :: coefd_PEM          ! (FC) d_k coefficients [SI]
-real, allocatable, dimension(:,:),   save :: alph_PEM           ! (FC) alpha_k coefficients [SI]
-real, allocatable, dimension(:,:),   save :: beta_PEM           ! beta_k coefficients [SI]
-real,                                save :: mu_PEM             ! mu coefficient [SI]
-real,                                save :: fluxgeo            ! Geothermal flux [W/m^2]
-real,                                save :: depth_breccia      ! Depth at which we have breccia [m]
-real,                                save :: depth_bedrock      ! Depth at which we have bedrock [m]
-integer,                             save :: index_breccia      ! last index of the depth grid before having breccia
-integer,                             save :: index_bedrock      ! last index of the depth grid before having bedrock
-logical                                   :: soil_pem           ! True by default, to run with the subsurface physic. Read in pem.def
-real,                                save :: ini_huge_h2oice    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
-logical,                             save :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
+real, allocatable, dimension(:,:,:) :: tsoil_PEM          ! sub-surface temperatures [K]
+real, allocatable, dimension(:,:)   :: mthermdiff_PEM     ! (FC) mid-layer thermal diffusivity [SI]
+real, allocatable, dimension(:,:)   :: thermdiff_PEM      ! (FC) inter-layer thermal diffusivity [SI]
+real, allocatable, dimension(:)     :: coefq_PEM          ! (FC) q_{k+1/2} coefficients [SI]
+real, allocatable, dimension(:,:)   :: coefd_PEM          ! (FC) d_k coefficients [SI]
+real, allocatable, dimension(:,:)   :: alph_PEM           ! (FC) alpha_k coefficients [SI]
+real, allocatable, dimension(:,:)   :: beta_PEM           ! beta_k coefficients [SI]
+real                                :: mu_PEM             ! mu coefficient [SI]
+real                                :: fluxgeo            ! Geothermal flux [W/m^2]
+real                                :: depth_breccia      ! Depth at which we have breccia [m]
+real                                :: depth_bedrock      ! Depth at which we have bedrock [m]
+integer                             :: index_breccia      ! last index of the depth grid before having breccia
+integer                             :: index_bedrock      ! last index of the depth grid before having bedrock
+logical                             :: soil_pem           ! True by default, to run with the subsurface physic. Read in pem.def
+real                                :: ini_huge_h2oice    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
+logical                             :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
 
 !=======================================================================
Index: trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3532)
@@ -24,6 +24,6 @@
 !!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do
 !!!          the ice transfer
-!!!          
-!!!          
+!!!
+!!!
 !!! Author: LL
 !!!
@@ -47,8 +47,8 @@
 real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
 real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
-! Local 
+! Local
 !------
 real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
-real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
+real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
 
 write(*,*) "Flow of CO2 glaciers"
@@ -67,6 +67,6 @@
 !!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do
 !!!          the ice transfer
-!!!          
-!!!          
+!!!
+!!!
 !!! Author: LL
 !!!
@@ -87,6 +87,6 @@
 real, dimension(ngrid,nslope), intent(inout) :: flag_h2oflow      ! flag to see if there is flow on the subgrid slopes
 real, dimension(ngrid),        intent(inout) :: flag_h2oflow_mesh ! same but within the mesh
-! Local 
-real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2  glacier before flow 
+! Local
+real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2  glacier before flow
 
 write(*,*) "Flow of H2O glaciers"
@@ -111,5 +111,5 @@
 !!!
 !!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow
-!!!          
+!!!
 !!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
 !!!
@@ -222,5 +222,5 @@
             endif ! co2ice > hmax
         endif ! iflat
-    enddo !islope 
+    enddo !islope
 enddo !ig
 
@@ -237,8 +237,8 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  
+
 use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2
 
-implicit none 
+implicit none
 
 ! arguments:
@@ -253,5 +253,5 @@
 
 ! OUTPUT
-real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged 
+real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged
 
 ! LOCAL
Index: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3532)
@@ -401,5 +401,5 @@
 ! Locals
 ! ------
-integer :: ik       ! Loop variables 
+integer :: ik       ! Loop variables
 integer :: indexice ! Index of the ice
 
@@ -412,5 +412,5 @@
     if(ice_depth < mlayer_PEM(0)) then
         indexice = 0.
-    else     
+    else
         do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
             if (mlayer_PEM(ik) <= ice_depth .and. mlayer_PEM(ik + 1) > ice_depth) then
Index: trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90	(revision 3532)
@@ -31,5 +31,5 @@
 logical       :: ok
 integer       :: cstat
-character(10) :: frmt
+character(50) :: info_frmt
 character(20) :: ich1, ich2, ich3, ich4, fch1, fch2, fch3
 
@@ -53,5 +53,6 @@
     ! Martian date, Number of Martians years done by the PEM run, Number of Martians years done by the chainded simulation, Code of the stopping criterion
     ! The conversion ratio from Planetary years to Earth years is given in the header of the file
-    write(1,*) year_bp_ini + i_myear, i_myear_leg, i_myear, stopPEM
+    info_frmt = '(f'//int2str(nb_digits(year_bp_ini + i_myear) + 5)//'.4,'//'f'//int2str(nb_digits(i_myear_leg) + 5)//'.4,'//'f'//int2str(nb_digits(i_myear) + 5)//'.4,'//'i0)'
+    write(1,trim(adjustl(info_frmt))) year_bp_ini + i_myear, i_myear_leg, i_myear, stopPEM
     close(1)
 else
Index: trunk/LMDZ.COMMON/libf/evolution/interpol_TI_PEM2PCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/interpol_TI_PEM2PCM_mod.F90	(revision 3527)
+++ 	(revision )
@@ -1,39 +1,0 @@
-MODULE interpol_TI_PEM2PCM_mod
-
-implicit none
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE interpol_TI_PEM2PCM(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PEM,TI_PCM)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!
-!!! Purpose: Transfer the thermal inertia from the PEM vertical grid to the PCM vertical grid
-!!!
-!!!
-!!! Author: LL
-!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-implicit none
-
-!======================================================================
-!  arguments
-!  ---------
-!  inputs:
-integer,                                 intent(in) :: ngrid     ! # of horizontal grid points
-integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
-integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
-integer,                                 intent(in) :: nsoil_PCM ! # of soil layers in the GCM
-real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM    ! Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
-
-real, dimension(ngrid,nsoil_PCM,nslope), intent(inout) :: TI_PCM ! Thermal inertia in the PCM vertical grid [J/m^2/K/s^{1/2}]
-
-!----- Code
-TI_PCM = TI_PEM(:,:nsoil_PCM,:)
-
-END SUBROUTINE interpol_TI_PEM2PCM
-
-END MODULE interpol_TI_PEM2PCM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/math_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/math_mod.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/math_mod.F90	(revision 3532)
@@ -3,5 +3,5 @@
 !!!
 !!! 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
 !!!
@@ -17,6 +17,6 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
-!!! Purpose: Compute the first derivative of a function y(z) on an irregular grid 
-!!!           
+!!! Purpose: Compute the first derivative of a function y(z) on an irregular grid
+!!!
 !!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
 !!!
@@ -59,6 +59,6 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
-!!! Purpose: Compute the second derivative of a function y(z) on an irregular grid 
-!!!           
+!!! Purpose: Compute the second derivative of a function y(z) on an irregular grid
+!!!
 !!! Author: N.S (raw copy/paste from MSIM)
 !!!
@@ -97,5 +97,5 @@
 !!!
 !!! Purpose: First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
-!!!           
+!!!
 !!! Author: N.S (raw copy/paste from MSIM)
 !!!
@@ -128,11 +128,11 @@
 !!!
 !!! Purpose:  Column integrates y on irregular grid
-!!!           
+!!!
 !!! Author: N.S (raw copy/paste from MSIM)
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  
+
 implicit none
-  
+
 integer, intent(in) :: nz, i1, i2
 real,    intent(in) :: y(nz), z(nz)
@@ -156,5 +156,5 @@
 !!!
 !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
-!!! 
+!!!
 !!! Author: LL
 !!!
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3532)
@@ -62,5 +62,4 @@
 use compute_tend_mod,           only: compute_tend
 use info_PEM_mod,               only: info_PEM
-use interpol_TI_PEM2PCM_mod,    only: interpol_TI_PEM2PCM
 use nb_time_step_PCM_mod,       only: nb_time_step_PCM
 use pemetat0_mod,               only: pemetat0
@@ -298,5 +297,5 @@
     call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
     if (cstat /= 0) then
-        call execute_command_line('qstat -f '//trim(job_id)//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt', cmdstat = cstat)
+        call execute_command_line('qstat -f '//trim(job_id)//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt',cmdstat = cstat)
         if (cstat > 0) then
             error stop 'pem: command execution failed!'
@@ -939,5 +938,5 @@
             do ig = 1,ngrid
                 do islope = 1,nslope
-                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2),ice_porefilling(ig,:,islope),porefill,ssi_depth)
+                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth)
                     icetable_depth(ig,islope) = ssi_depth
                     ice_porefilling(ig,:,islope) = porefill
@@ -961,5 +960,5 @@
             do ig = 1,ngrid
                 do islope = 1,nslope
-                    do l = 1,nsoilmx_PEM 
+                    do l = 1,nsoilmx_PEM
                         if (l == 1) then
                             totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
@@ -1008,5 +1007,5 @@
                 call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
                 call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
-            endif                       
+            endif
         endif
     enddo
@@ -1101,5 +1100,5 @@
 ! III_a.2 Tsoil update (for startfi)
 if (soil_pem) then
-    call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
+    inertiesoil = TI_PEM(:,:nsoilmx,:)
     tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen)
 #ifndef CPP_STD
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3532)
@@ -45,5 +45,5 @@
 real,                           intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
 real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1       ! surface temperature at the first year of PCM call [K]
-real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second  year of PCM call [K]
+real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second year of PCM call [K]
 real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
 real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
@@ -281,6 +281,5 @@
                 call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
             else
-! predictor corrector: restart from year 1 of the PCM and build the evolution of
-! tsoil at depth
+! predictor corrector: restart from year 1 of the PCM and build the evolution of tsoil at depth
                 tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
                 call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
@@ -289,5 +288,5 @@
                 call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
 
-                do iloop = nsoil_PCM+1,nsoil_PEM
+                do iloop = nsoil_PCM + 1,nsoil_PEM
                     tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
                 enddo
@@ -295,12 +294,8 @@
 
             do it = 1,timelen
-                do isoil = nsoil_PCM+1,nsoil_PEM
-                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
-                enddo
-            enddo
-            do isoil = nsoil_PCM+1,nsoil_PEM
-                do ig = 1,ngrid
-                    watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
-                enddo
+                tsoil_inst(:,nsoil_PCM + 1:nsoil_PEM,islope,it) = tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)
+            enddo
+            do isoil = nsoil_PCM + 1,nsoil_PEM
+                watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
             enddo
         enddo ! islope
@@ -314,5 +309,5 @@
                 write(*,*)'PEM settings: failed loading <ice_table_depth>'
                 write(*,*)'will reconstruct the values of the ice table given the current state'
-                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
+                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
                 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
                 do islope = 1,nslope
@@ -322,4 +317,9 @@
             write(*,*) 'PEMETAT0: ICE TABLE done'
         else if (icetable_dynamic) then
+            call get_field("ice_porefilling",ice_porefilling,found)
+            if (.not. found) then
+                write(*,*)'PEM settings: failed loading <ice_porefilling>'
+                ice_porefilling = 0.
+            endif
             call get_field("ice_table_depth",ice_table_depth,found)
             if (.not. found) then
@@ -331,9 +331,4 @@
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
                 enddo
-            endif
-            call get_field("ice_porefilling",ice_porefilling,found)
-            if (.not. found) then
-                write(*,*)'PEM settings: failed loading <ice_porefilling>'
-                ice_porefilling = 1.
             endif
             write(*,*) 'PEMETAT0: ICE TABLE done'
@@ -510,5 +505,5 @@
             write(*,*) 'PEMETAT0: Ice table done'
         else if (icetable_dynamic) then
-            ice_porefilling = 1.
+            ice_porefilling = 0.
             ice_table_depth = -9999.
             call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
@@ -518,5 +513,5 @@
             write(*,*) 'PEMETAT0: Ice table done'
         endif
-        
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !d) Regolith adsorbed
Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3532)
@@ -41,5 +41,5 @@
 
 write(*,*) "Update of the CO2 tendency from the current pressure"
-    
+
 ! Evolution of the water ice for each physical point
 do i = 1,ngrid
Index: trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90	(revision 3532)
@@ -15,5 +15,5 @@
 
 !=======================================================================
-!  Author: LL, based on work by  Ehouarn Millour (07/2006)
+!  Author: LL, based on work by Ehouarn Millour (07/2006)
 !
 !  Purpose: Read and/or initialise soil depths and properties
@@ -53,9 +53,9 @@
 alpha = 2
 do iloop = 0,nsoil_PCM - 1
-    mlayer_PEM(iloop) = lay1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
+    mlayer_PEM(iloop) = lay1*(1 + iloop**2.9*(1 - exp(-real(iloop)/20.)))
 enddo
 
-do iloop = 0, nsoil_PEM-1
-    if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM-1)) then
+do iloop = 0,nsoil_PEM - 1
+    if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then
         index_powerlaw = iloop
         exit
@@ -63,13 +63,13 @@
 enddo
 
-do iloop = nsoil_PCM, nsoil_PEM-1 
-    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop-nsoil_PCM)-0.5))
+do iloop = nsoil_PCM,nsoil_PEM - 1
+    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5))
 enddo
 
 ! Build layer_PEM()
-do iloop=1,nsoil_PEM-1
-layer_PEM(iloop)=(mlayer_PEM(iloop)+mlayer_PEM(iloop-1))/2.
+do iloop = 1,nsoil_PEM - 1
+layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2.
 enddo
-layer_PEM(nsoil_PEM)=2*mlayer_PEM(nsoil_PEM-1) - mlayer_PEM(nsoil_PEM-2)
+layer_PEM(nsoil_PEM) = 2*mlayer_PEM(nsoil_PEM - 1) - mlayer_PEM(nsoil_PEM - 2)
 
 
Index: trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90	(revision 3527)
+++ trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90	(revision 3532)
@@ -680,8 +680,8 @@
     call abort_physic("writediagsoilpem","firstname too short",1)
   endif
-  
+
   ! Set output sample rate
   isample = int(ecritpem) ! same as for diagpem outputs
-  
+
   ! Create output NetCDF file
   if (is_master) then
@@ -724,5 +724,5 @@
     enddo
    endif
-   
+
    ! write "header" of file (longitudes, latitudes, geopotential, ...)
    if (klon_glo>1) then ! general 3D case
@@ -733,8 +733,8 @@
 
   endif ! of if (is_master)
-  
+
   ! set zitau to -1 to be compatible with zitau incrementation step below
   zitau=-1
-  
+
 else
   ! If not an initialization call, simply open the NetCDF file
@@ -757,5 +757,5 @@
     ntime=ntime+1
     date = float(zitau + 1)
-    
+
     if (is_master) then
      ! Get NetCDF ID for "time"
@@ -769,5 +769,5 @@
      if (ierr.ne.NF_NOERR) then
       write(*,*)"writediagsoilpem: Failed writing date to time variable"
-      call abort_physic("writediagsoilpem","failed writing time",1) 
+      call abort_physic("writediagsoilpem","failed writing time",1)
      endif
     endif ! of if (is_master)
@@ -810,5 +810,5 @@
    data3_1d(1,1:nsoilmx_PEM)=px(1,1:nsoilmx_PEM)
   endif
-  
+
   ! B. Write (append) the variable to the NetCDF file
   if (is_master) then
@@ -828,5 +828,5 @@
     call def_var(nid,name,title,units,4,id,varid,ierr)
   endif ! of if (ierr.ne.NF_NOERR)
-  
+
   ! B.2. Prepare things to be able to write/append the variable
   corners(1)=1
@@ -834,5 +834,5 @@
   corners(3)=1
   corners(4)=ntime
-  
+
   if (klon_glo==1) then
     edges(1)=1
@@ -843,5 +843,5 @@
   edges(3)=nsoilmx_PEM
   edges(4)=1
-  
+
   ! B.3. Write the slab of data
 !#ifdef NC_DOUBLE
@@ -917,5 +917,5 @@
   corners(2)=1
   corners(3)=ntime
-  
+
   if (klon_glo==1) then
     edges(1)=1
@@ -925,5 +925,5 @@
   edges(2)=nbp_lat
   edges(3)=1
-  
+
   ! B.3. Write the slab of data
 !#ifdef NC_DOUBLE
@@ -966,5 +966,5 @@
   ! B.2. Prepare things to be able to write/append the variable
   corners(1)=ntime
-  
+
   edges(1)=1
 
