Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3159)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3161)
@@ -179,2 +179,9 @@
 - Addition of (CO2 and H2O) ice metamorphism: if the frost is above a given value, then the excess frost is transferred to the perennial ice.
 - Thresholds value for ice management can now be set in the "run_PEM.def".
+
+== 19/12/2023 == JBC
+PEM:
+- Addition of flags defined in the "run_PEM.def" to decide to do or not CO2 & H2O ice metamorphism: 'metam_co2ice' and 'metam_h2oice' (default is false).
+- The variations of infinite reservoirs ('watercap') during the PCM years are now taken into account to update H2O ice at the PEM initialization.
+- 'ini_h2o_bigreservoir' is renamed into 'ini_huge_h2oice'.
+- Some cleanings, in particular for the main program "pem.F90".
Index: trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 3159)
+++ trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 3161)
@@ -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(:),     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,                                save :: ini_h2o_bigreservoir ! Initial value for the large reservoir of water [kg/m^2]
-logical,                             save :: reg_thprop_dependp   ! thermal properites of the regolith vary with the surface pressure
+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
 
 !=======================================================================
Index: trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3159)
+++ trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3161)
@@ -25,10 +25,9 @@
 use time_evol_mod,         only: year_bp_ini, dt_pem, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
                           evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
-use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
+use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_huge_h2oice, depth_breccia, depth_bedrock, reg_thprop_dependp
 use adsorption_mod,        only: adsorption_pem
-use glaciers_mod,          only: co2_ice_flow, h2o_ice_flow
+use glaciers_mod,          only: h2oice_flow, co2ice_flow, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold, metam_h2oice, metam_co2ice
 use ice_table_mod,         only: icetable_equilibrium, icetable_dynamic
 use time_phylmdz_mod,      only: ecritphy
-use constants_marspem_mod, only: inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold
 
 implicit none
@@ -101,10 +100,4 @@
 call getin('adsorption_pem',adsorption_pem)
 
-co2_ice_flow = .true.
-call getin('co2_ice_flow',co2_ice_flow)
-
-h2o_ice_flow = .true.
-call getin('h2o_ice_flow',h2o_ice_flow)
-
 reg_thprop_dependp = .false.
 call getin('reg_thprop_dependp',reg_thprop_dependp)
@@ -158,16 +151,28 @@
 
 !#---------- Ice management parameters ----------#
-ini_h2o_bigreservoir = 1.e4
-call getin('ini_h2o_bigreservoir',ini_h2o_bigreservoir)
+ini_huge_h2oice = 1.e4
+call getin('ini_huge_h2oice',ini_huge_h2oice)
 
 inf_h2oice_threshold = 2.e3
 call getin('inf_h2oice_threshold',inf_h2oice_threshold)
 
+metam_h2oice = .false.
+call getin('metam_h2oice',metam_h2oice)
+
 metam_h2oice_threshold = 5.e-2
 call getin('metam_h2oice_threshold',metam_h2oice_threshold)
+
+h2oice_flow = .true.
+call getin('h2oice_flow',h2oice_flow)
+
+metam_co2ice = .false.
+call getin('metam_co2ice',metam_co2ice)
 
 metam_co2ice_threshold = 16.e3
 call getin('metam_co2ice_threshold',metam_co2ice_threshold)
 
+co2ice_flow = .true.
+call getin('co2ice_flow',co2ice_flow)
+
 END SUBROUTINE conf_pem
 
Index: trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90	(revision 3159)
+++ trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90	(revision 3161)
@@ -37,8 +37,3 @@
 real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
 
-! Thresholds for ice management
-real, save :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir
-real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice
-real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice
-
 END MODULE constants_marspem_mod
Index: trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def	(revision 3159)
+++ trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def	(revision 3161)
@@ -28,5 +28,5 @@
 #  Max_iter_pem=100000000
 
-# Acceptance rate of sublimating H2o ice surface change? Default = 0.2
+# Acceptance rate of sublimating H2O ice surface change? Default = 0.2
 # h2o_ice_crit=0.2
 
