Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3672)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3673)
@@ -616,2 +616,7 @@
 == 03/03/2025 == JBC
 Adjusting the anticipation time for more security regarding the job time limit detection.
+
+== 07/03/2025 == JBC
+- Correction for the update of average pressure at the end.
+- Addition of pore-filled ice information for 'stratum' in the layering algorithm.
+- Deactivating the soil temperature shifting because of bad values given by interpolation.
Index: trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py	(revision 3672)
+++ trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py	(revision 3673)
@@ -67,5 +67,5 @@
     stratif_h2oice = np.zeros((ngrid,nfile,nslope,max_nb_str))
     stratif_dust = np.zeros((ngrid,nfile,nslope,max_nb_str))
-    stratif_air = np.zeros((ngrid,nfile,nslope,max_nb_str))
+    stratif_pore = np.zeros((ngrid,nfile,nslope,max_nb_str))
     for var_name in datasets[0].variables:
         if 'top_elevation' in var_name:
@@ -97,10 +97,10 @@
                             stratif_dust[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
             print(f"Processed variable: {var_name}")
-        elif 'air_volfrac' in var_name:
+        elif 'pore_volfrac' in var_name:
             for i in range(0,ngrid):
                 for j in range(0,nfile):
                     for k in range(0,nslope):
                         if f'slope{k + 1:02d}' in var_name:
-                            stratif_air[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
+                            stratif_pore[i,j,k,:datasets[j].variables[var_name].shape[1]] = datasets[j].variables[var_name][0,:,i]
             print(f"Processed variable: {var_name}")
 
@@ -109,5 +109,5 @@
         ds.close()
 
-    stratif_data = [stratif_heights,stratif_co2ice,stratif_h2oice,stratif_dust,stratif_air]
+    stratif_data = [stratif_heights,stratif_co2ice,stratif_h2oice,stratif_dust,stratif_pore]
     return stratif_data, max_top_elevation, longitude, latitude
 
@@ -158,5 +158,5 @@
 
 ### Figures plotting
-subtitle = ['CO2 ice','H2O ice','Dust','Air']
+subtitle = ['CO2 ice','H2O ice','Dust','Pore']
 cmap = plt.get_cmap('viridis').copy()
 cmap.set_under('white')
Index: trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py	(revision 3672)
+++ trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py	(revision 3673)
@@ -43,5 +43,5 @@
 stratif_h2oice_volfrac = []
 stratif_dust_volfrac = []
-stratif_air_volfrac = []
+stratif_pore_volfrac = []
 for i in range(1,nslope + 1):
     stratif_thickness.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_thickness'][:])
@@ -50,10 +50,10 @@
     stratif_h2oice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_h2oice_volfrac'][:])
     stratif_dust_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_dust_volfrac'][:])
-    stratif_air_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_air_volfrac'][:])
+    stratif_pore_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_pore_volfrac'][:])
 
 ### Display the data
 igrid = igrid - 1
 islope = islope - 1
-labels = 'CO2 ice', 'H2O ice', 'Dust', 'Air'
+labels = 'CO2 ice', 'H2O ice', 'Dust', 'Pore'
 contents = np.zeros([4,len(stratif_top_elevation[islope][0,:,igrid]) + 1])
 height = np.zeros(len(stratif_top_elevation[islope][0,:,igrid]) + 1)
@@ -64,5 +64,5 @@
     contents[1,1 + i] = stratif_h2oice_volfrac[islope][0,i,igrid]
     contents[2,1 + i] = stratif_dust_volfrac[islope][0,i,igrid]
-    contents[3,1 + i] = stratif_air_volfrac[islope][0,i,igrid]
+    contents[3,1 + i] = stratif_pore_volfrac[islope][0,i,igrid]
 contents[:,0] = contents[:,1]
 
Index: trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3672)
+++ trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3673)
@@ -33,12 +33,13 @@
 ! Stratum = node of the linked list
 type :: stratum
