Context Navigation

pem.F90 @ 3790

Last change on this file since 3790 was 3789, checked in by jbclement, 3 weeks ago
PEM: Few corrections for the layering algorithm, in particular when a dust lag layer is created. JBC
File size: 71.4 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
28	PROGRAM pem
29
30	use phyetat0_mod, only: phyetat0
31	use phyredem, only: physdem0, physdem1
32	use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
33	use turb_mod, only: q2, wstar
34	use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h
35	use logic_mod, only: iflag_phys
36	use mod_const_mpi, only: COMM_LMDZ
37	use infotrac
38	use geometry_mod, only: latitude_deg
39	use conf_pem_mod, only: conf_pem
40	use pemredem, only: pemdem0, pemdem1
41	use 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
43	use stopping_crit_mod, only: stopping_crit_h2o_ice, stopping_crit_co2
44	use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
45	use evol_ice_mod, only: evol_co2_ice, evol_h2o_ice
46	use 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
50	use 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
53	use time_evol_mod, only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
54	use orbit_param_criterion_mod, only: orbit_param_criterion
55	use recomp_orb_param_mod, only: recomp_orb_param
56	use 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
58	use soil_thermalproperties_mod, only: update_soil_thermalproperties
59	use time_phylmdz_mod, only: daysec, dtphys
60	use abort_pem_mod, only: abort_pem
61	use soil_settings_PEM_mod, only: soil_settings_PEM
62	use compute_tend_mod, only: compute_tend
63	use info_PEM_mod, only: info_PEM
64	use get_timelen_PCM_mod, only: get_timelen_PCM
65	use pemetat0_mod, only: pemetat0
66	use read_data_PCM_mod, only: read_data_PCM
67	use recomp_tend_mod, only: recomp_tend_co2, recomp_tend_h2o
68	use compute_soiltemp_mod, only: compute_tsoil_pem, shift_tsoil2surf
69	use writediagpem_mod, only: writediagpem, writediagsoilpem
70	use co2condens_mod, only: CO2cond_ps
71	use 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
74	use dyn_ss_ice_m_mod, only: dyn_ss_ice_m
75	use version_info_mod, only: print_version_info
76	use 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
116	implicit none
117
118	include "dimensions.h"
119	include "paramet.h"
120	include "comgeom.h"
121	include "iniprint.h"
122
123	! Same variable names as in the PCM
124	integer, parameter :: ngridmx = 2 + (jjm - 1)*iim - 1/jjm
125	integer, parameter :: nlayer = llm ! Number of vertical layer
126	integer :: ngrid ! Number of physical grid points
127	integer :: nq ! Number of tracer
128	integer :: day_ini ! First day of the simulation
129	real :: pday ! Physical day
130	real :: time_phys ! Same as in PCM
131	real :: ptimestep ! Same as in PCM
132	real :: ztime_fin ! Same as in PCM
133
134	! Variables to read "start.nc"
135	character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM
136
137	! Dynamic variables
138	real, dimension(ip1jm,llm) :: vcov ! vents covariants
139	real, dimension(ip1jmp1,llm) :: ucov ! vents covariants
140	real, dimension(ip1jmp1,llm) :: teta ! Potential temperature
141	real, dimension(:,:,:), allocatable :: q ! champs advectes
142	real, dimension(:), allocatable :: pdyn ! pressure for the dynamic grid
143	real, dimension(:), allocatable :: ps_start ! surface pressure in the start file
144	real, dimension(:), allocatable :: ps_start0 ! surface pressure in the start file at the beginning
145	real, dimension(:), allocatable :: ps_avg ! (ngrid) Average surface pressure
146	real, dimension(:), allocatable :: ps_dev ! (ngrid x timelen) Surface pressure deviation
147	real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure
148	real, dimension(ip1jmp1,llm) :: masse ! Air mass
149	real, dimension(ip1jmp1) :: phis ! geopotentiel au sol
150	real :: time_0
151
152	! Variables to read starfi.nc
153	character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM
154	character(2) :: str2
155	integer :: ncid, status ! Variable for handling opening of files
156	integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files
157	integer :: apvarid, bpvarid ! Variable ID for Netcdf files
158
159	! Variables to read starfi.nc and write restartfi.nc
160	real, dimension(:), allocatable :: longitude ! Longitude read in startfi_name and written in restartfi
161	real, dimension(:), allocatable :: latitude ! Latitude read in startfi_name and written in restartfi
162	real, dimension(:), allocatable :: cell_area ! Cell_area read in startfi_name and written in restartfi
163	real :: total_surface ! Total surface of the planet
164
165	! Variables for h2o ice evolution
166	real, dimension(:,:), allocatable :: h2o_ice ! h2o ice in the PEM
167	real, dimension(:,:), allocatable :: d_h2oice ! physical point x slope field: Tendency of evolution of perennial h2o ice
168	real, dimension(:,:,:), allocatable :: min_h2o_ice ! Minima of h2o ice at each point for the PCM years [kg/m^2]
169	real :: h2oice_ini_surf ! Initial surface of sublimating h2o ice
170	logical, dimension(:,:), allocatable :: is_h2oice_sublim_ini ! Logical array to know if h2o ice is sublimating
171	real :: ps_avg_global_ini ! constant: Global average pressure at initialization [Pa]
172	real :: ps_avg_global_old ! constant: Global average pressure of previous time step
173	real :: ps_avg_global_new ! constant: Global average pressure of current time step
174	real, dimension(:,:), allocatable :: zplev_new ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
175	real, dimension(:,:), allocatable :: zplev_start0 ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
176	real, dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
177	real, dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
178	integer :: 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
179	real :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg]
180	real, 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]
181	integer :: timelen ! # time samples
182	real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1
183
184	! Variables for co2 ice evolution
185	real, dimension(:,:), allocatable :: co2_ice ! co2 ice in the PEM
186	real, dimension(:,:), allocatable :: d_co2ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
187	real, dimension(:,:), allocatable :: d_co2ice_ini ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
188	logical, dimension(:,:), allocatable :: is_co2ice_ini ! Was there co2 ice initially in the PEM?
189	real, dimension(:,:,:), allocatable :: min_co2_ice ! Minimum of co2 ice at each point for the first year [kg/m^2]
190	real :: co2ice_sublim_surf_ini ! Initial surface of sublimating co2 ice
191	logical, dimension(:,:), allocatable :: is_co2ice_sublim_ini ! Logical array to know if co2 ice is sublimating
192	real, dimension(:,:), allocatable :: vmr_co2_PCM ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
193	real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Grid points x Times co2 volume mixing ratio used in the PEM
194	real, 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]
195	real(kind = 16) :: totmass_co2ice, totmass_atmco2 ! Current total CO2 masses
196	real(kind = 16) :: totmass_co2ice_ini, totmass_atmco2_ini ! Initial total CO2 masses
197
198	! Variables for the evolution of layered layerings_map
199	type(layering), dimension(:,:), allocatable :: layerings_map ! Layering for each grid point and slope
200	type(ptrarray), dimension(:,:), allocatable :: current ! Current active stratum in the layering
201	logical, dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm
202	real, dimension(:,:), allocatable :: h2o_ice_depth_old ! Old depth of subsurface ice layer
203
204	! Variables for slopes
205	integer(kind = 1), dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow
206	integer(kind = 1), dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow
207
208	! Variables for surface and soil
209	real, dimension(:,:), allocatable :: tsurf_avg ! Grid points x Slope field: Average surface temperature [K]
210	real, dimension(:,:), allocatable :: tsurf_dev ! Grid points x Slope field: Surface temperature deviation [K]
211	real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Grid points x Slope field: Average surface temperature of first call of the PCM [K]
212	real, dimension(:,:,:), allocatable :: tsoil_avg ! Grid points x Soil x Slope field: Average Soil Temperature [K]
213	real, dimension(:,:), allocatable :: tsoil_avg_old ! Grid points x Soil field: Average Soil Temperature at the previous time step [K]
214	real, dimension(:,:,:), allocatable :: tsoil_dev ! Grid points x Soil x Slope field: Soil temperature deviation [K]
215	real, dimension(:,:,:,:), allocatable :: tsoil_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
216	real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K]
217	real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K]
218	real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
219	real, dimension(:,:), allocatable :: watersurf_density_avg ! Grid points x Slope: Average water surface density [kg/m^3]
220	real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
221	real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Grid points x Soil x Slopes: Average water soil density [kg/m^3]
222	real, dimension(:), allocatable :: delta_co2_adsorbed ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
223	real, dimension(:), allocatable :: delta_h2o_adsorbed ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
224	real(kind = 16) :: totmass_adsco2, totmass_adsco2_ini ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
225	real :: totmass_adsh2o ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
226	logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
227	real, dimension(:,:), allocatable :: icetable_thickness_old ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
228	real, dimension(:,:,:), allocatable :: ice_porefilling_old ! ngrid x nslope: Ice pore filling at the previous iteration [m]
229	real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
230	real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table
231	real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table
232	real, dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m]
233	real, dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m]
234	real, dimension(:,:), allocatable :: icetable_depth_old ! Old depth of the ice table
235
236	! Some variables for the PEM run
237	real, parameter :: year_step = 1 ! Timestep for the pem
238	real :: i_myear_leg ! Number of iteration
239	real :: n_myear_leg ! Maximum number of iterations before stopping
240	real :: i_myear ! Global number of Martian years of the chained simulations
241	real :: n_myear ! Maximum number of Martian years of the chained simulations
242	real :: timestep ! Timestep [s]
243	character(100) :: arg ! To read command-line arguments program was invoked
244	logical :: timewall ! Flag to use the time limit stopping criterion in case of a PEM job
245	integer(kind = 8) :: cr ! Number of clock ticks per second (count rate)
246	integer(kind = 8) :: c1, c2 ! Counts of processor clock
247	character(100) :: chtimelimit ! Time limit for the PEM job outputted by the SLURM command
248	real :: timelimit ! Time limit for the PEM job in seconds
249	real, 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)
250	integer :: cstat, days, hours, minutes, seconds
251	character(1) :: sep
252	character(8) :: date
253	character(10) :: time
254	character(5) :: zone
255	integer, dimension(8) :: values
256	character(128) :: dir = ' '
257	character(32) :: logname = ' '
258	character(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
299	integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
300	logical :: num_str
301
302	! CODE
303	write(,) ' * . . + . * . + . . . '
304	write(,) ' + _______ ________ ____ ____ * + '
305	write(,) ' + . * \|_ __ \\|_ __ \|\|_ \ / _\| . *'
306	write(,) ' . . \| \|__) \| \| \|_ \_\| \| \/ \| * * . '
307	write(,) ' . \| ___/ \| _\| _ \| \|\ /\| \| . . '
308	write(,) '. * * _\| \|_ _\| \|__/ \| _\| \|_\/_\| \|_ * '
309	write(,) ' + \|_____\| \|________\|\|_____\|\|_____\| + . '
310	write(,) ' . * . * . + * . + .'
311
312	! Elapsed time with system clock
313	call system_clock(count_rate = cr)
314	call system_clock(c1)
315	timewall = .true.
316	timelimit = 86400 ! 86400 seconds = 24 h by default
317	if (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 = days86400 + hours3600 + 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 = hours3600 + minutes60 + 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
368	else
369	timewall = .false.
370	endif
371
372	! Some user info
373	call date_and_time(date,time,zone,values)
374	call getcwd(dir) ! Current directory
375	call getlog(logname) ! User name
376	call hostnm(hostname) ! Machine/station name
377	write(,)
378	write(,) '******* PEM information *******'
379	write(,) '+ User : '//trim(logname)
380	write(,) '+ Machine : '//trim(hostname)
381	write(,) '+ Directory: '//trim(dir)
382	write(*,'(a,i2,a,i2,a,i4)') ' + Date : ',values(3),'/',values(2),'/',values(1)
383	write(*,'(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
395	day_ini = 0
396	time_phys = 0.
397	ngrid = ngridmx
398	A = (1./m_co2 - 1./m_noco2)
399	B = 1./m_noco2
400	year_day = 669
401	daysec = 88775.
402	timestep = year_daydaysecyear_step
403
404	!----------------------------- INITIALIZATION --------------------------
405	!------------------------
406	! I Initialization
407	! I_a Read the "run.def"
408	!------------------------
409	write(,)
410	write(,) '******* PEM initialization *******'
411	write(,) '> 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
444	call 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"
451	write(,) '> Reading "start.nc"'
452	allocate(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
476	deallocate(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
480	status = nf90_open(startfi_name,NF90_NOWRITE,ncid)
481	status = nf90_inq_varid(ncid,"longitude",lonvarid)
482	status = nf90_get_var(ncid,lonvarid,longitude)
483	status = nf90_inq_varid(ncid,"latitude",latvarid)
484	status = nf90_get_var(ncid,latvarid,latitude)
485	status = nf90_inq_varid(ncid,"area",areavarid)
486	status = nf90_get_var(ncid,areavarid,cell_area)
487	status = nf90_inq_varid(ncid,"soildepth",sdvarid)
488	status = nf90_get_var(ncid,sdvarid,mlayer)
489	status = 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
555	do 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
562	enddo
563
564	!------------------------
565	! I Initialization
566	! I_c Subslope parametrisation
567	!------------------------
568	! Define some slope statistics
569	iflat = 1
570	do islope = 2,nslope
571	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
572	enddo
573	write(,) 'Flat slope for islope = ',iflat
574	write(,) '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
581	call get_timelen_PCM("data_PCM_Y1.nc",timelen)
582
583	allocate(vmr_co2_PCM(ngrid,timelen),q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
584	allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid))
585	allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
586	allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope))
587	allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_timeseries(ngrid,nsoilmx,nslope,timelen))
588	allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
589
590	call 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
594	allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
595	ps_dev = ps_start0 - ps_avg
596	tsurf_dev = tsurf - tsurf_avg
597	tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:)
598
599	!------------------------
600	! I Initialization
601	! I_e Initialization of the PEM variables and soil
602	!------------------------
603	call end_comsoil_h_PEM
604	call ini_comsoil_h_PEM(ngrid,nslope)
605	call end_adsorption_h_PEM
606	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
607	call end_ice_table
608	call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
609
610	allocate(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))
611	if (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
622	endif !soil_pem
623	deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries)
624
625	!------------------------
626	! I Initialization
627	! I_f Compute tendencies
628	!------------------------
629	allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope))
630	call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
631	call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
632	d_co2ice_ini = d_co2ice
633	deallocate(min_co2_ice,min_h2o_ice)
634
635	!------------------------
636	! I Initialization
637	! I_g Compute global surface pressure
638	!------------------------
639	total_surface = sum(cell_area)
640	ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
641	ps_avg_global_new = ps_avg_global_ini
642	write(,) "Total surface of the planet =", total_surface
643	write(,) "Initial global average pressure =", ps_avg_global_ini
644
645	!------------------------
646	! I Initialization
647	! I_h Read the "startpem.nc"
648	!------------------------
649	write(,) '> Reading "startpem.nc"'
650	allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),layerings_map(ngrid,nslope))
651	allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid))
652	delta_h2o_icetablesublim = 0.
653	call 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)
656	deallocate(tsurf_avg_yr1)
657
658	! We save the places where h2o ice is sublimating
659	! We compute the surface of h2o ice sublimating
660	allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope),co2ice_disappeared(ngrid,nslope))
661	co2ice_sublim_surf_ini = 0.
662	h2oice_ini_surf = 0.
663	is_co2ice_sublim_ini = .false.
664	is_h2oice_sublim_ini = .false.
665	is_co2ice_ini = .false.
666	co2ice_disappeared = .false.
667	totmass_co2ice_ini = 0.
668	totmass_atmco2_ini = 0.
669	if (layering_algo) then
670	do ig = 1,ngrid
671	do islope = 1,nslope
672	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
673	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
674	enddo
675	enddo
676	! We put the sublimating tendency coming from subsurface ice into the overall tendency
677	where (zdqsdif_ssi_tot < 0.)
678	d_h2oice = zdqsdif_ssi_tot
679	end where
680	endif
681	do i = 1,ngrid
682	totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
683	do islope = 1,nslope
684	totmass_co2ice_ini = totmass_co2ice_ini + co2_ice(i,islope)cell_area(i)subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)
685	if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
686	if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
687	is_co2ice_sublim_ini(i,islope) = .true.
688	co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
689	endif
690	if (d_h2oice(i,islope) < 0.) then
691	if (h2o_ice(i,islope) > 0.) then
692	is_h2oice_sublim_ini(i,islope) = .true.
693	h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
694	else if (h2o_ice_depth(i,islope) > 0.) then
695	is_h2oice_sublim_ini(i,islope) = .true.
696	endif
697	endif
698	enddo
699	enddo
700	write(,) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini
701	write(,) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf
702
703	totmass_adsco2_ini = 0.
704	totmass_adsh2o = 0.
705	if (adsorption_pem) then
706	do ig = 1,ngrid
707	do islope = 1,nslope
708	do l = 1,nsoilmx_PEM - 1
709	if (l == 1) then
710	totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
711	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
712	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
713	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
714	else
715	totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
716	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
717	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
718	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
719	endif
720	enddo
721	enddo
722	enddo
723	totmass_adsco2 = totmass_adsco2_ini
724	write(,) "Tot mass of CO2 in the regolith =", totmass_adsco2
725	write(,) "Tot mass of H2O in the regolith =", totmass_adsh2o
726	endif ! adsorption
727
728	!------------------------
729	! I Initialization
730	! I_i Compute orbit criterion
731	!------------------------
732	#ifndef CPP_STD
733	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
734	#else
735	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
736	#endif
737
738	n_myear_leg = Max_iter_pem
739	if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
740
741	!-------------------------- END INITIALIZATION -------------------------
742
743	!-------------------------------- RUN ----------------------------------
744	!------------------------
745	! II Run
746	! II_a Update pressure, ice and tracers
747	!------------------------
748	write(,)
749	write(,) '******* PEM cycle *******'
750	i_myear_leg = 0
751	stopPEM = 0
752	if (layering_algo) then
753	allocate(h2o_ice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
754	new_str = .true.
755	new_lag = .true.
756	do islope = 1,nslope
757	do ig = 1,ngrid
758	current(ig,islope)%p => layerings_map(ig,islope)%top
759	enddo
760	enddo
761	endif
762
763	do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
764	! II.a.1. Compute updated global pressure
765	write(,'(a,f10.2)') ' *** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
766	write(,) "> Updating the surface pressure"
767	ps_avg_global_old = ps_avg_global_new
768	do i = 1,ngrid
769	do islope = 1,nslope
770	ps_avg_global_new = ps_avg_global_new - CO2cond_psgcell_area(i)d_co2ice(i,islope)subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/total_surface
771	enddo
772	enddo
773	if (adsorption_pem) then
774	do i = 1,ngrid
775	ps_avg_global_new = ps_avg_global_new - gcell_area(i)delta_co2_adsorbed(i)/total_surface
776	enddo
777	endif
778	ps_avg = ps_avg*ps_avg_global_new/ps_avg_global_old
779	write(,) 'Global average pressure old time step:',ps_avg_global_old
780	write(,) 'Global average pressure new time step:',ps_avg_global_new
781
782	! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
783	write(,) "> Updating the surface pressure timeseries for the new pressure"
784	allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
785	do l = 1,nlayer + 1
786	do ig = 1,ngrid
787	zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
788	enddo
789	enddo
790	ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old
791	write(,) "> Updating the pressure levels timeseries for the new pressure"
792	allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
793	do l = 1,nlayer + 1
794	do ig = 1,ngrid
795	zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
796	enddo
797	enddo
798
799	! II.a.3. Tracers timeseries
800	write(,) "> Updating the tracer VMR timeseries for the new pressure"
801	allocate(vmr_co2_PEM_phys(ngrid,timelen))
802	l = 1
803	do ig = 1,ngrid
804	do t = 1,timelen
805	! H2O
806	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))/ &
807	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
808	if (q_h2o_PEM_phys(ig,t) < 0) then
809	q_h2o_PEM_phys(ig,t) = 1.e-30
810	else if (q_h2o_PEM_phys(ig,t) > 1) then
811	q_h2o_PEM_phys(ig,t) = 1.
812	endif
813	! CO2
814	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))/ &
815	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
816	+ ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
817	- (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/ &
818	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
819	if (q_co2_PEM_phys(ig,t) < 0) then
820	q_co2_PEM_phys(ig,t) = 1.e-30
821	else if (q_co2_PEM_phys(ig,t) > 1) then
822	q_co2_PEM_phys(ig,t) = 1.
823	endif
824	mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
825	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
826	enddo
827	enddo
828	deallocate(zplev_timeseries_new,zplev_timeseries_old)
829
830	!------------------------
831	! II Run
832	! II_b Evolution of ice
833	!------------------------
834	allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
835	if (layering_algo) then
836	h2o_ice_depth_old = h2o_ice_depth
837	do islope = 1,nslope
838	do ig = 1,ngrid
839	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)
840	co2_ice(ig,islope) = 0.
841	h2o_ice(ig,islope) = 0.
842	h2o_ice_depth(ig,islope) = 0.
843	if (is_co2ice_str(layerings_map(ig,islope)%top)) then
844	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
845	else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
846	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
847	else
848	call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope))
849	endif
850	enddo
851	enddo
852	call print_layering(layerings_map(1,1))
853	else
854	zlag = 0.
855	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
856	call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
857	endif
858
859	!------------------------
860	! II Run
861	! II_c Flow of glaciers
862	!------------------------
863	allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
864	if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
865	ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
866	if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
867	if (layering_algo) then
868	do islope = 1,nslope
869	do ig = 1,ngrid
870	layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
871	layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
872	enddo
873	enddo
874	endif
875
876	!------------------------
877	! II Run
878	! II_d Update surface and soil temperatures
879	!------------------------
880	! II_d.1 Update Tsurf
881	write(,) "> Updating surface temperature"
882	do ig = 1,ngrid
883	do islope = 1,nslope
884	! CO2 ice disappeared so we look for the closest point without CO2 ice
885	if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then
886	co2ice_disappeared(ig,islope) = .true.
887	if (latitude_deg(ig) > 0.) then ! North hemisphere
888	outer1: do ig_loop = ig,ngrid ! Go towards equator
889	do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be equator-ward)
890	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
891	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
892	exit outer1
893	endif
894	enddo
895	do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be pole-ward)
896	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
897	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
898	exit outer1
899	endif
900	enddo
901	enddo outer1
902	else ! South hemisphere
903	outer2: do ig_loop = ig,1,-1 ! Go towards equator
904	do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be equator-ward)
905	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
906	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
907	exit outer2
908	endif
909	enddo
910	do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be pole-ward)
911	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
912	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
913	exit outer2
914	endif
915	enddo
916	enddo outer2
917	endif
918	else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as Tcond CO2
919	call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_ini,ps_avg_global_new,tsurf_avg)
920	endif
921	enddo
922	enddo
923
924	if (soil_pem) then
925	! II_d.2 Shifting soil temperature to surface
926	call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
927
928	! II_d.3 Update soil temperature
929	write(,)"> Updating soil temperature profile"
930	allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen))
931	tsoil_PEM_timeseries_old = tsoil_PEM_timeseries
932	do islope = 1,nslope
933	tsoil_avg_old = tsoil_PEM(:,:,islope)
934	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
935	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
936
937	do t = 1,timelen
938	do ig = 1,ngrid
939	do isoil = 1,nsoilmx_PEM
940	! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
941	tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
942	! Update of watersoil density
943	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)/(mugazr)
944	if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
945	enddo
946	enddo
947	enddo
948	enddo
949	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
950	deallocate(tsoil_avg_old)
951
952	! II_d.4 Update the ice table
953	allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope))
954	if (icetable_equilibrium) then
955	write(,) "> Updating ice table (equilibrium method)"
956	icetable_thickness_old = icetable_thickness
957	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
958	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
959	else if (icetable_dynamic) then
960	write(,) "> Updating ice table (dynamic method)"
961	ice_porefilling_old = ice_porefilling
962	icetable_depth_old = icetable_depth
963	allocate(porefill(nsoilmx_PEM))
964	do ig = 1,ngrid
965	do islope = 1,nslope
966	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)
967	icetable_depth(ig,islope) = ssi_depth
968	ice_porefilling(ig,:,islope) = porefill
969	enddo
970	enddo
971	deallocate(porefill)
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	endif
974	deallocate(icetable_thickness_old,ice_porefilling_old)
975
976	! II_d.5 Update the soil thermal properties
977	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)
978
979	! II_d.6 Update the mass of the regolith adsorbed
980	totmass_adsco2 = 0.
981	totmass_adsh2o = 0.
982	if (adsorption_pem) then
983	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
984	tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
985	h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
986	do ig = 1,ngrid
987	do islope = 1,nslope
988	do l = 1,nsoilmx_PEM
989	if (l == 1) then
990	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
991	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
992	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
993	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
994	else
995	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
996	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
997	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
998	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
999	endif
1000	enddo
1001	enddo
1002	enddo
1003	write(,) "Total mass of CO2 in the regolith =", totmass_adsco2
1004	write(,) "Total mass of H2O in the regolith =", totmass_adsh2o
1005	endif
1006	endif !soil_pem
1007	deallocate(zshift_surf,zlag)
1008
1009	!------------------------
1010	! II Run
1011	! II_e Outputs
1012	!------------------------
1013	call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
1014	do islope = 1,nslope
1015	write(str2(1:2),'(i2.2)') islope
1016	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
1017	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
1018	call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
1019	call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
1020	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope)))
1021	call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope)))
1022	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
1023	if (icetable_equilibrium) then
1024	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
1025	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope))
1026	else if (icetable_dynamic) then
1027	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
1028	endif
1029
1030	if (soil_pem) then
1031	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
1032	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
1033	if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
1034	if (adsorption_pem) then
1035	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope))
1036	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope))
1037	endif
1038	endif
1039	enddo
1040	deallocate(flag_co2flow,flag_h2oflow)
1041
1042	! Checking mass balance for CO2
1043	totmass_co2ice = 0.
1044	totmass_atmco2 = 0.
1045	do ig = 1,ngrid
1046	totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
1047	do islope = 1,nslope
1048	totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)cell_area(ig)subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
1049	enddo
1050	enddo
1051	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), ' %'
1052	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
1053	write(,) ' /!\ Warning: mass balance is not conseved!'
1054	write(,'(a,f8.3,a)') ' Atmospheric CO2 mass balance = ', 100.(totmass_atmco2 - totmass_atmco2_ini)/totmass_atmco2_ini, ' %'
1055	write(,'(a,f8.3,a)') ' CO2 ice mass balance = ', 100.(totmass_co2ice - totmass_co2ice_ini)/totmass_co2ice_ini, ' %'
1056	write(,'(a,f8.3,a)') ' Adsorbed CO2 mass balance = ', 100.(totmass_adsco2 - totmass_adsco2_ini)/totmass_adsco2_ini, ' %'
1057	endif
1058
1059	!------------------------
1060	! II Run
1061	! II_f Update the tendencies
1062	!------------------------
1063	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)
1064	write(,) "> Updating the H2O sub-surface ice tendency due to lag layer"
1065	if (layering_algo) then
1066	do ig = 1,ngrid
1067	do islope = 1,nslope
1068	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))
1069	enddo
1070	enddo
1071	!~ else
1072	!~ do ig = 1,ngrid
1073	!~ do islope = 1,nslope
1074	!~ 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))
1075	!~ enddo
1076	!~ enddo
1077	endif
1078	if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old)
1079	deallocate(vmr_co2_PEM_phys)
1080	write(,) "> Updating the H2O sub-surface ice depth"
1081
1082	!------------------------
1083	! II Run
1084	! II_g Checking the stopping criterion
1085	!------------------------
1086	write(,) "> Checking the stopping criteria"
1087	call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
1088	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)
1089	i_myear_leg = i_myear_leg + dt
1090	i_myear = i_myear + dt
1091	if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
1092	if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
1093	call system_clock(c2)
1094	if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
1095	if (stopPEM > 0) then
1096	select case (stopPEM)
1097	case(1)
1098	write(,'(a,i0,a)') " *** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
1099	case(2)
1100	write(,'(a,i0,a)') " *** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
1101	case(3)
1102	write(,'(a,i0,a)') " *** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
1103	case(4)
1104	write(,'(a,i0,a)') " *** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
1105	case(5)
1106	write(,'(a,i0)') " *** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
1107	case(6)
1108	write(,'(a,i0)') " *** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
1109	case(7)
1110	write(,'(a,i0)') " *** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
1111	case default
1112	write(,'(a,i0)') " *** STOPPING with unknown stopping criterion code: ", stopPEM
1113	end select
1114	exit
1115	else
1116	write(,'(a,f10.2,a)') ' *** The chained simulation has run for ',i_myear,' Martian years.'
1117	write(,) '**** The PEM can continue!'
1118	write(,) '****'
1119	endif
1120	enddo ! big time iteration loop of the pem
1121	deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed)
1122	deallocate(watersoil_density_PEM_avg,watersurf_density_avg)
1123	deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries)
1124	deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim)
1125	deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
1126	deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
1127	if (layering_algo) deallocate(h2o_ice_depth_old,new_str,new_lag,current)
1128	!------------------------------ END RUN --------------------------------
1129
1130	!------------------------------- OUTPUT --------------------------------
1131	!------------------------
1132	! III Output
1133	! III_a Update surface values for the PCM start files
1134	!------------------------
1135	write(,)
1136	write(,) '******* PEM finalization *******'
1137	! III_a.1 Ice update for start file
1138	write(,) '> Reconstructing perennial ice and frost for the PCM'
1139	watercap = 0.
1140	perennial_co2ice = co2_ice
1141	do ig = 1,ngrid
1142	! H2O ice metamorphism
1143	!if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_h2oice_threshold) then
1144	! h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1145	! qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
1146	!endif
1147
1148	! Is H2O ice still considered as an infinite reservoir for the PCM?
1149	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
1150	! There is enough ice to be considered as an infinite reservoir
1151	watercaptag(ig) = .true.
1152	else
1153	! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
1154	watercaptag(ig) = .false.
1155	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1156	h2o_ice(ig,:) = 0.
1157	endif
1158
1159	! CO2 ice metamorphism
1160	!if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_co2ice_threshold) then
1161	! perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1162	! qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
1163	!endif
1164	enddo
1165
1166	! III.a.3. Tsurf update for start file
1167	write(,) '> Reconstructing the surface temperature for the PCM'
1168	tsurf = tsurf_avg + tsurf_dev
1169	deallocate(tsurf_dev)
1170
1171	! III_a.4 Tsoil update for start file
1172	if (soil_pem) then
1173	write(,) '> Reconstructing the soil temperature profile for the PCM'
1174	inertiesoil = TI_PEM(:,:nsoilmx,:)
1175	! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value
1176	do isoil = 1,nsoilmx
1177	tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:)
1178	enddo
1179	tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev
1180	#ifndef CPP_STD
1181	flux_geo = fluxgeo
1182	#endif
1183	endif
1184	deallocate(tsurf_avg,tsoil_dev)
1185
1186	! III_a.5 Pressure update for start file
1187	write(,) '> Reconstructing the pressure for the PCM'
1188	allocate(ps_start(ngrid))
1189	! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop
1190	ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini
1191	deallocate(ps_avg,ps_dev)
1192
1193	! III_a.6 Tracers update for start file
1194	write(,) '> Reconstructing the tracer VMR for the PCM'
1195	allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
1196	do l = 1,nlayer + 1
1197	zplev_start0(:,l) = ap(l) + bp(l)*ps_start0
1198	zplev_new(:,l) = ap(l) + bp(l)*ps_start
1199	enddo
1200
1201	do nnq = 1,nqtot
1202	if (noms(nnq) /= "co2") then
1203	do l = 1,llm - 1
1204	do ig = 1,ngrid
1205	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))
1206	enddo
1207	q(:,llm,nnq) = q(:,llm - 1,nnq)
1208	enddo
1209	else
1210	do l = 1,llm - 1
1211	do ig = 1,ngrid
1212	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)) &
1213	+ ((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))
1214	enddo
1215	q(:,llm,nnq) = q(:,llm - 1,nnq)
1216	enddo
1217	endif
1218	enddo
1219	deallocate(zplev_start0)
1220
1221	! Conserving the tracers mass for start file
1222	do nnq = 1,nqtot
1223	do ig = 1,ngrid
1224	do l = 1,llm - 1
1225	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
1226	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1227	q(ig,l,nnq) = 1.
1228	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1229	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1230	endif
1231	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1232	enddo
1233	enddo
1234	enddo
1235	deallocate(zplev_new)
1236
1237	! III_a.7 Albedo update for start file
1238	write(,) '> Reconstructing the albedo for the PCM'
1239	do ig = 1,ngrid
1240	if (latitude(ig) < 0.) then
1241	icap = 2 ! Southern hemisphere
1242	else
1243	icap = 1 ! Northern hemisphere
1244	endif
1245	do islope = 1,nslope
1246	! Bare ground
1247	albedo(ig,:,islope) = albedodat(ig)
1248	emis(ig,islope) = emissiv
1249
1250	! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant
1251	! H2O ice
1252	if (h2o_ice(ig,islope) > 0.) then
1253	albedo(ig,:,islope) = albedo_h2o_cap
1254	emis(ig,islope) = 1.
1255	endif
1256	! CO2 ice
1257	if (co2_ice(ig,islope) > 0.) then
1258	albedo(ig,:,islope) = albedo_perennialco2(icap)
1259	emis(ig,islope) = emisice(icap)
1260	endif
1261	! H2O frost
1262	if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
1263	albedo(ig,:,islope) = albedo_h2o_frost
1264	emis(ig,islope) = 1.
1265	endif
1266	! CO2 frost
1267	if (qsurf(ig,igcm_co2,islope) > 0.) then
1268	albedo(ig,:,islope) = albedice(icap)
1269	emis(ig,islope) = emisice(icap)
1270	endif
1271	enddo
1272	enddo
1273
1274	! III_a.8 Orbital parameters update for start file
1275	write(,) '> Setting the new orbital parameters'
1276	if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
1277
1278	!------------------------
1279	! III Output
1280	! III_b Write "restart.nc" and "restartfi.nc"
1281	!------------------------
1282	! III_b.1 Write "restart.nc"
1283	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1284	pday = day_ini
1285	ztime_fin = time_phys
1286	#ifndef CPP_1D
1287	write(,) '> Writing "restart.nc"'
1288	! Correction on teta due to surface pressure changes
1289	allocate(pdyn(ip1jmp1))
1290	call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn)
1291	do i = 1,ip1jmp1
1292	teta(i,:) = teta(i,:)pdyn(i)*rcp
1293	enddo
1294	! Correction on atmospheric pressure
1295	allocate(p(ip1jmp1,nlayer + 1))
1296	call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start,pdyn)
1297	call pression(ip1jmp1,ap,bp,pdyn,p)
1298	! Correction on the mass of atmosphere
1299	call massdair(p,masse)
1300	call dynredem0("restart.nc",day_ini,phis)
1301	call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,pdyn)
1302	deallocate(ap,bp,p,pdyn)
1303	#else
1304	write(,) '> Writing "restart1D.txt"'
1305	call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1306	#endif
1307	deallocate(ps_start0,ps_start)
1308
1309	! III_b.2 Write the "restartfi.nc"
1310	write(,) '> Writing "restartfi.nc"'
1311	#ifndef CPP_STD
1312	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1313	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1314	inertiedat,def_slope,subslope_dist)
1315	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1316	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1317	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1318	wstar,watercap,perennial_co2ice)
1319	#else
1320	if (allocated(noms)) deallocate(noms)
1321	deallocate(qsurf,tsurf,tsoil,emis,watercap,watercaptag,albedo,inertiesoil)
1322	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1323	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
1324	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1325	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1326	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1327	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
1328	tsea_ice,sea_ice)
1329	#endif
1330
1331	!------------------------
1332	! III Output
1333	! III_c Write the "restartpem.nc"
1334	!------------------------
1335	write(,) '> Writing "restartpem.nc"'
1336	if (layering_algo) nb_str_max = get_nb_str_max(layerings_map,ngrid,nslope) ! Get the maximum number of "stratum" in the layerings_mapication (layerings)
1337	call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1338	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
1339	co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
1340
1341	call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
1342
1343	write(,)
1344	write(,) '**** PEM final information *****'
1345	write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
1346	write(,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myearconvert_years, " Earth years."
1347	write(,'(a,f16.4,a)') " + The reached date is now ", (year_bp_ini + i_myear)convert_years, " Earth years."
1348	write(,) "+ PEM: so far, so good!"
1349	write(,) '************************************'
1350
1351	if (layering_algo) then
1352	do islope = 1,nslope
1353	do i = 1,ngrid
1354	call del_layering(layerings_map(i,islope))
1355	enddo
1356	enddo
1357	endif
1358	deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
1359	deallocate(co2_ice,h2o_ice,layerings_map)
1360	!----------------------------- END OUTPUT ------------------------------
1361
1362	END PROGRAM pem

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format