@@ -47,12 +47,9 @@
 # adsorption_pem=.true.
 
-# Do you want the glaciers to flow along subslope inside a mesh? Default = .true.
-# co2glaciersflow=.true.
-
 # Do you want to modify the soil thermal properties with the pressure? Default = .false.
 # reg_thprop_dependp=.false.
 
 #---------- Layering parameters ----------#
-# Value of the geothermal flux? Default= 0.
+# Value of the geothermal flux? Default = 0.
 # fluxgeo=0.
 
@@ -70,16 +67,28 @@
 
 #---------- Ice management parameters ----------#
-# Amount of h2o ice to initialize the big reservoir if the variable is not present in "startfi_PEM.nc"? Default = 1.e4
-# ini_h2o_bigreservoir=1.e4
+# Amount of H2O ice to initialize the huge reservoir if the variable is not present in "startfi_PEM.nc"? Default = 1.e4
+# ini_huge_h2oice=1.e4
 
 # Threshold to consider the amount of H2O ice as an infinite reservoir? Default = 2.e3
 # inf_h2oice_threshold=2.e3
 
+# Do you want H2O frost to transform into perennial H2O ice? Default = .false.
+# metam_h2oice=.false.
+
 # Threshold to consider frost is becoming perennial H2O ice? Default = 5.e-2
 # metam_h2oice_threshold=5.e-2
+
+# Do you want the H2O ice to flow along subslope inside a cell? Default = .true.
+# h2oice_flow=.true.
+
+# Do you want CO2 frost to transform into perennial CO2 ice? Default = .false.
+# metam_co2ice=.false.
 
 # Threshold to consider frost is becoming perennial CO2 ice? Default = 16.e3
 # metam_co2ice_threshold=16.e3
 
+# Do you want the CO2 ice to flow along subslope inside a cell? Default = .true.
+# co2ice_flow=.true.
+
 # Some definitions for the physics, in file 'callphys.def'
 INCLUDEDEF=callphys.def
Index: trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3159)
+++ trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3161)
@@ -3,6 +3,15 @@
 implicit none
 
-logical :: co2_ice_flow ! True by default, to compute co2 ice flow. Read in run_PEM.def
-logical :: h2o_ice_flow ! True by default, to compute h2o ice flow. Read in run_PEM.def
+
+! Flags for ice management
+logical :: h2oice_flow ! True by default, to compute H2O ice flow. Read in "run_PEM.def"
+logical :: co2ice_flow ! True by default, to compute CO2 ice flow. Read in "run_PEM.def"
+logical :: metam_h2oice ! False by default, to compute H2O ice metamorphism. Read in "run_PEM.def"
+logical :: metam_co2ice ! False by default, to compute CO2 ice metamorphism. Read in "run_PEM.def"
+
+! Thresholds for ice management
+real, save :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir
+real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice
+real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice
 
 !=======================================================================
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3159)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3161)
@@ -1,12 +1,12 @@
 !------------------------
 ! I   Initialization
-!    I_a Read run.def
-!    I_b Read of start_evol.nc and starfi_evol.nc
+!    I_a Read the "run.def"
+!    I_b Read the "start_evol.nc" and "startfi_evol.nc"
 !    I_c Subslope parametrisation
-!    I_d Read PCM data and convert to the physical grid
+!    I_d Read the PCM data and convert them to the physical grid
 !    I_e Initialization of the PEM variable and soil
 !    I_f Compute tendencies
 !    I_g Save initial PCM situation
-!    I_h Read the startpem.nc
+!    I_h Read the "startpem.nc"
 !    I_i Compute orbit criterion
 
@@ -22,51 +22,51 @@
 ! III Output
 !    III_a Update surface value for the PCM start files
-!    III_b Write restart_evol.nc and restartfi_evol.nc
-!    III_c Write restartpem.nc
+!    III_b Write the "restart_evol.nc" and "restartfi_evol.nc"
+!    III_c Write the "restartpem.nc"
 !------------------------
 
 PROGRAM pem
 