-    real                   :: thickness      ! Layer thickness [m]
-    real                   :: top_elevation  ! Layer top_elevation (top height from the surface) [m]
-    real                   :: co2ice_volfrac ! CO2 ice volumetric fraction
-    real                   :: h2oice_volfrac ! H2O ice volumetric fraction
-    real                   :: dust_volfrac   ! Dust volumetric fraction
-    real                   :: air_volfrac    ! Air volumetric fraction inside pores
-    type(stratum), pointer :: up => null()   ! Upper stratum (next node)
-    type(stratum), pointer :: down => null() ! Lower stratum (previous node)
+    real                   :: thickness       ! Layer thickness [m]
+    real                   :: top_elevation   ! Layer top_elevation (top height from the surface) [m]
+    real                   :: co2ice_volfrac  ! CO2 ice volumetric fraction
+    real                   :: h2oice_volfrac  ! H2O ice volumetric fraction
+    real                   :: dust_volfrac    ! Dust volumetric fraction
+    real                   :: pore_volfrac    ! Pore volumetric fraction = pores
+    real                   :: poreice_volfrac ! Pore-filled ice volumetric fraction
+    type(stratum), pointer :: up => null()    ! Upper stratum (next node)
+    type(stratum), pointer :: down => null()  ! Lower stratum (previous node)
 end type stratum
 
@@ -85,10 +86,11 @@
 
 !---- Code
-write(*,'(a20,es13.6)') 'Thickness      : ', this%thickness
-write(*,'(a20,es13.6)') 'Top elevation  : ', this%top_elevation
-write(*,'(a20,es13.6)') 'CO2 ice (m3/m3): ', this%co2ice_volfrac
-write(*,'(a20,es13.6)') 'H2O ice (m3/m3): ', this%h2oice_volfrac
-write(*,'(a20,es13.6)') 'Dust (m3/m3)   : ', this%dust_volfrac
-write(*,'(a20,es13.6)') 'Air (m3/m3)    : ', this%air_volfrac
+write(*,'(a,es13.6)') 'Thickness              : ', this%thickness
+write(*,'(a,es13.6)') 'Top elevation          : ', this%top_elevation
+write(*,'(a,es13.6)') 'CO2 ice (m3/m3)        : ', this%co2ice_volfrac
+write(*,'(a,es13.6)') 'H2O ice (m3/m3)        : ', this%h2oice_volfrac
+write(*,'(a,es13.6)') 'Dust (m3/m3)           : ', this%dust_volfrac
+write(*,'(a,es13.6)') 'Porosity (m3/m3)       : ', this%pore_volfrac
+write(*,'(a,es13.6)') 'Pore-filled ice (m3/m3): ', this%poreice_volfrac
 
 END SUBROUTINE print_stratum
@@ -96,5 +98,5 @@
 !=======================================================================
 ! To add a stratum to the top of a layering
-SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)
+SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
 
 implicit none
@@ -102,5 +104,5 @@
 !---- Arguments
 type(layering), intent(inout)  :: this
-real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, air
+real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, pore, poreice
 
 !---- Local variables
@@ -116,5 +118,6 @@
 str%h2oice_volfrac = 0.
 str%dust_volfrac = 0.
-str%air_volfrac = 1.
+str%pore_volfrac = 1.
+str%poreice_volfrac = 0.
 if (present(thickness)) str%thickness = thickness
 if (present(top_elevation)) str%top_elevation = top_elevation
@@ -122,8 +125,9 @@
 if (present(h2oice)) str%h2oice_volfrac = h2oice
 if (present(dust)) str%dust_volfrac = dust
-if (present(air)) str%air_volfrac = air
+if (present(pore)) str%pore_volfrac = pore
+if (present(poreice)) str%poreice_volfrac = poreice
 
 ! Verification of volume fraction
-if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) then
+if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%pore_volfrac)) > tol) then
     call print_stratum(str)
     error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!'
@@ -147,5 +151,5 @@
 !=======================================================================
 ! To insert a stratum after another one in a layering
-SUBROUTINE insert_stratum(this,str_prev,thickness,top_elevation,co2ice,h2oice,dust,air)
+SUBROUTINE insert_stratum(this,str_prev,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
 
 implicit none
@@ -154,5 +158,5 @@
 type(layering),         intent(inout) :: this
 type(stratum), pointer, intent(inout) :: str_prev
-real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, air
+real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, pore, poreice
 
 !---- Local variables
@@ -161,5 +165,5 @@
 !---- Code
 if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
-    call add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)
+    call add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice)
 else
     ! Creation of the stratum
@@ -171,5 +175,6 @@
     str%h2oice_volfrac = 0.
     str%dust_volfrac = 0.
-    str%air_volfrac = 1.
+    str%pore_volfrac = 1.
+    str%poreice_volfrac = 0.
     if (present(thickness)) str%thickness = thickness
     if (present(top_elevation)) str%top_elevation = top_elevation
@@ -177,8 +182,9 @@
     if (present(h2oice)) str%h2oice_volfrac = h2oice
     if (present(dust)) str%dust_volfrac = dust
