Context Navigation

pem.F90 @ 3616

Last change on this file since 3616 was 3616, checked in by jbclement, 10 days ago

PEM:

Improvement of the loops to update surface temperature in case of CO2 ice is disappearing.
Improvement of messages printed on the terminal.

JBC

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

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

Download in other formats:

Original Format