Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3136)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3143)
@@ -145,2 +145,7 @@
 == 22/11/2023 == JBC
 Update of files in the deftank: addition of variables in the xml definition files, inclusion of "callphys.def" into "run_PEM.def" and minor typo in "launch_pem.sh".
+
+== 29/11/2023 == JBC
+'Watercap' has been removed from the water ice evolution since 'water_reservoir' does already the job + Some cleanings to simplify the code.
+
+/!\ Commit for the PEM management of h2o ice before a rework of ice management in the PEM!
Index: trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 3136)
+++ trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 3143)
@@ -1,69 +1,76 @@
-module comsoil_h_PEM
+MODULE comsoil_h_PEM
 
 implicit none
-  integer, parameter :: nsoilmx_PEM = 68                 ! number of layers in the PEM
-  real,save,allocatable,dimension(:) :: layer_PEM        ! soil layer depths   [m]
-  real,save,allocatable,dimension(:) :: mlayer_PEM       ! soil mid-layer depths [m]
-  real,save,allocatable,dimension(:,:,:) :: TI_PEM       ! soil thermal inertia [SI]
-  real,save,allocatable,dimension(:,:) :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI]
-  ! variables (FC: built in firstcall in soil.F)
-  REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:)              ! sub-surface temperatures [K]
-  real,save,allocatable :: mthermdiff_PEM(:,:)           ! (FC) mid-layer thermal diffusivity [SI]
-  real,save,allocatable :: thermdiff_PEM(:,:)            ! (FC) inter-layer thermal diffusivity [SI]
-  real,save,allocatable :: coefq_PEM(:)                  ! (FC) q_{k+1/2} coefficients [SI]
-  real,save,allocatable :: coefd_PEM(:,:)                ! (FC) d_k coefficients [SI]
-  real,save,allocatable :: alph_PEM(:,:)                 ! (FC) alpha_k coefficients [SI]
-  real,save,allocatable :: 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,allocatable,dimension(:) :: water_reservoir  ! Large reserve of water   [kg/m^2]
-  real,save :: water_reservoir_nom                       ! Nominal value for the large reserve of water   [kg/m^2]
-  logical, save :: reg_thprop_dependp                    ! thermal properites of the regolith vary with the surface pressure
 
+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]
+! 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, allocatable, dimension(:),     save :: water_reservoir     ! Large reserve of water [kg/m^2]
+real,                                save :: water_reservoir_nom ! Nominal value for the large reserve of water [kg/m^2]
+logical,                             save :: reg_thprop_dependp  ! thermal properites of the regolith vary with the surface pressure
+
+!=======================================================================
 contains
+!=======================================================================
 
-  subroutine ini_comsoil_h_PEM(ngrid,nslope)
-  
-  implicit none
-  integer,intent(in) :: ngrid ! number of atmospheric columns
-  integer,intent(in) :: nslope ! number of slope within a mesh 
+SUBROUTINE ini_comsoil_h_PEM(ngrid,nslope)
 
-    allocate(layer_PEM(nsoilmx_PEM)) 
-    allocate(mlayer_PEM(0:nsoilmx_PEM-1)) 
-    allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
-    allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope)) 
-    allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1))
-    allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1))
-    allocate(coefq_PEM(0:nsoilmx_PEM-1))
-    allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
-    allocate(alph_PEM(ngrid,nsoilmx_PEM-1))
-    allocate(beta_PEM(ngrid,nsoilmx_PEM-1))
-    allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) 
-    allocate(water_reservoir(ngrid))
-  end subroutine ini_comsoil_h_PEM
+implicit none
 
+integer, intent(in) :: ngrid  ! number of atmospheric columns
+integer, intent(in) :: nslope ! number of slope within a mesh
 
-  subroutine end_comsoil_h_PEM
+allocate(layer_PEM(nsoilmx_PEM))
+allocate(mlayer_PEM(0:nsoilmx_PEM-1))
+allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
+allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
+allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1))
+allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1))
+allocate(coefq_PEM(0:nsoilmx_PEM-1))
+allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
+allocate(alph_PEM(ngrid,nsoilmx_PEM-1))
+allocate(beta_PEM(ngrid,nsoilmx_PEM-1))
+allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
+allocate(water_reservoir(ngrid))
 
