Context Navigation

pem.F90 @ 3851

Last change on this file since 3851 was 3851, checked in by jbclement, 4 days ago

PEM:

Making the computation of CO2 mass balance more robust, especially regarding 'CO2cond_ps'.
Small correction about the dust tendency management for the layering algorithm.
Small improvement of the visualization in "visu_evol_layering.py".
File "log_launchPEM.txt" renamed into "launchPEM.log".

JBC

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

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

Download in other formats:

Original Format