Context Navigation

pem.F90 @ 3836

Last change on this file since 3836 was 3836, checked in by jbclement, 8 days ago

COMMON:
Rework related to the command-line parsing:

Replace ad-hoc argument parsing with a unified 'parse_args' subroutine, allowing easier extension to other programs and options across models;;
Use of '--version' (with ab optional output file) to print compilation/version details;
Addition of 'job_timelimit_mod' module to handle SLURM/PBS job time-limit via '--jobid' (currently only used in the PEM), allowing easier extension to other programs;
Replace manual SSO handling with 'parse_args' for the Mars start2archive;
Clean-up related legacy code in the programs supporting the version option.

JBC

File size: 68.6 KB

Line
1	!------------------------
2	! I Initialization
3	! I_a Read the "run.def"
4	! I_b Read the "start.nc" and "startfi.nc"
5	! I_c Subslope parametrisation
6	! I_d Read the PCM data and convert them to the physical grid
7	! I_e Initialization of the PEM variable and soil
8	! I_f Compute tendencies
9	! I_g Compute global surface pressure
10	! I_h Read the "startpem.nc"
11	! I_i Compute orbit criterion
12
13	! II Run
14	! II_a Update pressure, ice and tracers
15	! II_b Evolution of ice
16	! II_c Flow of glaciers
17	! II_d Update surface and soil temperatures
18	! II_e Outputs
19	! II_f Update the tendencies
20	! II_g Checking the stopping criterion
21
22	! III Output
23	! III_a Update surface values for the PCM start files
24	! III_b Write the "restart.nc" and "restartfi.nc"
25	! III_c Write the "restartpem.nc"
26	!------------------------
27
28	PROGRAM pem
29
30	use phyetat0_mod, only: phyetat0
31	use phyredem, only: physdem0, physdem1
32	use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
33	use turb_mod, only: q2, wstar
34	use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h
35	use logic_mod, only: iflag_phys
36	use mod_const_mpi, only: COMM_LMDZ
37	use infotrac
38	use geometry_mod, only: latitude_deg
39	use conf_pem_mod, only: conf_pem
40	use pemredem, only: pemdem0, pemdem1
41	use glaciers_mod, only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
42	metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
43	use stopping_crit_mod, only: stopping_crit_h2o_ice, stopping_crit_co2
44	use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
45	use evol_ice_mod, only: evol_co2_ice, evol_h2o_ice
46	use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
47	TI_PEM, & ! Soil thermal inertia
48	tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
49	fluxgeo ! Geothermal flux for the PEM and PCM
50	use adsorption_mod, only: regolith_adsorption, adsorption_pem, & ! Bool to check if adsorption, main subroutine
51	ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
52	co2_adsorbed_phys, h2o_adsorbed_phys ! Mass of co2 and h2O adsorbed
53	use time_evol_mod, only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
54	use orbit_param_criterion_mod, only: orbit_param_criterion
55	use recomp_orb_param_mod, only: recomp_orb_param
56	use ice_table_mod, only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, &
57	ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi
58	use soil_thermalproperties_mod, only: update_soil_thermalproperties
59	use time_phylmdz_mod, only: daysec, dtphys
60	use abort_pem_mod, only: abort_pem
61	use soil_settings_PEM_mod, only: soil_settings_PEM
62	use compute_tend_mod, only: compute_tend
63	use info_PEM_mod, only: info_PEM
64	use get_timelen_PCM_mod, only: get_timelen_PCM
65	use pemetat0_mod, only: pemetat0
66	use read_data_PCM_mod, only: read_data_PCM
67	use recomp_tend_mod, only: recomp_tend_co2, recomp_tend_h2o
68	use compute_soiltemp_mod, only: compute_tsoil_pem, shift_tsoil2surf
69	use writediagpem_mod, only: writediagpem, writediagsoilpem
70	use co2condens_mod, only: CO2cond_ps
71	use layering_mod, only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, &
72	ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str, &
73	print_layering
74	use dyn_ss_ice_m_mod, only: dyn_ss_ice_m
75	use parse_args_mod, only: parse_args
76	use job_timelimit_mod, only: timelimit, antetime, timewall
77	use paleoclimate_mod, only: h2o_ice_depth, zdqsdif_ssi_tot
78
79	#ifndef CPP_STD
80	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil
81	use surfdat_h, only: tsurf, qsurf, emis, emissiv, emisice, ini_surfdat_h, &
82	albedodat, albedice, albedo_h2o_frost, albedo_h2o_cap, &
83	zmea, zstd, zsig, zgam, zthe, frost_albedo_threshold, &
84	watercap, watercaptag, perennial_co2ice, albedo_perennialco2
85	use dimradmars_mod, only: totcloudfrac, albedo
86	use dust_param_mod, only: tauscaling
87	use tracer_mod, only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
88	use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
89	use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
90	use surfini_mod, only: surfini
91	use comcstfi_h, only: mugaz
92	#else
93	use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
94	use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, &
95	PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
96	ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
97	use aerosol_mod, only: iniaerosol
98	use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit
99	use comcstfi_mod, only: pi, rad, g, r, cpp, rcp, mugaz
100	#endif
101
102	#ifndef CPP_1D
103	use comconst_mod, only: pi, rad, g, r, cpp, rcp => kappa
104	use iniphysiq_mod, only: iniphysiq
105	use control_mod, only: iphysiq, day_step, nsplit_phys
106	#else
107	use comcstfi_h, only: pi, rad, g, r, cpp, rcp
108	use time_phylmdz_mod, only: iphysiq, steps_per_sol
109	use regular_lonlat_mod, only: init_regular_lonlat
110	use physics_distribution_mod, only: init_physics_distribution
111	use mod_grid_phy_lmdz, only: regular_lonlat
112	use init_testphys1d_mod, only: init_testphys1d
113	use comvert_mod, only: ap, bp
114	use writerestart1D_mod, only: writerestart1D
115	#endif
116
117	implicit none
118
119	include "dimensions.h"
120	include "paramet.h"
121	include "comgeom.h"
122	include "iniprint.h"
123
124	! Same variable names as in the PCM
125	integer, parameter :: ngridmx = 2 + (jjm - 1)*iim - 1/jjm
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) Average 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	real(kind = 16) :: totmass_co2ice, totmass_atmco2 ! Current total CO2 masses
197	real(kind = 16) :: totmass_co2ice_ini, totmass_atmco2_ini ! Initial total CO2 masses
198
199	! Variables for the evolution of layered layerings_map
200	type(layering), dimension(:,:), allocatable :: layerings_map ! Layering for each grid point and slope
201	type(ptrarray), dimension(:,:), allocatable :: current ! Current active stratum in the layering
202	logical, dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm
203	real, dimension(:,:), allocatable :: h2o_ice_depth_old ! Old depth of subsurface ice layer
204
205	! Variables for slopes
206	integer(kind = 1), dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow
207	integer(kind = 1), dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow
208
209	! Variables for surface and soil
210	real, dimension(:,:), allocatable :: tsurf_avg ! Grid points x Slope field: Average surface temperature [K]
211	real, dimension(:,:), allocatable :: tsurf_dev ! Grid points x Slope field: Surface temperature deviation [K]
212	real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Grid points x Slope field: Average surface temperature of first call of the PCM [K]
213	real, dimension(:,:,:), allocatable :: tsoil_avg ! Grid points x Soil x Slope field: Average Soil Temperature [K]
214	real, dimension(:,:), allocatable :: tsoil_avg_old ! Grid points x Soil field: Average Soil Temperature at the previous time step [K]
215	real, dimension(:,:,:), allocatable :: tsoil_dev ! Grid points x Soil x Slope field: Soil temperature deviation [K]
216	real, dimension(:,:,:,:), allocatable :: tsoil_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
217	real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K]
218	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]
219	real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
220	real, dimension(:,:), allocatable :: watersurf_density_avg ! Grid points x Slope: Average water surface density [kg/m^3]
221	real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
222	real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Grid points x Soil x Slopes: Average water soil density [kg/m^3]
223	real, dimension(:), allocatable :: delta_co2_adsorbed ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
224	real, dimension(:), allocatable :: delta_h2o_adsorbed ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
225	real(kind = 16) :: totmass_adsco2, totmass_adsco2_ini ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
226	real :: totmass_adsh2o ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
227	logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
228	real, dimension(:,:), allocatable :: icetable_thickness_old ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
229	real, dimension(:,:,:), allocatable :: ice_porefilling_old ! ngrid x nslope: Ice pore filling at the previous iteration [m]
230	real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
231	real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table
232	real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table
233	real, dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m]
234	real, dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m]
235	real, dimension(:,:), allocatable :: icetable_depth_old ! Old depth of the ice table
236
237	! Some variables for the PEM run
238	real, parameter :: year_step = 1 ! Timestep for the pem
239	real :: i_myear_leg ! Number of iteration
240	real :: n_myear_leg ! Maximum number of iterations before stopping
241	real :: i_myear ! Global number of Martian years of the chained simulations
242	real :: n_myear ! Maximum number of Martian years of the chained simulations
243	real :: timestep ! Timestep [s]
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(8) :: date
247	character(10) :: time
248	character(5) :: zone
249	integer, dimension(8) :: values
250	character(128) :: dir = ' '
251	character(32) :: logname = ' '
252	character(32) :: hostname = ' '
253
254	#ifdef CPP_STD
255	real :: frost_albedo_threshold = 0.05 ! Frost albedo threeshold to convert fresh frost to old ice
256	real :: albedo_h2o_frost ! Albedo of h2o frost
257	real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
258	real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
259	real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
260	real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
261	real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
262	real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case
263	real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case
264	real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case
265	real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case
266	real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
267	logical, dimension(:), allocatable :: watercaptag ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
268	real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case
269	real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case
270	#endif
271
272	#ifdef CPP_1D
273	integer :: nsplit_phys
274	integer, parameter :: jjm_value = jjm - 1
275	integer :: day_step
276
277	! Dummy variables to use the subroutine 'init_testphys1d'
278	logical :: therestart1D, therestartfi, prescribed_h2ovap
279	integer :: ndt, day0
280	real :: ptif, pks, day, gru, grv, h2ovap_relax_time
281	real, dimension(:), allocatable :: zqsat
282	real, dimension(:,:,:), allocatable :: dq, dqdyn
283	real, dimension(nlayer) :: play, w, q_prescribed_h2o_vap
284	real, dimension(nlayer + 1) :: plev
285	#else
286	integer, parameter :: jjm_value = jjm
287	real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
288	real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
289	real, dimension(:,:), allocatable :: p ! Grid points x Atmosphere: pressure to recompute and write in restart (ip1jmp1,llmp1)
290	#endif
291
292	! Loop variables
293	integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
294	logical :: num_str
295
296	! CODE
297	write(,) ' * . . + . * . + . . . '
298	write(,) ' + _______ ________ ____ ____ * + '
299	write(,) ' + . * \|_ __ \\|_ __ \|\|_ \ / _\| . *'
300	write(,) ' . . \| \|__) \| \| \|_ \_\| \| \/ \| * * . '
301	write(,) ' . \| ___/ \| _\| _ \| \|\ /\| \| . . '
302	write(,) '. * * _\| \|_ _\| \|__/ \| _\| \|_\/_\| \|_ * '
303	write(,) ' + \|_____\| \|________\|\|_____\|\|_____\| + . '
304	write(,) ' . * . * . + * . + .'
305
306	! Elapsed time with system clock
307	call system_clock(count_rate = cr)
308	call system_clock(c1)
309	call parse_args()
310
311	! Some user info
312	call date_and_time(date,time,zone,values)
313	call getcwd(dir) ! Current directory
314	call getlog(logname) ! User name
315	call hostnm(hostname) ! Machine/station name
316	write(,)
317	write(,) '******* PEM information *******'
318	write(,) '+ User : '//trim(logname)
319	write(,) '+ Machine : '//trim(hostname)
320	write(,) '+ Directory: '//trim(dir)
321	write(*,'(a,i2,a,i2,a,i4)') ' + Date : ',values(3),'/',values(2),'/',values(1)
322	write(*,'(a,i2,a,i2,a,i2,a)') ' + Time : ',values(5),':',values(6),':',values(7)
323
324	! Parallel variables
325	#ifndef CPP_STD
326	is_sequential = .true.
327	is_parallel = .false.
328	is_mpi_root = .true.
329	is_omp_root = .true.
330	is_master = .true.
331	#endif
332
333	! Some constants
334	day_ini = 0
335	time_phys = 0.
336	ngrid = ngridmx
337	A = (1./m_co2 - 1./m_noco2)
338	B = 1./m_noco2
339	year_day = 669
340	daysec = 88775.
341	timestep = year_daydaysecyear_step
342
343	!----------------------------- INITIALIZATION --------------------------
344	!------------------------
345	! I Initialization
346	! I_a Read the "run.def"
347	!------------------------
348	write(,)
349	write(,) '******* PEM initialization *******'
350	write(,) '> Reading "run.def" (PEM)'
351	#ifndef CPP_1D
352	dtphys = daysec/48. ! Dummy value (overwritten in phyetat0)
353	call conf_gcm(99,.true.)
354	call infotrac_init
355	nq = nqtot
356	allocate(q(ip1jmp1,llm,nqtot))
357	allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
358	#else
359	allocate(q(1,llm,nqtot),pdyn(1))
360	allocate(longitude(1),latitude(1),cell_area(1))
361
362	therestart1D = .false. ! Default value
363	inquire(file = 'start1D.txt',exist = therestart1D)
364	if (.not. therestart1D) then
365	write(,) 'There is no "start1D.txt" file!'
366	error stop 'Initialization cannot be done for the 1D PEM.'
367	endif
368	therestartfi = .false. ! Default value
369	inquire(file = 'startfi.nc',exist = therestartfi)
370	if (.not. therestartfi) then
371	write(,) 'There is no "startfi.nc" file!'
372	error stop 'Initialization cannot be done for the 1D PEM.'
373	endif
374
375	write(,) '> Reading "start1D.txt" and "startfi.nc"'
376	call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, &
377	time_0,pdyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
378	play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap)
379	nsplit_phys = 1
380	day_step = steps_per_sol
381	#endif
382
383	call conf_pem(i_myear,n_myear)
384
385	!------------------------
386	! I Initialization
387	! I_b Read of the "start.nc" and "starfi.nc"
388	!------------------------
389	! I_b.1 Read "start.nc"
390	write(,) '> Reading "start.nc"'
391	allocate(ps_start0(ngrid))
392	#ifndef CPP_1D
393	allocate(pdyn(ip1jmp1))
394	call dynetat0(start_name,vcov,ucov,teta,q,masse,pdyn,phis,time_0)
395
396	call gr_dyn_fi(1,iip1,jjp1,ngridmx,pdyn,ps_start0)
397
398	call iniconst ! Initialization of dynamical constants (comconst_mod)
399	call inigeom ! Initialization of geometry
400
401	allocate(ap(nlayer + 1))
402	allocate(bp(nlayer + 1))
403	status = nf90_open(start_name,NF90_NOWRITE,ncid)
404	status = nf90_inq_varid(ncid,"ap",apvarid)
405	status = nf90_get_var(ncid,apvarid,ap)
406	status = nf90_inq_varid(ncid,"bp",bpvarid)
407	status = nf90_get_var(ncid,bpvarid,bp)
408	status = nf90_close(ncid)
409
410	! Initialization of physics constants and variables (comcstfi_h)
411	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)
412	#else
413	ps_start0(1) = pdyn(1)
414	#endif
415	deallocate(pdyn)
416
417	! In the PCM, these values are given to the physic by the dynamic.
418	! Here we simply read them in the "startfi.nc" file
419	status = nf90_open(startfi_name,NF90_NOWRITE,ncid)
420	status = nf90_inq_varid(ncid,"longitude",lonvarid)
421	status = nf90_get_var(ncid,lonvarid,longitude)
422	status = nf90_inq_varid(ncid,"latitude",latvarid)
423	status = nf90_get_var(ncid,latvarid,latitude)
424	status = nf90_inq_varid(ncid,"area",areavarid)
425	status = nf90_get_var(ncid,areavarid,cell_area)
426	status = nf90_inq_varid(ncid,"soildepth",sdvarid)
427	status = nf90_get_var(ncid,sdvarid,mlayer)
428	status = nf90_close(ncid)
429
430	! I_b.2 Read the "startfi.nc"
431	! First we read the initial state (starfi.nc)
432	#ifndef CPP_STD
433	write(,) '> Reading "startfi.nc"'
434	call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
435	tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, &
436	watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
437
438	! Remove unphysical values of surface tracer
439	where (qsurf < 0.) qsurf = 0.
440
441	call surfini(ngrid,nslope,qsurf)
442	#else
443	call phys_state_var_init(nq)
444	if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
445	call initracer(ngrid,nq)
446	call iniaerosol()
447	allocate(tsurf_read_generic(ngrid))
448	allocate(qsurf_read_generic(ngrid,nq))
449	allocate(tsoil_read_generic(ngrid,nsoilmx))
450	allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
451	allocate(emis_read_generic(ngrid))
452	allocate(albedo_read_generic(ngrid,2))
453	allocate(qsurf(ngrid,nq,1))
454	allocate(tsurf(ngrid,1))
455	allocate(tsoil(ngrid,nsoilmx,1))
456	allocate(emis(ngrid,1))
457	allocate(watercap(ngrid,1))
458	allocate(watercaptag(ngrid))
459	allocate(albedo(ngrid,2,1))
460	allocate(inertiesoil(ngrid,nsoilmx,1))
461	call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
462	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
463	qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, &
464	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
465	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
466
467	nslope = 1
468	call ini_comslope_h(ngrid,1)
469
470	qsurf(:,:,1) = qsurf_read_generic
471	tsurf(:,1) = tsurf_read_generic
472	tsoil(:,:,1) = tsoil_read_generic
473	emis(:,1) = emis_read_generic
474	watercap(:,1) = 0.
475	watercaptag(:) = .false.
476	albedo(:,1,1) = albedo_read_generic(:,1)
477	albedo(:,2,1) = albedo_read_generic(:,2)
478	inertiesoil(:,:,1) = inertiedat
479
480	if (nslope == 1) then
481	def_slope(1) = 0
482	def_slope(2) = 0
483	def_slope_mean = 0
484	subslope_dist(:,1) = 1.
485	endif
486
487	! Remove unphysical values of surface tracer
488	qsurf(:,:,1) = qsurf_read_generic
489	where (qsurf < 0.) qsurf = 0.
490
491	deallocate(tsurf_read_generic,qsurf_read_generic,qsoil_read_generic,emis_read_generic)
492	#endif
493
494	do nnq = 1,nqtot ! Why not using ini_tracer?
495	if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
496	if (noms(nnq) == "h2o_vap") then
497	igcm_h2o_vap = nnq
498	mmol(igcm_h2o_vap) = 18.
499	endif
500	if (noms(nnq) == "co2") igcm_co2 = nnq
501	enddo
502
503	!------------------------
504	! I Initialization
505	! I_c Subslope parametrisation
506	!------------------------
507	! Define some slope statistics
508	iflat = 1
509	do islope = 2,nslope
510	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
511	enddo
512	write(,) 'Flat slope for islope = ',iflat
513	write(,) 'Corresponding criterium = ',def_slope_mean(iflat)
514
515	!------------------------
516	! I Initialization
517	! I_d Read the PCM data and convert them to the physical grid
518	!------------------------
519	! 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
520	call get_timelen_PCM("data_PCM_Y1.nc",timelen)
521
522	allocate(vmr_co2_PCM(ngrid,timelen),q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
523	allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid))
524	allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
525	allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope))
526	allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_timeseries(ngrid,nsoilmx,nslope,timelen))
527	allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
528
529	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, &
530	tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2_PEM_phys,q_h2o_PEM_phys,watersurf_density_avg,watersoil_density_timeseries)
531
532	! Compute the deviation from the average
533	allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
534	ps_dev = ps_start0 - ps_avg
535	tsurf_dev = tsurf - tsurf_avg
536	tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:)
537
538	!------------------------
539	! I Initialization
540	! I_e Initialization of the PEM variables and soil
541	!------------------------
542	call end_comsoil_h_PEM
543	call ini_comsoil_h_PEM(ngrid,nslope)
544	call end_adsorption_h_PEM
545	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
546	call end_ice_table
547	call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
548
549	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))
550	if (soil_pem) then
551	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
552	tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
553	watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
554	tsoil_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_timeseries
555	do l = nsoilmx + 1,nsoilmx_PEM
556	tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
557	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
558	tsoil_PEM_timeseries(:,l,:,:) = tsoil_timeseries(:,nsoilmx,:,:)
559	enddo
560	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
561	endif !soil_pem
562	deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries)
563
564	!------------------------
565	! I Initialization
566	! I_f Compute tendencies
567	!------------------------
568	allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope))
569	call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
570	call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
571	d_co2ice_ini = d_co2ice
572	deallocate(min_co2_ice,min_h2o_ice)
573
574	!------------------------
575	! I Initialization
576	! I_g Compute global surface pressure
577	!------------------------
578	total_surface = sum(cell_area)
579	ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
580	ps_avg_global_new = ps_avg_global_ini
581	write(,) "Total surface of the planet =", total_surface
582	write(,) "Initial global average pressure =", ps_avg_global_ini
583
584	!------------------------
585	! I Initialization
586	! I_h Read the "startpem.nc"
587	!------------------------
588	write(,) '> Reading "startpem.nc"'
589	allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),layerings_map(ngrid,nslope))
590	allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid))
591	delta_h2o_icetablesublim = 0.
592	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
593	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, &
594	watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,layerings_map)
595	deallocate(tsurf_avg_yr1)
596
597	! We save the places where h2o ice is sublimating
598	! We compute the surface of h2o ice sublimating
599	allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope),co2ice_disappeared(ngrid,nslope))
600	co2ice_sublim_surf_ini = 0.
601	h2oice_ini_surf = 0.
602	is_co2ice_sublim_ini = .false.
603	is_h2oice_sublim_ini = .false.
604	is_co2ice_ini = .false.
605	co2ice_disappeared = .false.
606	totmass_co2ice_ini = 0.
607	totmass_atmco2_ini = 0.
608	if (layering_algo) then
609	do ig = 1,ngrid
610	do islope = 1,nslope
611	if (is_co2ice_str(layerings_map(ig,islope)%top)) then
612	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
613	else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
614	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
615	else
616	call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
617	endif
618	enddo
619	enddo
620	! We put the sublimating tendency coming from subsurface ice into the overall tendency
621	where (zdqsdif_ssi_tot < 0.)
622	d_h2oice = zdqsdif_ssi_tot
623	end where
624	endif
625	do i = 1,ngrid
626	totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
627	do islope = 1,nslope
628	totmass_co2ice_ini = totmass_co2ice_ini + co2_ice(i,islope)cell_area(i)subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)
629	if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
630	if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
631	is_co2ice_sublim_ini(i,islope) = .true.
632	co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
633	endif
634	if (d_h2oice(i,islope) < 0.) then
635	if (h2o_ice(i,islope) > 0.) then
636	is_h2oice_sublim_ini(i,islope) = .true.
637	h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
638	else if (h2o_ice_depth(i,islope) > 0.) then
639	is_h2oice_sublim_ini(i,islope) = .true.
640	endif
641	endif
642	enddo
643	enddo
644	write(,) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini
645	write(,) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf
646
647	totmass_adsco2_ini = 0.
648	totmass_adsh2o = 0.
649	if (adsorption_pem) then
650	do ig = 1,ngrid
651	do islope = 1,nslope
652	do l = 1,nsoilmx_PEM - 1
653	if (l == 1) then
654	totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
655	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
656	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
657	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
658	else
659	totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
660	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
661	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
662	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
663	endif
664	enddo
665	enddo
666	enddo
667	totmass_adsco2 = totmass_adsco2_ini
668	write(,) "Tot mass of CO2 in the regolith =", totmass_adsco2
669	write(,) "Tot mass of H2O in the regolith =", totmass_adsh2o
670	endif ! adsorption
671
672	!------------------------
673	! I Initialization
674	! I_i Compute orbit criterion
675	!------------------------
676	#ifndef CPP_STD
677	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
678	#else
679	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
680	#endif
681
682	n_myear_leg = Max_iter_pem
683	if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
684
685	!-------------------------- END INITIALIZATION -------------------------
686
687	!-------------------------------- RUN ----------------------------------
688	!------------------------
689	! II Run
690	! II_a Update pressure, ice and tracers
691	!------------------------
692	write(,)
693	write(,) '******* PEM cycle *******'
694	i_myear_leg = 0
695	stopPEM = 0
696	if (layering_algo) then
697	allocate(h2o_ice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
698	new_str = .true.
699	new_lag = .true.
700	do islope = 1,nslope
701	do ig = 1,ngrid
702	current(ig,islope)%p => layerings_map(ig,islope)%top
703	enddo
704	enddo
705	endif
706
707	do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
708	! II.a.1. Compute updated global pressure
709	write(,'(a,f10.2)') ' *** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
710	write(,) "> Updating the surface pressure"
711	ps_avg_global_old = ps_avg_global_new
712	do i = 1,ngrid
713	do islope = 1,nslope
714	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
715	enddo
716	enddo
717	if (adsorption_pem) then
718	do i = 1,ngrid
719	ps_avg_global_new = ps_avg_global_new - gcell_area(i)delta_co2_adsorbed(i)/total_surface
720	enddo
721	endif
722	ps_avg = ps_avg*ps_avg_global_new/ps_avg_global_old
723	write(,) 'Global average pressure old time step:',ps_avg_global_old
724	write(,) 'Global average pressure new time step:',ps_avg_global_new
725
726	! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
727	write(,) "> Updating the surface pressure timeseries for the new pressure"
728	allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
729	do l = 1,nlayer + 1
730	do ig = 1,ngrid
731	zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
732	enddo
733	enddo
734	ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old
735	write(,) "> Updating the pressure levels timeseries for the new pressure"
736	allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
737	do l = 1,nlayer + 1
738	do ig = 1,ngrid
739	zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
740	enddo
741	enddo
742
743	! II.a.3. Tracers timeseries
744	write(,) "> Updating the tracer VMR timeseries for the new pressure"
745	allocate(vmr_co2_PEM_phys(ngrid,timelen))
746	l = 1
747	do ig = 1,ngrid
748	do t = 1,timelen
749	! H2O
750	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))/ &
751	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
752	if (q_h2o_PEM_phys(ig,t) < 0) then
753	q_h2o_PEM_phys(ig,t) = 1.e-30
754	else if (q_h2o_PEM_phys(ig,t) > 1) then
755	q_h2o_PEM_phys(ig,t) = 1.
756	endif
757	! CO2
758	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))/ &
759	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
760	+ ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
761	- (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/ &
762	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
763	if (q_co2_PEM_phys(ig,t) < 0) then
764	q_co2_PEM_phys(ig,t) = 1.e-30
765	else if (q_co2_PEM_phys(ig,t) > 1) then
766	q_co2_PEM_phys(ig,t) = 1.
767	endif
768	mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
769	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
770	enddo
771	enddo
772	deallocate(zplev_timeseries_new,zplev_timeseries_old)
773
774	!------------------------
775	! II Run
776	! II_b Evolution of ice
777	!------------------------
778	allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
779	if (layering_algo) then
780	h2o_ice_depth_old = h2o_ice_depth
781	do islope = 1,nslope
782	do ig = 1,ngrid
783	call make_layering(layerings_map(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p)
784	co2_ice(ig,islope) = 0.
785	h2o_ice(ig,islope) = 0.
786	h2o_ice_depth(ig,islope) = 0.
787	if (is_co2ice_str(layerings_map(ig,islope)%top)) then
788	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
789	else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
790	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
791	else
792	call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
793	endif
794	enddo
795	enddo
796	else
797	zlag = 0.
798	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
799	call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
800	endif
801
802	!------------------------
803	! II Run
804	! II_c Flow of glaciers
805	!------------------------
806	allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
807	if (co2ice_flow .and. nslope > 1) then
808	call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
809	ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
810	if (layering_algo) then
811	do islope = 1,nslope
812	do ig = 1,ngrid
813	layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
814	enddo
815	enddo
816	endif
817	endif
818	if (h2oice_flow .and. nslope > 1) then
819	call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
820	if (layering_algo) then
821	do islope = 1,nslope
822	do ig = 1,ngrid
823	!~ layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
824	enddo
825	enddo
826	endif
827	endif
828
829	!------------------------
830	! II Run
831	! II_d Update surface and soil temperatures
832	!------------------------
833	! II_d.1 Update Tsurf
834	write(,) "> Updating surface temperature"
835	do ig = 1,ngrid
836	do islope = 1,nslope
837	! CO2 ice disappeared so we look for the closest point without CO2 ice
838	if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then
839	co2ice_disappeared(ig,islope) = .true.
840	if (latitude_deg(ig) > 0.) then ! North hemisphere
841	outer1: do ig_loop = ig,ngrid ! Go towards equator
842	do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be equator-ward)
843	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
844	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
845	exit outer1
846	endif
847	enddo
848	do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be pole-ward)
849	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
850	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
851	exit outer1
852	endif
853	enddo
854	enddo outer1
855	else ! South hemisphere
856	outer2: do ig_loop = ig,1,-1 ! Go towards equator
857	do islope_loop = islope + 1,nslope ! Go over the slopes (forward 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 outer2
861	endif
862	enddo
863	do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward 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 outer2
867	endif
868	enddo
869	enddo outer2
870	endif
871	else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as Tcond CO2
872	call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_ini,ps_avg_global_new,tsurf_avg)
873	endif
874	enddo
875	enddo
876
877	if (soil_pem) then
878	! II_d.2 Shifting soil temperature to surface
879	call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
880
881	! II_d.3 Update soil temperature
882	write(,)"> Updating soil temperature profile"
883	allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen))
884	tsoil_PEM_timeseries_old = tsoil_PEM_timeseries
885	do islope = 1,nslope
886	tsoil_avg_old = tsoil_PEM(:,:,islope)
887	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
888	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
889
890	do t = 1,timelen
891	do ig = 1,ngrid
892	do isoil = 1,nsoilmx_PEM
893	! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
894	tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
895	! Update of watersoil density
896	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)
897	if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
898	enddo
899	enddo
900	enddo
901	enddo
902	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
903	deallocate(tsoil_avg_old)
904
905	! II_d.4 Update the ice table
906	allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope))
907	if (icetable_equilibrium) then
908	write(,) "> Updating ice table (equilibrium method)"
909	icetable_thickness_old = icetable_thickness
910	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
911	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
912	else if (icetable_dynamic) then
913	write(,) "> Updating ice table (dynamic method)"
914	ice_porefilling_old = ice_porefilling
915	icetable_depth_old = icetable_depth
916	allocate(porefill(nsoilmx_PEM))
917	do ig = 1,ngrid
918	do islope = 1,nslope
919	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)
920	icetable_depth(ig,islope) = ssi_depth
921	ice_porefilling(ig,:,islope) = porefill
922	enddo
923	enddo
924	deallocate(porefill)
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	endif
927	deallocate(icetable_thickness_old,ice_porefilling_old)
928
929	! II_d.5 Update the soil thermal properties
930	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)
931
932	! II_d.6 Update the mass of the regolith adsorbed
933	totmass_adsco2 = 0.
934	totmass_adsh2o = 0.
935	if (adsorption_pem) then
936	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
937	tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
938	h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
939	do ig = 1,ngrid
940	do islope = 1,nslope
941	do l = 1,nsoilmx_PEM
942	if (l == 1) then
943	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
944	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
945	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
946	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
947	else
948	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
949	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
950	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
951	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
952	endif
953	enddo
954	enddo
955	enddo
956	write(,) "Total mass of CO2 in the regolith =", totmass_adsco2
957	write(,) "Total mass of H2O in the regolith =", totmass_adsh2o
958	endif
959	endif !soil_pem
960	deallocate(zshift_surf,zlag)
961
962	!------------------------
963	! II Run
964	! II_e Outputs
965	!------------------------
966	call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
967	do islope = 1,nslope
968	write(str2(1:2),'(i2.2)') islope
969	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
970	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
971	call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
972	call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
973	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope)))
974	call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope)))
975	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
976	if (icetable_equilibrium) then
977	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
978	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope))
979	else if (icetable_dynamic) then
980	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
981	endif
982
983	if (soil_pem) then
984	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
985	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
986	if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
987	if (adsorption_pem) then
988	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope))
989	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope))
990	endif
991	endif
992	enddo
993	deallocate(flag_co2flow,flag_h2oflow)
994
995	! Checking mass balance for CO2
996	totmass_co2ice = 0.
997	totmass_atmco2 = 0.
998	do ig = 1,ngrid
999	totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
1000	do islope = 1,nslope
1001	totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)cell_area(ig)subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
1002	enddo
1003	enddo
1004	write(,'(a,f8.3,a)') " > Relative total CO2 mass balance = ", 100.(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini), ' %'
1005	if ((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini) > 0.01) then
1006	write(,) ' /!\ Warning: mass balance is not conseved!'
1007	write(,'(a,f8.3,a)') ' Atmospheric CO2 mass balance = ', 100.(totmass_atmco2 - totmass_atmco2_ini)/totmass_atmco2_ini, ' %'
1008	write(,'(a,f8.3,a)') ' CO2 ice mass balance = ', 100.(totmass_co2ice - totmass_co2ice_ini)/totmass_co2ice_ini, ' %'
1009	write(,'(a,f8.3,a)') ' Adsorbed CO2 mass balance = ', 100.(totmass_adsco2 - totmass_adsco2_ini)/totmass_adsco2_ini, ' %'
1010	endif
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	if (layering_algo) then
1019	do ig = 1,ngrid
1020	do islope = 1,nslope
1021	if (is_h2oice_sublim_ini(ig,islope) .and. h2o_ice_depth(ig,islope) > 0.) call recomp_tend_h2o(h2o_ice_depth_old(ig,islope),h2o_ice_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope))
1022	enddo
1023	enddo
1024	!~ else
1025	!~ do ig = 1,ngrid
1026	!~ do islope = 1,nslope
1027	!~ 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))
1028	!~ enddo
1029	!~ enddo
1030	endif
1031	if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old)
1032	deallocate(vmr_co2_PEM_phys)
1033
1034	!------------------------
1035	! II Run
1036	! II_g Checking the stopping criterion
1037	!------------------------
1038	write(,) "> Checking the stopping criteria"
1039	call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
1040	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)
1041	i_myear_leg = i_myear_leg + dt
1042	i_myear = i_myear + dt
1043	if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
1044	if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
1045	call system_clock(c2)
1046	if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
1047	if (stopPEM > 0) then
1048	select case (stopPEM)
1049	case(1)
1050	write(,'(a,i0,a)') " *** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
1051	case(2)
1052	write(,'(a,i0,a)') " *** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
1053	case(3)
1054	write(,'(a,i0,a)') " *** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
1055	case(4)
1056	write(,'(a,i0,a)') " *** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
1057	case(5)
1058	write(,'(a,i0)') " *** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
1059	case(6)
1060	write(,'(a,i0)') " *** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
1061	case(7)
1062	write(,'(a,i0)') " *** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
1063	case default
1064	write(,'(a,i0)') " *** STOPPING with unknown stopping criterion code: ", stopPEM
1065	end select
1066	exit
1067	else
1068	write(,'(a,f10.2,a)') ' *** The chained simulation has run for ',i_myear,' Martian years.'
1069	write(,) '**** The PEM can continue!'
1070	write(,) '****'
1071	endif
1072	enddo ! big time iteration loop of the pem
1073	deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed)
1074	deallocate(watersoil_density_PEM_avg,watersurf_density_avg)
1075	deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries)
1076	deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim)
1077	deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
1078	deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
1079	if (layering_algo) deallocate(h2o_ice_depth_old,new_str,new_lag,current)
1080	!------------------------------ END RUN --------------------------------
1081
1082	!------------------------------- OUTPUT --------------------------------
1083	!------------------------
1084	! III Output
1085	! III_a Update surface values for the PCM start files
1086	!------------------------
1087	write(,)
1088	write(,) '******* PEM finalization *******'
1089	! III_a.1 Ice update for start file
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	! 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.3. 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.4 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.5 Pressure update for start file
1139	write(,) '> Reconstructing the pressure for the PCM'
1140	allocate(ps_start(ngrid))
1141	! The pressure deviation is rescaled as well 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.6 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.7 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.8 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)*rcp
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,pdyn)
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(layerings_map,ngrid,nslope) ! Get the maximum number of "stratum" in the layerings_mapication (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,layerings_map)
1292
1293	call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
1294
1295	write(,)
1296	write(,) '**** PEM final information *****'
1297	write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
1298	write(,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myearconvert_years, " Earth years."
1299	write(,'(a,f16.4,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(layerings_map(i,islope))
1307	enddo
1308	enddo
1309	endif
1310	deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
1311	deallocate(co2_ice,h2o_ice,layerings_map)
1312	!----------------------------- END OUTPUT ------------------------------
1313
1314	END PROGRAM pem

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

Download in other formats:

Original Format