-  implicit none
+END SUBROUTINE ini_comsoil_h_PEM
 
-    if (allocated(layer_PEM)) deallocate(layer_PEM)
-    if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
-    if (allocated(TI_PEM)) deallocate(TI_PEM)
-    if (allocated(tsoil_PEM)) deallocate(tsoil_PEM)
-    if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM)
-    if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM)
-    if (allocated(coefq_PEM)) deallocate(coefq_PEM) 
-    if (allocated(coefd_PEM)) deallocate(coefd_PEM)
-    if (allocated(alph_PEM)) deallocate(alph_PEM)
-    if (allocated(beta_PEM)) deallocate(beta_PEM)
-    if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
-    if (allocated(water_reservoir)) deallocate(water_reservoir)
-  end subroutine end_comsoil_h_PEM
+!=======================================================================
 
-end module comsoil_h_PEM
+SUBROUTINE end_comsoil_h_PEM
+
+implicit none
+
+if (allocated(layer_PEM)) deallocate(layer_PEM)
+if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
+if (allocated(TI_PEM)) deallocate(TI_PEM)
+if (allocated(tsoil_PEM)) deallocate(tsoil_PEM)
+if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM)
+if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM)
+if (allocated(coefq_PEM)) deallocate(coefq_PEM)
+if (allocated(coefd_PEM)) deallocate(coefd_PEM)
+if (allocated(alph_PEM)) deallocate(alph_PEM)
+if (allocated(beta_PEM)) deallocate(beta_PEM)
+if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
+if (allocated(water_reservoir)) deallocate(water_reservoir)
+
+END SUBROUTINE end_comsoil_h_PEM
+
+END MODULE comsoil_h_PEM
Index: trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90	(revision 3136)
+++ trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90	(revision 3143)
@@ -38,7 +38,6 @@
 
 ! Conversion H2O/CO2 frost to perennial frost and vice versa
-real, parameter :: threshold_water_frost2perennial = 1000. !~ 1 m
-real, parameter :: threshold_co2_frost2perennial = 16000.   !~ 10 m
+real, parameter :: threshold_h2o_frost2perennial = 1000.  !~ 1 m
+real, parameter :: threshold_co2_frost2perennial = 16000. !~ 10 m
 
 END MODULE constants_marspem_mod
-
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 3136)
+++ trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 3143)
@@ -3,5 +3,7 @@
 implicit none
 
+! ----------------------------------------------------------------------
 contains
+! ----------------------------------------------------------------------
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -44,21 +46,28 @@
 STOPPING = .false.
 
-! Computation of the present surface of water ice sublimating
+! Computation of the present surface of h2o ice
 present_surf = 0.
 do i = 1,ngrid
     do islope = 1,nslope
-        !if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
-        if (initial_h2o_ice(i,islope) > 0.5) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
+        if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
     enddo
 enddo
 
 ! Check of the criterion
-if (present_surf < ini_surf*(1. - water_ice_criterion) .or. present_surf > ini_surf*(1. + water_ice_criterion)) then
+if (present_surf < ini_surf*(1. - water_ice_criterion)) then
     STOPPING = .true.
     write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
+    write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", present_surf < ini_surf*(1. - water_ice_criterion)
     write(*,*) "Current surface of water ice sublimating =", present_surf
     write(*,*) "Initial surface of water ice sublimating =", ini_surf
     write(*,*) "Percentage of change accepted =", water_ice_criterion*100
-    write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", (present_surf < ini_surf*(1. - water_ice_criterion))
+endif
+if (present_surf > ini_surf*(1. + water_ice_criterion)) then
+    STOPPING = .true.
+    write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
+    write(*,*) "present_surf > ini_surf*(1. + water_ice_criterion)", present_surf > ini_surf*(1. + water_ice_criterion)
+    write(*,*) "Current surface of water ice sublimating =", present_surf
+    write(*,*) "Initial surface of water ice sublimating =", ini_surf
+    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
 endif
 
@@ -105,5 +114,5 @@
 STOPPING_ps = .false.
 
-! Computation of the actual surface
+! Computation of the present surface of co2 ice
 present_surf = 0.
 do i = 1,ngrid
@@ -114,22 +123,38 @@
 
 ! Check of the criterion