-use phyetat0_mod,                only: phyetat0
-use phyredem,                    only: physdem0, physdem1
-use netcdf,                      only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
-use turb_mod,                    only: q2, wstar
-use comslope_mod,                only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
-use logic_mod,                   only: iflag_phys
-use mod_const_mpi,               only: COMM_LMDZ
+use phyetat0_mod,               only: phyetat0
+use phyredem,                   only: physdem0, physdem1
+use netcdf,                     only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
+use turb_mod,                   only: q2, wstar
+use comslope_mod,               only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
+use logic_mod,                  only: iflag_phys
+use mod_const_mpi,              only: COMM_LMDZ
 use infotrac
-use geometry_mod,                only: latitude_deg
-use conf_pem_mod,                only: conf_pem
-use pemredem,                    only: pemdem0, pemdem1
-use glaciers_mod,                only: flow_co2glaciers, flow_h2oglaciers, co2_ice_flow, h2o_ice_flow
-use stopping_crit_mod,           only: stopping_crit_h2o_ice, stopping_crit_co2
-use constants_marspem_mod,       only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
-                                       m_noco2, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold
-use evol_ice_mod,                only: evol_co2_ice, evol_h2o_ice
-use comsoil_h_PEM,               only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
-                                       TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
-                                       tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
-                                       fluxgeo                             ! Geothermal flux for the PEM and PCM
-use adsorption_mod,              only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
-                                       ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
-                                       co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
-use time_evol_mod,               only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
-use orbit_param_criterion_mod,   only: orbit_param_criterion
-use recomp_orb_param_mod,        only: recomp_orb_param
-use ice_table_mod,               only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
-                                       ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
-use soil_thermalproperties_mod,  only: update_soil_thermalproperties
-use time_phylmdz_mod,            only: daysec, dtphys, day_end
-use abort_pem_mod,               only: abort_pem
-use soil_settings_PEM_mod,       only: soil_settings_PEM
-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
-use read_data_PCM_mod,           only: read_data_PCM
-use recomp_tend_co2_slope_mod,   only: recomp_tend_co2_slope
-use soil_pem_compute_mod,        only: soil_pem_compute
-use writediagpem_mod,            only: writediagpem
+use geometry_mod,               only: latitude_deg
+use conf_pem_mod,               only: conf_pem
+use pemredem,                   only: pemdem0, pemdem1
+use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
+                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice
+use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2
+use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
+use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
+use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
+                                      TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
+                                      tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
+                                      fluxgeo                             ! Geothermal flux for the PEM and PCM
+use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
+                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
+                                      co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
+use time_evol_mod,              only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
+use orbit_param_criterion_mod,  only: orbit_param_criterion
+use recomp_orb_param_mod,       only: recomp_orb_param
+use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
+                                      ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
+use soil_thermalproperties_mod, only: update_soil_thermalproperties
+use time_phylmdz_mod,           only: daysec, dtphys, day_end
+use abort_pem_mod,              only: abort_pem
+use soil_settings_PEM_mod,      only: soil_settings_PEM
+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
+use read_data_PCM_mod,          only: read_data_PCM
+use recomp_tend_co2_slope_mod,  only: recomp_tend_co2_slope
+use soil_pem_compute_mod,       only: soil_pem_compute
+use writediagpem_mod,           only: writediagpem
 
 #ifndef CPP_STD
@@ -291,5 +291,5 @@
 !------------------------
 ! I   Initialization
-!    I_a Read run.def
+!    I_a Read the "run.def"
 !------------------------
 #ifndef CPP_1D
@@ -328,7 +328,7 @@
 !------------------------
 ! I   Initialization
-!    I_b Read of start_evol.nc and starfi_evol.nc
-!------------------------
-! I_b.1 Read start_evol.nc
+!    I_b Read of the "start_evol.nc" and starfi_evol.nc
+!------------------------
+! I_b.1 Read "start_evol.nc"
 allocate(ps_start_PCM(ngrid))
 #ifndef CPP_1D
@@ -355,5 +355,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
+! Here we simply read them in the "startfi_evol.nc" file
 status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
 
