Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3027)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3028)
@@ -1,11 +1,9 @@
-
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
+!------------------------
+! I   Initialization
+!    I_a READ run.def
 !    I_b READ of start_evol.nc and starfi_evol.nc
 !    I_c Subslope parametrisation
-!    I_d READ GCM data and convert to the physical grid 
-!    I_e Initialisation of the PEM variable and soil 
+!    I_d READ GCM data and convert to the physical grid
+!    I_e Initialization of the PEM variable and soil
 !    I_f Compute tendencies & Save initial situation
 !    I_g Save initial PCM situation
@@ -14,1404 +12,1158 @@
 
 ! II  Run
-!    II_a update pressure, ice and tracers 
+!    II_a Update pressure, ice and tracers
 !    II_b Evolution of the ice
-!    II_c CO2 glaciers flows 
+!    II_c CO2 & H2O glaciers flows
 !    II_d Update surface and soil temperatures
-!    II_e Update the tendencies 
+!    II_e Update the tendencies
 !    II_f Checking the stopping criterion
 
 ! III Output
 !    III_a Update surface value for the PCM start files
-!    III_b Write start and starfi.nc
-!    III_c Write start_pem
-
+!    III_b Write restart_evol.nc and restartfi_evol.nc
+!    III_c Write restartfi_PEM.nc
 !------------------------
 
 PROGRAM pem
 
-! 1: Modules needed for reading and writting startfi:
-      use phyetat0_mod, only: phyetat0
-      use phyredem,     only: physdem0, physdem1
-      use netcdf,       only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
-                         nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
-                         nf90_inquire_dimension,nf90_close
-      ! netcdf is needed to read info like lat and lon in the physiq file
-      use turb_mod,     only: q2, wstar
-! 1a: Modules specific from the Marsian physiq
+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 comconst_mod,               only: rad, g, cpp, pi
+use infotrac
+use geometry_mod,               only: latitude_deg
+use conf_pem_mod,               only: conf_pem
+use pemredem,                   only: pemdem0, pemdem1
+use glaciers_mod,               only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
+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_frost2perenial, threshold_co2_frost2perenial
+use evol_co2_ice_s_mod,         only: evol_co2_ice_s
+use evol_h2o_ice_s_mod,         only: evol_h2o_ice_s
+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 GCM
+                                      water_reservoir                     ! Water ressources
+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 temps_mod_evol,             only: dt_pem, evol_orbit_pem, Max_iter_pem
+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
+use soil_thermalproperties_mod, only: update_soil_thermalproperties
+
 #ifndef CPP_STD
-      use comsoil_h,    only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa,inertiesoil
-      use surfdat_h,    only: tsurf, emis,&
-                           qsurf,watercap, ini_surfdat_h, &
-                           albedodat, zmea, zstd, zsig, zgam, zthe, &
-                           hmons, summit, base,albedo_h2o_frost, &
-                           frost_albedo_threshold,emissiv,watercaptag,perenial_co2ice
-      use dimradmars_mod, only: totcloudfrac, albedo
-      use dust_param_mod, only: tauscaling
-      use tracer_mod,     only: noms,igcm_h2o_ice,igcm_co2,mmol,igcm_h2o_vap ! tracer names and molar masses
-#else 
-! 1b: Modules specific from the Generic physiq
-      use comsoil_h,    only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
-      use surfdat_h,    only: albedodat, zmea, zstd, zsig, zgam, zthe, &
-                           emissiv
-      use tracer_h,     only: noms,igcm_h2o_ice,igcm_co2 ! tracer names
-      use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, &
-                                    PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE,          &
-                                    ALBEDO_BAREGROUND,ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
-      use aerosol_mod,  only : iniaerosol
+    use comsoil_h,          only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil
+    use surfdat_h,          only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
+                                  albedodat, zmea, zstd, zsig, zgam, zthe,     &
+                                  hmons, summit, base,albedo_h2o_frost,        &
+                                  frost_albedo_threshold, emissiv, watercaptag, perenial_co2ice
+    use dimradmars_mod,     only: totcloudfrac, albedo
+    use dust_param_mod,     only: tauscaling
+    use tracer_mod,         only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
+    use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
+    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit
+    use comcstfi_h,         only: r, mugaz
+#else
+    use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
+    use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT,   &
+                                  PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
+                                  ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
+    use aerosol_mod,  only: iniaerosol
+    use planete_mod,  only: apoastr, periastr, year_day, peri_day, obliquit
+    use comcstfi_mod, only: r, mugaz
 #endif
-! 1c: Modules specific from the 1d marsian physiq
+
 #ifndef CPP_1D
-      USE iniphysiq_mod, ONLY: iniphysiq
-      USE control_mod,   ONLY: iphysiq, day_step,nsplit_phys
+    use iniphysiq_mod,            only: iniphysiq
+    use control_mod,              only: iphysiq, day_step, nsplit_phys
 #else
-      use time_phylmdz_mod, only: daysec, iphysiq,day_step
+    use time_phylmdz_mod,         only: daysec, iphysiq, day_step, dtphys
 #endif
 
+#ifdef CPP_1D
+    use regular_lonlat_mod,       only: init_regular_lonlat
+    use physics_distribution_mod, only: init_physics_distribution
+    use mod_grid_phy_lmdz,        only: regular_lonlat
+    use init_phys_1d_mod,         only: init_phys_1d
+#endif
+
+IMPLICIT NONE
+
+include "dimensions.h"
+include "paramet.h"
+include "comgeom.h"
+include "iniprint.h"
+
+integer ngridmx
+parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
+
+! Same variable names as in the GCM
+integer :: ngrid     ! Number of physical grid points
+integer :: nlayer    ! Number of vertical layer
+integer :: nq        ! Number of tracer
+integer :: day_ini   ! First day of the simulation
+real    :: pday      ! Physical day
+real    :: time_phys ! Same as GCM
+real    :: ptimestep ! Same as GCM
+real    :: ztime_fin ! Same as GCM
+
+! Variables to read start.nc
+character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM
+
+! Dynamic variables
+real                                :: vcov(ip1jm,llm), ucov(ip1jmp1,llm) ! vents covariants
+real                                :: teta(ip1jmp1,llm)                  ! temperature potentielle 
+real, dimension(:,:,:), allocatable :: q                                  ! champs advectes
+real                                :: ps(ip1jmp1)                        ! 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                                :: masse(ip1jmp1,llm)                 ! masse d'air
+real                                :: phis(ip1jmp1)                      ! geopotentiel au sol
+real                                :: time_0
+
+! 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
+
+! 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 :: 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
+real, dimension(:), allocatable :: cell_area     ! Cell_area read in FILE_NAME and written in restartfi 
+real                            :: Total_surface ! Total surface of the planet
+
+! Variables for h2o_ice evolution
+real                                :: ini_surf_h2o          ! Initial surface of sublimating h2o ice
+real                                :: ini_surf_co2          ! Initial surface of sublimating co2 ice
+real                                :: global_ave_press_GCM  ! constant: global average pressure retrieved in the GCM [Pa]
+real                                :: global_ave_press_old  ! constant: Global average pressure of initial/previous time step
+real                                :: global_ave_press_new  ! constant: Global average pressure of current time step
+real, dimension(:,:),   allocatable ::  zplev_new            ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
+real, dimension(:,:),   allocatable ::  zplev_gcm            ! same but retrieved from the gcm [kg/m^2]
+real, dimension(:,:,:), allocatable ::  zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
+real, dimension(:,:,:), allocatable :: zplev_old_timeseries  ! same but with the time series, for oldest time step
+logical                             :: STOPPING_water        ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
+logical                             :: STOPPING_1_water      ! Logical : is there still water ice to sublimate?
+logical                             :: STOPPING_co2          ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
+logical                             :: STOPPING_pressure     ! Logical : is the criterion (% of change in the surface pressure) reached?
+integer                             :: criterion_stop        ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
+real, save                          :: A , B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
+real, allocatable                   :: vmr_co2_gcm(:,:)      ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
+real, allocatable                   :: vmr_co2_pem_phys(:,:) ! Physics x Times  co2 volume mixing ratio used in the PEM
+real, allocatable                   :: q_co2_PEM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]
+real, allocatable                   :: q_h2o_PEM_phys(:,:)   ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
+integer                             :: timelen               ! # time samples
+real                                :: ave                   ! intermediate varibale to compute average
+real, allocatable                   :: p(:,:)                ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
+real                                :: extra_mass            ! Intermediate variables Extra mass of a tracer if it is greater than 1
+
+! Variables for slopes
+real, dimension(:,:),   allocatable :: min_co2_ice_1          ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
+real, dimension(:,:),   allocatable :: min_co2_ice_2          ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
+real, dimension(:,:),   allocatable :: min_h2o_ice_1          ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
+real, dimension(:,:),   allocatable :: min_h2o_ice_2          ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
+real, dimension(:,:,:), allocatable :: co2_ice_GCM            ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
+real, dimension(:,:),   allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice
+real, dimension(:,:),   allocatable :: initial_h2o_ice        ! physical point field : Logical array indicating if there is water ice at initial state
+real, dimension(:,:),   allocatable :: initial_co2_ice        ! physical point field : Logical array indicating if there is co2 ice at initial state
+real, dimension(:,:),   allocatable :: tendencies_co2_ice     ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
+real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
+real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
+real, dimension(:,:),   allocatable :: flag_co2flow(:,:)      !(ngrid,nslope): Flag where there is a CO2 glacier flow
+real, dimension(:),     allocatable :: flag_co2flow_mesh(:)   !(ngrid)       : Flag where there is a CO2 glacier flow
+real, dimension(:,:),   allocatable :: flag_h2oflow(:,:)      !(ngrid,nslope): Flag where there is a H2O glacier flow
+real, dimension(:),     allocatable :: flag_h2oflow_mesh(:)   !(ngrid)       : Flag where there is a H2O glacier flow
+
+! Variables for surface and soil
+real, allocatable :: tsurf_ave(:,:)              ! Physic x SLOPE field : Averaged Surface Temperature [K]
+real, allocatable :: tsoil_ave(:,:,:)            ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
+real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
+real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
+real, allocatable :: tsoil_GCM_timeseries(:,:,:,:)      !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
+real, allocatable :: tsurf_ave_yr1(:,:)                 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
+real, allocatable :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
+real, allocatable :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
+real, allocatable :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
+real, allocatable :: watersoil_density_timeseries(:,:,:,:)     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
+real, allocatable :: watersurf_density_ave(:,:)                ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
+real, allocatable :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
+real, allocatable :: watersoil_density_PEM_ave(:,:,:)          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
+real, allocatable :: Tsurfave_before_saved(:,:)                ! Surface temperature saved from previous time step [K]
+real, allocatable :: delta_co2_adsorbded(:)                    ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
+real, allocatable :: delta_h2o_adsorbded(:)                    ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
+real              :: totmassco2_adsorbded                      ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
+real              :: totmassh2o_adsorbded                      ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
+logical           :: bool_sublim                               ! logical to check if there is sublimation or not
+
+! Some variables for the PEM run
+real, parameter :: year_step = 1 ! timestep for the pem
+integer         :: year_iter     ! number of iteration
+integer         :: year_iter_max ! maximum number of iterations before stopping
+real            :: timestep      ! timestep [s]
+real            :: watercap_sum  ! total mass of water cap [kg/m^2] 
+real            :: water_sum     ! total mass of water in the mesh [kg/m^2] 
+
+#ifdef CPP_STD
+    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,    allocatable :: tsurf_read_generic(:)         ! Temporary variable to do the subslope transfert dimensiion when reading form generic
+    real,    allocatable :: qsurf_read_generic(:,:)       ! Temporary variable to do the subslope transfert dimensiion when reading form generic 
+    real,    allocatable :: tsoil_read_generic(:,:)       ! Temporary variable to do the subslope transfert dimensiion when reading form generic
+    real,    allocatable :: emis_read_generic(:)          ! Temporary variable to do the subslope transfert dimensiion when reading form generic
+    real,    allocatable :: albedo_read_generic(:,:)      ! Temporary variable to do the subslope transfert dimensiion when reading form generic
+    real,    allocatable :: tsurf(:,:)                    ! Subslope variable, only needed in the GENERIC case
+    real,    allocatable :: qsurf(:,:,:)                  ! Subslope variable, only needed in the GENERIC case
+    real,    allocatable :: tsoil(:,:,:)                  ! Subslope variable, only needed in the GENERIC case
+    real,    allocatable :: emis(:,:)                     ! Subslope variable, only needed in the GENERIC case
+    real,    allocatable :: watercap(:,:)                 ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
+    logical, allocatable :: WATERCAPTAG(:)                ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
+    real,    allocatable :: albedo(:,:,:)                 ! Subslope variable, only needed in the GENERIC case
+    real,    allocatable :: inertiesoil(:,:,:)            ! Subslope variable, only needed in the GENERIC case
+#endif
+
+#ifdef CPP_1D
+    integer            :: nsplit_phys
+    integer, parameter :: jjm_value = jjm - 1
+    integer            :: ierr
+#else
+    integer, parameter :: jjm_value = jjm
+#endif
+
+! Loop variables
+integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil
+
+! Parallel variables
 #ifndef CPP_STD