-if (present_surf < ini_surf*(1. - co2_ice_criterion) .or. present_surf > ini_surf*(1. + co2_ice_criterion)) then
+if (present_surf < ini_surf*(1. - co2_ice_criterion)) then
     STOPPING_ice = .true.
-    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold"
+    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reached the threshold"
+    write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", present_surf < ini_surf*(1. - co2_ice_criterion)
     write(*,*) "Current surface of co2 ice sublimating =", present_surf
     write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
     write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
-    write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", (present_surf < ini_surf*(1. - co2_ice_criterion))
+endif
+if (present_surf > ini_surf*(1. + co2_ice_criterion)) then
+    STOPPING_ice = .true.
+    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold"
+    write(*,*) "present_surf > ini_surf*(1. + co2_ice_criterion)", present_surf > ini_surf*(1. + co2_ice_criterion)
+    write(*,*) "Current surface of co2 ice sublimating =", present_surf
+    write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
+    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
 endif
 
 if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false.
 
-if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion) .or. global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then
+if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) then
     STOPPING_ps = .true.
     write(*,*) "Reason of stopping: the global pressure reach the threshold"
+    write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)
     write(*,*) "Current global pressure =", global_ave_press_new
     write(*,*) "PCM global pressure =", global_ave_press_GCM
     write(*,*) "Percentage of change accepted =", ps_criterion*100.
-    write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion))
+endif
+if (global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then
+    STOPPING_ps = .true.
+    write(*,*) "Reason of stopping: the global pressure reach the threshold"
+    write(*,*) "global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)", global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)
+    write(*,*) "Current global pressure =", global_ave_press_new
+    write(*,*) "PCM global pressure =", global_ave_press_GCM
+    write(*,*) "Percentage of change accepted =", ps_criterion*100.
 endif
 
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3136)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3143)
@@ -1,10 +1,10 @@
 !------------------------
 ! I   Initialization
-!    I_a READ run.def
-!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_a Read run.def
+!    I_b Read of start_evol.nc and starfi_evol.nc
 !    I_c Subslope parametrisation
-!    I_d READ PCM data and convert to the physical grid
+!    I_d Read PCM data and convert to the physical grid
 !    I_e Initialization of the PEM variable and soil
-!    I_f Compute tendencies & Save initial situation
+!    I_f Compute tendencies
 !    I_g Save initial PCM situation
 !    I_h Read the startpem.nc
@@ -42,5 +42,5 @@
 use criterion_pem_stop_mod,       only: criterion_waterice_stop, criterion_co2_stop
 use constants_marspem_mod,        only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
-                                        m_noco2, threshold_water_frost2perennial, threshold_co2_frost2perennial
+                                        m_noco2, threshold_h2o_frost2perennial, threshold_co2_frost2perennial
 use evol_co2_ice_s_mod,           only: evol_co2_ice_s
 use evol_h2o_ice_s_mod,           only: evol_h2o_ice_s
@@ -86,4 +86,5 @@
     use paleoclimate_mod,   only: albedo_perennialco2
     use comcstfi_h,         only: pi, rad, g, cpp, mugaz, r
+    use surfini_mod,        only: surfini
 #else
     use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
@@ -131,5 +132,5 @@
 
 ! Variables to read start.nc
-character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM
+character(*), parameter :: start_name = "start_evol.nc" ! Name of the file used to initialize the PEM
 
 ! Dynamic variables
@@ -138,7 +139,7 @@
 real, dimension(ip1jmp1,llm)        :: teta          ! temperature potentielle
 real, dimension(:,:,:), allocatable :: q             ! champs advectes
-real, dimension(ip1jmp1)            :: ps            ! pression  au sol
+real, dimension(ip1jmp1)            :: ps            ! pression au sol
 real, dimension(:),     allocatable :: ps_start_GCM  ! (ngrid) pression  au sol
-real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) ! pression  au sol instantannées
+real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées
 real, dimension(ip1jmp1,llm)        :: masse         ! masse d'air
 real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
@@ -146,15 +147,15 @@
 
 ! Variables to read starfi.nc
