Context Navigation

pem.F90 @ 3939

Last change on this file since 3939 was 3938, checked in by jbclement, 3 months ago

PEM:

Correction of the process to balance the H2O flux from and into the atmosphere accross reservoirs: (i) computation of H2O amount going from/in the atmosphere is corrected (missing dt and wrong condition); (ii) computation of the scaling factor is corrected because we need to balance all the icetable/adsorbed/surface ice flux but it affects only sublimating/condensing surface ice flux; (iii) process of balancing is modified to be applied to every flux by scaling instead of applying it only to positive or negative flux by trimming; (iv) balance is now fully processed by the tendency before evolving the ice. This is done through dedicated subroutines.
Correction of ice conversion between the layering algorithm and PEM ice variables (density factor was missing).
Addition of H2O flux balance and related stopping criterion for the layering algorithm.
Few updates for files in the deftank to be in line with the wiki Planeto pages.

JBC

File size: 67.4 KB

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

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

Download in other formats:

Original Format