-      use mod_phys_lmdz_para, only: is_parallel, is_sequential, &
-                                   is_mpi_root, is_omp_root,    &
-                                   is_master
-      use planete_h,     only: aphelie, periheli, year_day, peri_day, &
-                          obliquit
+    is_sequential = .true.
+    is_parallel = .false.
+    is_mpi_root = .true.
+    is_omp_root = .true.
+    is_master = .true.
+#endif
+
+day_ini = 0    ! test
+time_phys = 0. ! test
+
+! Some constants
+ngrid = ngridmx
+nlayer = llm
+
+A = (1/m_co2 - 1/m_noco2)
+B = 1/m_noco2
+
+year_day = 669
+daysec = 88775.
+timestep = year_day*daysec/year_step
+
+!----------------------------- INITIALIZATION --------------------------
+!------------------------
+! I   Initialization
+!    I_a READ run.def
+!------------------------
+#ifndef CPP_1D
+    dtphys = 0
+    call conf_gcm(99,.true.)
+    call infotrac_init
+    nq = nqtot
+    allocate(q(ip1jmp1,llm,nqtot))
 #else
-      use planete_mod,   only: apoastr, periastr, year_day, peri_day, &
-                          obliquit
-#endif
-
-      USE comslope_mod,  ONLY: nslope,def_slope,def_slope_mean, &
-                           subslope_dist,iflat,                &
-                           major_slope,ini_comslope_h
-
-#ifndef CPP_STD
-      USE comcstfi_h, only: r, mugaz
-#else
-      USE comcstfi_mod, only: r, mugaz
-#endif
-
-      USE logic_mod,        ONLY: iflag_phys
-      USE mod_const_mpi,    ONLY: COMM_LMDZ
-      use time_phylmdz_mod, only: daysec,dtphys
-      USE comconst_mod,     ONLY: rad,g,cpp,pi
-      USE infotrac
-      USE geometry_mod,     only: latitude_deg 
-
-      use conf_pem_mod,     only: conf_pem
-      use pemredem,         only:  pemdem0,pemdem1
-      use glaciers_mod,        only: co2glaciers_evol,h2oglaciers_evol,co2glaciersflow,h2oglaciersflow
-      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_frost2perenial,threshold_co2_frost2perenial
-      use evol_co2_ice_s_mod, only: evol_co2_ice_s
-      use evol_h2o_ice_s_mod, only: evol_h2o_ice_s
-      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 GCM
-                               water_reservoir                                         ! Water ressources
-      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 temps_mod_evol,            ONLY: dt_pem, evol_orbit_pem, Max_iter_pem
-      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
-      use soil_thermalproperties_mod, only: update_soil_thermalproperties
-
-#ifdef CPP_1D
-      use regular_lonlat_mod,       only: init_regular_lonlat
-      use physics_distribution_mod, only: init_physics_distribution
-      use mod_grid_phy_lmdz,        only : regular_lonlat
-      use init_phys_1d_mod,         only : init_phys_1d
-#endif
-
-
-  IMPLICIT NONE
-
-  include "dimensions.h"
-  include "paramet.h"
-
-  INTEGER ngridmx
-  PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
-
-  include "comgeom.h"
-  include "iniprint.h"
-! Same variable's name as in the GCM
-
-  INTEGER :: ngrid      !Number of physical grid points
-  INTEGER :: nlayer     !Number of vertical layer
-  INTEGER :: nq         !Number of tracer
-  INTEGER :: day_ini    !First day of the simulation
-  REAL :: pday          !Physical day
-  REAL :: time_phys     !Same as GCM
-  REAL :: ptimestep     !Same as GCM
-  REAL :: ztime_fin     !Same as GCM
-
-! Variable for reading start.nc
-      character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
-  !   variables dynamiques
-  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
-  REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
-  REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
-  REAL ps(ip1jmp1)                       ! 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 masse(ip1jmp1,llm)                ! masse d'air
-  REAL phis(ip1jmp1)                     ! geopotentiel au sol
-  REAL time_0
-
-! Variable for reading 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
-
-! Variable for reading starfi.nc and writting 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 :: 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
-      REAL, dimension(:),allocatable  :: cell_area                 !Cell_area read in FILE_NAME and written in restartfi 
-      REAL :: Total_surface                                        !Total surface of the planet
-
-! Variable for h2o_ice evolution
-      REAL :: ini_surf_h2o                                         ! Initial surface of sublimating h2o ice
-      REAL :: ini_surf_co2                                         ! Initial surface of sublimating co2 ice
-
-      REAL :: global_ave_press_GCM           ! constant: global average pressure retrieved in the GCM [Pa]
-      REAL :: global_ave_press_old           ! constant: Global average pressure of initial/previous time step
-      REAL :: global_ave_press_new           ! constant: Global average pressure of current time step
-
-      REAL , dimension(:,:), allocatable ::  zplev_new ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
-      REAL , dimension(:,:), allocatable ::  zplev_gcm ! same but retrieved from the gcm [kg/m^2]
-      REAL , dimension(:,:,:), allocatable ::  zplev_new_timeseries  ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
-      REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries   ! same but with the time series, for oldest time step
-
-      LOGICAL :: STOPPING_water   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
-      LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
-      LOGICAL :: STOPPING_co2     ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
-      LOGICAL :: STOPPING_pressure! Logical : is the criterion (% of change in the surface pressure) reached?
-      INTEGER :: criterion_stop   ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
-
-      real,save ::     A , B, mmean              ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
-      real ,allocatable :: vmr_co2_gcm(:,:)      ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
-      real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times  co2 volume mixing ratio used in the PEM
-      real ,allocatable :: q_co2_PEM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]
-      REAL ,allocatable :: q_h2o_PEM_phys(:,:)   ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
-      integer :: timelen                         ! # time samples
-      REAL :: ave                                ! intermediate varibale to compute average
-
-      REAL, ALLOCATABLE :: p(:,:)  ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
-      REAL :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
-
-
-!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
-      REAL, dimension(:,:),  allocatable :: min_co2_ice_1           ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
-      REAL, dimension(:,:),  allocatable :: min_co2_ice_2           ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
-      REAL, dimension(:,:),  allocatable :: min_h2o_ice_1           ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
-      REAL, dimension(:,:),  allocatable :: min_h2o_ice_2           ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
-      REAL, dimension(:,:,:),allocatable :: co2_ice_GCM             ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
-      REAL, dimension(:,:),  allocatable :: initial_co2_ice_sublim  ! physical point field : Logical array indicating sublimating point of co2 ice
-      REAL, dimension(:,:),  allocatable :: initial_h2o_ice         ! physical point field : Logical array indicating if there is water ice at initial state
-      REAL, dimension(:,:),  allocatable :: initial_co2_ice         ! physical point field : Logical array indicating if there is co2 ice at initial state
-      REAL, dimension(:,:),  allocatable :: tendencies_co2_ice      ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
-      REAL, dimension(:,:),  allocatable :: tendencies_co2_ice_ini  ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
-      REAL, dimension(:,:),  allocatable :: tendencies_h2o_ice      ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
-      REAL, dimension(:,:),  allocatable :: flag_co2flow(:,:)       !(ngrid,nslope): Flag where there is a CO2 glacier flow
-      REAL, dimension(:),    allocatable :: flag_co2flow_mesh(:)    !(ngrid)       : Flag where there is a CO2 glacier flow
-      REAL, dimension(:,:),  allocatable :: flag_h2oflow(:,:)       !(ngrid,nslope): Flag where there is a H2O glacier flow
-      REAL, dimension(:),    allocatable :: flag_h2oflow_mesh(:)    !(ngrid)       : Flag where there is a H2O glacier flow
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL 
-
-     REAL, ALLOCATABLE :: tsurf_ave(:,:)       ! Physic x SLOPE field : Averaged Surface Temperature [K]
-     REAL, ALLOCATABLE :: tsoil_ave(:,:,:)     ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
-     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:)    ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
-
-     REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-     REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:)      !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-     REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:)                 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
-
-     REAL,ALLOCATABLE  :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
-     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
-     REAL,ALLOCATABLE  :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
-     REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:)     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
-     REAL,ALLOCATABLE  :: watersurf_density_ave(:,:)                ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
-     REAL,ALLOCATABLE  :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
-     REAL,ALLOCATABLE  :: watersoil_density_PEM_ave(:,:,:)          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
-     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)                ! Surface temperature saved from previous time step [K]
-     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)                    ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
-     REAL, ALLOCATABLE :: delta_h2o_adsorbded(:)                    ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
-     REAL :: totmassco2_adsorbded                                   ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
-     REAL :: totmassh2o_adsorbded                                   ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
-     LOGICAL :: bool_sublim                                         ! logical to check if there is sublimation or not
-      
-
-!! Some parameters for the PEM run
-    REAL, PARAMETER :: year_step = 1     ! timestep for the pem 
-    INTEGER ::  year_iter                ! number of iteration
-    INTEGER :: year_iter_max             ! maximum number of iterations before stopping
-    REAL :: timestep                     ! timestep [s]
-    REAL :: watercap_sum                 ! total mass of water cap [kg/m^2] 
-    REAL :: water_sum                    ! total mass of water in the mesh [kg/m^2] 
-
-#ifdef CPP_STD
-     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,allocatable :: tsurf_read_generic(:)       ! Temporary variable to do the subslope transfert dimensiion when reading form generic
-     REAL,allocatable :: qsurf_read_generic(:,:)     ! Temporary variable to do the subslope transfert dimensiion when reading form generic 
-     REAL,allocatable :: tsoil_read_generic(:,:)     ! Temporary variable to do the subslope transfert dimensiion when reading form generic
-     REAL,allocatable :: emis_read_generic(:)        ! Temporary variable to do the subslope transfert dimensiion when reading form generic
-     REAL,allocatable :: albedo_read_generic(:,:)    ! Temporary variable to do the subslope transfert dimensiion when reading form generic
-     REAL,allocatable :: tsurf(:,:)                  ! Subslope variable, only needed in the GENERIC case
-     REAL,allocatable :: qsurf(:,:,:)                ! Subslope variable, only needed in the GENERIC case
-     REAL,allocatable :: tsoil(:,:,:)                ! Subslope variable, only needed in the GENERIC case
-     REAL,allocatable :: emis(:,:)                   ! Subslope variable, only needed in the GENERIC case
-     REAL,allocatable :: watercap(:,:)               ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
-     LOGICAL, allocatable :: WATERCAPTAG(:)          ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
-     REAL,allocatable :: albedo(:,:,:)               ! Subslope variable, only needed in the GENERIC case
-     REAL,allocatable :: inertiesoil(:,:,:)          ! Subslope variable, only needed in the GENERIC case
-#endif
-
-#ifdef CPP_1D
-  INTEGER :: nsplit_phys
-  integer,parameter:: jjm_value=jjm-1
-  integer :: ierr
-#else
-  integer,parameter:: jjm_value=jjm
-#endif
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! Loop variable
-     INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil 
-#ifndef CPP_STD
-! Parallel variables
-      is_sequential=.true.
-      is_parallel=.false.
-      is_mpi_root=.true.
-      is_omp_root=.true.
-      is_master=.true.
-#endif
-
-      day_ini=0    !test
-      time_phys=0. !test
-
-! Some constants
-
-      ngrid=ngridmx
-      nlayer=llm
-   
-      A =(1/m_co2 - 1/m_noco2)
-      B=1/m_noco2
-
-      year_day=669
-      daysec=88775.
-      timestep=year_day*daysec/year_step
-
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
-
-!------------------------
-
-!----------------------------READ run.def ---------------------
-
-#ifndef CPP_1D
-      dtphys=0
-      CALL conf_gcm( 99, .TRUE. )
-      call infotrac_init
-      nq=nqtot
-     allocate(q(ip1jmp1,llm,nqtot))
-#else
-      ! load tracer names from file 'traceur.def'
-        open(90,file='traceur.def',status='old',form='formatted',&
-            iostat=ierr)
-        if (ierr.ne.0) then
-          write(*,*) 'Cannot find required file "traceur.def"'
-          write(*,*) ' If you want to run with tracers, I need it'
-          write(*,*) ' ... might as well stop here ...'
-          stop
-        else
-          write(*,*) "pem1d: Reading file traceur.def"
-          ! read number of tracers:
-          read(90,*,iostat=ierr) nq
-          print *, "nq",nq
-          nqtot=nq ! set value of nqtot (in infotrac module) as nq
-          if (ierr.ne.0) then
+    ! load tracer names from file 'traceur.def'
+    open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
+    if (ierr /= 0) then
+        write(*,*) 'Cannot find required file "traceur.def"'
+        write(*,*) ' If you want to run with tracers, I need it'
+        write(*,*) ' ... might as well stop here ...'
+        stop
+    else
+        write(*,*) "pem1d: Reading file traceur.def"
+        ! read number of tracers:
+        read(90,*,iostat = ierr) nq
+        write(*,*) "nq",nq
+        nqtot = nq ! set value of nqtot (in infotrac module) as nq
+        if (ierr /= 0) then
             write(*,*) "pem1d: error reading number of tracers"
             write(*,*) "   (first line of traceur.def) "
             stop