-    if (present(air)) str%air_volfrac = air
+    if (present(pore)) str%pore_volfrac = pore
+    if (present(poreice)) str%poreice_volfrac = poreice
 
     ! Verification of volume fraction
-    if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) then
+    if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%pore_volfrac)) > tol) then
         call print_stratum(str)
         error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!'
@@ -326,5 +332,5 @@
 i = this%nb_str
 do while (associated(current))
-    write(*,'(a8,i4)') 'Stratum ', i
+    write(*,'(a,i4)') 'Stratum ', i
     call print_stratum(current)
     write(*,*)
@@ -386,5 +392,6 @@
             stratif_array(ig,islope,k,4) = current%h2oice_volfrac
             stratif_array(ig,islope,k,5) = current%dust_volfrac
-            stratif_array(ig,islope,k,6) = current%air_volfrac
+            stratif_array(ig,islope,k,6) = current%pore_volfrac
+            stratif_array(ig,islope,k,7) = current%poreice_volfrac
             current => current%up
             k = k + 1
@@ -412,13 +419,14 @@
 do islope = 1,nslope
     do ig = 1,ngrid
-        stratif(ig,islope)%bottom%thickness      = stratif_array(ig,islope,1,1)
-        stratif(ig,islope)%bottom%top_elevation  = stratif_array(ig,islope,1,2)
-        stratif(ig,islope)%bottom%co2ice_volfrac = stratif_array(ig,islope,1,3)
-        stratif(ig,islope)%bottom%h2oice_volfrac = stratif_array(ig,islope,1,4)
-        stratif(ig,islope)%bottom%dust_volfrac   = stratif_array(ig,islope,1,5)
-        stratif(ig,islope)%bottom%air_volfrac    = stratif_array(ig,islope,1,6)
+        stratif(ig,islope)%bottom%thickness       = stratif_array(ig,islope,1,1)
+        stratif(ig,islope)%bottom%top_elevation   = stratif_array(ig,islope,1,2)
+        stratif(ig,islope)%bottom%co2ice_volfrac  = stratif_array(ig,islope,1,3)
+        stratif(ig,islope)%bottom%h2oice_volfrac  = stratif_array(ig,islope,1,4)
+        stratif(ig,islope)%bottom%dust_volfrac    = stratif_array(ig,islope,1,5)
+        stratif(ig,islope)%bottom%pore_volfrac    = stratif_array(ig,islope,1,6)
+        stratif(ig,islope)%bottom%poreice_volfrac = stratif_array(ig,islope,1,7)
         do k = 2,size(stratif_array,3)
             if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit
-            call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6))
+            call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6),stratif_array(ig,islope,k,7))
         enddo
     enddo
@@ -490,5 +498,5 @@
     do while (h2lift_tot > 0. .and. associated(current1))
         ! Is the considered stratum a dust lag?
-        if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we can lift dust
+        if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we can lift dust
             h_dust_old = current1%dust_volfrac*current1%thickness
             ! How much dust can be lifted in the considered dust lag?
@@ -524,5 +532,5 @@
         if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
-                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity)
+                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity,0.)
                  new_str = .false.
              else
@@ -540,5 +548,5 @@
         if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
-                 call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness))
+                 call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness),0.)
                  new_str = .false.
              else
@@ -563,5 +571,5 @@
             if (associated(current1%up)) then ! If there is an upper stratum
                 currentloc => current1%up
-                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
                     h_toplag = h_toplag + currentloc%thickness
                     if (.not. associated(currentloc%up)) exit
@@ -578,5 +586,5 @@
             else if (h_co2ice_old < eps) then ! There is nothing to sublimate
                 ! Is the considered stratum a dust lag?
-                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
+                if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
                     current1 => current1%down
                     new_lag1 = .true.
@@ -600,5 +608,5 @@
         if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
-                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness))
+                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness),0.)
                  new_str = .false.
              else
@@ -623,5 +631,5 @@
             if (associated(current1%up)) then ! If there is an upper stratum
                 currentloc => current1%up
-                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
                     h_toplag = h_toplag + currentloc%thickness
                     if (.not. associated(currentloc%up)) exit
@@ -638,5 +646,5 @@
             else if (h_co2ice_old < eps) then ! There is nothing to sublimate
                 ! Is the considered stratum a dust lag?
-                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
+                if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
                     current1 => current1%down
                     new_lag1 = .true.
