Context Navigation

pem.F90 @ 3620

Last change on this file since 3620 was 3620, checked in by jbclement, 28 hours ago
PEM: Few small adjustments in the code. JBC
File size: 67.8 KB

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

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

Download in other formats:

Original Format