Context Navigation

pem.F90 @ 3910

Last change on this file since 3910 was 3907, checked in by jbclement, 5 months ago
PEM: Correction for the update of surface temperature when ice is disappearing: now it is really set as the surface temperature of the nearest bare soil point because we check the lon-lat grid and not the reduced grid. JBC
File size: 66.4 KB

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

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

Download in other formats:

Original Format