-character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsing the PEM
-character(2)                   :: str2
-integer                        :: ncid, varid, status                      ! Variable for handling opening of files
-integer                        :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files
-integer                        :: lonvarid, latvarid, areavarid, sdvarid   ! Variable ID for Netcdf files
-integer                        :: apvarid, bpvarid                         ! Variable ID for Netcdf files
+character(*), parameter :: startfi_name = "startfi_evol.nc" ! Name of the file used to initialize the PEM
+character(2)            :: str2
+integer                 :: ncid, varid, status                      ! Variable for handling opening of files
+integer                 :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files
+integer                 :: lonvarid, latvarid, areavarid, sdvarid   ! Variable ID for Netcdf files
+integer                 :: apvarid, bpvarid                         ! Variable ID for Netcdf files
 
 ! Variables to read starfi.nc and write restartfi.nc
-real, dimension(:), allocatable :: longitude     ! Longitude read in FILE_NAME and written in restartfi
-real, dimension(:), allocatable :: latitude      ! Latitude read in FILE_NAME and written in restartfi
-real, dimension(:), allocatable :: cell_area     ! Cell_area read in FILE_NAME and written in restartfi
+real, dimension(:), allocatable :: longitude     ! Longitude read in startfi_name and written in restartfi
+real, dimension(:), allocatable :: latitude      ! Latitude read in startfi_name and written in restartfi
+real, dimension(:), allocatable :: cell_area     ! Cell_area read in startfi_name and written in restartfi
 real                            :: Total_surface ! Total surface of the planet
 
@@ -187,6 +188,4 @@
 real, dimension(:,:), allocatable :: co2frost_ini         ! Initial amount of co2 frost (at the start of the PEM run)
 real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run)
-
-
 
 ! Variables for slopes
@@ -243,9 +242,9 @@
     real                                :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
     real                                :: albedo_h2o_frost              ! albedo of h2o frost
-    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
-    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
-    real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
-    real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimensiion when reading form generic
-    real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimensiion when reading form generic
+    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
+    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
+    real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
+    real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimension when reading form generic
+    real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimension when reading form generic
     real, dimension(:,:),   allocatable :: tsurf                         ! Subslope variable, only needed in the GENERIC case
     real, dimension(:,:,:), allocatable :: qsurf                         ! Subslope variable, only needed in the GENERIC case
@@ -259,7 +258,7 @@
 
 #ifdef CPP_1D
-    integer                           :: nsplit_phys
-    integer, parameter                :: jjm_value = jjm - 1
-    integer                           :: ierr
+    integer            :: nsplit_phys
+    integer, parameter :: jjm_value = jjm - 1
+    integer            :: ierr
 
     ! Dummy variables to use the subroutine 'init_testphys1d'
@@ -272,7 +271,7 @@
     real, dimension(nlayer + 1)         :: plev
 #else
-    integer, parameter                :: jjm_value = jjm
-    real, dimension(:), allocatable   :: ap ! Coefficient ap read in FILE_NAME_start and written in restart
-    real, dimension(:), allocatable   :: bp ! Coefficient bp read in FILE_NAME_start and written in restart
+    integer, parameter              :: jjm_value = jjm
+    real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
+    real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
 #endif
 
@@ -302,5 +301,5 @@
 !------------------------
 ! I   Initialization
-!    I_a READ run.def
+!    I_a Read run.def
 !------------------------
 #ifndef CPP_1D
@@ -315,5 +314,5 @@
     allocate(longitude(1),latitude(1),cell_area(1))
 
-    therestart1D = .false.
+    therestart1D = .false. ! Default value
     inquire(file = 'start1D_evol.txt',exist = therestart1D)
     if (.not. therestart1D) then
@@ -321,5 +320,5 @@
         error stop 'Initialization cannot be done for the 1D PEM.'
     endif
-    therestartfi = .false.
+    therestartfi = .false. ! Default value
     inquire(file = 'startfi_evol.nc',exist = therestartfi)
     if (.not. therestartfi) then
@@ -339,10 +338,10 @@
 !------------------------
 ! I   Initialization
-!    I_b READ of start_evol.nc and starfi_evol.nc
-!------------------------
-! I_b.1 READ start_evol.nc
+!    I_b Read of start_evol.nc and starfi_evol.nc
+!------------------------
+! I_b.1 Read start_evol.nc
 allocate(ps_start_GCM(ngrid))
 #ifndef CPP_1D
-    call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0)
+    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
 
     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
@@ -353,5 +352,5 @@
     allocate(ap(nlayer + 1))
     allocate(bp(nlayer + 1))