-          endif
-          if (nq<1) then
+        endif
+        if (nq < 1) then
             write(*,*) "pem1d: error number of tracers"
             write(*,*) "is nq=",nq," but must be >=1!"
             stop
-          endif
         endif
-     nq=nqtot
-     allocate(q(ip1jmp1,llm,nqtot))
-     allocate(ap(nlayer+1))
-     allocate(bp(nlayer+1))
-     call init_phys_1d(llm,nqtot,vcov,ucov, &
-                    teta,q,ps, time_0,ap,bp)   
-     pi=2.E+0*asin(1.E+0) 
-     g=3.72 
-     nsplit_phys=1
+    endif
+    nq = nqtot
+    allocate(q(ip1jmp1,llm,nqtot))
+    allocate(ap(nlayer + 1))
+    allocate(bp(nlayer + 1))
+    call init_phys_1d(llm,nqtot,vcov,ucov,teta,q,ps,time_0,ap,bp)
+    pi = 2.*asin(1.) 
+    g = 3.72 
+    nsplit_phys = 1
 #endif
-      CALL conf_pem     
- 
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
+
+call conf_pem
+
+!------------------------
+! I   Initialization
 !    I_b READ of start_evol.nc and starfi_evol.nc
-
-!------------------------
-
-!----------------------------Initialisation : READ some constant of startfi_evol.nc --------------------- 
-
-!----------------------------READ start.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)
-#endif
-
-     allocate(ps_start_GCM(ngrid))
-#ifndef CPP_1D
-     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
-
-     CALL iniconst !new
-     CALL inigeom
-
-     allocate(ap(nlayer+1))
-     allocate(bp(nlayer+1))
-     status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
-     status = nf90_inq_varid(ncid, "ap", apvarid)
-     status = nf90_get_var(ncid, apvarid, ap)
-     status = nf90_inq_varid(ncid, "bp", bpvarid)
-     status = nf90_get_var(ncid, bpvarid, bp)
-     status =nf90_close(ncid)
-
+    call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0)
+
+    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
+
+    call iniconst !new
+    call inigeom
+
+    allocate(ap(nlayer + 1))
+    allocate(bp(nlayer + 1))
+    status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid)
+    status = nf90_inq_varid(ncid,"ap",apvarid)
+    status = nf90_get_var(ncid,apvarid,ap)
+    status = nf90_inq_varid(ncid,"bp",bpvarid)
+    status = nf90_get_var(ncid,bpvarid,bp)
+    status = nf90_close(ncid)
+
+    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, &
+                   rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
 #else
-
-     ps_start_GCM(1)=ps(1)
-
-#endif
-
-
-#ifndef CPP_1D
-     CALL iniphysiq(iim,jjm,llm, &
-          (jjm-1)*iim+2,comm_lmdz, &
-          daysec,day_ini,dtphys/nsplit_phys, &
-          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
-          iflag_phys)
+    ps_start_GCM(1) = ps(1)
 #endif
 
 ! In the gcm, 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)
-
-      allocate(longitude(ngrid))
-      allocate(latitude(ngrid))
-      allocate(cell_area(ngrid))
-
-      status = nf90_inq_varid(ncid, "longitude", lonvarid)
-      status = nf90_get_var(ncid, lonvarid, longitude)
-
-      status = nf90_inq_varid(ncid, "latitude", latvarid)
-      status = nf90_get_var(ncid, latvarid, latitude)
-
-      status = nf90_inq_varid(ncid, "area", areavarid)
-      status = nf90_get_var(ncid, areavarid, cell_area)
-
-      status = nf90_inq_varid(ncid, "soildepth", sdvarid)
-      status = nf90_get_var(ncid, sdvarid, mlayer)
-
-      status =nf90_close(ncid)
-
-!----------------------------READ startfi.nc --------------------- 
-
+status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
+
+allocate(longitude(ngrid))
+allocate(latitude(ngrid))
+allocate(cell_area(ngrid))
+
+status = nf90_inq_varid(ncid,"longitude",lonvarid)
+status = nf90_get_var(ncid,lonvarid,longitude)
+
+status = nf90_inq_varid(ncid,"latitude",latvarid)
+status = nf90_get_var(ncid,latvarid,latitude)
+
+status = nf90_inq_varid(ncid,"area",areavarid)
+status = nf90_get_var(ncid,areavarid,cell_area)
+
+status = nf90_inq_varid(ncid,"soildepth",sdvarid)
+status = nf90_get_var(ncid,sdvarid,mlayer)
+
+status = nf90_close(ncid)
+
+! 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,   &
-              day_ini,time_phys,         &
-              tsurf,tsoil,albedo,emis,   &
-              q2,qsurf,tauscaling,totcloudfrac,wstar,     &
-              watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
+    call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,day_ini,time_phys,tsurf, &
+                  tsoil,albedo,emis,q2,qsurf,tauscaling,totcloudfrac,wstar,      &
+                  watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
 
 ! Remove unphysical values of surface tracer