@@ -667,5 +675,5 @@
             if (associated(current2%up)) then ! If there is an upper stratum
                 currentloc => current2%up
-                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol)
                     h_toplag = h_toplag + currentloc%thickness
                     if (.not. associated(currentloc%up)) exit
@@ -682,5 +690,5 @@
             else if (h_h2oice_old < eps) then ! There is nothing to sublimate
                 ! Is the considered stratum a dust lag?
-                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
+                if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum
                     current1 => current1%down
                     new_lag1 = .true.
@@ -743,5 +751,6 @@
     str%h2oice_volfrac = h_h2oice_old/str%thickness
     str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
-    str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
+    str%pore_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
+    str%poreice_volfrac = 0.
     ! Correct the value of top_elevation for the upper strata
     current => str%up
@@ -801,5 +810,6 @@
     str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
     str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
-    str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
+    str%pore_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
+    str%poreice_volfrac = 0.
     ! Correct the value of top_elevation for the upper strata
     current => str%up
@@ -838,6 +848,6 @@
 !---- Code
 ! Update of properties in the eroding dust lag
-str%thickness = str%thickness - h2lift/(1. - str%air_volfrac)
-str%top_elevation = str%top_elevation - h2lift/(1. - str%air_volfrac)
+str%thickness = str%thickness - h2lift/(1. - str%pore_volfrac)
+str%top_elevation = str%top_elevation - h2lift/(1. - str%pore_volfrac)
 ! Remove the eroding dust lag if there is no dust anymore
 if (str%thickness < tol) call rm_stratum(this,str)
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3672)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3673)
@@ -145,5 +145,5 @@
 real, dimension(:),     allocatable :: ps_start      ! surface pressure in the start file
 real, dimension(:),     allocatable :: ps_start0     ! surface pressure in the start file at the beginning
-real, dimension(:),     allocatable :: ps_avg        ! (ngrid) Averaged surface pressure
+real, dimension(:),     allocatable :: ps_avg        ! (ngrid) Average surface pressure
 real, dimension(:),     allocatable :: ps_dev        ! (ngrid x timelen) Surface pressure deviation
 real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure
@@ -206,9 +206,9 @@
 
 ! Variables for surface and soil
-real, dimension(:,:),     allocatable :: tsurf_avg                        ! Grid points x Slope field: Averaged surface temperature [K]
+real, dimension(:,:),     allocatable :: tsurf_avg                        ! Grid points x Slope field: Average surface temperature [K]
 real, dimension(:,:),     allocatable :: tsurf_dev                        ! Grid points x Slope field: Surface temperature deviation [K]
-real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Grid points x Slope field: Averaged surface temperature of first call of the PCM [K]
-real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Grid points x Soil x Slope field: Averaged Soil Temperature [K]
-real, dimension(:,:),     allocatable :: tsoil_avg_old                    ! Grid points x Soil field: Averaged Soil Temperature at the previous time step [K]
+real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Grid points x Slope field: Average surface temperature of first call of the PCM [K]
+real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Grid points x Soil x Slope field: Average Soil Temperature [K]
+real, dimension(:,:),     allocatable :: tsoil_avg_old                    ! Grid points x Soil field: Average Soil Temperature at the previous time step [K]
 real, dimension(:,:,:),   allocatable :: tsoil_dev                        ! Grid points x Soil x Slope field: Soil temperature deviation [K]
 real, dimension(:,:,:,:), allocatable :: tsoil_timeseries                 ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
@@ -216,7 +216,7 @@
 real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old         ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K]
 real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
-real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Grid points x Slope: Averaged  water surface density [kg/m^3]
+real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Grid points x Slope: Average water surface density [kg/m^3]
 real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
-real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Averaged water soil density [kg/m^3]
+real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Average water soil density [kg/m^3]
 real, dimension(:),       allocatable :: delta_co2_adsorbed               ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
 real, dimension(:),       allocatable :: delta_h2o_adsorbed               ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
@@ -374,5 +374,5 @@
 call hostnm(hostname) ! Machine/station name
 write(*,*)
-write(*,*) '######### PEM information #########'
+write(*,*) '********* PEM information *********'
 write(*,*) '+ User     : '//trim(logname)
 write(*,*) '+ Machine  : '//trim(hostname)
@@ -406,5 +406,5 @@
 !------------------------
 write(*,*)
-write(*,*) '######### PEM initialization #########'
+write(*,*) '********* PEM initialization *********'
 write(*,*) '> Reading "run.def" (PEM)'
 #ifndef CPP_1D
