Context Navigation

pem.F90 @ 3707

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

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

Download in other formats:

Original Format