-    status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid)
+    status = nf90_open(start_name,NF90_NOWRITE,ncid)
     status = nf90_inq_varid(ncid,"ap",apvarid)
     status = nf90_get_var(ncid,apvarid,ap)
@@ -368,5 +367,5 @@
 ! In the PCM, these values are given to the physic by the dynamic.
 ! Here we simply read them in the startfi_evol.nc file
-status = nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
+status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
 
 status = nf90_inq_varid(ncid,"longitude",lonvarid)
@@ -384,8 +383,8 @@
 status = nf90_close(ncid)
 
-! I_b.2 READ startfi_evol.nc
+! I_b.2 Read startfi_evol.nc
 ! First we read the initial state (starfi.nc)
 #ifndef CPP_STD
-    call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
+    call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
                   tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,       &
                   watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
@@ -394,5 +393,5 @@
     where (qsurf < 0.) qsurf = 0.
 
-    call surfini(ngrid,qsurf)
+    call surfini(ngrid,nslope,qsurf)
 #else
     call phys_state_var_init(nq)
@@ -414,5 +413,5 @@
     allocate(albedo(ngrid,2,1))
     allocate(inertiesoil(ngrid,nsoilmx,1))
-    call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
+    call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
                   tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,            &
                   qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,     &
@@ -449,5 +448,5 @@
     if (noms(nnq) == "h2o_vap") then
         igcm_h2o_vap = nnq
-        mmol(igcm_h2o_vap)=18.
+        mmol(igcm_h2o_vap) = 18.
     endif
     if (noms(nnq) == "co2") igcm_co2 = nnq
@@ -480,5 +479,5 @@
 !------------------------
 ! I   Initialization
-!    I_d READ PCM data and convert to the physical grid
+!    I_d Read PCM data and convert to the physical grid
 !------------------------
 ! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value
@@ -552,5 +551,5 @@
 !------------------------
 ! I   Initialization
-!    I_f Compute tendencies & Save initial situation
+!    I_f Compute tendencies
 !------------------------
 allocate(tendencies_co2_ice(ngrid,nslope))
@@ -603,5 +602,5 @@
 write(*,*) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2
 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o
-write(*,*) "Total surface of the planet=", Total_surface
+write(*,*) "Total surface of the planet =", Total_surface
 allocate(zplev_gcm(ngrid,nlayer + 1))
 
@@ -617,5 +616,5 @@
 global_ave_press_GCM = global_ave_press_old
 global_ave_press_new = global_ave_press_old
-write(*,*) "Initial global average pressure=", global_ave_press_GCM
+write(*,*) "Initial global average pressure =", global_ave_press_GCM
 
 !------------------------
@@ -638,5 +637,5 @@
 do ig = 1,ngrid
     do islope = 1,nslope
-        qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
+        qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
         qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope)
     enddo
@@ -657,6 +656,6 @@
     enddo
 
-    write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
-    write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
+    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
+    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
 endif ! adsorption
 deallocate(tsurf_ave_yr1)