-     DO i=1,ngrid
-       DO nnq=1,nqtot
-         DO islope=1,nslope
-           if(qsurf(i,nnq,islope).LT.0) then
-             qsurf(i,nnq,islope)=0.
-           endif
-         enddo
-       enddo
-     enddo
-   
-     call surfini(ngrid,qsurf)
-
+    do i = 1,ngrid
+        do nnq = 1,nqtot
+            do islope = 1,nslope
+                if (qsurf(i,nnq,islope) < 0) qsurf(i,nnq,islope) = 0.
+            enddo
+        enddo
+    enddo
+
+    call surfini(ngrid,qsurf)
 #else
-        call phys_state_var_init(nq)
-         IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
-         call initracer(ngrid,nq)
-         call iniaerosol()
-         allocate(tsurf_read_generic(ngrid))
-         allocate(qsurf_read_generic(ngrid,nq))
-         allocate(tsoil_read_generic(ngrid,nsoilmx))
-         allocate(emis_read_generic(ngrid))
-         allocate(tsurf(ngrid,1))
-         allocate(qsurf(ngrid,nq,1))
-         allocate(tsoil(ngrid,nsoilmx,1))
-         allocate(emis(ngrid,1))
-         allocate(watercap(ngrid,1))
-         allocate(watercaptag(ngrid))
-         allocate(albedo_read_generic(ngrid,2))
-         allocate(albedo(ngrid,2,1))
-         allocate(inertiesoil(ngrid,nsoilmx,1))
-         call phyetat0(.true.,                                 &
-                       ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,      &
-                       day_ini,time_phys,tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,qsurf_read_generic,   &
-                       cloudfrac,totcloudfrac,hice,                   &
-                       rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
-         call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 
-
-         nslope=1
-         call ini_comslope_h(ngrid,1)
-
-         qsurf(:,:,1)=qsurf_read_generic(:,:)
-         tsurf(:,1)=tsurf_read_generic(:)
-         tsoil(:,:,1)=tsoil_read_generic(:,:)
-         emis(:,1)=emis_read_generic(:)
-         watercap(:,1)=0.
-         watercaptag(:)=.false.
-         albedo(:,1,1)=albedo_read_generic(:,1)
-         albedo(:,2,1)=albedo_read_generic(:,2)
-         inertiesoil(:,:,1)=inertiedat(:,:)
-
-if (nslope.eq.1) then
- def_slope(1) = 0
- def_slope(2) = 0
- def_slope_mean=0
- subslope_dist(:,1) = 1.
-endif
+    call phys_state_var_init(nq)
+    if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
+    call initracer(ngrid,nq)
+    call iniaerosol()
+    allocate(tsurf_read_generic(ngrid))
+    allocate(qsurf_read_generic(ngrid,nq))
+    allocate(tsoil_read_generic(ngrid,nsoilmx))
+    allocate(emis_read_generic(ngrid))
+    allocate(tsurf(ngrid,1))
+    allocate(qsurf(ngrid,nq,1))
+    allocate(tsoil(ngrid,nsoilmx,1))
+    allocate(emis(ngrid,1))
+    allocate(watercap(ngrid,1))
+    allocate(watercaptag(ngrid))
+    allocate(albedo_read_generic(ngrid,2))
+    allocate(albedo(ngrid,2,1))
+    allocate(inertiesoil(ngrid,nsoilmx,1))
+    call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,day_ini,time_phys, &
+                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,     &
+                  qsurf_read_generic,cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, &
+                  tslab, tsea_ice,sea_ice)
+    call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 
+
+    nslope = 1
+    call ini_comslope_h(ngrid,1)
+
+    qsurf(:,:,1) = qsurf_read_generic(:,:)
+    tsurf(:,1) = tsurf_read_generic(:)
+    tsoil(:,:,1) = tsoil_read_generic(:,:)
+    emis(:,1) = emis_read_generic(:)
+    watercap(:,1) = 0.
+    watercaptag(:) = .false.
+    albedo(:,1,1) = albedo_read_generic(:,1)
+    albedo(:,2,1) = albedo_read_generic(:,2)
+    inertiesoil(:,:,1) = inertiedat(:,:)
+
+    if (nslope == 1) then
+        def_slope(1) = 0
+        def_slope(2) = 0
+        def_slope_mean = 0
+        subslope_dist(:,1) = 1.
+    endif
 
 ! Remove unphysical values of surface tracer
-     DO i=1,ngrid
-       DO nnq=1,nqtot
-         qsurf(i,nnq,1)=qsurf_read_generic(i,nnq)
-         if(qsurf(i,nnq,1).LT.0) then
-           qsurf(i,nnq,1)=0.
-         endif
-       enddo
-     enddo
+    do i = 1,ngrid
+        do nnq = 1,nqtot
+            qsurf(i,nnq,1) = qsurf_read_generic(i,nnq)
+            if (qsurf(i,nnq,1) < 0) qsurf(i,nnq,1) = 0.
+        enddo
+    enddo
 #endif
 
-     DO nnq=1,nqtot  ! Why not using ini_tracer ?
-       if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
-       if(noms(nnq).eq."h2o_vap") then
-          igcm_h2o_vap = nnq
-          mmol(igcm_h2o_vap)=18.
-       endif
-       if(noms(nnq).eq."co2") igcm_co2 = nnq
-     ENDDO 
-     r= 8.314511E+0 *1000.E+0/mugaz
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
-!    I_b READ of start_evol.nc and starfi_evol.nc
+do nnq = 1,nqtot  ! Why not using ini_tracer ?
+    if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
+    if (noms(nnq) == "h2o_vap") then
+        igcm_h2o_vap = nnq
+        mmol(igcm_h2o_vap)=18.
+    endif
+    if (noms(nnq) == "co2") igcm_co2 = nnq
+enddo 
+r = 8.314511E+0*1000.E+0/mugaz
+
+!------------------------
+! I   Initialization
 !    I_c Subslope parametrisation
-
-!------------------------
-
-!----------------------------Subslope parametrisation definition --------------------- 
-
-!     Define some slope statistics
-     iflat=1
-     DO islope=2,nslope
-       IF(abs(def_slope_mean(islope)).lt. &
-         abs(def_slope_mean(iflat))) THEN
-         iflat = islope
-       ENDIF     
-     ENDDO
-
-     PRINT*,'Flat slope for islope = ',iflat
-     PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
-
-     allocate(flag_co2flow(ngrid,nslope))
-     allocate(flag_co2flow_mesh(ngrid))
-     allocate(flag_h2oflow(ngrid,nslope))
-     allocate(flag_h2oflow_mesh(ngrid))
-
-     flag_co2flow(:,:) = 0     
-     flag_co2flow_mesh(:) = 0
-     flag_h2oflow(:,:) = 0     
-     flag_h2oflow_mesh(:) = 0
-!---------------------------- READ GCM data --------------------- 
-
-! I   Initialisation
-!    I_a READ run.def 
-!    I_b READ of start_evol.nc and starfi_evol.nc
-!    I_c Subslope parametrisation
-!    I_d READ GCM data
-
-!------------------------
-
+!------------------------
+! Define some slope statistics
+iflat = 1
+do islope = 2,nslope
+    if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
+enddo
+
+write(*,*) 'Flat slope for islope = ',iflat
+write(*,*) 'corresponding criterium = ',def_slope_mean(iflat)
+
+allocate(flag_co2flow(ngrid,nslope))
+allocate(flag_co2flow_mesh(ngrid))
+allocate(flag_h2oflow(ngrid,nslope))
+allocate(flag_h2oflow_mesh(ngrid))
+
+flag_co2flow(:,:) = 0     
+flag_co2flow_mesh(:) = 0
+flag_h2oflow(:,:) = 0     
+flag_h2oflow_mesh(:) = 0
+
+!------------------------
+! I   Initialization
+!    I_d READ GCM 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 GCM run, saving only the minimum value
-
-     call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
-
-     allocate(tsoil_ave(ngrid,nsoilmx,nslope))
-     allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
-
-     allocate(vmr_co2_gcm(ngrid,timelen))
-     allocate(ps_timeseries(ngrid,timelen))
-     allocate(min_co2_ice_1(ngrid,nslope))
-     allocate(min_h2o_ice_1(ngrid,nslope))
-     allocate(min_co2_ice_2(ngrid,nslope))
-     allocate(min_h2o_ice_2(ngrid,nslope))
-     allocate(tsurf_ave_yr1(ngrid,nslope))
-     allocate(tsurf_ave(ngrid,nslope))
-
-     allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
-     allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
-     allocate(q_co2_PEM_phys(ngrid,timelen))
-     allocate(q_h2o_PEM_phys(ngrid,timelen))
-     allocate(co2_ice_GCM(ngrid,nslope,timelen))
-     allocate(watersurf_density_ave(ngrid,nslope))
-     allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
-
-     allocate(Tsurfave_before_saved(ngrid,nslope))
-     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
-     allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
-     allocate(delta_co2_adsorbded(ngrid))
-     allocate(delta_h2o_adsorbded(ngrid))
-     allocate(vmr_co2_pem_phys(ngrid,timelen))
-
-     print *, "Downloading data Y1..."
-
-     call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1,&   
-                       tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &      
-                       watersurf_density_ave,watersoil_density_timeseries)
+call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
+
+allocate(tsoil_ave(ngrid,nsoilmx,nslope))
+allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
+allocate(vmr_co2_gcm(ngrid,timelen))
+allocate(ps_timeseries(ngrid,timelen))
+allocate(min_co2_ice_1(ngrid,nslope))
+allocate(min_h2o_ice_1(ngrid,nslope))
+allocate(min_co2_ice_2(ngrid,nslope))
+allocate(min_h2o_ice_2(ngrid,nslope))
+allocate(tsurf_ave_yr1(ngrid,nslope))
+allocate(tsurf_ave(ngrid,nslope))
+allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
+allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
+allocate(q_co2_PEM_phys(ngrid,timelen))
+allocate(q_h2o_PEM_phys(ngrid,timelen))
+allocate(co2_ice_GCM(ngrid,nslope,timelen))
+allocate(watersurf_density_ave(ngrid,nslope))
+allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
+allocate(Tsurfave_before_saved(ngrid,nslope))
+allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
+allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
+allocate(delta_co2_adsorbded(ngrid))
+allocate(delta_h2o_adsorbded(ngrid))
+allocate(vmr_co2_pem_phys(ngrid,timelen))
+
+write(*,*) "Downloading data Y1..."
+call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
+                   tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,           &
+                   co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
+write(*,*) "Downloading data Y1 done"
 
 ! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