@@ -638,6 +638,6 @@
 ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
 ps_avg_global_new = ps_avg_global_ini
-write(*,*) "Total surface of the planet      =", total_surface
-write(*,*) "Initial global averaged pressure =", ps_avg_global_ini
+write(*,*) "Total surface of the planet     =", total_surface
+write(*,*) "Initial global average pressure =", ps_avg_global_ini
 
 !------------------------
@@ -732,5 +732,5 @@
 !------------------------
 write(*,*)
-write(*,*) '######### PEM cycle #########'
+write(*,*) '********* PEM cycle *********'
 i_myear_leg = 0
 stopPEM = 0
@@ -750,5 +750,5 @@
 do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
 ! II.a.1. Compute updated global pressure
-    write(*,'(a,f10.2)') ' #### Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
+    write(*,'(a,f10.2)') ' **** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
     write(*,*) "> Updating the surface pressure"
     ps_avg_global_old = ps_avg_global_new
@@ -775,7 +775,5 @@
         enddo
     enddo
-    do ig = 1,ngrid
-        ps_timeseries(ig,:) = ps_timeseries(ig,:)*ps_avg_global_new/ps_avg_global_old
-    enddo
+    ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old
     write(*,*) "> Updating the pressure levels timeseries for the new pressure"
     allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
@@ -893,5 +891,5 @@
     if (soil_pem) then
 ! II_d.2 Shifting soil temperature to surface
-        call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
+        !call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
 
 ! II_d.3 Update soil temperature
@@ -939,5 +937,5 @@
             enddo
             deallocate(porefill)
-            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
         endif
         deallocate(icetable_thickness_old,ice_porefilling_old)
@@ -1040,27 +1038,27 @@
         select case (stopPEM)
             case(1)
-                write(*,'(a,i0,a)') " #### STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
+                write(*,'(a,i0,a)') " **** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
             case(2)
-                write(*,'(a,i0,a)') " #### STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
+                write(*,'(a,i0,a)') " **** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
             case(3)
-                write(*,'(a,i0,a)') " #### STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
+                write(*,'(a,i0,a)') " **** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
             case(4)
-                write(*,'(a,i0,a)') " #### STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
+                write(*,'(a,i0,a)') " **** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
             case(5)
-                write(*,'(a,i0)') " #### STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
+                write(*,'(a,i0)') " **** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
             case(6)
-                write(*,'(a,i0)') " #### STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
+                write(*,'(a,i0)') " **** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
             case(7)
-                write(*,'(a,i0)') " #### STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
+                write(*,'(a,i0)') " **** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
             case(8)
-                write(*,'(a,i0)') " #### STOPPING because the layering algorithm met an hasty end: ", stopPEM
+                write(*,'(a,i0)') " **** STOPPING because the layering algorithm met an hasty end: ", stopPEM
             case default
-                write(*,'(a,i0)') " #### STOPPING with unknown stopping criterion code: ", stopPEM
+                write(*,'(a,i0)') " **** STOPPING with unknown stopping criterion code: ", stopPEM
         end select
         exit
     else
-        write(*,'(a,f10.2,a)') ' #### The chained simulation has run for ',i_myear,' Martian years.'
-        write(*,*) '#### The PEM can continue!'
-        write(*,*) '####'
+        write(*,'(a,f10.2,a)') ' **** The chained simulation has run for ',i_myear,' Martian years.'
+        write(*,*) '**** The PEM can continue!'
+        write(*,*) '****'
     endif
 enddo ! big time iteration loop of the pem
@@ -1081,5 +1079,5 @@
 ! III_a.1 Ice update for start file
 write(*,*)
-write(*,*) '######### PEM finalization #########'
+write(*,*) '********* PEM finalization *********'
 write(*,*) '> Reconstructing perennial ice and frost for the PCM'
 watercap = 0.
@@ -1133,6 +1131,5 @@
 write(*,*) '> Reconstructing the pressure for the PCM'
 allocate(ps_start(ngrid))
-! The pressure deviation is rescaled to avoid disproportionate oscillations in case of huge average pressure drop
-ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini
+ps_start = (ps_avg + ps_dev)*ps_avg_global_new/ps_avg_global_ini ! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop
 deallocate(ps_avg,ps_dev)
 
@@ -1288,10 +1285,10 @@
 
 write(*,*)
