source: trunk/LMDZ.COMMON/libf/evolution/pem.F90 @ 3820

Last change on this file since 3820 was 3809, checked in by jbclement, 3 weeks ago

PEM:

  • Correction to detect whether the stratum below the dust lag is made of ice or not.
  • Correction for the initialization of surface ice and 'h2o_ice_depth' according to the layerings map.
  • Correction to handle flow of glaciers with the layering algorithm.
  • Simplification of "recomp_tend_h2o", making it more robust.
  • Making the removing of a stratum more robust.
  • Plotting the orbital variations over the simulated time in "visu_evol_layering.py".
  • Lowering the value of dust tendency.
  • Addition of "_co2" in the variables name for the CO2 flux correction.

JBC

File size: 71.8 KB
Line 
1!------------------------
2! I   Initialization
3!    I_a Read the "run.def"
4!    I_b Read the "start.nc" and "startfi.nc"
5!    I_c Subslope parametrisation
6!    I_d Read the PCM data and convert them to the physical grid
7!    I_e Initialization of the PEM variable and soil
8!    I_f Compute tendencies
9!    I_g Compute global surface pressure
10!    I_h Read the "startpem.nc"
11!    I_i Compute orbit criterion
12
13! II  Run
14!    II_a Update pressure, ice and tracers
15!    II_b Evolution of ice
16!    II_c Flow of glaciers
17!    II_d Update surface and soil temperatures
18!    II_e Outputs
19!    II_f Update the tendencies
20!    II_g Checking the stopping criterion
21
22! III Output
23!    III_a Update surface values for the PCM start files
24!    III_b Write the "restart.nc" and "restartfi.nc"
25!    III_c Write the "restartpem.nc"
26!------------------------
27
28PROGRAM pem
29
30use phyetat0_mod,               only: phyetat0
31use phyredem,                   only: physdem0, physdem1
32use netcdf,                     only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
33use turb_mod,                   only: q2, wstar
34use comslope_mod,               only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h
35use logic_mod,                  only: iflag_phys
36use mod_const_mpi,              only: COMM_LMDZ
37use infotrac
38use geometry_mod,               only: latitude_deg
39use conf_pem_mod,               only: conf_pem
40use pemredem,                   only: pemdem0, pemdem1
41use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
42                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
43use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2
44use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
45use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
46use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
47                                      TI_PEM,               & ! Soil thermal inertia
48                                      tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
49                                      fluxgeo                 ! Geothermal flux for the PEM and PCM
50use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
51                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
52                                      co2_adsorbed_phys, h2o_adsorbed_phys          ! Mass of co2 and h2O adsorbed
53use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
54use orbit_param_criterion_mod,  only: orbit_param_criterion
55use recomp_orb_param_mod,       only: recomp_orb_param
56use ice_table_mod,              only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, &
57                                      ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi
58use soil_thermalproperties_mod, only: update_soil_thermalproperties
59use time_phylmdz_mod,           only: daysec, dtphys
60use abort_pem_mod,              only: abort_pem
61use soil_settings_PEM_mod,      only: soil_settings_PEM
62use compute_tend_mod,           only: compute_tend
63use info_PEM_mod,               only: info_PEM
64use get_timelen_PCM_mod,        only: get_timelen_PCM
65use pemetat0_mod,               only: pemetat0
66use read_data_PCM_mod,          only: read_data_PCM
67use recomp_tend_mod,            only: recomp_tend_co2, recomp_tend_h2o
68use compute_soiltemp_mod,       only: compute_tsoil_pem, shift_tsoil2surf
69use writediagpem_mod,           only: writediagpem, writediagsoilpem
70use co2condens_mod,             only: CO2cond_ps
71use layering_mod,               only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering,            &
72                                      ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str, &
73                                      print_layering
74use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
75use version_info_mod,           only: print_version_info
76use paleoclimate_mod,           only: h2o_ice_depth, zdqsdif_ssi_tot
77
78#ifndef CPP_STD
79    use comsoil_h,          only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil
80    use surfdat_h,          only: tsurf, qsurf, emis, emissiv, emisice, ini_surfdat_h,   &
81                                  albedodat, albedice, albedo_h2o_frost, albedo_h2o_cap, &
82                                  zmea, zstd, zsig, zgam, zthe, frost_albedo_threshold,  &
83                                  watercap, watercaptag, perennial_co2ice, albedo_perennialco2
84    use dimradmars_mod,     only: totcloudfrac, albedo
85    use dust_param_mod,     only: tauscaling
86    use tracer_mod,         only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
87    use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
88    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
89    use surfini_mod,        only: surfini
90    use comcstfi_h,         only: mugaz
91#else
92    use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
93    use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT,   &
94                                  PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
95                                  ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
96    use aerosol_mod,        only: iniaerosol
97    use planete_mod,        only: apoastr, periastr, year_day, peri_day, obliquit
98    use comcstfi_mod,       only: pi, rad, g, r, cpp, rcp, mugaz
99#endif
100
101#ifndef CPP_1D
102    use comconst_mod,             only: pi, rad, g, r, cpp, rcp => kappa
103    use iniphysiq_mod,            only: iniphysiq
104    use control_mod,              only: iphysiq, day_step, nsplit_phys
105#else
106    use comcstfi_h,               only: pi, rad, g, r, cpp, rcp
107    use time_phylmdz_mod,         only: iphysiq, steps_per_sol
108    use regular_lonlat_mod,       only: init_regular_lonlat
109    use physics_distribution_mod, only: init_physics_distribution
110    use mod_grid_phy_lmdz,        only: regular_lonlat
111    use init_testphys1d_mod,      only: init_testphys1d
112    use comvert_mod,              only: ap, bp
113    use writerestart1D_mod,       only: writerestart1D
114#endif
115
116implicit none
117
118include "dimensions.h"
119include "paramet.h"
120include "comgeom.h"
121include "iniprint.h"
122
123! Same variable names as in the PCM
124integer, parameter :: ngridmx = 2 + (jjm - 1)*iim - 1/jjm
125integer, parameter :: nlayer = llm ! Number of vertical layer
126integer            :: ngrid        ! Number of physical grid points
127integer            :: nq           ! Number of tracer
128integer            :: day_ini      ! First day of the simulation
129real               :: pday         ! Physical day
130real               :: time_phys    ! Same as in PCM
131real               :: ptimestep    ! Same as in PCM
132real               :: ztime_fin    ! Same as in PCM
133
134! Variables to read "start.nc"
135character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM
136
137! Dynamic variables
138real, dimension(ip1jm,llm)          :: vcov          ! vents covariants
139real, dimension(ip1jmp1,llm)        :: ucov          ! vents covariants
140real, dimension(ip1jmp1,llm)        :: teta          ! Potential temperature
141real, dimension(:,:,:), allocatable :: q             ! champs advectes
142real, dimension(:),     allocatable :: pdyn          ! pressure for the dynamic grid
143real, dimension(:),     allocatable :: ps_start      ! surface pressure in the start file
144real, dimension(:),     allocatable :: ps_start0     ! surface pressure in the start file at the beginning
145real, dimension(:),     allocatable :: ps_avg        ! (ngrid) Average surface pressure
146real, dimension(:),     allocatable :: ps_dev        ! (ngrid x timelen) Surface pressure deviation
147real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure
148real, dimension(ip1jmp1,llm)        :: masse         ! Air mass
149real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
150real                                :: time_0
151
152! Variables to read starfi.nc
153character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM
154character(2)            :: str2
155integer                 :: ncid, status                           ! Variable for handling opening of files
156integer                 :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files
157integer                 :: apvarid, bpvarid                       ! Variable ID for Netcdf files
158
159! Variables to read starfi.nc and write restartfi.nc
160real, dimension(:), allocatable :: longitude     ! Longitude read in startfi_name and written in restartfi
161real, dimension(:), allocatable :: latitude      ! Latitude read in startfi_name and written in restartfi
162real, dimension(:), allocatable :: cell_area     ! Cell_area read in startfi_name and written in restartfi
163real                            :: total_surface ! Total surface of the planet
164
165! Variables for h2o ice evolution
166real, dimension(:,:),    allocatable  :: h2o_ice              ! h2o ice in the PEM
167real, dimension(:,:),    allocatable  :: d_h2oice             ! physical point x slope field: Tendency of evolution of perennial h2o ice
168real, dimension(:,:,:),  allocatable  :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
169real                                  :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
170logical, dimension(:,:), allocatable  :: is_h2oice_sublim_ini ! Logical array to know if h2o ice is sublimating
171real                                  :: ps_avg_global_ini    ! constant: Global average pressure at initialization [Pa]
172real                                  :: ps_avg_global_old    ! constant: Global average pressure of previous time step
173real                                  :: ps_avg_global_new    ! constant: Global average pressure of current time step
174real,   dimension(:,:),   allocatable :: zplev_new            ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
175real,   dimension(:,:),   allocatable :: zplev_start0         ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
176real,   dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
177real,   dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
178integer                               :: stopPEM              ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
179real                                  :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
180real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
181integer                               :: timelen              ! # time samples
182real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
183
184! Variables for co2 ice evolution
185real,    dimension(:,:), allocatable :: co2_ice                ! co2 ice in the PEM
186real,    dimension(:,:), allocatable :: d_co2ice               ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
187real,    dimension(:,:), allocatable :: d_co2ice_ini           ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
188logical, dimension(:,:), allocatable :: is_co2ice_ini          ! Was there co2 ice initially in the PEM?
189real,  dimension(:,:,:), allocatable :: min_co2_ice            ! Minimum of co2 ice at each point for the first year [kg/m^2]
190real                                 :: co2ice_sublim_surf_ini ! Initial surface of sublimating co2 ice
191logical, dimension(:,:), allocatable :: is_co2ice_sublim_ini   ! Logical array to know if co2 ice is sublimating
192real,    dimension(:,:), allocatable :: vmr_co2_PCM            ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
193real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys       ! Grid points x Times co2 volume mixing ratio used in the PEM
194real,    dimension(:,:), allocatable :: q_co2_PEM_phys         ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
195real(kind = 16)                      :: totmass_co2ice, totmass_atmco2         ! Current total CO2 masses
196real(kind = 16)                      :: totmass_co2ice_ini, totmass_atmco2_ini ! Initial total CO2 masses
197
198! Variables for the evolution of layered layerings_map
199type(layering), dimension(:,:), allocatable :: layerings_map     ! Layering for each grid point and slope
200type(ptrarray), dimension(:,:), allocatable :: current           ! Current active stratum in the layering
201logical,        dimension(:,:), allocatable :: new_str, new_lag  ! Flags for the layering algorithm
202real,           dimension(:,:), allocatable :: h2o_ice_depth_old ! Old depth of subsurface ice layer
203
204! Variables for slopes
205integer(kind = 1), dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow
206integer(kind = 1), dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow
207
208! Variables for surface and soil
209real, dimension(:,:),     allocatable :: tsurf_avg                          ! Grid points x Slope field: Average surface temperature [K]
210real, dimension(:,:),     allocatable :: tsurf_dev                          ! Grid points x Slope field: Surface temperature deviation [K]
211real, dimension(:,:),     allocatable :: tsurf_avg_yr1                      ! Grid points x Slope field: Average surface temperature of first call of the PCM [K]
212real, dimension(:,:,:),   allocatable :: tsoil_avg                          ! Grid points x Soil x Slope field: Average Soil Temperature [K]
213real, dimension(:,:),     allocatable :: tsoil_avg_old                      ! Grid points x Soil field: Average Soil Temperature at the previous time step [K]
214real, dimension(:,:,:),   allocatable :: tsoil_dev                          ! Grid points x Soil x Slope field: Soil temperature deviation [K]
215real, dimension(:,:,:,:), allocatable :: tsoil_timeseries                   ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
216real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries               ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K]
217real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old           ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K]
218real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries       ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
219real, dimension(:,:),     allocatable :: watersurf_density_avg              ! Grid points x Slope: Average water surface density [kg/m^3]
220real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries   ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
221real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg          ! Grid points x Soil x Slopes: Average water soil density [kg/m^3]
222real, dimension(:),       allocatable :: delta_co2_adsorbed                 ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
223real, dimension(:),       allocatable :: delta_h2o_adsorbed                 ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
224real(kind = 16)                       :: totmass_adsco2, totmass_adsco2_ini ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
225real                                  :: totmass_adsh2o                     ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
226logical, dimension(:,:),  allocatable :: co2ice_disappeared                 ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
227real, dimension(:,:),     allocatable :: icetable_thickness_old             ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
228real, dimension(:,:,:),   allocatable :: ice_porefilling_old                ! ngrid x nslope: Ice pore filling at the previous iteration [m]
229real, dimension(:),       allocatable :: delta_h2o_icetablesublim           ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
230real, dimension(:),       allocatable :: porefill                           ! Pore filling (output) to compute the dynamic ice table
231real                                  :: ssi_depth                          ! Ice table depth (output) to compute the dynamic ice table
232real, dimension(:,:),     allocatable :: zshift_surf                        ! Elevation shift for the surface [m]
233real, dimension(:,:),     allocatable :: zlag                               ! Newly built lag thickness [m]
234real, dimension(:,:),     allocatable :: icetable_depth_old                 ! Old depth of the ice table
235
236! Some variables for the PEM run
237real, parameter       :: year_step = 1   ! Timestep for the pem
238real                  :: i_myear_leg     ! Number of iteration
239real                  :: n_myear_leg     ! Maximum number of iterations before stopping
240real                  :: i_myear         ! Global number of Martian years of the chained simulations
241real                  :: n_myear         ! Maximum number of Martian years of the chained simulations
242real                  :: timestep        ! Timestep [s]
243character(100)        :: arg             ! To read command-line arguments program was invoked
244logical               :: timewall        ! Flag to use the time limit stopping criterion in case of a PEM job
245integer(kind = 8)     :: cr              ! Number of clock ticks per second (count rate)
246integer(kind = 8)     :: c1, c2          ! Counts of processor clock
247character(100)        :: chtimelimit     ! Time limit for the PEM job outputted by the SLURM command
248real                  :: timelimit       ! Time limit for the PEM job in seconds
249real, parameter       :: antetime = 3600 ! Anticipation time to prevent reaching the job time limit: 3600 s by default (it should cover the computing time of the reshaping tool)
250integer               :: cstat, days, hours, minutes, seconds
251character(1)          :: sep
252character(8)          :: date
253character(10)         :: time
254character(5)          :: zone
255integer, dimension(8) :: values
256character(128)        :: dir = ' '
257character(32)         :: logname = ' '
258character(32)         :: hostname = ' '
259
260#ifdef CPP_STD
261    real                                :: frost_albedo_threshold = 0.05 ! Frost albedo threeshold to convert fresh frost to old ice
262    real                                :: albedo_h2o_frost              ! Albedo of h2o frost
263    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
264    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
265    real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
266    real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimension when reading form generic
267    real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimension when reading form generic
268    real, dimension(:,:),   allocatable :: tsurf                         ! Subslope variable, only needed in the GENERIC case
269    real, dimension(:,:,:), allocatable :: qsurf                         ! Subslope variable, only needed in the GENERIC case
270    real, dimension(:,:,:), allocatable :: tsoil                         ! Subslope variable, only needed in the GENERIC case
271    real, dimension(:,:),   allocatable :: emis                          ! Subslope variable, only needed in the GENERIC case
272    real, dimension(:,:),   allocatable :: watercap                      ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
273    logical, dimension(:),  allocatable :: watercaptag                   ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
274    real, dimension(:,:,:), allocatable :: albedo                        ! Subslope variable, only needed in the GENERIC case
275    real, dimension(:,:,:), allocatable :: inertiesoil                   ! Subslope variable, only needed in the GENERIC case
276#endif
277
278#ifdef CPP_1D
279    integer            :: nsplit_phys
280    integer, parameter :: jjm_value = jjm - 1
281    integer            :: day_step
282
283    ! Dummy variables to use the subroutine 'init_testphys1d'
284    logical                             :: therestart1D, therestartfi, prescribed_h2ovap
285    integer                             :: ndt, day0
286    real                                :: ptif, pks, day, gru, grv, h2ovap_relax_time
287    real, dimension(:),     allocatable :: zqsat
288    real, dimension(:,:,:), allocatable :: dq, dqdyn
289    real, dimension(nlayer)             :: play, w, q_prescribed_h2o_vap
290    real, dimension(nlayer + 1)         :: plev
291#else
292    integer, parameter                :: jjm_value = jjm
293    real, dimension(:),   allocatable :: ap ! Coefficient ap read in start_name and written in restart
294    real, dimension(:),   allocatable :: bp ! Coefficient bp read in start_name and written in restart
295    real, dimension(:,:), allocatable :: p  ! Grid points x Atmosphere: pressure to recompute and write in restart (ip1jmp1,llmp1)
296#endif
297
298! Loop variables
299integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
300logical :: num_str
301
302! CODE
303write(*,*) '  *    .          .   +     .    *        .  +      .    .       .      '
304write(*,*) '           +         _______  ________  ____    ____      *           + '
305write(*,*) ' +   .        *     |_   __ \|_   __  ||_   \  /   _|          .       *'
306write(*,*) '          .     .     | |__) | | |_ \_|  |   \/   |  *        *      .  '
307write(*,*) '       .              |  ___/  |  _| _   | |\  /| |      .        .     '
308write(*,*) '.  *          *      _| |_    _| |__/ | _| |_\/_| |_                  * '
309write(*,*) '            +       |_____|  |________||_____||_____|   +     .         '
310write(*,*) '  .      *          .   *       .   +       *          .        +      .'
311
312! Elapsed time with system clock
313call system_clock(count_rate = cr)
314call system_clock(c1)
315timewall = .true.
316timelimit = 86400 ! 86400 seconds = 24 h by default
317if (command_argument_count() > 0) then ! Get the number of command-line arguments
318    call get_command_argument(1,arg) ! Read the argument given to the program
319    num_str = .true.
320    do i = 1,len_trim(arg)
321        if (arg(i:i) < '0' .or. arg(i:i) > '9') then
322            num_str = .false.
323            exit
324        endif
325    enddo
326
327    if (num_str) then ! This is a numeric sting so we considerer this is the job id
328           ! Execute the system command
329            call execute_command_line('squeue -j '//trim(adjustl(arg))//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
330            if (cstat /= 0) then
331                call execute_command_line('qstat -f '//trim(adjustl(arg))//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt',cmdstat = cstat)
332                if (cstat > 0) then
333                    error stop 'pem: command execution failed!'
334                else if (cstat < 0) then
335                    error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!'
336                endif
337            endif
338            ! Read the output
339            open(1,file = 'tmp_cmdout.txt',status = 'old')
340            read(1,'(a)') chtimelimit
341            close(1)
342            chtimelimit = trim(adjustl(chtimelimit))
343            call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat)
344            if (cstat > 0) then
345                error stop 'pem: command execution failed!'
346            else if (cstat < 0) then
347                error stop 'pem: command execution not supported!'
348            endif
349            if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS"
350                read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds
351                timelimit = days*86400 + hours*3600 + minutes*60 + seconds
352            else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS"
353                read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds
354                timelimit = hours*3600 + minutes*60 + seconds
355            else ! 'chtimelimit' format is "MM:SS"
356                read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds
357                timelimit = minutes*60 + seconds
358            endif
359    else ! Arg is not a numeric string
360        select case (trim(adjustl(arg)))
361            case('version') ! Handle command‐line argument "version"
362                call print_version_info()
363                stop
364            case default
365                error stop "The argument given to the program is unknown!"
366        end select
367    endif
368else
369    timewall = .false.
370endif
371
372! Some user info
373call date_and_time(date,time,zone,values)
374call getcwd(dir)      ! Current directory
375call getlog(logname)  ! User name
376call hostnm(hostname) ! Machine/station name
377write(*,*)
378write(*,*) '********* PEM information *********'
379write(*,*) '+ User     : '//trim(logname)
380write(*,*) '+ Machine  : '//trim(hostname)
381write(*,*) '+ Directory: '//trim(dir)
382write(*,'(a,i2,a,i2,a,i4)') ' + Date     : ',values(3),'/',values(2),'/',values(1)
383write(*,'(a,i2,a,i2,a,i2,a)') ' + Time     : ',values(5),':',values(6),':',values(7)
384
385! Parallel variables
386#ifndef CPP_STD
387    is_sequential = .true.
388    is_parallel = .false.
389    is_mpi_root = .true.
390    is_omp_root = .true.
391    is_master = .true.
392#endif
393
394! Some constants
395day_ini = 0
396time_phys = 0.
397ngrid = ngridmx
398A = (1./m_co2 - 1./m_noco2)
399B = 1./m_noco2
400year_day = 669
401daysec = 88775.
402timestep = year_day*daysec*year_step
403
404!----------------------------- INITIALIZATION --------------------------
405!------------------------
406! I   Initialization
407!    I_a Read the "run.def"
408!------------------------
409write(*,*)
410write(*,*) '********* PEM initialization *********'
411write(*,*) '> Reading "run.def" (PEM)'
412#ifndef CPP_1D
413    dtphys = daysec/48. ! Dummy value (overwritten in phyetat0)
414    call conf_gcm(99,.true.)
415    call infotrac_init
416    nq = nqtot
417    allocate(q(ip1jmp1,llm,nqtot))
418    allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
419#else
420    allocate(q(1,llm,nqtot),pdyn(1))
421    allocate(longitude(1),latitude(1),cell_area(1))
422
423    therestart1D = .false. ! Default value
424    inquire(file = 'start1D.txt',exist = therestart1D)
425    if (.not. therestart1D) then
426        write(*,*) 'There is no "start1D.txt" file!'
427        error stop 'Initialization cannot be done for the 1D PEM.'
428    endif
429    therestartfi = .false. ! Default value
430    inquire(file = 'startfi.nc',exist = therestartfi)
431    if (.not. therestartfi) then
432        write(*,*) 'There is no "startfi.nc" file!'
433        error stop 'Initialization cannot be done for the 1D PEM.'
434    endif
435
436    write(*,*) '> Reading "start1D.txt" and "startfi.nc"'
437    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,         &
438                         time_0,pdyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
439                         play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap)
440    nsplit_phys = 1
441    day_step = steps_per_sol
442#endif
443
444call conf_pem(i_myear,n_myear)
445
446!------------------------
447! I   Initialization
448!    I_b Read of the "start.nc" and "starfi.nc"
449!------------------------
450! I_b.1 Read "start.nc"
451write(*,*) '> Reading "start.nc"'
452allocate(ps_start0(ngrid))
453#ifndef CPP_1D
454    allocate(pdyn(ip1jmp1))
455    call dynetat0(start_name,vcov,ucov,teta,q,masse,pdyn,phis,time_0)
456
457    call gr_dyn_fi(1,iip1,jjp1,ngridmx,pdyn,ps_start0)
458
459    call iniconst ! Initialization of dynamical constants (comconst_mod)
460    call inigeom ! Initialization of geometry
461
462    allocate(ap(nlayer + 1))
463    allocate(bp(nlayer + 1))
464    status = nf90_open(start_name,NF90_NOWRITE,ncid)
465    status = nf90_inq_varid(ncid,"ap",apvarid)
466    status = nf90_get_var(ncid,apvarid,ap)
467    status = nf90_inq_varid(ncid,"bp",bpvarid)
468    status = nf90_get_var(ncid,bpvarid,bp)
469    status = nf90_close(ncid)
470
471    ! Initialization of physics constants and variables (comcstfi_h)
472    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)
473#else
474    ps_start0(1) = pdyn(1)
475#endif
476deallocate(pdyn)
477
478! In the PCM, these values are given to the physic by the dynamic.
479! Here we simply read them in the "startfi.nc" file
480status = nf90_open(startfi_name,NF90_NOWRITE,ncid)
481status = nf90_inq_varid(ncid,"longitude",lonvarid)
482status = nf90_get_var(ncid,lonvarid,longitude)
483status = nf90_inq_varid(ncid,"latitude",latvarid)
484status = nf90_get_var(ncid,latvarid,latitude)
485status = nf90_inq_varid(ncid,"area",areavarid)
486status = nf90_get_var(ncid,areavarid,cell_area)
487status = nf90_inq_varid(ncid,"soildepth",sdvarid)
488status = nf90_get_var(ncid,sdvarid,mlayer)
489status = nf90_close(ncid)
490
491! I_b.2 Read the "startfi.nc"
492! First we read the initial state (starfi.nc)
493#ifndef CPP_STD
494    write(*,*) '> Reading "startfi.nc"'
495    call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
496                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,          &
497                  watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
498
499    ! Remove unphysical values of surface tracer
500    where (qsurf < 0.) qsurf = 0.
501
502    call surfini(ngrid,nslope,qsurf)
503#else
504    call phys_state_var_init(nq)
505    if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
506    call initracer(ngrid,nq)
507    call iniaerosol()
508    allocate(tsurf_read_generic(ngrid))
509    allocate(qsurf_read_generic(ngrid,nq))
510    allocate(tsoil_read_generic(ngrid,nsoilmx))
511    allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
512    allocate(emis_read_generic(ngrid))
513    allocate(albedo_read_generic(ngrid,2))
514    allocate(qsurf(ngrid,nq,1))
515    allocate(tsurf(ngrid,1))
516    allocate(tsoil(ngrid,nsoilmx,1))
517    allocate(emis(ngrid,1))
518    allocate(watercap(ngrid,1))
519    allocate(watercaptag(ngrid))
520    allocate(albedo(ngrid,2,1))
521    allocate(inertiesoil(ngrid,nsoilmx,1))
522    call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
523                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,               &
524                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,        &
525                  rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
526    call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
527
528    nslope = 1
529    call ini_comslope_h(ngrid,1)
530
531    qsurf(:,:,1) = qsurf_read_generic
532    tsurf(:,1) = tsurf_read_generic
533    tsoil(:,:,1) = tsoil_read_generic
534    emis(:,1) = emis_read_generic
535    watercap(:,1) = 0.
536    watercaptag(:) = .false.
537    albedo(:,1,1) = albedo_read_generic(:,1)
538    albedo(:,2,1) = albedo_read_generic(:,2)
539    inertiesoil(:,:,1) = inertiedat
540
541    if (nslope == 1) then
542        def_slope(1) = 0
543        def_slope(2) = 0
544        def_slope_mean = 0
545        subslope_dist(:,1) = 1.
546    endif
547
548    ! Remove unphysical values of surface tracer
549    qsurf(:,:,1) = qsurf_read_generic
550    where (qsurf < 0.) qsurf = 0.
551
552    deallocate(tsurf_read_generic,qsurf_read_generic,qsoil_read_generic,emis_read_generic)
553#endif
554
555do nnq = 1,nqtot  ! Why not using ini_tracer?
556    if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
557    if (noms(nnq) == "h2o_vap") then
558        igcm_h2o_vap = nnq
559        mmol(igcm_h2o_vap) = 18.
560    endif
561    if (noms(nnq) == "co2") igcm_co2 = nnq
562enddo
563
564!------------------------
565! I   Initialization
566!    I_c Subslope parametrisation
567!------------------------
568! Define some slope statistics
569iflat = 1
570do islope = 2,nslope
571    if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
572enddo
573write(*,*) 'Flat slope for islope = ',iflat
574write(*,*) 'Corresponding criterium = ',def_slope_mean(iflat)
575
576!------------------------
577! I   Initialization
578!    I_d Read the PCM data and convert them to the physical grid
579!------------------------
580! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value
581call get_timelen_PCM("data_PCM_Y1.nc",timelen)
582
583allocate(vmr_co2_PCM(ngrid,timelen),q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
584allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid))
585allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
586allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope))
587allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_timeseries(ngrid,nsoilmx,nslope,timelen))
588allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
589
590call read_data_PCM("data_PCM_Y1.nc","data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,ps_avg,tsurf_avg_yr1,tsurf_avg, &
591                   tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2_PEM_phys,q_h2o_PEM_phys,watersurf_density_avg,watersoil_density_timeseries)
592
593! Compute the deviation from the average
594allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
595ps_dev = ps_start0 - ps_avg
596tsurf_dev = tsurf - tsurf_avg
597tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:)
598
599!------------------------
600! I   Initialization
601!    I_e Initialization of the PEM variables and soil
602!------------------------
603call end_comsoil_h_PEM
604call ini_comsoil_h_PEM(ngrid,nslope)
605call end_adsorption_h_PEM
606call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
607call end_ice_table
608call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
609
610allocate(tsoil_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen),watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope),watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
611if (soil_pem) then
612    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
613    tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
614    watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
615    tsoil_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_timeseries
616    do l = nsoilmx + 1,nsoilmx_PEM
617        tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
618        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
619        tsoil_PEM_timeseries(:,l,:,:) = tsoil_timeseries(:,nsoilmx,:,:)
620    enddo
621    watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
622endif !soil_pem
623deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries)
624
625!------------------------
626! I   Initialization
627!    I_f Compute tendencies
628!------------------------
629allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope))
630call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
631call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
632d_co2ice_ini = d_co2ice
633deallocate(min_co2_ice,min_h2o_ice)
634
635!------------------------
636! I   Initialization
637!    I_g Compute global surface pressure
638!------------------------
639total_surface = sum(cell_area)
640ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
641ps_avg_global_new = ps_avg_global_ini
642write(*,*) "Total surface of the planet     =", total_surface
643write(*,*) "Initial global average pressure =", ps_avg_global_ini
644
645!------------------------
646! I   Initialization
647!    I_h Read the "startpem.nc"
648!------------------------
649write(*,*) '> Reading "startpem.nc"'
650allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),layerings_map(ngrid,nslope))
651allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid))
652delta_h2o_icetablesublim = 0.
653call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
654              tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice,            &
655              watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,layerings_map)
656deallocate(tsurf_avg_yr1)
657
658! We save the places where h2o ice is sublimating
659! We compute the surface of h2o ice sublimating
660allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope),co2ice_disappeared(ngrid,nslope))
661co2ice_sublim_surf_ini = 0.
662h2oice_ini_surf = 0.
663is_co2ice_sublim_ini = .false.
664is_h2oice_sublim_ini = .false.
665is_co2ice_ini = .false.
666co2ice_disappeared = .false.
667totmass_co2ice_ini = 0.
668totmass_atmco2_ini = 0.
669if (layering_algo) then
670    do ig = 1,ngrid
671        do islope = 1,nslope
672            if (is_co2ice_str(layerings_map(ig,islope)%top)) then
673                co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
674            else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
675                h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
676            else
677                call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
678            endif
679        enddo
680    enddo
681    ! We put the sublimating tendency coming from subsurface ice into the overall tendency
682    where (zdqsdif_ssi_tot < 0.)
683        d_h2oice = zdqsdif_ssi_tot
684    end where
685endif
686do i = 1,ngrid
687    totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
688    do islope = 1,nslope
689        totmass_co2ice_ini = totmass_co2ice_ini + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)
690        if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
691        if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
692            is_co2ice_sublim_ini(i,islope) = .true.
693            co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
694        endif
695        if (d_h2oice(i,islope) < 0.) then
696            if (h2o_ice(i,islope) > 0.) then
697                is_h2oice_sublim_ini(i,islope) = .true.
698                h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
699            else if (h2o_ice_depth(i,islope) > 0.) then
700                is_h2oice_sublim_ini(i,islope) = .true.
701            endif
702        endif
703    enddo
704enddo
705write(*,*) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini
706write(*,*) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf
707
708totmass_adsco2_ini = 0.
709totmass_adsh2o = 0.
710if (adsorption_pem) then
711    do ig = 1,ngrid
712        do islope = 1,nslope
713            do l = 1,nsoilmx_PEM - 1
714                if (l == 1) then
715                   totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
716                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
717                   totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
718                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
719                else
720                   totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
721                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
722                   totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
723                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
724                endif
725            enddo
726        enddo
727    enddo
728    totmass_adsco2 = totmass_adsco2_ini
729    write(*,*) "Tot mass of CO2 in the regolith =", totmass_adsco2
730    write(*,*) "Tot mass of H2O in the regolith =", totmass_adsh2o
731endif ! adsorption
732
733!------------------------
734! I   Initialization
735!    I_i Compute orbit criterion
736!------------------------
737#ifndef CPP_STD
738    call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
739#else
740    call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
741#endif
742
743n_myear_leg = Max_iter_pem
744if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
745
746!-------------------------- END INITIALIZATION -------------------------
747
748!-------------------------------- RUN ----------------------------------
749!------------------------
750! II  Run
751!    II_a Update pressure, ice and tracers
752!------------------------
753write(*,*)
754write(*,*) '********* PEM cycle *********'
755i_myear_leg = 0
756stopPEM = 0
757if (layering_algo) then
758    allocate(h2o_ice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
759    new_str = .true.
760    new_lag = .true.
761    do islope = 1,nslope
762        do ig = 1,ngrid
763            current(ig,islope)%p => layerings_map(ig,islope)%top
764        enddo
765    enddo
766endif
767
768do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
769! II.a.1. Compute updated global pressure
770    write(*,'(a,f10.2)') ' **** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
771    write(*,*) "> Updating the surface pressure"
772    ps_avg_global_old = ps_avg_global_new
773    do i = 1,ngrid
774        do islope = 1,nslope
775            ps_avg_global_new = ps_avg_global_new - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/total_surface
776        enddo
777    enddo
778    if (adsorption_pem) then
779        do i = 1,ngrid
780            ps_avg_global_new = ps_avg_global_new - g*cell_area(i)*delta_co2_adsorbed(i)/total_surface
781        enddo
782    endif
783    ps_avg = ps_avg*ps_avg_global_new/ps_avg_global_old
784    write(*,*) 'Global average pressure old time step:',ps_avg_global_old
785    write(*,*) 'Global average pressure new time step:',ps_avg_global_new
786
787! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
788    write(*,*) "> Updating the surface pressure timeseries for the new pressure"
789    allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
790    do l = 1,nlayer + 1
791        do ig = 1,ngrid
792            zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
793        enddo
794    enddo
795    ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old
796    write(*,*) "> Updating the pressure levels timeseries for the new pressure"
797    allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
798    do l = 1,nlayer + 1
799        do ig = 1,ngrid
800            zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
801        enddo
802    enddo
803
804! II.a.3. Tracers timeseries
805    write(*,*) "> Updating the tracer VMR timeseries for the new pressure"
806    allocate(vmr_co2_PEM_phys(ngrid,timelen))
807    l = 1
808    do ig = 1,ngrid
809        do t = 1,timelen
810            ! H2O
811            q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ &
812                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
813            if (q_h2o_PEM_phys(ig,t) < 0) then
814                q_h2o_PEM_phys(ig,t) = 1.e-30
815            else if (q_h2o_PEM_phys(ig,t) > 1) then
816                q_h2o_PEM_phys(ig,t) = 1.
817            endif
818            ! CO2
819            q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ &
820                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))                       &
821                                + ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))                       &
822                                -  (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/                     &
823                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
824            if (q_co2_PEM_phys(ig,t) < 0) then
825                q_co2_PEM_phys(ig,t) = 1.e-30
826            else if (q_co2_PEM_phys(ig,t) > 1) then
827                q_co2_PEM_phys(ig,t) = 1.
828            endif
829            mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
830            vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
831        enddo
832    enddo
833    deallocate(zplev_timeseries_new,zplev_timeseries_old)
834
835!------------------------
836! II  Run
837!    II_b Evolution of ice
838!------------------------
839    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
840    if (layering_algo) then
841        h2o_ice_depth_old = h2o_ice_depth
842        do islope = 1,nslope
843            do ig = 1,ngrid
844                call make_layering(layerings_map(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p)
845                co2_ice(ig,islope) = 0.
846                h2o_ice(ig,islope) = 0.
847                h2o_ice_depth(ig,islope) = 0.
848                if (is_co2ice_str(layerings_map(ig,islope)%top)) then
849                    co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
850                else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
851                    h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
852                else
853                    call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
854                endif
855            enddo
856        enddo
857    else
858        zlag = 0.
859        call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
860        call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
861    endif
862
863!------------------------
864! II  Run
865!    II_c Flow of glaciers
866!------------------------
867    allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
868    if (co2ice_flow .and. nslope > 1) then
869        call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
870                              ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
871        if (layering_algo) then
872            do islope = 1,nslope
873                do ig = 1,ngrid
874                    layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
875                enddo
876            enddo
877        endif
878    endif
879    if (h2oice_flow .and. nslope > 1) then
880        call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
881        if (layering_algo) then
882            do islope = 1,nslope
883                do ig = 1,ngrid
884!~                     layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
885                enddo
886            enddo
887        endif
888    endif
889
890!------------------------
891! II  Run
892!    II_d Update surface and soil temperatures
893!------------------------
894! II_d.1 Update Tsurf
895    write(*,*) "> Updating surface temperature"
896    do ig = 1,ngrid
897        do islope = 1,nslope
898            ! CO2 ice disappeared so we look for the closest point without CO2 ice
899            if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then
900                co2ice_disappeared(ig,islope) = .true.
901                if (latitude_deg(ig) > 0.) then ! North hemisphere
902                    outer1: do ig_loop = ig,ngrid ! Go towards equator
903                        do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be equator-ward)
904                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
905                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
906                                exit outer1
907                            endif
908                        enddo
909                        do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be pole-ward)
910                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
911                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
912                                exit outer1
913                            endif
914                        enddo
915                    enddo outer1
916                else ! South hemisphere
917                    outer2: do ig_loop = ig,1,-1 ! Go towards equator
918                        do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be equator-ward)
919                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
920                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
921                                exit outer2
922                            endif
923                        enddo
924                        do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be pole-ward)
925                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
926                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
927                                exit outer2
928                            endif
929                        enddo
930                    enddo outer2
931                endif
932            else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as Tcond CO2
933                call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_ini,ps_avg_global_new,tsurf_avg)
934            endif
935        enddo
936    enddo
937
938    if (soil_pem) then
939! II_d.2 Shifting soil temperature to surface
940        call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
941
942! II_d.3 Update soil temperature
943        write(*,*)"> Updating soil temperature profile"
944        allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen))
945        tsoil_PEM_timeseries_old = tsoil_PEM_timeseries
946        do islope = 1,nslope
947            tsoil_avg_old = tsoil_PEM(:,:,islope)
948            call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
949            call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
950
951            do t = 1,timelen
952                do ig = 1,ngrid
953                    do isoil = 1,nsoilmx_PEM
954                        ! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
955                        tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
956                        ! Update of watersoil density
957                        watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/tsoil_PEM_timeseries(ig,isoil,islope,t) + alpha_clap_h2o)/tsoil_PEM_timeseries(ig,isoil,islope,t)*mmol(igcm_h2o_vap)/(mugaz*r)
958                        if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
959                    enddo
960                enddo
961            enddo
962        enddo
963        watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
964        deallocate(tsoil_avg_old)
965
966! II_d.4 Update the ice table
967        allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope))
968        if (icetable_equilibrium) then
969            write(*,*) "> Updating ice table (equilibrium method)"
970            icetable_thickness_old = icetable_thickness
971            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
972            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
973        else if (icetable_dynamic) then
974            write(*,*) "> Updating ice table (dynamic method)"
975            ice_porefilling_old = ice_porefilling
976            icetable_depth_old = icetable_depth
977            allocate(porefill(nsoilmx_PEM))
978            do ig = 1,ngrid
979                do islope = 1,nslope
980                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps_avg(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth)
981                    icetable_depth(ig,islope) = ssi_depth
982                    ice_porefilling(ig,:,islope) = porefill
983                enddo
984            enddo
985            deallocate(porefill)
986            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
987        endif
988        deallocate(icetable_thickness_old,ice_porefilling_old)
989
990! II_d.5 Update the soil thermal properties
991        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,ps_avg_global_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
992
993! II_d.6 Update the mass of the regolith adsorbed
994        totmass_adsco2 = 0.
995        totmass_adsh2o = 0.
996        if (adsorption_pem) then
997            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
998                                     tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,       &
999                                     h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
1000            do ig = 1,ngrid
1001                do islope = 1,nslope
1002                    do l = 1,nsoilmx_PEM
1003                        if (l == 1) then
1004                            totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
1005                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
1006                            totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
1007                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
1008                        else
1009                            totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
1010                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
1011                            totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
1012                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
1013                        endif
1014                    enddo
1015                enddo
1016            enddo
1017            write(*,*) "Total mass of CO2 in the regolith =", totmass_adsco2
1018            write(*,*) "Total mass of H2O in the regolith =", totmass_adsh2o
1019        endif
1020    endif !soil_pem
1021    deallocate(zshift_surf,zlag)
1022
1023!------------------------
1024! II  Run
1025!    II_e Outputs
1026!------------------------
1027    call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
1028    do islope = 1,nslope
1029        write(str2(1:2),'(i2.2)') islope
1030        call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
1031        call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
1032        call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
1033        call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
1034        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope)))
1035        call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope)))
1036        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
1037        if (icetable_equilibrium) then
1038            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
1039            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope))
1040        else if (icetable_dynamic) then
1041            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
1042        endif
1043
1044        if (soil_pem) then
1045            call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
1046            call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
1047            if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
1048            if (adsorption_pem) then
1049                call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope))
1050                call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope))
1051            endif
1052        endif
1053    enddo
1054    deallocate(flag_co2flow,flag_h2oflow)
1055
1056    ! Checking mass balance for CO2
1057    totmass_co2ice = 0.
1058    totmass_atmco2 = 0.
1059    do ig = 1,ngrid
1060        totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
1061        do islope = 1,nslope
1062            totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)*cell_area(ig)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
1063        enddo
1064    enddo
1065    write(*,'(a,f8.3,a)') " > Relative total CO2 mass balance = ", 100.*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini), ' %'
1066    if ((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini) > 0.01) then
1067        write(*,*) '  /!\ Warning: mass balance is not conseved!'
1068        write(*,'(a,f8.3,a)') '       Atmospheric CO2 mass balance = ', 100.*(totmass_atmco2 - totmass_atmco2_ini)/totmass_atmco2_ini, ' %'
1069        write(*,'(a,f8.3,a)') '       CO2 ice mass balance         = ', 100.*(totmass_co2ice - totmass_co2ice_ini)/totmass_co2ice_ini, ' %'
1070        write(*,'(a,f8.3,a)') '       Adsorbed CO2 mass balance    = ', 100.*(totmass_adsco2 - totmass_adsco2_ini)/totmass_adsco2_ini, ' %'
1071    endif
1072
1073!------------------------
1074! II  Run
1075!    II_f Update the tendencies
1076!------------------------
1077    call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new)
1078    write(*,*) "> Updating the H2O sub-surface ice tendency due to lag layer"
1079    if (layering_algo) then
1080        do ig = 1,ngrid
1081            do islope = 1,nslope
1082                if (is_h2oice_sublim_ini(ig,islope) .and. h2o_ice_depth(ig,islope) > 0.) call recomp_tend_h2o(h2o_ice_depth_old(ig,islope),h2o_ice_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope))
1083            enddo
1084        enddo
1085!~     else
1086!~         do ig = 1,ngrid
1087!~             do islope = 1,nslope
1088!~                 call recomp_tend_h2o(icetable_depth_old(ig,islope),icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope))
1089!~             enddo
1090!~         enddo
1091    endif
1092    if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old)
1093    deallocate(vmr_co2_PEM_phys)
1094
1095!------------------------
1096! II  Run
1097!    II_g Checking the stopping criterion
1098!------------------------
1099    write(*,*) "> Checking the stopping criteria"
1100    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
1101    call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope)
1102    i_myear_leg = i_myear_leg + dt
1103    i_myear = i_myear + dt
1104    if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
1105    if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
1106    call system_clock(c2)
1107    if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
1108    if (stopPEM > 0) then
1109        select case (stopPEM)
1110            case(1)
1111                write(*,'(a,i0,a)') " **** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
1112            case(2)
1113                write(*,'(a,i0,a)') " **** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
1114            case(3)
1115                write(*,'(a,i0,a)') " **** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
1116            case(4)
1117                write(*,'(a,i0,a)') " **** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
1118            case(5)
1119                write(*,'(a,i0)') " **** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
1120            case(6)
1121                write(*,'(a,i0)') " **** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
1122            case(7)
1123                write(*,'(a,i0)') " **** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
1124            case default
1125                write(*,'(a,i0)') " **** STOPPING with unknown stopping criterion code: ", stopPEM
1126        end select
1127        exit
1128    else
1129        write(*,'(a,f10.2,a)') ' **** The chained simulation has run for ',i_myear,' Martian years.'
1130        write(*,*) '**** The PEM can continue!'
1131        write(*,*) '****'
1132    endif
1133enddo ! big time iteration loop of the pem
1134deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed)
1135deallocate(watersoil_density_PEM_avg,watersurf_density_avg)
1136deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries)
1137deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim)
1138deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
1139deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
1140if (layering_algo) deallocate(h2o_ice_depth_old,new_str,new_lag,current)
1141!------------------------------ END RUN --------------------------------
1142
1143!------------------------------- OUTPUT --------------------------------
1144!------------------------
1145! III Output
1146!    III_a Update surface values for the PCM start files
1147!------------------------
1148write(*,*)
1149write(*,*) '********* PEM finalization *********'
1150! III_a.1 Ice update for start file
1151write(*,*) '> Reconstructing perennial ice and frost for the PCM'
1152watercap = 0.
1153perennial_co2ice = co2_ice
1154do ig = 1,ngrid
1155    ! H2O ice metamorphism
1156    !if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
1157    !    h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1158    !    qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
1159    !endif
1160
1161    ! Is H2O ice still considered as an infinite reservoir for the PCM?
1162    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
1163        ! There is enough ice to be considered as an infinite reservoir
1164        watercaptag(ig) = .true.
1165    else
1166        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
1167        watercaptag(ig) = .false.
1168        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1169        h2o_ice(ig,:) = 0.
1170    endif
1171
1172    ! CO2 ice metamorphism
1173    !if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
1174    !    perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1175    !    qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
1176    !endif
1177enddo
1178
1179! III.a.3. Tsurf update for start file
1180write(*,*) '> Reconstructing the surface temperature for the PCM'
1181tsurf = tsurf_avg + tsurf_dev
1182deallocate(tsurf_dev)
1183
1184! III_a.4 Tsoil update for start file
1185if (soil_pem) then
1186    write(*,*) '> Reconstructing the soil temperature profile for the PCM'
1187    inertiesoil = TI_PEM(:,:nsoilmx,:)
1188    ! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value
1189    do isoil = 1,nsoilmx
1190        tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:)
1191    enddo
1192    tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev
1193#ifndef CPP_STD
1194    flux_geo = fluxgeo
1195#endif
1196endif
1197deallocate(tsurf_avg,tsoil_dev)
1198
1199! III_a.5 Pressure update for start file
1200write(*,*) '> Reconstructing the pressure for the PCM'
1201allocate(ps_start(ngrid))
1202! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop
1203ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini
1204deallocate(ps_avg,ps_dev)
1205
1206! III_a.6 Tracers update for start file
1207write(*,*) '> Reconstructing the tracer VMR for the PCM'
1208allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
1209do l = 1,nlayer + 1
1210    zplev_start0(:,l) = ap(l) + bp(l)*ps_start0
1211    zplev_new(:,l) = ap(l) + bp(l)*ps_start
1212enddo
1213
1214do nnq = 1,nqtot
1215    if (noms(nnq) /= "co2") then
1216        do l = 1,llm - 1
1217            do ig = 1,ngrid
1218                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
1219            enddo
1220            q(:,llm,nnq) = q(:,llm - 1,nnq)
1221        enddo
1222    else
1223        do l = 1,llm - 1
1224            do ig = 1,ngrid
1225                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
1226                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_start0(ig,l) - zplev_start0(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
1227            enddo
1228            q(:,llm,nnq) = q(:,llm - 1,nnq)
1229        enddo
1230    endif
1231enddo
1232deallocate(zplev_start0)
1233
1234! Conserving the tracers mass for start file
1235do nnq = 1,nqtot
1236    do ig = 1,ngrid
1237        do l = 1,llm - 1
1238            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
1239                extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1240                q(ig,l,nnq) = 1.
1241                q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1242                write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1243            endif
1244            if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1245        enddo
1246    enddo
1247enddo
1248deallocate(zplev_new)
1249
1250! III_a.7 Albedo update for start file
1251write(*,*) '> Reconstructing the albedo for the PCM'
1252do ig = 1,ngrid
1253    if (latitude(ig) < 0.) then
1254        icap = 2 ! Southern hemisphere
1255    else
1256        icap = 1 ! Northern hemisphere
1257    endif
1258    do islope = 1,nslope
1259        ! Bare ground
1260        albedo(ig,:,islope) = albedodat(ig)
1261        emis(ig,islope) = emissiv
1262
1263        ! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant
1264        ! H2O ice
1265        if (h2o_ice(ig,islope) > 0.) then
1266            albedo(ig,:,islope) = albedo_h2o_cap
1267            emis(ig,islope) = 1.
1268        endif
1269        ! CO2 ice
1270        if (co2_ice(ig,islope) > 0.) then
1271            albedo(ig,:,islope) = albedo_perennialco2(icap)
1272            emis(ig,islope) = emisice(icap)
1273        endif
1274        ! H2O frost
1275        if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
1276            albedo(ig,:,islope) = albedo_h2o_frost
1277            emis(ig,islope) = 1.
1278        endif
1279        ! CO2 frost
1280        if (qsurf(ig,igcm_co2,islope) > 0.) then
1281            albedo(ig,:,islope) = albedice(icap)
1282            emis(ig,islope) = emisice(icap)
1283        endif
1284    enddo
1285enddo
1286
1287! III_a.8 Orbital parameters update for start file
1288write(*,*) '> Setting the new orbital parameters'
1289if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
1290
1291!------------------------
1292! III Output
1293!    III_b Write "restart.nc" and "restartfi.nc"
1294!------------------------
1295! III_b.1 Write "restart.nc"
1296ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1297pday = day_ini
1298ztime_fin = time_phys
1299#ifndef CPP_1D
1300    write(*,*) '> Writing "restart.nc"'
1301    ! Correction on teta due to surface pressure changes
1302    allocate(pdyn(ip1jmp1))
1303    call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn)
1304    do i = 1,ip1jmp1
1305        teta(i,:) = teta(i,:)*pdyn(i)**rcp
1306    enddo
1307    ! Correction on atmospheric pressure
1308    allocate(p(ip1jmp1,nlayer + 1))
1309    call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start,pdyn)
1310    call pression(ip1jmp1,ap,bp,pdyn,p)
1311    ! Correction on the mass of atmosphere
1312    call massdair(p,masse)
1313    call dynredem0("restart.nc",day_ini,phis)
1314    call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,pdyn)
1315    deallocate(ap,bp,p,pdyn)
1316#else
1317    write(*,*) '> Writing "restart1D.txt"'
1318    call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1319#endif
1320deallocate(ps_start0,ps_start)
1321
1322! III_b.2 Write the "restartfi.nc"
1323write(*,*) '> Writing "restartfi.nc"'
1324#ifndef CPP_STD
1325    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1326                  nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1327                  inertiedat,def_slope,subslope_dist)
1328    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,      &
1329                  ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,        &
1330                  albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1331                  wstar,watercap,perennial_co2ice)
1332#else
1333    if (allocated(noms)) deallocate(noms)
1334    deallocate(qsurf,tsurf,tsoil,emis,watercap,watercaptag,albedo,inertiesoil)
1335    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1336                  nlayer,nq,ptimestep,pday,time_phys,cell_area,    &
1337                  albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1338    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,       &
1339                  ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1340                  cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab,   &
1341                  tsea_ice,sea_ice)
1342#endif
1343
1344!------------------------
1345! III Output
1346!    III_c Write the "restartpem.nc"
1347!------------------------
1348write(*,*) '> Writing "restartpem.nc"'
1349if (layering_algo) nb_str_max = get_nb_str_max(layerings_map,ngrid,nslope) ! Get the maximum number of "stratum" in the layerings_mapication (layerings)
1350call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1351call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
1352             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
1353
1354call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
1355
1356write(*,*)
1357write(*,*) '****** PEM final information *******'
1358write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
1359write(*,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years."
1360write(*,'(a,f16.4,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years."
1361write(*,*) "+ PEM: so far, so good!"
1362write(*,*) '************************************'
1363
1364if (layering_algo) then
1365    do islope = 1,nslope
1366        do i = 1,ngrid
1367            call del_layering(layerings_map(i,islope))
1368        enddo
1369    enddo
1370endif
1371deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
1372deallocate(co2_ice,h2o_ice,layerings_map)
1373!----------------------------- END OUTPUT ------------------------------
1374
1375END PROGRAM pem
Note: See TracBrowser for help on using the repository browser.