-
-     print *, "Downloading data Y1 done"
-     print *, "Downloading data Y2"
-
-     call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
-                  tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &      
-                  watersurf_density_ave,watersoil_density_timeseries)
-
-     print *, "Downloading data Y2 done"
-
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
-!    I_b READ of start_evol.nc and starfi_evol.nc
-!    I_c Subslope parametrisation
-!    I_d READ GCM data and convert to the physical grid 
-!    I_e Initialisation of the PEM variable and soil 
-
-!------------------------
-
-!---------------------------- Initialisation of the PEM soil and values --------------------- 
-
-   call end_comsoil_h_PEM
-   call ini_comsoil_h_PEM(ngrid,nslope)
-   call end_adsorption_h_PEM
-   call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
-   call end_ice_table_porefilling
-   call ini_ice_table_porefilling(ngrid,nslope)
-
-   if(soil_pem) then
-      call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
-      DO l=1,nsoilmx
+write(*,*) "Downloading data Y2"
+call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
+                   tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,              &
+                   co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
+write(*,*) "Downloading data Y2 done"
+
+!------------------------
+! I   Initialization
+!    I_e Initialization of the PEM variables and soil
+!------------------------
+call end_comsoil_h_PEM
+call ini_comsoil_h_PEM(ngrid,nslope)
+call end_adsorption_h_PEM
+call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
+call end_ice_table_porefilling
+call ini_ice_table_porefilling(ngrid,nslope)
+
+if (soil_pem) then
+    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
+    do l=1,nsoilmx
         tsoil_PEM(:,l,:)=tsoil_ave(:,l,:)
         tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:)
         watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:)
-      ENDDO
-      DO l=nsoilmx+1,nsoilmx_PEM
+    enddo
+    do l=nsoilmx+1,nsoilmx_PEM
         tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:)
         watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:)
-      ENDDO
-      watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
-  endif !soil_pem
-      deallocate(tsoil_ave)
-      deallocate(tsoil_GCM_timeseries)
-
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
-!    I_b READ of start_evol.nc and starfi_evol.nc
-!    I_c Subslope parametrisation
-!    I_d READ GCM data and convert to the physical grid 
-!    I_e Initialisation of the PEM variable and soil 
+    enddo
+    watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
+endif !soil_pem
+deallocate(tsoil_ave,tsoil_GCM_timeseries)
+
+!------------------------
+! I   Initialization
 !    I_f Compute tendencies & Save initial situation
-
-!----- Compute tendencies from the PCM run
-
-     allocate(tendencies_co2_ice(ngrid,nslope))
-     allocate(tendencies_co2_ice_ini(ngrid,nslope))
-     allocate(tendencies_h2o_ice(ngrid,nslope))
-
-!  Compute the tendencies of the evolution of ice over the years
-
-      call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,&
-             min_co2_ice_2,tendencies_co2_ice)
-
-      tendencies_co2_ice_ini(:,:)=tendencies_co2_ice(:,:)
-
-      call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,&
-             min_h2o_ice_2,tendencies_h2o_ice)
-
-     deallocate(min_co2_ice_1)
-     deallocate(min_co2_ice_2)
-     deallocate(min_h2o_ice_1)
-     deallocate(min_h2o_ice_2)
-
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
-!    I_b READ of start_evol.nc and starfi_evol.nc
-!    I_c Subslope parametrisation
-!    I_d READ GCM data and convert to the physical grid 
-!    I_e Initialisation of the PEM variable and soil 
-!    I_f Compute tendencies & Save initial situation
+!------------------------
+allocate(tendencies_co2_ice(ngrid,nslope))
+allocate(tendencies_co2_ice_ini(ngrid,nslope))
+allocate(tendencies_h2o_ice(ngrid,nslope))
+
+! Compute the tendencies of the evolution of ice over the years
+call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice)
+tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:)
+call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice)
+
+deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2)
+
+!------------------------
+! I   Initialization
 !    I_g Save initial PCM situation
-
-!---------------------------- Save initial PCM situation --------------------- 
-
-     allocate(initial_co2_ice_sublim(ngrid,nslope))
-     allocate(initial_co2_ice(ngrid,nslope))
-     allocate(initial_h2o_ice(ngrid,nslope))
+!------------------------
+allocate(initial_co2_ice_sublim(ngrid,nslope))
+allocate(initial_co2_ice(ngrid,nslope))
+allocate(initial_h2o_ice(ngrid,nslope))
 
 ! We save the places where water ice is sublimating
 ! We compute the surface of water ice sublimating
-  ini_surf_co2=0.
-  ini_surf_h2o=0.
-  Total_surface=0.
-  do i=1,ngrid
-    Total_surface=Total_surface+cell_area(i)
-    do islope=1,nslope
-      if (tendencies_co2_ice(i,islope).LT.0) then
-         initial_co2_ice_sublim(i,islope)=1.
-         ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
-      else
-         initial_co2_ice_sublim(i,islope)=0.         
-      endif
-      if (qsurf(i,igcm_co2,islope).GT.0) then
-         initial_co2_ice(i,islope)=1.
-      else
-         initial_co2_ice(i,islope)=0.         
-      endif
-      if (tendencies_h2o_ice(i,islope).LT.0) then
-         initial_h2o_ice(i,islope)=1.
-         ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
-      else
-         initial_h2o_ice(i,islope)=0.         
-      endif
-    enddo
-  enddo
-
-     print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
-     print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
-     print *, "Total surface of the planet=", Total_surface
-    
-     allocate(zplev_gcm(ngrid,nlayer+1))
-
-     DO l=1,nlayer+1
-       DO ig=1,ngrid
-         zplev_gcm(ig,l) = ap(l)  + bp(l)*ps_start_GCM(ig)
-       ENDDO
-     ENDDO
-
-     global_ave_press_old=0.
-     do i=1,ngrid
-       global_ave_press_old=global_ave_press_old+cell_area(i)*ps_start_GCM(i)/Total_surface
-     enddo
-
-     global_ave_press_GCM=global_ave_press_old
-     global_ave_press_new=global_ave_press_old
-     print *, "Initial global average pressure=", global_ave_press_GCM
-
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
-!    I_b READ of start_evol.nc and starfi_evol.nc
-!    I_c Subslope parametrisation
-!    I_d READ GCM data and convert to the physical grid 
-!    I_e Initialisation of the PEM variable and soil 
-!    I_f Compute tendencies & Save initial situation
-!    I_g Save initial PCM situation
+ini_surf_co2 = 0.
+ini_surf_h2o = 0.
+Total_surface = 0.
+do i = 1,ngrid
+    Total_surface = Total_surface+cell_area(i)
+    do islope = 1,nslope
+        if (tendencies_co2_ice(i,islope) < 0) then
+            initial_co2_ice_sublim(i,islope) = 1.
+            ini_surf_co2 = ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
+        else
+            initial_co2_ice_sublim(i,islope) = 0.
+        endif
+        if (qsurf(i,igcm_co2,islope) > 0) then
+            initial_co2_ice(i,islope) = 1.
+        else
+            initial_co2_ice(i,islope) = 0.
+        endif
+        if (tendencies_h2o_ice(i,islope) < 0) then
+            initial_h2o_ice(i,islope) = 1.
+            ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
+        else
+            initial_h2o_ice(i,islope) = 0.
+        endif
+    enddo
+enddo
+
+write(*,*) "Total initial surface of co2ice 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
+allocate(zplev_gcm(ngrid,nlayer + 1))
+
+do l = 1,nlayer + 1
+    do ig = 1,ngrid
+        zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
+    enddo
+enddo
+
+global_ave_press_old = 0.
+do i = 1,ngrid
+       global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
+enddo
+
+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
+
+!------------------------
+! I   Initialization
 !    I_h Read the PEMstart