-write(*,*) '######### PEM finalization #########'
+write(*,*) '********* PEM finalization *********'
 write(*,'(a,f10.2,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
 write(*,'(a,f10.2,a,f10.2,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years."
 write(*,'(a,f10.2,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years."
 write(*,*) "+ PEM: so far, so good!"
-write(*,*) '####################################'
+write(*,*) '************************************'
 
 if (layering_algo) then
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3672)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3673)
@@ -65,12 +65,12 @@
 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
 ! local
-real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem           ! soil temperature saved in the start [K]
-real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM              ! soil thermal inertia saved in the start [SI]
-logical                                 :: found, found2            ! check if variables are found in the start
-integer                                 :: iloop, ig, islope, isoil ! index for loops
-real                                    :: delta                    ! Depth of the interface regolith-breccia, breccia -bedrock [m]
-character(2)                            :: num                      ! intermediate string to read PEM start sloped variables
-logical                                 :: startpem_file            ! boolean to check if we read the startfile or not
-real, dimension(:,:,:,:), allocatable   :: stratif_array            ! Array for stratification (layerings)
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem    ! soil temperature saved in the start [K]
+real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM       ! soil thermal inertia saved in the start [SI]
+logical                                 :: found, found2     ! check if variables are found in the start
+integer                                 :: iloop, ig, islope ! index for loops
+real                                    :: delta             ! Depth of the interface regolith-breccia, breccia -bedrock [m]
+character(2)                            :: num               ! intermediate string to read PEM start sloped variables
+logical                                 :: startpem_file     ! boolean to check if we read the startfile or not
+real, dimension(:,:,:,:), allocatable   :: stratif_array     ! Array for stratification (layerings)
 
 #ifdef CPP_STD
@@ -143,5 +143,5 @@
             nb_str_max = int(inquire_dimension_length('nb_str_max'))
         endif
-        allocate(stratif_array(ngrid,nslope,nb_str_max,6))
+        allocate(stratif_array(ngrid,nslope,nb_str_max,7))
         stratif_array = 0.
         do islope = 1,nslope
@@ -156,5 +156,7 @@
             call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
             found = found .or. found2
-            call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
+            call get_field('stratif_slope'//num//'_pore_volfrac',stratif_array(:,islope,:,6),found2)
+            found = found .or. found2
+            call get_field('stratif_slope'//num//'_icepore_volfrac',stratif_array(:,islope,:,7),found2)
             found = found .or. found2
             if (.not. found) then
@@ -267,12 +269,8 @@
                 call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
 
-                do iloop = nsoil_PCM + 1,nsoil_PEM
-                    tsoil_PEM(:,iloop,islope) = tsoil_startpem(:,iloop,islope)
-                enddo
+                tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) = tsoil_startpem(:,nsoil_PCM + 1:nsoil_PEM,islope)
             endif !found
 
-            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
+            watersoil_avg(:,nsoil_PCM + 1:nsoil_PEM,islope) = exp(beta_clap_h2o/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) + alpha_clap_h2o)/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
         enddo ! islope
         write(*,*) 'PEMETAT0: TSOIL done'
@@ -431,9 +429,5 @@
 
 ! First raw initialization
-            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
-            enddo
+            watersoil_avg(:,nsoil_PCM + 1:nsoil_PEM,islope) = exp(beta_clap_h2o/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) + alpha_clap_h2o)/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
         enddo !islope
         write(*,*) 'PEMETAT0: TSOIL done'
Index: trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 3672)
+++ trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 3673)
@@ -99,5 +99,5 @@
 
 if (layering_algo) then
-    allocate(stratif_array(ngrid,nslope,nb_str_max,6))
+    allocate(stratif_array(ngrid,nslope,nb_str_max,7))
     call stratif2array(stratif,ngrid,nslope,stratif_array)
     do islope = 1,nslope
@@ -108,5 +108,6 @@
         call put_field('stratif_slope'//num//'_h2oice_volfrac','Layering H2O ice volume fraction',stratif_array(:,islope,:,4),Year)
         call put_field('stratif_slope'//num//'_dust_volfrac','Layering dust volume fraction',stratif_array(:,islope,:,5),Year)
-        call put_field('stratif_slope'//num//'_air_volfrac','Layering air volume fraction',stratif_array(:,islope,:,6),Year)
+        call put_field('stratif_slope'//num//'_pore_volfrac','Layering pore volume fraction',stratif_array(:,islope,:,6),Year)
+        call put_field('stratif_slope'//num//'_icepore_volfrac','Layering ice pore volume fraction',stratif_array(:,islope,:,7),Year)
     enddo
     deallocate(stratif_array)
