Context Navigation

pem.F90 @ 3842

Last change on this file since 3842 was 3842, checked in by jbclement, 44 hours ago

PEM:

Handling the situation in the layering algorithm when CO2 ice sublimation and H2O ice condensation happen simultaneously.
Correction of the incorporation of the sub-surface sublimation tendency into the overall tendency to be more robust.
Revision of the layering initialization in the case where there is no "startpem.nc" file.

JBC

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

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

Download in other formats:

Original Format