-
-!---------------------------- Read the PEMstart --------------------- 
-
-      call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth,porefillingice_thickness, &
-      tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,&
-      tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,&
-      watersurf_density_ave,watersoil_density_PEM_ave, &
-      co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
+!------------------------
+call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
+              porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,          &
+              tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),                      &
+              qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave,             &
+              co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
+
+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.)
+    enddo
+enddo
+
+if (adsorption_pem) then
+    totmassco2_adsorbded = 0.
+    totmassh2o_adsorbded = 0.
+    do ig = 1,ngrid
+        do islope = 1, nslope
+            do l = 1,nsoilmx_PEM - 1
+                totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
+                totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
+            enddo
+        enddo
+    enddo
+
+    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) 
+
+!------------------------
+! I   Initialization
+!    I_i Compute orbit criterion
+!------------------------
+#ifndef CPP_STD
+     call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
+#else
+     call iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
+#endif
+
+if (evol_orbit_pem) then
+    call orbit_param_criterion(year_iter_max)
+else
+    year_iter_max = Max_iter_pem
+endif
+!-------------------------- END INITIALIZATION -------------------------
+
+!-------------------------------- RUN ----------------------------------
+!------------------------
+! II  Run
+!    II_a Update pressure, ice and tracers 
+!------------------------
+year_iter = 0
+do while (year_iter < year_iter_max)
+
+! II.a.1. Compute updated global pressure
+    write(*,*) "Recomputing the new pressure..."
+
+    do i = 1,ngrid
+        do islope = 1,nslope
+            global_ave_press_new = global_ave_press_new-g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
+        enddo
+    enddo
+    write(*,*) 'Global average pressure old time step',global_ave_press_old
+    call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
+     
+    if (adsorption_pem) then
+        do i = 1,ngrid
+            global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
+        enddo 
+        write(*,*) 'Global average pressure old time step',global_ave_press_old
+        write(*,*) 'Global average pressure new time step',global_ave_press_new
+    endif
+
+! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
+    allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
+    write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
+    do l = 1,nlayer + 1
+        do ig = 1,ngrid
+            zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
+        enddo
+    enddo
+
+! II.a.3. Surface pressure timeseries
+    write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..."
+    do ig = 1,ngrid
+        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
+    enddo
+
+! II.a.4. New pressure levels timeseries
+    allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
+    write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
+    do l = 1,nlayer + 1
+        do ig = 1,ngrid
+            zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
+        enddo
+    enddo
+
+! II.a.5. Tracers timeseries
+    write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..."
+
+    l = 1
+    do ig = 1,ngrid
+        do t = 1,timelen
+            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.
+        enddo
+    enddo
+
+    do ig = 1,ngrid
+        do t = 1, timelen
+            q_co2_PEM_phys(ig,t) = q_co2_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))  &
+                                + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,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_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
+                q_co2_PEM_phys(ig,t) = 1.
+            endif
+            mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
+            vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
+        enddo
+    enddo
+
+    deallocate(zplev_new_timeseries,zplev_old_timeseries)
+
+!------------------------
+! II  Run
+!    II_b Evolution of the ice
+!------------------------
+    write(*,*) "Evolution of h2o ice"
+    call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)
+
+    write(*,*) "Evolution of co2 ice"
+    call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
+
+    do islope=1, nslope
+        write(str2(1:2),'(i2.2)') islope
+        call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
+        call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
+        call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
+        call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
+    enddo
+
+!------------------------
+! II  Run
+!    II_c CO2 & H2O glaciers flows
+!------------------------
+    write(*,*) "CO2 glacier flows"
+
+    if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
+                                               global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
   
-      do ig = 1,ngrid
+    write(*,*) "H2O glacier flows"
+
+    if (h2oglaciersflow) call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
+
+    do islope = 1,nslope
+        write(str2(1:2),'(i2.2)') islope
+        call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
+        call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
+        call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
+    enddo
+
+!------------------------
+! II  Run
+!    II_d Update surface and soil temperatures
+!------------------------
+! II_d.1 Update Tsurf
+    write(*,*) "Updating the new Tsurf"
+    bool_sublim = .false.
+    Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
+    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.)
-        enddo
-      enddo
-
-    if(adsorption_pem) then
-     totmassco2_adsorbded = 0.
-     totmassh2o_adsorbded = 0.
-     do ig = 1,ngrid
-      do islope =1, nslope
-       do l = 1,nsoilmx_PEM - 1
-          totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
-          subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
-          cell_area(ig)
-          totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
-          subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
-          cell_area(ig)
-       enddo
-      enddo
-     enddo
-
-     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) 
-
-!------------------------
-
-! I   Initialisation
-!    I_a READ run.def 
-!    I_b READ of start_evol.nc and starfi_evol.nc
-!    I_c Subslope parametrisation
-!    I_d READ GCM data and convert to the physical grid 
-!    I_e Initialisation of the PEM variable and soil 
-!    I_f Compute tendencies & Save initial situation
-!    I_g Save initial PCM situation
-!    I_h Read the PEMstar
-!    I_i Compute orbit criterion
-
-#ifndef CPP_STD
-     CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
-#else
-     CALL iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
-#endif
-
-     if (evol_orbit_pem) then
-       call orbit_param_criterion(year_iter_max)
-     else
-       year_iter_max = Max_iter_pem
-     endif
-
-!--------------------------- END INITIALISATION ---------------------
-
-!---------------------------- RUN ---------------------
-
-!------------------------
-
+            if (initial_co2_ice(ig,islope) > 0.5 .and. qsurf(ig,igcm_co2,islope) < 1.e-10) then !co2ice disappeared, look for closest point without co2ice
+                if (latitude_deg(ig) > 0) then
+                    do ig_loop = ig,ngrid
+                        do islope_loop = islope,iflat,-1
+                            if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
+                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
+                                bool_sublim = .true.
+                                exit
+                            endif
+                        enddo
+                        if (bool_sublim) exit
+                    enddo
+                else
+                    do ig_loop = ig,1,-1
+                        do islope_loop = islope,iflat
+                            if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
+                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
+                                bool_sublim = .true.
+                                exit
+                            endif
+                        enddo
+                        if (bool_sublim) exit
+                    enddo
+                endif
+                initial_co2_ice(ig,islope) = 0
+                if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
+                    albedo(ig,1,islope) = albedo_h2o_frost
+                    albedo(ig,2,islope) = albedo_h2o_frost
+                else
+                    albedo(ig,1,islope) = albedodat(ig)
+                    albedo(ig,2,islope) = albedodat(ig) 
+                    emis(ig,islope) = emissiv
+                endif
+            else if ((qsurf(ig,igcm_co2,islope) > 1.e-3) .and. (tendencies_co2_ice(ig,islope) > 1.e-10)) then !Put tsurf as tcond co2
+                ave = 0.
+                do t = 1,timelen
+                    if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
+                        ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
+                    else
+                        ave = ave + tsurf_GCM_timeseries(ig,islope,t)
+                    endif
+                enddo
+                tsurf_ave(ig,islope) = ave/timelen
+            endif
+        enddo
+    enddo
+
+    do t = 1,timelen
+        tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:)) 
+    enddo
+    ! for the start
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
+        enddo
+    enddo
+
+    do islope = 1,nslope
+        write(str2(1:2),'(i2.2)') islope
+        call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
+    enddo
+
+    if (soil_pem) then
+
+! II_d.2 Update soil temperature
+        allocate(TI_locslope(ngrid,nsoilmx_PEM))
+        allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
+        allocate(Tsurf_locslope(ngrid))
+        write(*,*)"Updating soil temperature"
+
+        ! Soil averaged
+        do islope = 1,nslope
+            TI_locslope(:,:) = TI_PEM(:,:,islope)
+            do t = 1,timelen
+                Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
+                Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
+                call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
+                call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
+                tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
+                do ig = 1,ngrid
+                    do isoil = 1,nsoilmx_PEM
+                        watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r)
+                        if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) 
+                    enddo
+                enddo
+            enddo
+        enddo
+        tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
+        watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
+
+        write(*,*) "Update of soil temperature done"
+
+        deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
+        write(*,*) "Compute ice table"
+
+! II_d.3 Update the ice table
+        call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
+        write(*,*) "Update soil propreties"
+
+! II_d.4 Update the soil thermal properties
+        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
+        if (adsorption_pem) then
+            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
+                                     qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, &
+                                     q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
+
+            totmassco2_adsorbded = 0.
+            totmassh2o_adsorbded = 0.
+            do ig = 1,ngrid
+                do islope =1, nslope
+                    do l = 1,nsoilmx_PEM - 1
+                        totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
+                        cell_area(ig)
+                        totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
+                        cell_area(ig)
+                    enddo
+                enddo
+            enddo
+            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
+            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
+        endif
+    endif !soil_pem
+
+!------------------------
 ! II  Run
