Context Navigation

pem.F90 @ 3921

Last change on this file since 3921 was 3912, checked in by jbclement, 5 months ago
Mars PCM 1D: Prescription of the atmospheric water ice profile in 1D. The boolean 'ctrl_h2oice' activates the option, the file "profile_ref_h2ovap" defines the prescribed profile and the value 'relaxtime_h2oice', if positive, activates the relaxation and gives the constant of time. JBC
File size: 66.5 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	use dyn_ss_ice_m_mod, only: dyn_ss_ice_m
74	use parse_args_mod, only: parse_args
75	use job_timelimit_mod, only: timelimit, antetime, timewall
76	use paleoclimate_mod, only: h2o_ice_depth, zdqsdif_ssi_tot
77	use surf_temp, only: update_tsurf_nearest_baresoil
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, ctrl_h2ovap, ctrl_h2oice
279	integer :: ndt, day0
280	real :: ptif, pks, day, gru, grv, relaxtime_h2ovap, relaxtime_h2oice
281	real, dimension(:), allocatable :: zqsat
282	real, dimension(:,:,:), allocatable :: dq, dqdyn
283	real, dimension(nlayer) :: play, w, qref_h2ovap, qref_h2oice
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	real :: totmass_ini
295	logical :: num_str
296
297	! CODE
298	! Elapsed time with system clock
299	call system_clock(count_rate = cr)
300	call system_clock(c1)
301
302	! Parse command-line options
303	call parse_args()
304
305	! Header
306	write(,) ' * . . + . * . + . . . '
307	write(,) ' + _______ ________ ____ ____ * + '
308	write(,) ' + . * \|_ __ \\|_ __ \|\|_ \ / _\| . *'
309	write(,) ' . . \| \|__) \| \| \|_ \_\| \| \/ \| * * . '
310	write(,) ' . \| ___/ \| _\| _ \| \|\ /\| \| . . '
311	write(,) '. * * _\| \|_ _\| \|__/ \| _\| \|_\/_\| \|_ * '
312	write(,) ' + \|_____\| \|________\|\|_____\|\|_____\| + . '
313	write(,) ' . * . * . + * . + .'
314
315	! Some user info
316	call date_and_time(date,time,zone,values)
317	call getcwd(dir) ! Current directory
318	call getlog(logname) ! User name
319	call hostnm(hostname) ! Machine/station name
320	write(,)
321	write(,) '******* PEM information *******'
322	write(,) '+ User : '//trim(logname)
323	write(,) '+ Machine : '//trim(hostname)
324	write(,) '+ Directory: '//trim(dir)
325	write(*,'(a,i2,a,i2,a,i4)') ' + Date : ',values(3),'/',values(2),'/',values(1)
326	write(*,'(a,i2,a,i2,a,i2,a)') ' + Time : ',values(5),':',values(6),':',values(7)
327
328	! Parallel variables
329	#ifndef CPP_STD
330	is_sequential = .true.
331	is_parallel = .false.
332	is_mpi_root = .true.
333	is_omp_root = .true.
334	is_master = .true.
335	#endif
336
337	! Some constants
338	day_ini = 0
339	time_phys = 0.
340	ngrid = ngridmx
341	A = (1./m_co2 - 1./m_noco2)
342	B = 1./m_noco2
343	year_day = 669
344	daysec = 88775.
345	timestep = year_daydaysecyear_step
346
347	!----------------------------- INITIALIZATION --------------------------
348	!------------------------
349	! I Initialization
350	! I_a Read the "run.def"
351	!------------------------
352	write(,)
353	write(,) '******* PEM initialization *******'
354	write(,) '> Reading "run.def" (PEM)'
355	#ifndef CPP_1D
356	dtphys = daysec/48. ! Dummy value (overwritten in phyetat0)
357	call conf_gcm(99,.true.)
358	call infotrac_init
359	nq = nqtot
360	allocate(q(ip1jmp1,llm,nqtot))
361	allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
362	#else
363	allocate(q(1,llm,nqtot),pdyn(1))
364	allocate(longitude(1),latitude(1),cell_area(1))
365
366	therestart1D = .false. ! Default value
367	inquire(file = 'start1D.txt',exist = therestart1D)
368	if (.not. therestart1D) then
369	write(,) 'There is no "start1D.txt" file!'
370	error stop 'Initialization cannot be done for the 1D PEM.'
371	endif
372	therestartfi = .false. ! Default value
373	inquire(file = 'startfi.nc',exist = therestartfi)
374	if (.not. therestartfi) then
375	write(,) 'There is no "startfi.nc" file!'
376	error stop 'Initialization cannot be done for the 1D PEM.'
377	endif
378
379	write(,) '> Reading "start1D.txt" and "startfi.nc"'
380	call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, &
381	time_0,pdyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
382	play,plev,latitude,longitude,cell_area, &
383	ctrl_h2ovap,relaxtime_h2ovap,qref_h2ovap,ctrl_h2oice,relaxtime_h2oice,qref_h2oice)
384	nsplit_phys = 1
385	day_step = steps_per_sol
386	#endif
387
388	call conf_pem(i_myear,n_myear)
389
390	!------------------------
391	! I Initialization
392	! I_b Read of the "start.nc" and "starfi.nc"
393	!------------------------
394	! I_b.1 Read "start.nc"
395	write(,) '> Reading "start.nc"'
396	allocate(ps_start0(ngrid))
397	#ifndef CPP_1D
398	allocate(pdyn(ip1jmp1))
399	call dynetat0(start_name,vcov,ucov,teta,q,masse,pdyn,phis,time_0)
400
401	call gr_dyn_fi(1,iip1,jjp1,ngridmx,pdyn,ps_start0)
402
403	call iniconst ! Initialization of dynamical constants (comconst_mod)
404	call inigeom ! Initialization of geometry
405
406	allocate(ap(nlayer + 1))
407	allocate(bp(nlayer + 1))
408	status = nf90_open(start_name,NF90_NOWRITE,ncid)
409	status = nf90_inq_varid(ncid,"ap",apvarid)
410	status = nf90_get_var(ncid,apvarid,ap)
411	status = nf90_inq_varid(ncid,"bp",bpvarid)
412	status = nf90_get_var(ncid,bpvarid,bp)
413	status = nf90_close(ncid)
414
415	! Initialization of physics constants and variables (comcstfi_h)
416	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)
417	#else
418	ps_start0(1) = pdyn(1)
419	#endif
420	deallocate(pdyn)
421
422	! In the PCM, these values are given to the physic by the dynamic.
423	! Here we simply read them in the "startfi.nc" file
424	status = nf90_open(startfi_name,NF90_NOWRITE,ncid)
425	status = nf90_inq_varid(ncid,"longitude",lonvarid)
426	status = nf90_get_var(ncid,lonvarid,longitude)
427	status = nf90_inq_varid(ncid,"latitude",latvarid)
428	status = nf90_get_var(ncid,latvarid,latitude)
429	status = nf90_inq_varid(ncid,"area",areavarid)
430	status = nf90_get_var(ncid,areavarid,cell_area)
431	status = nf90_inq_varid(ncid,"soildepth",sdvarid)
432	status = nf90_get_var(ncid,sdvarid,mlayer)
433	status = nf90_close(ncid)
434
435	! I_b.2 Read the "startfi.nc"
436	! First we read the initial state (starfi.nc)
437	#ifndef CPP_STD
438	write(,) '> Reading "startfi.nc"'
439	call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
440	tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, &
441	watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
442
443	! Remove unphysical values of surface tracer
444	where (qsurf < 0.) qsurf = 0.
445
446	call surfini(ngrid,nslope,qsurf)
447	#else
448	call phys_state_var_init(nq)
449	if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
450	call initracer(ngrid,nq)
451	call iniaerosol()
452	allocate(tsurf_read_generic(ngrid))
453	allocate(qsurf_read_generic(ngrid,nq))
454	allocate(tsoil_read_generic(ngrid,nsoilmx))
455	allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
456	allocate(emis_read_generic(ngrid))
457	allocate(albedo_read_generic(ngrid,2))
458	allocate(qsurf(ngrid,nq,1))
459	allocate(tsurf(ngrid,1))
460	allocate(tsoil(ngrid,nsoilmx,1))
461	allocate(emis(ngrid,1))
462	allocate(watercap(ngrid,1))
463	allocate(watercaptag(ngrid))
464	allocate(albedo(ngrid,2,1))
465	allocate(inertiesoil(ngrid,nsoilmx,1))
466	call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
467	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
468	qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, &
469	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
470	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
471
472	nslope = 1
473	call ini_comslope_h(ngrid,1)
474
475	qsurf(:,:,1) = qsurf_read_generic
476	tsurf(:,1) = tsurf_read_generic
477	tsoil(:,:,1) = tsoil_read_generic
478	emis(:,1) = emis_read_generic
479	watercap(:,1) = 0.
480	watercaptag(:) = .false.
481	albedo(:,1,1) = albedo_read_generic(:,1)
482	albedo(:,2,1) = albedo_read_generic(:,2)
483	inertiesoil(:,:,1) = inertiedat
484
485	if (nslope == 1) then
486	def_slope(1) = 0
487	def_slope(2) = 0
488	def_slope_mean = 0
489	subslope_dist(:,1) = 1.
490	endif
491
492	! Remove unphysical values of surface tracer
493	qsurf(:,:,1) = qsurf_read_generic
494	where (qsurf < 0.) qsurf = 0.
495
496	deallocate(tsurf_read_generic,qsurf_read_generic,qsoil_read_generic,emis_read_generic)
497	#endif
498
499	do nnq = 1,nqtot ! Why not using ini_tracer?
500	if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
501	if (noms(nnq) == "h2o_vap") then
502	igcm_h2o_vap = nnq
503	mmol(igcm_h2o_vap) = 18.
504	endif
505	if (noms(nnq) == "co2") igcm_co2 = nnq
506	enddo
507
508	!------------------------
509	! I Initialization
510	! I_c Subslope parametrisation
511	!------------------------
512	! Define some slope statistics
513	iflat = 1
514	do islope = 2,nslope
515	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
516	enddo
517	write(,) 'Flat slope for islope = ',iflat
518	write(,) 'Corresponding criterium = ',def_slope_mean(iflat)
519
520	!------------------------
521	! I Initialization
522	! I_d Read the PCM data and convert them to the physical grid
523	!------------------------
524	! 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
525	call get_timelen_PCM("data_PCM_Y1.nc",timelen)
526
527	allocate(vmr_co2_PCM(ngrid,timelen),q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
528	allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid))
529	allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
530	allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope))
531	allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_timeseries(ngrid,nsoilmx,nslope,timelen))
532	allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
533
534	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, &
535	tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2_PEM_phys,q_h2o_PEM_phys,watersurf_density_avg,watersoil_density_timeseries)
536
537	! Compute the deviation from the average
538	allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
539	ps_dev = ps_start0 - ps_avg
540	tsurf_dev = tsurf - tsurf_avg
541	tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:)
542
543	!------------------------
544	! I Initialization
545	! I_e Initialization of the PEM variables and soil
546	!------------------------
547	call end_comsoil_h_PEM
548	call ini_comsoil_h_PEM(ngrid,nslope)
549	call end_adsorption_h_PEM
550	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
551	call end_ice_table
552	call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
553
554	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))
555	if (soil_pem) then
556	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
557	tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
558	watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
559	tsoil_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_timeseries
560	do l = nsoilmx + 1,nsoilmx_PEM
561	tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
562	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
563	tsoil_PEM_timeseries(:,l,:,:) = tsoil_timeseries(:,nsoilmx,:,:)
564	enddo
565	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
566	endif !soil_pem
567	deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries)
568
569	!------------------------
570	! I Initialization
571	! I_f Compute tendencies
572	!------------------------
573	allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope))
574	call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
575	call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
576	d_co2ice_ini = d_co2ice
577	deallocate(min_co2_ice,min_h2o_ice)
578
579	!------------------------
580	! I Initialization
581	! I_g Compute global surface pressure
582	!------------------------
583	total_surface = sum(cell_area)
584	ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
585	ps_avg_global_new = ps_avg_global_ini
586	write(,) "Total surface of the planet =", total_surface
587	write(,) "Initial global average pressure =", ps_avg_global_ini
588
589	!------------------------
590	! I Initialization
591	! I_h Read the "startpem.nc"
592	!------------------------
593	write(,) '> Reading "startpem.nc"'
594	allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),layerings_map(ngrid,nslope))
595	allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid))
596	delta_h2o_icetablesublim = 0.
597	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
598	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, &
599	watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,layerings_map)
600	deallocate(tsurf_avg_yr1)
601
602	! We save the places where h2o ice is sublimating
603	! We compute the surface of h2o ice sublimating
604	allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope),co2ice_disappeared(ngrid,nslope))
605	co2ice_sublim_surf_ini = 0.
606	h2oice_ini_surf = 0.
607	is_co2ice_sublim_ini = .false.
608	is_h2oice_sublim_ini = .false.
609	is_co2ice_ini = .false.
610	co2ice_disappeared = .false.
611	totmass_co2ice_ini = 0.
612	totmass_atmco2_ini = 0.
613	if (layering_algo) then
614	do ig = 1,ngrid
615	do islope = 1,nslope
616	if (is_co2ice_str(layerings_map(ig,islope)%top)) then
617	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
618	else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
619	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
620	else
621	call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
622	endif
623	enddo
624	enddo
625	! We put the sublimating tendency coming from subsurface ice into the overall tendency
626	where (h2o_ice_depth > 0. .and. zdqsdif_ssi_tot < -1.e-10) d_h2oice = zdqsdif_ssi_tot
627	endif
628	do i = 1,ngrid
629	totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
630	do islope = 1,nslope
631	totmass_co2ice_ini = totmass_co2ice_ini + co2_ice(i,islope)cell_area(i)subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)
632	if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
633	if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
634	is_co2ice_sublim_ini(i,islope) = .true.
635	co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
636	endif
637	if (d_h2oice(i,islope) < 0.) then
638	if (h2o_ice(i,islope) > 0.) then
639	is_h2oice_sublim_ini(i,islope) = .true.
640	h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
641	else if (h2o_ice_depth(i,islope) > 0.) then
642	is_h2oice_sublim_ini(i,islope) = .true.
643	endif
644	endif
645	enddo
646	enddo
647	write(,) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini
648	write(,) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf
649
650	totmass_adsco2_ini = 0.
651	totmass_adsh2o = 0.
652	if (adsorption_pem) then
653	do ig = 1,ngrid
654	do islope = 1,nslope
655	do l = 1,nsoilmx_PEM - 1
656	if (l == 1) then
657	totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
658	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
659	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
660	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
661	else
662	totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
663	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
664	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
665	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
666	endif
667	enddo
668	enddo
669	enddo
670	totmass_adsco2 = totmass_adsco2_ini
671	write(,) "Tot mass of CO2 in the regolith =", totmass_adsco2
672	write(,) "Tot mass of H2O in the regolith =", totmass_adsh2o
673	endif ! adsorption
674
675	!------------------------
676	! I Initialization
677	! I_i Compute orbit criterion
678	!------------------------
679	#ifndef CPP_STD
680	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
681	#else
682	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
683	#endif
684
685	n_myear_leg = Max_iter_pem
686	if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
687
688	!-------------------------- END INITIALIZATION -------------------------
689
690	!-------------------------------- RUN ----------------------------------
691	!------------------------
692	! II Run
693	! II_a Update pressure, ice and tracers
694	!------------------------
695	write(,)
696	write(,) '******* PEM cycle *******'
697	i_myear_leg = 0
698	stopPEM = 0
699	if (layering_algo) then
700	allocate(h2o_ice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
701	new_str = .true.
702	new_lag = .true.
703	do islope = 1,nslope
704	do ig = 1,ngrid
705	current(ig,islope)%p => layerings_map(ig,islope)%top
706	enddo
707	enddo
708	endif
709
710	do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
711	! II.a.1. Compute updated global pressure
712	write(,'(a,f10.2)') ' *** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
713	write(,) "> Updating the surface pressure"
714	ps_avg_global_old = ps_avg_global_new
715	do i = 1,ngrid
716	do islope = 1,nslope
717	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
718	enddo
719	enddo
720	if (adsorption_pem) then
721	do i = 1,ngrid
722	ps_avg_global_new = ps_avg_global_new - gcell_area(i)delta_co2_adsorbed(i)/total_surface
723	enddo
724	endif
725	ps_avg = ps_avg*ps_avg_global_new/ps_avg_global_old
726	write(,) 'Global average pressure old time step:',ps_avg_global_old
727	write(,) 'Global average pressure new time step:',ps_avg_global_new
728
729	! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
730	write(,) "> Updating the surface pressure timeseries for the new pressure"
731	allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
732	do l = 1,nlayer + 1
733	do ig = 1,ngrid
734	zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
735	enddo
736	enddo
737	ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old
738	write(,) "> Updating the pressure levels timeseries for the new pressure"
739	allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
740	do l = 1,nlayer + 1
741	do ig = 1,ngrid
742	zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
743	enddo
744	enddo
745
746	! II.a.3. Tracers timeseries
747	write(,) "> Updating the tracer VMR timeseries for the new pressure"
748	allocate(vmr_co2_PEM_phys(ngrid,timelen))
749	l = 1
750	do ig = 1,ngrid
751	do t = 1,timelen
752	! H2O
753	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))/ &
754	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
755	if (q_h2o_PEM_phys(ig,t) < 0) then
756	q_h2o_PEM_phys(ig,t) = 1.e-30
757	else if (q_h2o_PEM_phys(ig,t) > 1) then
758	q_h2o_PEM_phys(ig,t) = 1.
759	endif
760	! CO2
761	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))/ &
762	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
763	+ ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
764	- (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	if (q_co2_PEM_phys(ig,t) < 0) then
767	q_co2_PEM_phys(ig,t) = 1.e-30
768	else if (q_co2_PEM_phys(ig,t) > 1) then
769	q_co2_PEM_phys(ig,t) = 1.
770	endif
771	mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
772	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
773	enddo
774	enddo
775	deallocate(zplev_timeseries_new,zplev_timeseries_old)
776
777	!------------------------
778	! II Run
779	! II_b Evolution of ice
780	!------------------------
781	allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
782	if (layering_algo) then
783	h2o_ice_depth_old = h2o_ice_depth
784	do islope = 1,nslope
785	do ig = 1,ngrid
786	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)
787	!call print_layering(layerings_map(ig,islope))
788	co2_ice(ig,islope) = 0.
789	h2o_ice(ig,islope) = 0.
790	h2o_ice_depth(ig,islope) = 0.
791	if (is_co2ice_str(layerings_map(ig,islope)%top)) then
792	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
793	else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
794	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
795	else
796	call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
797	endif
798	enddo
799	enddo
800	else
801	zlag = 0.
802	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
803	call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
804	endif
805
806	!------------------------
807	! II Run
808	! II_c Flow of glaciers
809	!------------------------
810	allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
811	if (co2ice_flow .and. nslope > 1) then
812	call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
813	ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
814	if (layering_algo) then
815	do islope = 1,nslope
816	do ig = 1,ngrid
817	layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
818	enddo
819	enddo
820	endif
821	endif
822	if (h2oice_flow .and. nslope > 1) then
823	call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
824	if (layering_algo) then
825	do islope = 1,nslope
826	do ig = 1,ngrid
827	!~ layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
828	enddo
829	enddo
830	endif
831	endif
832
833	!------------------------
834	! II Run
835	! II_d Update surface and soil temperatures
836	!------------------------
837	! II_d.1 Update Tsurf
838	call update_tsurf_nearest_baresoil(ngrid,nslope,iim,jjm,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
839
840	if (soil_pem) then
841	! II_d.2 Shifting soil temperature to surface
842	call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
843
844	! II_d.3 Update soil temperature
845	write(,)"> Updating soil temperature profile"
846	allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen))
847	tsoil_PEM_timeseries_old = tsoil_PEM_timeseries
848	do islope = 1,nslope
849	tsoil_avg_old = tsoil_PEM(:,:,islope)
850	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
851	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
852
853	do t = 1,timelen
854	do ig = 1,ngrid
855	do isoil = 1,nsoilmx_PEM
856	! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
857	tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
858	! Update of watersoil density
859	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)
860	if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
861	enddo
862	enddo
863	enddo
864	enddo
865	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
866	deallocate(tsoil_avg_old)
867
868	! II_d.4 Update the ice table
869	allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope))
870	if (icetable_equilibrium) then
871	write(,) "> Updating ice table (equilibrium method)"
872	icetable_thickness_old = icetable_thickness
873	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
874	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
875	else if (icetable_dynamic) then
876	write(,) "> Updating ice table (dynamic method)"
877	ice_porefilling_old = ice_porefilling
878	icetable_depth_old = icetable_depth
879	allocate(porefill(nsoilmx_PEM))
880	do ig = 1,ngrid
881	do islope = 1,nslope
882	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)
883	icetable_depth(ig,islope) = ssi_depth
884	ice_porefilling(ig,:,islope) = porefill
885	enddo
886	enddo
887	deallocate(porefill)
888	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
889	endif
890	deallocate(icetable_thickness_old,ice_porefilling_old)
891
892	! II_d.5 Update the soil thermal properties
893	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)
894
895	! II_d.6 Update the mass of the regolith adsorbed
896	totmass_adsco2 = 0.
897	totmass_adsh2o = 0.
898	if (adsorption_pem) then
899	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
900	tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
901	h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
902	do ig = 1,ngrid
903	do islope = 1,nslope
904	do l = 1,nsoilmx_PEM
905	if (l == 1) then
906	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
907	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
908	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
909	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
910	else
911	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
912	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
913	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
914	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
915	endif
916	enddo
917	enddo
918	enddo
919	write(,) "Total mass of CO2 in the regolith =", totmass_adsco2
920	write(,) "Total mass of H2O in the regolith =", totmass_adsh2o
921	endif
922	endif !soil_pem
923	deallocate(zshift_surf,zlag)
924
925	!------------------------
926	! II Run
927	! II_e Outputs
928	!------------------------
929	call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
930	do islope = 1,nslope
931	write(str2(1:2),'(i2.2)') islope
932	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
933	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
934	call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
935	call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
936	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope)))
937	call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope)))
938	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
939	if (icetable_equilibrium) then
940	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
941	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope))
942	else if (icetable_dynamic) then
943	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
944	endif
945
946	if (soil_pem) then
947	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
948	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
949	if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
950	if (adsorption_pem) then
951	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope))
952	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope))
953	endif
954	endif
955	enddo
956	deallocate(flag_co2flow,flag_h2oflow)
957
958	! Checking mass balance for CO2
959	if (abs(CO2cond_ps - 1.) < 1.e-10) then
960	totmass_co2ice = 0.
961	totmass_atmco2 = 0.
962	do ig = 1,ngrid
963	totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
964	do islope = 1,nslope
965	totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)cell_area(ig)subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
966	enddo
967	enddo
968	totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,1.e-10)
969	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, ' %'
970	if ((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini > 0.01) then
971	write(,) ' /!\ Warning: mass balance is not conseved!'
972	totmass_ini = max(totmass_atmco2_ini,1.e-10)
973	write(,'(a,f8.3,a)') ' Atmospheric CO2 mass balance = ', 100.(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini, ' %'
974	totmass_ini = max(totmass_co2ice_ini,1.e-10)
975	write(,'(a,f8.3,a)') ' CO2 ice mass balance = ', 100.(totmass_co2ice - totmass_co2ice_ini)/totmass_ini, ' %'
976	totmass_ini = max(totmass_adsco2_ini,1.e-10)
977	write(,'(a,f8.3,a)') ' Adsorbed CO2 mass balance = ', 100.(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini, ' %'
978	endif
979	endif
980
981	!------------------------
982	! II Run
983	! II_f Update the tendencies
984	!------------------------
985	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)
986	write(,) "> Updating the H2O sub-surface ice tendency due to lag layer"
987	if (layering_algo) then
988	do ig = 1,ngrid
989	do islope = 1,nslope
990	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))
991	enddo
992	enddo
993	!~ else
994	!~ do ig = 1,ngrid
995	!~ do islope = 1,nslope
996	!~ 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))
997	!~ enddo
998	!~ enddo
999	endif
1000	if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old)
1001	deallocate(vmr_co2_PEM_phys)
1002
1003	!------------------------
1004	! II Run
1005	! II_g Checking the stopping criterion
1006	!------------------------
1007	write(,) "> Checking the stopping criteria"
1008	call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
1009	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)
1010	i_myear_leg = i_myear_leg + dt
1011	i_myear = i_myear + dt
1012	if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
1013	if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
1014	call system_clock(c2)
1015	if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
1016	if (stopPEM > 0) then
1017	select case (stopPEM)
1018	case(1)
1019	write(,'(a,i0,a)') " *** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
1020	case(2)
1021	write(,'(a,i0,a)') " *** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
1022	case(3)
1023	write(,'(a,i0,a)') " *** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
1024	case(4)
1025	write(,'(a,i0,a)') " *** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
1026	case(5)
1027	write(,'(a,i0)') " *** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
1028	case(6)
1029	write(,'(a,i0)') " *** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
1030	case(7)
1031	write(,'(a,i0)') " *** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
1032	case default
1033	write(,'(a,i0)') " *** STOPPING with unknown stopping criterion code: ", stopPEM
1034	end select
1035	exit
1036	else
1037	write(,'(a,f10.2,a)') ' *** The chained simulation has run for ',i_myear,' Martian years.'
1038	write(,) '**** The PEM can continue!'
1039	write(,) '****'
1040	endif
1041	enddo ! big time iteration loop of the pem
1042	deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed)
1043	deallocate(watersoil_density_PEM_avg,watersurf_density_avg)
1044	deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries)
1045	deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim)
1046	deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
1047	deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
1048	if (layering_algo) deallocate(h2o_ice_depth_old,new_str,new_lag,current)
1049	!------------------------------ END RUN --------------------------------
1050
1051	!------------------------------- OUTPUT --------------------------------
1052	!------------------------
1053	! III Output
1054	! III_a Update surface values for the PCM start files
1055	!------------------------
1056	write(,)
1057	write(,) '******* PEM finalization *******'
1058	! III_a.1 Ice update for start file
1059	write(,) '> Reconstructing perennial ice and frost for the PCM'
1060	watercap = 0.
1061	perennial_co2ice = co2_ice
1062	do ig = 1,ngrid
1063	! H2O ice metamorphism
1064	!if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_h2oice_threshold) then
1065	! h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1066	! qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
1067	!endif
1068
1069	! Is H2O ice still considered as an infinite reservoir for the PCM?
1070	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
1071	! There is enough ice to be considered as an infinite reservoir
1072	watercaptag(ig) = .true.
1073	else
1074	! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
1075	watercaptag(ig) = .false.
1076	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1077	h2o_ice(ig,:) = 0.
1078	endif
1079
1080	! CO2 ice metamorphism
1081	!if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_co2ice_threshold) then
1082	! perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1083	! qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
1084	!endif
1085	enddo
1086
1087	! III.a.3. Tsurf update for start file
1088	write(,) '> Reconstructing the surface temperature for the PCM'
1089	tsurf = tsurf_avg + tsurf_dev
1090	deallocate(tsurf_dev)
1091
1092	! III_a.4 Tsoil update for start file
1093	if (soil_pem) then
1094	write(,) '> Reconstructing the soil temperature profile for the PCM'
1095	inertiesoil = TI_PEM(:,:nsoilmx,:)
1096	! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value
1097	do isoil = 1,nsoilmx
1098	tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:)
1099	enddo
1100	tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev
1101	#ifndef CPP_STD
1102	flux_geo = fluxgeo
1103	#endif
1104	endif
1105	deallocate(tsurf_avg,tsoil_dev)
1106
1107	! III_a.5 Pressure update for start file
1108	write(,) '> Reconstructing the pressure for the PCM'
1109	allocate(ps_start(ngrid))
1110	! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop
1111	ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini
1112	deallocate(ps_avg,ps_dev)
1113
1114	! III_a.6 Tracers update for start file
1115	write(,) '> Reconstructing the tracer VMR for the PCM'
1116	allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
1117	do l = 1,nlayer + 1
1118	zplev_start0(:,l) = ap(l) + bp(l)*ps_start0
1119	zplev_new(:,l) = ap(l) + bp(l)*ps_start
1120	enddo
1121
1122	do nnq = 1,nqtot
1123	if (noms(nnq) /= "co2") then
1124	do l = 1,llm - 1
1125	do ig = 1,ngrid
1126	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))
1127	enddo
1128	q(:,llm,nnq) = q(:,llm - 1,nnq)
1129	enddo
1130	else
1131	do l = 1,llm - 1
1132	do ig = 1,ngrid
1133	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)) &
1134	+ ((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))
1135	enddo
1136	q(:,llm,nnq) = q(:,llm - 1,nnq)
1137	enddo
1138	endif
1139	enddo
1140	deallocate(zplev_start0)
1141
1142	! Conserving the tracers mass for start file
1143	do nnq = 1,nqtot
1144	do ig = 1,ngrid
1145	do l = 1,llm - 1
1146	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
1147	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1148	q(ig,l,nnq) = 1.
1149	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1150	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1151	endif
1152	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1153	enddo
1154	enddo
1155	enddo
1156	deallocate(zplev_new)
1157
1158	! III_a.7 Albedo update for start file
1159	write(,) '> Reconstructing the albedo for the PCM'
1160	do ig = 1,ngrid
1161	if (latitude(ig) < 0.) then
1162	icap = 2 ! Southern hemisphere
1163	else
1164	icap = 1 ! Northern hemisphere
1165	endif
1166	do islope = 1,nslope
1167	! Bare ground
1168	albedo(ig,:,islope) = albedodat(ig)
1169	emis(ig,islope) = emissiv
1170
1171	! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant
1172	! H2O ice
1173	if (h2o_ice(ig,islope) > 0.) then
1174	albedo(ig,:,islope) = albedo_h2o_cap
1175	emis(ig,islope) = 1.
1176	endif
1177	! CO2 ice
1178	if (co2_ice(ig,islope) > 0.) then
1179	albedo(ig,:,islope) = albedo_perennialco2(icap)
1180	emis(ig,islope) = emisice(icap)
1181	endif
1182	! H2O frost
1183	if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
1184	albedo(ig,:,islope) = albedo_h2o_frost
1185	emis(ig,islope) = 1.
1186	endif
1187	! CO2 frost
1188	if (qsurf(ig,igcm_co2,islope) > 0.) then
1189	albedo(ig,:,islope) = albedice(icap)
1190	emis(ig,islope) = emisice(icap)
1191	endif
1192	enddo
1193	enddo
1194
1195	! III_a.8 Orbital parameters update for start file
1196	write(,) '> Setting the new orbital parameters'
1197	if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
1198
1199	!------------------------
1200	! III Output
1201	! III_b Write "restart.nc" and "restartfi.nc"
1202	!------------------------
1203	! III_b.1 Write "restart.nc"
1204	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1205	pday = day_ini
1206	ztime_fin = time_phys
1207	#ifndef CPP_1D
1208	write(,) '> Writing "restart.nc"'
1209	! Correction on teta due to surface pressure changes
1210	allocate(pdyn(ip1jmp1))
1211	call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn)
1212	do i = 1,ip1jmp1
1213	teta(i,:) = teta(i,:)pdyn(i)*rcp
1214	enddo
1215	! Correction on atmospheric pressure
1216	allocate(p(ip1jmp1,nlayer + 1))
1217	call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start,pdyn)
1218	call pression(ip1jmp1,ap,bp,pdyn,p)
1219	! Correction on the mass of atmosphere
1220	call massdair(p,masse)
1221	call dynredem0("restart.nc",day_ini,phis)
1222	call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,pdyn)
1223	deallocate(ap,bp,p,pdyn)
1224	#else
1225	write(,) '> Writing "restart1D.txt"'
1226	call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1227	#endif
1228	deallocate(ps_start0,ps_start)
1229
1230	! III_b.2 Write the "restartfi.nc"
1231	write(,) '> Writing "restartfi.nc"'
1232	#ifndef CPP_STD
1233	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1234	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1235	inertiedat,def_slope,subslope_dist)
1236	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1237	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1238	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1239	wstar,watercap,perennial_co2ice)
1240	#else
1241	if (allocated(noms)) deallocate(noms)
1242	deallocate(qsurf,tsurf,tsoil,emis,watercap,watercaptag,albedo,inertiesoil)
1243	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1244	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
1245	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1246	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1247	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1248	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
1249	tsea_ice,sea_ice)
1250	#endif
1251
1252	!------------------------
1253	! III Output
1254	! III_c Write the "restartpem.nc"
1255	!------------------------
1256	write(,) '> Writing "restartpem.nc"'
1257	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)
1258	call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1259	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
1260	co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
1261
1262	call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
1263
1264	write(,)
1265	write(,) '**** PEM final information *****'
1266	write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
1267	write(,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myearconvert_years, " Earth years."
1268	write(,'(a,f16.4,a)') " + The reached date is now ", (year_bp_ini + i_myear)convert_years, " Earth years."
1269	write(,) "+ PEM: so far, so good!"
1270	write(,) '************************************'
1271
1272	if (layering_algo) then
1273	do islope = 1,nslope
1274	do i = 1,ngrid
1275	call del_layering(layerings_map(i,islope))
1276	enddo
1277	enddo
1278	endif
1279	deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
1280	deallocate(co2_ice,h2o_ice,layerings_map)
1281	!----------------------------- END OUTPUT ------------------------------
1282
1283	END PROGRAM pem

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

Download in other formats:

Original Format