@@ -372,5 +372,5 @@
 status = nf90_close(ncid)
 
-! I_b.2 Read startfi_evol.nc
+! I_b.2 Read the "startfi_evol.nc"
 ! First we read the initial state (starfi.nc)
 #ifndef CPP_STD
@@ -468,5 +468,5 @@
 !------------------------
 ! I   Initialization
-!    I_d Read PCM data and convert to the physical grid
+!    I_d Read the PCM data and convert them 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
@@ -579,5 +579,5 @@
 !------------------------
 ! I   Initialization
-!    I_h Read the startpem.nc
+!    I_h Read the "startpem.nc"
 !------------------------
 allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
@@ -655,5 +655,5 @@
     write(*,*) 'Global average pressure new time step',global_avg_press_new
 
-! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
+! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
     allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
     write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
@@ -725,7 +725,7 @@
 !    II_c Flow of glaciers
 !------------------------
-    if (co2_ice_flow) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
+    if (co2ice_flow) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
                                             global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
-    if (h2o_ice_flow) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
+    if (h2oice_flow) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
 
 !------------------------
@@ -944,5 +944,5 @@
 do ig = 1,ngrid
     ! H2O ice metamorphism
-    if (sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
+    if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
         h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
         qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
@@ -961,5 +961,5 @@
 
     ! CO2 ice metamorphism
-    if (sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
+    if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
         perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
         qsurf(ig,igcm_co2,:) = metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
@@ -1022,7 +1022,7 @@
 !------------------------
 ! III Output
-!    III_b Write restart_evol.nc and restartfi_evol.nc
-!------------------------
-! III_b.1 Write restart_evol.nc
+!    III_b Write "restart_evol.nc" and "restartfi_evol.nc"
+!------------------------
+! III_b.1 Write "restart_evol.nc"
 ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
 pday = day_ini
@@ -1041,5 +1041,5 @@
 #endif
 
-! III_b.2 Write restartfi_evol.nc
+! III_b.2 Write the "restartfi_evol.nc"
 #ifndef CPP_STD
     call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
@@ -1063,5 +1063,5 @@
 !------------------------
 ! III Output
-!    III_c Write restartpem.nc
+!    III_c Write the "restartpem.nc"
 !------------------------
 call pemdem0("restartpem.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3159)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3161)
@@ -13,5 +13,5 @@
 
 use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
-use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, index_breccia, index_bedrock
+use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
 use comsoil_h,                  only: volcapa, inertiedat
 use adsorption_mod,             only: regolith_adsorption, adsorption_pem
@@ -23,10 +23,11 @@
 use soil_pem_compute_mod,       only: soil_pem_compute
 use soil_pem_ini_mod,           only: soil_pem_ini
+use comslope_mod,               only: def_slope_mean, subslope_dist
 
 #ifndef CPP_STD
-    use comcstfi_h,   only: r, mugaz
-    use surfdat_h,    only: watercaptag, perennial_co2ice
+    use comcstfi_h,   only: r, mugaz, pi
+    use surfdat_h,    only: watercaptag, watercap, perennial_co2ice
 #else
-    use comcstfi_mod, only: r, mugaz
+    use comcstfi_mod, only: r, mugaz, pi
 #endif
 
@@ -120,8 +121,11 @@
         write(*,*)'Pemetat0: failed loading <h2o_ice>'
         write(*,*)'will reconstruct the values from watercaptag'
-        write(*,*)'with default value ''ini_h2o_bigreservoir'''
+        write(*,*)'with default value ''ini_huge_h2oice'''
         do ig = 1,ngrid
-            if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
-        enddo
+            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
+        enddo
+    else
+        ! The variations of infinite reservoirs during the PCM years are taken into account
+        h2o_ice = h2o_ice + watercap
     endif
 
@@ -319,7 +323,7 @@
     ! h2o ice
     h2o_ice = 0.
-    write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_h2o_bigreservoir''.'
+    write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_huge_h2oice''.'
         do ig = 1,ngrid
-            if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
+            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
         enddo
 