@@ -735,6 +734,9 @@
             q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
                                    (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
-            if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30
-            if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1.
+            if (q_h2o_PEM_phys(ig,t) < 0) then
+                q_h2o_PEM_phys(ig,t) = 1.e-30
+            else if (q_h2o_PEM_phys(ig,t) > 1) then
+                q_h2o_PEM_phys(ig,t) = 1.
+            endif
         enddo
     enddo
@@ -749,5 +751,5 @@
             if (q_co2_PEM_phys(ig,t) < 0) then
                 q_co2_PEM_phys(ig,t) = 1.e-30
-            elseif (q_co2_PEM_phys(ig,t) > 1) then
+            else if (q_co2_PEM_phys(ig,t) > 1) then
                 q_co2_PEM_phys(ig,t) = 1.
             endif
@@ -898,5 +900,5 @@
         call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
 
-! II_d.5 Update the mass of the regolith adsorbded
+! II_d.5 Update the mass of the regolith adsorbed
         if (adsorption_pem) then
             call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
@@ -985,5 +987,5 @@
         exit
     else
-        write(*,*) "We continue!"
+        write(*,*) "The PEM can continue!"
         write(*,*) "Number of iterations of the PEM: year_iter =", year_iter
         write(*,*) "Number of simulated Martian years: i_myear =", i_myear
@@ -1007,9 +1009,11 @@
         watercap_sum = 0.
         do islope = 1,nslope
-            if (qsurf(ig,igcm_h2o_ice,islope) > (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy
-                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost
+            ! water_reservoir and watercap have not changed since PCM call: here we check if we have accumulate frost or not.
+            if (qsurf(ig,igcm_h2o_ice,islope) > water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) then
+                ! 1st case: more ice than initialy. Put back ancien frost
+                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
             else
-!               2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
-                watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
+                ! 2nd case: ice sublimated. Then let's put qsurf = 0 and add the difference in watercap
+                watercap(ig,islope) = qsurf(ig,igcm_h2o_ice,islope)
                 qsurf(ig,igcm_h2o_ice,islope) = 0.
             endif
@@ -1027,13 +1031,13 @@
     enddo
     if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
-        if (water_sum > threshold_water_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
+        if (water_sum > threshold_h2o_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
             watercaptag(ig) = .true.
-            water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
+            water_reservoir(ig) = water_reservoir(ig) + threshold_h2o_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
             do islope = 1,nslope
-                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.)
+                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_h2o_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.)
             enddo
         endif
     else ! let's check that the infinite source of water disapear
-        if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perennial) then
+        if ((water_sum + water_reservoir(ig)) < threshold_h2o_frost2perennial) then
             watercaptag(ig) = .false.
             do islope = 1,nslope
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3136)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3143)
@@ -104,5 +104,5 @@
 !0. Check if the start_PEM exist.
 
-inquire(file=filename,exist =  startpem_file)
+inquire(file = filename,exist = startpem_file)
 
 write(*,*)'Is start PEM?',startpem_file
@@ -312,15 +312,10 @@
 !. 5 water reservoir
 #ifndef CPP_STD
+   water_reservoir = 0.
    call get_field("water_reservoir",water_reservoir,found)
-   if(.not.found) then
-      write(*,*)'Pemetat0: failed loading <water_reservoir>'
-      write(*,*)'will reconstruct the values from watercaptag'
-      do ig=1,ngrid
-        if(watercaptag(ig)) then
-           water_reservoir(ig)=water_reservoir_nom
-        else
-           water_reservoir(ig)=0.
-        endif
-      enddo
+   if (.not. found) then
+       write(*,*)'Pemetat0: failed loading <water_reservoir>'
+       write(*,*)'will reconstruct the values from watercaptag'
+       where (watercaptag) water_reservoir = water_reservoir_nom
    endif
 #endif
@@ -481,11 +476,6 @@
 !. e) water reservoir
 #ifndef CPP_STD
-    do ig = 1,ngrid
-        if (watercaptag(ig)) then
-            water_reservoir(ig) = water_reservoir_nom
-        else
-            water_reservoir(ig) = 0.
-        endif
-    enddo
+    water_reservoir = 0.
+    where (watercaptag) water_reservoir = water_reservoir_nom
 #endif
 
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3136)
+++ trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3143)
@@ -104,4 +104,7 @@
     write(*,*) "Data for h2o_ice_s downloaded!"
 
+    call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:))
+    write(*,*) "Data for tsurf downloaded!"
+
 #ifndef CPP_STD
     call get_var3("watercap",watercap(:,:,1,:))
@@ -110,10 +113,5 @@
     call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
     write(*,*) "Data for perennial_co2ice downloaded!"
-#endif
-
-    call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:))
-    write(*,*) "Data for tsurf downloaded!"
-
-#ifndef CPP_STD
+
     if (soil_pem) then
         call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
@@ -140,4 +138,10 @@
     write(*,*) "Data for h2o_ice_s downloaded!"
 
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
+    enddo
+    write(*,*) "Data for tsurf downloaded!"
+
 #ifndef CPP_STD
     do islope = 1,nslope
@@ -152,13 +156,5 @@
     enddo
     write(*,*) "Data for perennial_co2ice downloaded!"
-#endif
-
-    do islope = 1,nslope
-        write(num,'(i2.2)') islope
-        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
-    enddo
-    write(*,*) "Data for tsurf downloaded!"
-
-#ifndef CPP_STD
+
     if (soil_pem) then
         do islope = 1,nslope