-!    II_a update pressure,ice and tracers 
-
-!------------------------
-     year_iter=0
-     do while (year_iter.LT.year_iter_max)
-
-! II.a.1. Compute updated global pressure
-     print *, "Recomputing the new pressure..."
-
-     do i=1,ngrid
-       do islope=1,nslope
-           global_ave_press_new=global_ave_press_new-g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
-      enddo
-     enddo
-     print *, 'Global average pressure old time step',global_ave_press_old
-     call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
-     
-     if(adsorption_pem) then
-      do i=1,ngrid
-        global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
-       enddo 
-       print *, 'Global average pressure old time step',global_ave_press_old
-       print *, 'Global average pressure new time step',global_ave_press_new
-     endif
-
-! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
-     allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
-     print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
-     DO l=1,nlayer+1
-       DO ig=1,ngrid
-         zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_timeseries(ig,:)
-       ENDDO
-     ENDDO
-
-! II.a.3. Surface pressure timeseries
-     print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
-     do ig = 1,ngrid
-       ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
-     enddo
-
-! II.a.4. New pressure levels timeseries
-     allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
-     print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
-     do l=1,nlayer+1
-       do ig=1,ngrid
-         zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_timeseries(ig,:)
-       enddo
-     enddo
-
-! II.a.5. Tracers timeseries
-      print *, "Recomputing of tracer VMR timeseries for the new pressure..."
-
-     l=1
-     DO ig=1,ngrid
-       DO t = 1, timelen
-         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).LT.0) then
-           q_h2o_PEM_phys(ig,t)=1E-30
-         endif
-         if(q_h2o_PEM_phys(ig,t).GT.1) then
-           q_h2o_PEM_phys(ig,t)=1.
-         endif
-       enddo
-     enddo
-
-     DO ig=1,ngrid
-       DO t = 1, timelen
-         q_co2_PEM_phys(ig,t)=q_co2_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))  &
-                                + (   (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,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_co2_PEM_phys(ig,t).LT.0) then
-              q_co2_PEM_phys(ig,t)=1E-30
-         elseif (q_co2_PEM_phys(ig,t).GT.1) then
-              q_co2_PEM_phys(ig,t)=1.
-         endif
-         mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
-         vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
-       ENDDO
-     ENDDO
-   
-     deallocate(zplev_new_timeseries)
-     deallocate(zplev_old_timeseries)
-
+!    II_e Update the tendencies
+!------------------------
+    write(*,*) "Adaptation of the new co2 tendencies to the current pressure"
+    call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
+                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
+
+!------------------------
 ! II  Run
-!    II_a update pressure, ice and tracers 
-!    II_b Evolution of the ice
-
-! II.b. Evolution of the ice
-     print *, "Evolution of h2o ice"
-     call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)
-     
-     print *, "Evolution of co2 ice"
-     call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
-
-     DO islope=1, nslope
-       write(str2(1:2),'(i2.2)') islope
-       call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
-       call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
-       call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
-       call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
-     ENDDO
-
-!------------------------
-
-! II  Run
-!    II_a update pressure, ice and tracers 
-!    II_b Evolution of the ice
-!    II_c CO2  & H2O glaciers flows 
-
-!------------------------
-
-      print *, "CO2 glacier flows"
-
-      if (co2glaciersflow) then
-       call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,&
-                         global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
-      endif
-
-      print *, "H2O glacier flows"
-
-      if (h2oglaciersflow) then
-        call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
-      endif
-     DO islope=1, nslope
-       write(str2(1:2),'(i2.2)') islope
-       call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
-       call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
-       call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
-     ENDDO
-
-!------------------------
-
-! II  Run
-!    II_a update pressure, ice and tracers 
-!    II_b Evolution of the ice
-!    II_c CO2 glaciers flows 
-!    II_d Update surface and soil temperatures
-
-!------------------------
-     
-! II_d.1 Update Tsurf
-
-      print *, "Updating the new Tsurf"
-      bool_sublim=.false.
-      Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
-      DO ig = 1,ngrid
-        DO islope = 1,nslope
-          if(initial_co2_ice(ig,islope).gt.0.5 .and. qsurf(ig,igcm_co2,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice
-            if(latitude_deg(ig).gt.0) then
-              do ig_loop=ig,ngrid
-                DO islope_loop=islope,iflat,-1
-                  if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then
-                    tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
-                    bool_sublim=.true.
-                    exit
-                  endif
-                enddo
-                if (bool_sublim.eqv. .true.) then
-                  exit
-                endif
-              enddo
-            else
-              do ig_loop=ig,1,-1
-                DO islope_loop=islope,iflat
-                  if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then
-                    tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
-                    bool_sublim=.true.
-                    exit
-                  endif
-                enddo
-                if (bool_sublim.eqv. .true.) then
-                  exit
-                endif
-              enddo
-            endif
-            initial_co2_ice(ig,islope)=0
-            if ((qsurf(ig,igcm_co2,islope).lt.1e-10).and. (qsurf(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then   
-              albedo(ig,1,islope) = albedo_h2o_frost
-              albedo(ig,2,islope) = albedo_h2o_frost
-            else
-              albedo(ig,1,islope) = albedodat(ig)
-              albedo(ig,2,islope) = albedodat(ig) 
-              emis(ig,islope) = emissiv         
-            endif
-          elseif( (qsurf(ig,igcm_co2,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
-            ave=0.
-            do t=1,timelen
-              if(co2_ice_GCM(ig,islope,t).gt.1e-3) then
-                ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
-              else
-                ave = ave + tsurf_GCM_timeseries(ig,islope,t)
-              endif
-            enddo
-            tsurf_ave(ig,islope)=ave/timelen
-          endif
-        enddo
-      enddo
-
-      do t = 1,timelen
-        tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) +( tsurf_ave(:,:) -Tsurfave_before_saved(:,:)) 
-      enddo
-      ! for the start
-      do ig = 1,ngrid
-        do islope = 1,nslope
-          tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope))
-        enddo
-      enddo
-
-     DO islope=1, nslope
-       write(str2(1:2),'(i2.2)') islope
-       call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
-     ENDDO
-
-    if(soil_pem) then
-
-! II_d.2 Update soil temperature
-
-  allocate(TI_locslope(ngrid,nsoilmx_PEM))
-  allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
-  allocate(Tsurf_locslope(ngrid))
-  print *,"Updating soil temperature"
-
-! Soil averaged
-  do islope = 1,nslope
-     TI_locslope(:,:) = TI_PEM(:,:,islope)
-     do t = 1,timelen
-       Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
-       Tsurf_locslope(:) =  tsurf_GCM_timeseries(:,islope,t)
-       call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
-       call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
-       tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
-       do ig = 1,ngrid
-         do isoil = 1,nsoilmx_PEM
-          watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r)
-          if(isnan(Tsoil_locslope(ig,isoil))) then
-            call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1) 
-          endif
-         enddo
-       enddo
-     enddo
-  enddo
-     tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
-     watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
-
-  print *, "Update of soil temperature done"
-
-  deallocate(TI_locslope)
-  deallocate(Tsoil_locslope)
-  deallocate(Tsurf_locslope)
-  write(*,*) "Compute ice table"
-
-! II_d.3 Update the ice table
-     call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
-     print *, "Update soil propreties"
-
-! II_d.4 Update the soil thermal properties
-      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
-     if(adsorption_pem) then
-       call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:), &
-                             tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
-                             h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
-  
-          
-      totmassco2_adsorbded = 0.
-      totmassh2o_adsorbded = 0.
-      do ig = 1,ngrid
-       do islope =1, nslope
-        do l = 1,nsoilmx_PEM - 1
-           totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
-           subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
-           cell_area(ig)
-           totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
-           subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
-           cell_area(ig)
-        enddo
-       enddo
-      enddo
-      write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
-      write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded                
-      endif
-  endif !soil_pem
-
-!------------------------
-
-! II  Run
-!    II_a update pressure, ice and tracers 
-!    II_b Evolution of the ice
-!    II_c CO2 glaciers flows 
-!    II_d Update surface and soil temperatures
-!    II_e Update the tendencies 
-
-!------------------------
-
-      print *, "Adaptation of the new co2 tendencies to the current pressure"
-      call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,&     
-                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
-
-!------------------------
-
-! II  Run
-!    II_a update pressure, ice and tracers 
-!    II_b Evolution of the ice
-!    II_c CO2 glaciers flows 
-!    II_d Update surface and soil temperatures
-!    II_e Update the tendencies 
 !    II_f Checking the stopping criterion
-
-!------------------------
-      call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
-
-      call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
-                                   initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
-
-      year_iter=year_iter+dt_pem     
-
-      print *, "Checking all the stopping criterion."
-      if (STOPPING_water) then
-        print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
-        criterion_stop=1
-      endif
-      if (STOPPING_1_water) then
-        print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
-        criterion_stop=1
-      endif
-      if (STOPPING_co2) then
-        print *, "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
-        criterion_stop=2
-      endif
-      if (STOPPING_pressure) then
-        print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
-        criterion_stop=3
-      endif
-      if (year_iter.ge.year_iter_max) then
-        print *, "STOPPING because maximum number of iterations reached"
-        criterion_stop=4
-      endif
-
-     if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
+!------------------------
+    call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
+
+    call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
+                            initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
+
+    year_iter = year_iter + dt_pem
+
+    write(*,*) "Checking all the stopping criterion."
+    if (STOPPING_water) then
+        write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
+        criterion_stop = 1
+    endif
+    if (STOPPING_1_water) then
+        write(*,*) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
+        criterion_stop = 1
+    endif
+    if (STOPPING_co2) then
+        write(*,*) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
+        criterion_stop = 2
+    endif
+    if (STOPPING_pressure) then
+        write(*,*) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
+        criterion_stop = 3
+    endif
+    if (year_iter >= year_iter_max) then
+        write(*,*) "STOPPING because maximum number of iterations reached"
+        criterion_stop = 4
+    endif
+
+    if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
         exit
-      else
-        print *, "We continue!"
-        print *, "Number of iteration of the PEM : year_iter=", year_iter
-      endif
-
-      global_ave_press_old=global_ave_press_new
-
-      enddo  ! big time iteration loop of the pem
-
-
-!---------------------------- END RUN PEM ---------------------
-
-!---------------------------- OUTPUT ---------------------
-
-!------------------------
-
+    else
+        write(*,*) "We continue!"
+        write(*,*) "Number of iteration of the PEM : year_iter=", year_iter
+    endif
+
+    global_ave_press_old=global_ave_press_new
+
+enddo  ! big time iteration loop of the pem
+!------------------------------ END RUN --------------------------------
+
+!------------------------------- OUTPUT --------------------------------
+!------------------------
 ! III Output
 !    III_a Update surface value for the PCM start files
-
-!------------------------
-
+!------------------------
 ! III_a.1 Ice update (for startfi)
 
 ! H2O ice
-     DO ig=1,ngrid
-       if(watercaptag(ig)) then
-         watercap_sum=0.
-         DO islope=1,nslope
-           if(qsurf(ig,igcm_h2o_ice,islope).GT. (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
-           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.))
-             qsurf(ig,igcm_h2o_ice,islope)=0.
+do ig = 1,ngrid
+    if (watercaptag(ig)) then
+        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
+            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.))
+                qsurf(ig,igcm_h2o_ice,islope)=0.
+            endif
+            watercap_sum = watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+            watercap(ig,islope) = 0.
+        enddo
+        water_reservoir(ig) = water_reservoir(ig) + watercap_sum
+    endif
+enddo
+
+do ig = 1,ngrid
+    water_sum = 0.
+    do islope = 1,nslope
+        water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+    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_frost2perenial) 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_frost2perenial/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_frost2perenial/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_frost2perenial) then
+            watercaptag(ig) = .false.
+            do islope = 1,nslope
+                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
+            enddo
+            water_reservoir(ig) = 0.
+        endif
+    endif
+enddo
+
+! CO2 ice
+do ig = 1,ngrid
+    do islope = 1,nslope
+        if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then
+            perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
+            qsurf(ig,igcm_co2,islope)  = 0.5*qsurf(ig,igcm_co2,islope)
+        endif
+    enddo
+enddo
+
+! III_a.2  Tsoil update (for startfi)
+if (soil_pem) then
+    call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
+    tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)      
+endif
+
+! III_a.4 Pressure (for start)
+do i = 1,ip1jmp1
+    ps(i) = ps(i)*global_ave_press_new/global_ave_press_GCM
+enddo
+
+do i = 1,ngrid
+    ps_start_GCM(i) = ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM
+enddo
+
+! III_a.5 Tracer (for start)
+allocate(zplev_new(ngrid,nlayer + 1))
+
+do l = 1,nlayer + 1
+    do ig = 1,ngrid
+        zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
+    enddo
+enddo
+
+do nnq = 1,nqtot
+    if (noms(nnq) /= "co2") then
+        do l = 1,llm - 1
+            do ig = 1,ngrid
+                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
+            enddo
+            q(:,llm,nnq) = q(:,llm - 1,nnq)
+        enddo
+    else
+        do l = 1,llm - 1
+            do ig = 1,ngrid
+                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
+                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ &
+                              (zplev_new(ig,l) - zplev_new(ig,l + 1))
+            enddo
+            q(:,llm,nnq) = q(:,llm - 1,nnq)
+        enddo
+    endif
+enddo
+
+! Conserving the tracers mass for GCM start files
+do nnq = 1,nqtot
+    do ig = 1,ngrid
+        do l = 1,llm - 1
+            if (q(ig,l,nnq) > 1 .and. (noms(nnq) /= "dust_number") .and. (noms(nnq) /= "ccn_number") .and. (noms(nnq) /= "stormdust_number") .and. (noms(nnq) /= "topdust_number")) then
+                extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
+                q(ig,l,nnq)=1.
+                q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
+                write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
            endif
