Context Navigation

pem.F90 @ 3742

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

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

Download in other formats:

Original Format