-           watercap_sum=watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-           watercap(ig,islope)=0.
-         enddo
-           water_reservoir(ig)=water_reservoir(ig)+watercap_sum
-       endif
-     enddo
-
-     DO ig=1,ngrid
-       water_sum = 0.
-       DO islope=1,nslope
-         water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-       ENDDO
-       if(watercaptag(ig).eqv..false.) then ! let's check if we have an 'infinite' source of water that has been forming. 
-         if(water_sum.gt.threshold_water_frost2perenial) 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_frost2perenial/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_frost2perenial/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)).lt.threshold_water_frost2perenial) then
-           watercaptag(ig)=.false.
-           DO islope = 1,nslope
-              qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
-           ENDDO
-           water_reservoir(ig) = 0.
-         endif
-       endif
-     enddo
-
-! CO2 ice
-
-     DO ig=1,ngrid
-       DO islope=1,nslope
-         if(qsurf(ig,igcm_co2,islope).gt.threshold_co2_frost2perenial) then
-           perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
-           qsurf(ig,igcm_co2,islope)  = 0.5*qsurf(ig,igcm_co2,islope)
-         endif
-       ENDDO
-     ENDDO
-       
-
-! III_a.2  Tsoil update (for startfi)
-
-   if(soil_pem) then
-     call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
-     tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)      
-   endif !soil_pem
-
-! III_a.4 Pressure (for start)
-     do i=1,ip1jmp1
-       ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
-     enddo
-
-     do i = 1,ngrid
-       ps_start_GCM(i)=ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM
-     enddo
-
-! III_a.5 Tracer (for start)
-     allocate(zplev_new(ngrid,nlayer+1))
-
-     do l=1,nlayer+1
-       do ig=1,ngrid
-         zplev_new(ig,l) = ap(l)  + bp(l)*ps_start_GCM(ig)
-       enddo
-     enddo
-
-     DO nnq=1,nqtot
-       if (noms(nnq).NE."co2") then
-         DO l=1,llm-1
-           DO ig=1,ngrid
-             q(ig,l,nnq)=q(ig,l,nnq)*(zplev_gcm(ig,l)-zplev_gcm(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))
-           ENDDO
-           q(:,llm,nnq)=q(:,llm-1,nnq)
-         ENDDO
-       else
-         DO l=1,llm-1
-           DO ig=1,ngrid
-             q(ig,l,nnq)=q(ig,l,nnq)*(zplev_gcm(ig,l)-zplev_gcm(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))  &
-                             + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
-                                   (zplev_gcm(ig,l)-zplev_gcm(ig,l+1))     )  / &
-                                   (zplev_new(ig,l)-zplev_new(ig,l+1))
-           ENDDO
-           q(:,llm,nnq)=q(:,llm-1,nnq)
-         ENDDO
-       endif
-     ENDDO
-
-! Conserving the tracers's mass for GCM start files
-     DO nnq=1,nqtot
-       DO ig=1,ngrid
-         DO l=1,llm-1
-           if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") .and. (noms(nnq).NE."stormdust_number") .and. (noms(nnq).NE."topdust_number")) then
-             extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
-             q(ig,l,nnq)=1.
-             q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
-             write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number"
-           endif
-           if(q(ig,l,nnq).LT.0) then
-             q(ig,l,nnq)=1E-30
-           endif
-         ENDDO
-       enddo
-     enddo
-    
-!------------------------
-     if(evol_orbit_pem) then
-      call recomp_orb_param(year_iter)
-     endif
-
+            if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
+        enddo
+    enddo
+enddo
+
+if (evol_orbit_pem) call recomp_orb_param(year_iter)
+
+!------------------------
 ! III Output
-!    III_a Update surface value for the PCM start files
-!    III_b Write start and starfi.nc
-
-!------------------------
-
-! III_b.1 WRITE restart.nc 
-
-      ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
-      pday=day_ini
-      ztime_fin=0.
-
-     allocate(p(ip1jmp1,nlayer+1))
+!    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
+pday = day_ini
+ztime_fin = 0.
+
+allocate(p(ip1jmp1,nlayer + 1))
 #ifndef CPP_1D
-     CALL pression (ip1jmp1,ap,bp,ps,p)
-     CALL massdair(p,masse)
-
-      CALL dynredem0("restart_evol.nc", day_ini, phis)
-
-      CALL dynredem1("restart_evol.nc",   &
-      		time_0,vcov,ucov,teta,q,masse,ps)
-      print *, "restart_evol.nc has been written"
-
+    call pression (ip1jmp1,ap,bp,ps,p)
+    call massdair(p,masse)
+    call dynredem0("restart_evol.nc", day_ini, phis)
+    call dynredem1("restart_evol.nc",   &
+    time_0,vcov,ucov,teta,q,masse,ps)
+    write(*,*) "restart_evol.nc has been written"
 #else
-     DO nnq = 1, nqtot
-      	call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf)
-     ENDDO
+    do nnq = 1, nqtot
+        call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf)
+    enddo
 #endif
 
-! III_b.2 WRITE restartfi.nc 
+! III_b.2 Write restartfi_evol.nc
 #ifndef CPP_STD
-      call physdem0("restartfi_evol.nc",longitude,latitude, &
-                        nsoilmx,ngrid,nlayer,nq,              &
-                        ptimestep,pday,0.,cell_area,          &
-                        albedodat,inertiedat,def_slope,   &
-                        subslope_dist)
-
-     call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
-                     ptimestep,ztime_fin,                        &
-                     tsurf,tsoil,inertiesoil,albedo,             &
-                     emis,q2,qsurf,tauscaling,totcloudfrac,wstar,&
-                     watercap,perenial_co2ice)
+    call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
+                  nlayer,nq,ptimestep,pday,0.,cell_area,albedodat,      &
+                  inertiedat,def_slope,subslope_dist)
+
+    call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
+                  ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,  &
+                  albedo,emis,q2,qsurf,tauscaling,totcloudfrac, &
+                  wstar,watercap,perenial_co2ice)
 #else
-     call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
-                         ptimestep,pday,time_phys,cell_area,          &
-                         albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
-
-     call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
-                          ptimestep,ztime_fin,                    &
-                          tsurf,tsoil,emis,q2,qsurf,         &
-                          cloudfrac,totcloudfrac,hice,            &
-                          rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+    call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
+                  nlayer,nq,ptimestep,pday,time_phys,cell_area,         &
+                  albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
+
+    call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,   &
+                  ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf, &
+                  cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,   &
+                  tslab,tsea_ice,sea_ice)
 #endif
-
-      print *, "restartfi_evol.nc has been written"
-!------------------------
-
+write(*,*) "restartfi_evol.nc has been written"
+
+!------------------------
 ! III Output
-!    III_a Update surface value for the PCM start files
-!    III_b Write start and starfi.nc
-!    III_c Write start_pem
-
-!------------------------
-        call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
-                         float(day_ini),0.,nslope,def_slope,subslope_dist)
-
-
-        call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
-                      tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, &
-                      co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
-
-        call info_run_PEM(year_iter, criterion_stop)
-
-      print *, "restartfi_PEM.nc has been written"
-      print *, "The PEM had run for ", year_iter, " years."
-      print *, "LL & RV : So far so good"
-
-     deallocate(vmr_co2_gcm)
-     deallocate(ps_timeseries)
-     deallocate(tsurf_GCM_timeseries)
-     deallocate(q_co2_PEM_phys)
-     deallocate(q_h2o_PEM_phys)
-     deallocate(co2_ice_GCM)
-     deallocate(watersurf_density_ave)
-     deallocate(watersoil_density_timeseries)
-     deallocate(Tsurfave_before_saved)
-     deallocate(tsoil_phys_PEM_timeseries)
-     deallocate(watersoil_density_PEM_timeseries)
-     deallocate(watersoil_density_PEM_ave)
-     deallocate(delta_co2_adsorbded)
-     deallocate(delta_h2o_adsorbded)
-     deallocate(vmr_co2_pem_phys)
+!    III_c Write restartfi_PEM.nc
+!------------------------
+call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
+             float(day_ini),0.,nslope,def_slope,subslope_dist)
+
+
+call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope ,           &
+             tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, &
+             co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
+
+call info_run_PEM(year_iter,criterion_stop)
+
+write(*,*) "restartfi_PEM.nc has been written"
+write(*,*) "The PEM had run for ", year_iter, " years."
+write(*,*) "LL & RV : So far so good"
+
+deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
+deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
+deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
+deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys)
+!----------------------------- END OUTPUT ------------------------------
 
 END PROGRAM pem
+
Index: trunk/LMDZ.MARS/changelog.txt
===================================================================
--- trunk/LMDZ.MARS/changelog.txt	(revision 3027)
+++ trunk/LMDZ.MARS/changelog.txt	(revision 3028)
@@ -4169,2 +4169,6 @@
 == 09/08/2023 == JBC
 Small fixes to be able to run the Mars PCM 1D without "water" + Improvements/addition of scripts in deftank/pem to run the PEM 1D model according to Laskar orbital parameters.
+
+== 14/08/2023 == JBC
+Related to commit r3026: improvement of error message in initracer.F (now it gives correctly the only identified tracers) + one small correction to run PCM 1D without water.
+Big cleaning of main program pem.F90 (indentation, declarations, comments, simplification of conditions/loops, etc).
