Context Navigation

pem.F90 @ 3525

Last change on this file since 3525 was 3525, checked in by jbclement, 4 days ago
PEM: Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings. JBC
File size: 61.8 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 Save the initial situation
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 value 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
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_adsorbded_phys, h2o_adsorbded_phys ! Mass of co2 and h2O adsorbded
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 interpol_TI_PEM2PCM_mod, only: interpol_TI_PEM2PCM
65	use nb_time_step_PCM_mod, only: nb_time_step_PCM
66	use pemetat0_mod, only: pemetat0
67	use read_data_PCM_mod, only: read_data_PCM
68	use recomp_tend_co2_slope_mod, only: recomp_tend_co2_slope
69	use compute_soiltemp_mod, only: compute_tsoil_pem
70	use writediagpem_mod, only: writediagpem, writediagsoilpem
71	use co2condens_mod, only: CO2cond_ps
72	use layering_mod, only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
73
74	#ifndef CPP_STD
75	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil
76	use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
77	albedodat, zmea, zstd, zsig, zgam, zthe, &
78	albedo_h2o_frost,frost_albedo_threshold, &
79	emissiv, watercaptag, perennial_co2ice, emisice, albedice
80	use dimradmars_mod, only: totcloudfrac, albedo
81	use dust_param_mod, only: tauscaling
82	use tracer_mod, only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
83	use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
84	use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
85	use comcstfi_h, only: pi, rad, g, mugaz, r
86	use surfini_mod, only: surfini
87	use comconst_mod, only: kappa, cpp
88	#else
89	use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
90	use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, &
91	PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
92	ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
93	use aerosol_mod, only: iniaerosol
94	use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit
95	use comcstfi_mod, only: pi, rad, g, mugaz, r
96	#endif
97
98	#ifndef CPP_1D
99	use iniphysiq_mod, only: iniphysiq
100	use control_mod, only: iphysiq, day_step, nsplit_phys
101	#else
102	use time_phylmdz_mod, only: iphysiq, steps_per_sol
103	use regular_lonlat_mod, only: init_regular_lonlat
104	use physics_distribution_mod, only: init_physics_distribution
105	use mod_grid_phy_lmdz, only: regular_lonlat
106	use init_testphys1d_mod, only: init_testphys1d
107	use comvert_mod, only: ap, bp
108	use writerestart1D_mod, only: writerestart1D
109	#endif
110
111	implicit none
112
113	include "dimensions.h"
114	include "paramet.h"
115	include "comgeom.h"
116	include "iniprint.h"
117	include "callkeys.h"
118
119	integer ngridmx
120	parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
121
122	! Same variable names as in the PCM
123	integer, parameter :: nlayer = llm ! Number of vertical layer
124	integer :: ngrid ! Number of physical grid points
125	integer :: nq ! Number of tracer
126	integer :: day_ini ! First day of the simulation
127	real :: pday ! Physical day
128	real :: time_phys ! Same as in PCM
129	real :: ptimestep ! Same as in PCM
130	real :: ztime_fin ! Same as in PCM
131
132	! Variables to read start.nc
133	character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM
134
135	! Dynamic variables
136	real, dimension(ip1jm,llm) :: vcov ! vents covariants
137	real, dimension(ip1jmp1,llm) :: ucov ! vents covariants
138	real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle
139	real, dimension(:,:,:), allocatable :: q ! champs advectes
140	real, dimension(ip1jmp1) :: ps ! pression au sol
141	real, dimension(ip1jmp1) :: ps0 ! pression au sol initiale
142	real, dimension(:), allocatable :: ps_start_PCM ! (ngrid) surface pressure
143	real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure
144	real, dimension(ip1jmp1,llm) :: masse ! masse d'air
145	real, dimension(ip1jmp1) :: phis ! geopotentiel au sol
146	real :: time_0
147
148	! Variables to read starfi.nc
149	character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM
150	character(2) :: str2
151	integer :: ncid, status ! Variable for handling opening of files
152	integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files
153	integer :: apvarid, bpvarid ! Variable ID for Netcdf files
154
155	! Variables to read starfi.nc and write restartfi.nc
156	real, dimension(:), allocatable :: longitude ! Longitude read in startfi_name and written in restartfi
157	real, dimension(:), allocatable :: latitude ! Latitude read in startfi_name and written in restartfi
158	real, dimension(:), allocatable :: cell_area ! Cell_area read in startfi_name and written in restartfi
159	real :: Total_surface ! Total surface of the planet
160
161	! Variables for h2o_ice evolution
162	real, dimension(:,:), allocatable :: h2o_ice ! h2o ice in the PEM
163	real, dimension(:,:,:), allocatable :: min_h2o_ice ! Minima of h2o ice at each point for the PCM years [kg/m^2]
164	real :: h2oice_ini_surf ! Initial surface of sublimating h2o ice
165	logical, dimension(:,:), allocatable :: ini_h2oice_sublim ! Logical array to know if h2o ice is sublimating
166	real :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
167	real :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
168	real :: global_avg_press_new ! constant: Global average pressure of current time step
169	real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
170	real, dimension(:,:), allocatable :: zplev_PCM ! same but retrieved from the PCM [kg/m^2]
171	real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
172	real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
173	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
174	real, save :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg]
175	real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
176	integer :: timelen ! # time samples
177	real :: ave ! intermediate varibale to compute average
178	real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
179	real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1
180
181	! Variables for co2_ice evolution
182	real, dimension(:,:), allocatable :: co2_ice ! co2 ice in the PEM
183	logical, dimension(:,:), allocatable :: is_co2ice_ini ! Was there co2 ice initially in the PEM?
184	real, dimension(:,:,:), allocatable :: min_co2_ice ! Minimum of co2 ice at each point for the first year [kg/m^2]
185	real :: co2ice_ini_surf ! Initial surface of sublimating co2 ice
186	logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating
187	real, dimension(:,:), allocatable :: vmr_co2_PCM ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
188	real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM
189	real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
190
191	! Variables for stratification (layering) evolution
192	type(layering), dimension(:,:), allocatable :: stratif ! Layering (linked list of "stratum") for each grid point and slope
193	type(ptrarray), dimension(:,:), allocatable :: current1, current2 ! Current active stratum in the layering
194	logical, dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm
195
196	! Variables for slopes
197	real, dimension(:,:), allocatable :: d_co2ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
198	real, dimension(:,:), allocatable :: d_co2ice_ini ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
199	real, dimension(:,:), allocatable :: d_h2oice ! physical point x slope field: Tendency of evolution of perennial h2o ice
200	real, dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow
201	real, dimension(:), allocatable :: flag_co2flow_mesh ! (ngrid) : Flag where there is a CO2 glacier flow
202	real, dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow
203	real, dimension(:), allocatable :: flag_h2oflow_mesh ! (ngrid) : Flag where there is a H2O glacier flow
204
205	! Variables for surface and soil
206	real, dimension(:,:), allocatable :: tsurf_avg ! Physic x SLOPE field: Averaged Surface Temperature [K]
207	real, dimension(:,:,:), allocatable :: tsoil_avg ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
208	real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
209	real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K]
210	real, dimension(:,:,:,:), allocatable :: tsoil_anom ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K]
211	real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
212	real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: Intermediate when computing Tsoil [K]
213	real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
214	real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
215	real, dimension(:,:), allocatable :: watersurf_density_avg ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
216	real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
217	real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
218	real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K]
219	real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
220	real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
221	real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
222	real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
223	logical :: bool_sublim ! logical to check if there is sublimation or not
224	logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
225	real, dimension(:,:), allocatable :: icetable_thickness_old ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
226	real, dimension(:,:,:), allocatable :: ice_porefilling_old ! ngrid x nslope: Ice pore filling at the previous iteration [m]
227	real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
228	real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table
229	real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table
230
231	! Some variables for the PEM run
232	real, parameter :: year_step = 1 ! Timestep for the pem
233	real :: i_myear_leg ! Number of iteration
234	real :: n_myear_leg ! Maximum number of iterations before stopping
235	real :: i_myear ! Global number of Martian years of the chained simulations
236	real :: n_myear ! Maximum number of Martian years of the chained simulations
237	real :: timestep ! Timestep [s]
238	character(20) :: job_id ! Job id provided as argument passed on the command line when the program was invoked
239	logical :: timewall ! Flag to use the time limit stopping criterion in case of a PEM job
240	integer(kind=8) :: cr ! Number of clock ticks per second (count rate)
241	integer(kind=8) :: c1, c2 ! Counts of processor clock
242	character(100) :: chtimelimit ! Time limit for the PEM job outputted by the SLURM command
243	real :: timelimit ! Time limit for the PEM job in seconds
244	real, parameter :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default
245	integer :: cstat, days, hours, minutes, seconds
246	character(1) :: sep
247
248	#ifdef CPP_STD
249	real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
250	real :: albedo_h2o_frost ! albedo of h2o frost
251	real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
252	real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
253	real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
254	real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
255	real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
256	real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case
257	real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case
258	real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case
259	real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case
260	real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
261	logical, dimension(:), allocatable :: watercaptag ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
262	real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case
263	real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case
264	#endif
265
266	#ifdef CPP_1D
267	integer :: nsplit_phys
268	integer, parameter :: jjm_value = jjm - 1
269	integer :: day_step
270
271	! Dummy variables to use the subroutine 'init_testphys1d'
272	logical :: therestart1D, therestartfi
273	integer :: ndt, day0
274	real :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau
275	real, dimension(:), allocatable :: zqsat
276	real, dimension(:,:,:), allocatable :: dq, dqdyn
277	real, dimension(nlayer) :: play, w
278	real, dimension(nlayer + 1) :: plev
279	#else
280	integer, parameter :: jjm_value = jjm
281	real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
282	real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
283	#endif
284
285	! Loop variables
286	integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
287
288	! Elapsed time with system clock
289	call system_clock(count_rate = cr)
290	call system_clock(c1)
291	timewall = .true.
292	timelimit = 86400 ! 86400 seconds = 24 h by default
293	if (command_argument_count() > 0) then
294	! Read the job id passed as argument to the program
295	call get_command_argument(1,job_id)
296	! Execute the system command
297	call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
298	if (cstat /= 0) then
299	call execute_command_line('qstat -f '//trim(job_id)//' \| grep "Walltime" \| awk ''{print $3}'' > tmp_cmdout.txt', cmdstat = cstat)
300	if (cstat > 0) then
301	error stop 'pem: command execution failed!'
302	else if (cstat < 0) then
303	error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!'
304	endif
305	endif
306	! Read the output
307	open(1,file = 'tmp_cmdout.txt',status = 'old')
308	read(1,'(a)') chtimelimit
309	close(1)
310	chtimelimit = trim(chtimelimit)
311	call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat)
312	if (cstat > 0) then
313	error stop 'pem: command execution failed!'
314	else if (cstat < 0) then
315	error stop 'pem: command execution not supported!'
316	endif
317	if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS"
318	read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds
319	timelimit = days86400 + hours3600 + minutes*60 + seconds
320	else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS"
321	read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds
322	timelimit = hours3600 + minutes60 + seconds
323	else ! 'chtimelimit' format is "MM:SS"
324	read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds
325	timelimit = minutes*60 + seconds
326	endif
327	else
328	timewall = .false.
329	endif
330
331	! Parallel variables
332	#ifndef CPP_STD
333	is_sequential = .true.
334	is_parallel = .false.
335	is_mpi_root = .true.
336	is_omp_root = .true.
337	is_master = .true.
338	#endif
339
340	! Some constants
341	day_ini = 0 ! test
342	time_phys = 0. ! test
343	ngrid = ngridmx
344	A = (1/m_co2 - 1/m_noco2)
345	B = 1/m_noco2
346	year_day = 669
347	daysec = 88775.
348	timestep = year_daydaysecyear_step
349
350	!----------------------------- INITIALIZATION --------------------------
351	!------------------------
352	! I Initialization
353	! I_a Read the "run.def"
354	!------------------------
355	#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))
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	call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, &
380	time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
381	play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
382	ps(2) = ps(1)
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_evol.nc
392	!------------------------
393	! I_b.1 Read "start.nc"
394	allocate(ps_start_PCM(ngrid))
395	#ifndef CPP_1D
396	call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
397
398	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
399
400	call iniconst !new
401	call inigeom
402
403	allocate(ap(nlayer + 1))
404	allocate(bp(nlayer + 1))
405	status = nf90_open(start_name,NF90_NOWRITE,ncid)
406	status = nf90_inq_varid(ncid,"ap",apvarid)
407	status = nf90_get_var(ncid,apvarid,ap)
408	status = nf90_inq_varid(ncid,"bp",bpvarid)
409	status = nf90_get_var(ncid,bpvarid,bp)
410	status = nf90_close(ncid)
411
412	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)
413	#else
414	ps_start_PCM(1) = ps(1)
415	#endif
416
417	! Save initial surface pressure
418	ps0 = ps
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
424	status = nf90_inq_varid(ncid,"longitude",lonvarid)
425	status = nf90_get_var(ncid,lonvarid,longitude)
426
427	status = nf90_inq_varid(ncid,"latitude",latvarid)
428	status = nf90_get_var(ncid,latvarid,latitude)
429
430	status = nf90_inq_varid(ncid,"area",areavarid)
431	status = nf90_get_var(ncid,areavarid,cell_area)
432
433	status = nf90_inq_varid(ncid,"soildepth",sdvarid)
434	status = nf90_get_var(ncid,sdvarid,mlayer)
435
436	status = nf90_close(ncid)
437
438	! I_b.2 Read the "startfi.nc"
439	! First we read the initial state (starfi.nc)
440	#ifndef CPP_STD
441	call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
442	tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, &
443	watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
444
445	! Remove unphysical values of surface tracer
446	where (qsurf < 0.) qsurf = 0.
447
448	call surfini(ngrid,nslope,qsurf)
449	#else
450	call phys_state_var_init(nq)
451	if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
452	call initracer(ngrid,nq)
453	call iniaerosol()
454	allocate(tsurf_read_generic(ngrid))
455	allocate(qsurf_read_generic(ngrid,nq))
456	allocate(tsoil_read_generic(ngrid,nsoilmx))
457	allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
458	allocate(emis_read_generic(ngrid))
459	allocate(tsurf(ngrid,1))
460	allocate(qsurf(ngrid,nq,1))
461	allocate(tsoil(ngrid,nsoilmx,1))
462	allocate(emis(ngrid,1))
463	allocate(watercap(ngrid,1))
464	allocate(watercaptag(ngrid))
465	allocate(albedo_read_generic(ngrid,2))
466	allocate(albedo(ngrid,2,1))
467	allocate(inertiesoil(ngrid,nsoilmx,1))
468	call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
469	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
470	qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, &
471	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
472	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
473
474	nslope = 1
475	call ini_comslope_h(ngrid,1)
476
477	qsurf(:,:,1) = qsurf_read_generic
478	tsurf(:,1) = tsurf_read_generic
479	tsoil(:,:,1) = tsoil_read_generic
480	emis(:,1) = emis_read_generic
481	watercap(:,1) = 0.
482	watercaptag(:) = .false.
483	albedo(:,1,1) = albedo_read_generic(:,1)
484	albedo(:,2,1) = albedo_read_generic(:,2)
485	inertiesoil(:,:,1) = inertiedat
486
487	if (nslope == 1) then
488	def_slope(1) = 0
489	def_slope(2) = 0
490	def_slope_mean = 0
491	subslope_dist(:,1) = 1.
492	endif
493
494	! Remove unphysical values of surface tracer
495	qsurf(:,:,1) = qsurf_read_generic
496	where (qsurf < 0.) qsurf = 0.
497	#endif
498
499	do nnq = 1,nqtot ! Why not using ini_tracer ?
500	if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
501	if (noms(nnq) == "h2o_vap") then
502	igcm_h2o_vap = nnq
503	mmol(igcm_h2o_vap) = 18.
504	endif
505	if (noms(nnq) == "co2") igcm_co2 = nnq
506	enddo
507	r = 8.314511*1000./mugaz
508
509	!------------------------
510	! I Initialization
511	! I_c Subslope parametrisation
512	!------------------------
513	! Define some slope statistics
514	iflat = 1
515	do islope = 2,nslope
516	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
517	enddo
518
519	write(,) 'Flat slope for islope = ',iflat
520	write(,) 'corresponding criterium = ',def_slope_mean(iflat)
521
522	allocate(flag_co2flow(ngrid,nslope))
523	allocate(flag_co2flow_mesh(ngrid))
524	allocate(flag_h2oflow(ngrid,nslope))
525	allocate(flag_h2oflow_mesh(ngrid))
526
527	flag_co2flow = 0
528	flag_co2flow_mesh = 0
529	flag_h2oflow = 0
530	flag_h2oflow_mesh = 0
531
532	!------------------------
533	! I Initialization
534	! I_d Read the PCM data and convert them to the physical grid
535	!------------------------
536	! 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
537	call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
538
539	allocate(tsoil_avg(ngrid,nsoilmx,nslope))
540	allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope))
541	allocate(vmr_co2_PCM(ngrid,timelen))
542	allocate(ps_timeseries(ngrid,timelen))
543	allocate(min_co2_ice(ngrid,nslope,2))
544	allocate(min_h2o_ice(ngrid,nslope,2))
545	allocate(tsurf_avg_yr1(ngrid,nslope))
546	allocate(tsurf_avg(ngrid,nslope))
547	allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
548	allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen))
549	allocate(q_co2_PEM_phys(ngrid,timelen))
550	allocate(q_h2o_PEM_phys(ngrid,timelen))
551	allocate(watersurf_density_avg(ngrid,nslope))
552	allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
553	allocate(Tsurfavg_before_saved(ngrid,nslope))
554	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
555	allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
556	allocate(delta_co2_adsorbded(ngrid))
557	allocate(co2ice_disappeared(ngrid,nslope))
558	allocate(icetable_thickness_old(ngrid,nslope))
559	allocate(ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
560	allocate(delta_h2o_icetablesublim(ngrid))
561	allocate(delta_h2o_adsorbded(ngrid))
562	allocate(vmr_co2_PEM_phys(ngrid,timelen))
563
564	write(,) "Downloading data Y1..."
565	call read_data_PCM("data_PCM_Y1.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), &
566	tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, &
567	watersurf_density_avg,watersoil_density_timeseries)
568	write(,) "Downloading data Y1 done!"
569
570	! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value
571	write(,) "Downloading data Y2..."
572	call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,2),min_h2o_ice(:,:,2), &
573	tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, &
574	watersurf_density_avg,watersoil_density_timeseries)
575	write(,) "Downloading data Y2 done!"
576
577	!------------------------
578	! I Initialization
579	! I_e Initialization of the PEM variables and soil
580	!------------------------
581	call end_comsoil_h_PEM
582	call ini_comsoil_h_PEM(ngrid,nslope)
583	call end_adsorption_h_PEM
584	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
585	call end_ice_table
586	call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
587
588	if (soil_pem) then
589	do t = 1,timelen
590	tsoil_anom(:,:,:,t) = tsoil_anom(:,:,:,t) - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - <Tsoil> to recompute properly tsoil in the restart
591	enddo
592	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
593	tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
594	tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_anom
595	watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
596	do l = nsoilmx + 1,nsoilmx_PEM
597	tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
598	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
599	enddo
600	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
601	endif !soil_pem
602	deallocate(tsoil_avg)
603
604	!------------------------
605	! I Initialization
606	! I_f Compute tendencies
607	!------------------------
608	allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope))
609	allocate(d_co2ice_ini(ngrid,nslope))
610
611	! Compute the tendencies of the evolution of ice over the years
612	call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
613	call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
614	d_co2ice_ini = d_co2ice
615	deallocate(min_co2_ice,min_h2o_ice)
616
617	!------------------------
618	! I Initialization
619	! I_g Save the initial situation
620	!------------------------
621	allocate(zplev_PCM(ngrid,nlayer + 1))
622	Total_surface = 0.
623	do ig = 1,ngrid
624	Total_surface = Total_surface + cell_area(ig)
625	zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
626	enddo
627	global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
628	global_avg_press_PCM = global_avg_press_old
629	global_avg_press_new = global_avg_press_old
630	write(,) "Total surface of the planet =", Total_surface
631	write(,) "Initial global average pressure =", global_avg_press_PCM
632
633	!------------------------
634	! I Initialization
635	! I_h Read the "startpem.nc"
636	!------------------------
637	allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
638
639	allocate(stratif(ngrid,nslope))
640	if (layering_algo) then
641	do islope = 1,nslope
642	do i = 1,ngrid
643	call ini_layering(stratif(i,islope))
644	enddo
645	enddo
646	endif
647
648	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
649	icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys, &
650	ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice, &
651	global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys, &
652	delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
653
654	! We save the places where h2o ice is sublimating
655	! We compute the surface of h2o ice sublimating
656	allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
657	co2ice_ini_surf = 0.
658	h2oice_ini_surf = 0.
659	ini_co2ice_sublim = .false.
660	ini_h2oice_sublim = .false.
661	is_co2ice_ini = .false.
662	co2ice_disappeared = .false.
663	do i = 1,ngrid
664	do islope = 1,nslope
665	if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
666	if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
667	ini_co2ice_sublim(i,islope) = .true.
668	co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
669	endif
670	if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
671	ini_h2oice_sublim(i,islope) = .true.
672	h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
673	endif
674	enddo
675	enddo
676	write(,) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf
677	write(,) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
678
679	delta_h2o_icetablesublim = 0.
680
681	if (adsorption_pem) then
682	totmassco2_adsorbded = 0.
683	totmassh2o_adsorbded = 0.
684	do ig = 1,ngrid
685	do islope = 1,nslope
686	do l = 1,nsoilmx_PEM - 1
687	if (l==1) then
688	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
689	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
690	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
691	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
692	else
693	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
694	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
695	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
696	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
697	endif
698	enddo
699	enddo
700	enddo
701	write(,) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
702	write(,) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
703	endif ! adsorption
704	deallocate(tsurf_avg_yr1)
705
706	!------------------------
707	! I Initialization
708	! I_i Compute orbit criterion
709	!------------------------
710	#ifndef CPP_STD
711	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
712	#else
713	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
714	#endif
715
716	n_myear_leg = Max_iter_pem
717	if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
718
719	!-------------------------- END INITIALIZATION -------------------------
720
721	!-------------------------------- RUN ----------------------------------
722	!------------------------
723	! II Run
724	! II_a Update pressure, ice and tracers
725	!------------------------
726	i_myear_leg = 0
727	stopPEM = 0
728	if (layering_algo) then
729	allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
730	new_str = .true.
731	new_lag1 = .true.
732	new_lag2 = .true.
733	do islope = 1,nslope
734	do ig = 1,ngrid
735	current1(ig,islope)%p => stratif(ig,islope)%top
736	current2(ig,islope)%p => stratif(ig,islope)%top
737	enddo
738	enddo
739	endif
740
741	do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
742	! II.a.1. Compute updated global pressure
743	write(,) "Recomputing the new pressure..."
744	do i = 1,ngrid
745	do islope = 1,nslope
746	global_avg_press_new = global_avg_press_new - CO2cond_psgcell_area(i)d_co2ice(i,islope)subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
747	enddo
748	enddo
749
750	if (adsorption_pem) then
751	do i = 1,ngrid
752	global_avg_press_new = global_avg_press_new - gcell_area(i)delta_co2_adsorbded(i)/Total_surface
753	enddo
754	endif
755	write(,) 'Global average pressure old time step',global_avg_press_old
756	write(,) 'Global average pressure new time step',global_avg_press_new
757
758	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
759	allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
760	write(,) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
761	do l = 1,nlayer + 1
762	do ig = 1,ngrid
763	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
764	enddo
765	enddo
766
767	! II.a.3. Surface pressure timeseries
768	write(,) "Recomputing the surface pressure timeserie adapted to the new pressure..."
769	do ig = 1,ngrid
770	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
771	enddo
772
773	! II.a.4. New pressure levels timeseries
774	allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
775	write(,) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
776	do l = 1,nlayer + 1
777	do ig = 1,ngrid
778	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
779	enddo
780	enddo
781
782	! II.a.5. Tracers timeseries
783	write(,) "Recomputing of tracer VMR timeseries for the new pressure..."
784
785	l = 1
786	do ig = 1,ngrid
787	do t = 1,timelen
788	q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
789	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
790	if (q_h2o_PEM_phys(ig,t) < 0) then
791	q_h2o_PEM_phys(ig,t) = 1.e-30
792	else if (q_h2o_PEM_phys(ig,t) > 1) then
793	q_h2o_PEM_phys(ig,t) = 1.
794	endif
795	enddo
796	enddo
797
798	do ig = 1,ngrid
799	do t = 1,timelen
800	q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
801	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
802	+ ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
803	- (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
804	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
805	if (q_co2_PEM_phys(ig,t) < 0) then
806	q_co2_PEM_phys(ig,t) = 1.e-30
807	else if (q_co2_PEM_phys(ig,t) > 1) then
808	q_co2_PEM_phys(ig,t) = 1.
809	endif
810	mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
811	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
812	enddo
813	enddo
814
815	deallocate(zplev_new_timeseries,zplev_old_timeseries)
816
817	!------------------------
818	! II Run
819	! II_b Evolution of ice
820	!------------------------
821	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
822	call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
823	if (layering_algo) then
824	do islope = 1,nslope
825	do ig = 1,ngrid
826	call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
827	enddo
828	enddo
829	endif
830
831	!------------------------
832	! II Run
833	! II_c Flow of glaciers
834	!------------------------
835	if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
836	global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
837	if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
838
839	!------------------------
840	! II Run
841	! II_d Update surface and soil temperatures
842	!------------------------
843	! II_d.1 Update Tsurf
844	write(,) "Updating the new Tsurf"
845	bool_sublim = .false.
846	Tsurfavg_before_saved = tsurf_avg
847	do ig = 1,ngrid
848	do islope = 1,nslope
849	if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then ! co2 ice disappeared, look for closest point without co2ice
850	co2ice_disappeared(ig,islope) = .true.
851	if (latitude_deg(ig) > 0) then
852	do ig_loop = ig,ngrid
853	do islope_loop = islope,iflat,-1
854	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
855	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
856	bool_sublim = .true.
857	exit
858	endif
859	enddo
860	if (bool_sublim) exit
861	enddo
862	else
863	do ig_loop = ig,1,-1
864	do islope_loop = islope,iflat
865	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
866	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
867	bool_sublim = .true.
868	exit
869	endif
870	enddo
871	if (bool_sublim) exit
872	enddo
873	endif
874	if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
875	albedo(ig,1,islope) = albedo_h2o_frost
876	albedo(ig,2,islope) = albedo_h2o_frost
877	else
878	albedo(ig,1,islope) = albedodat(ig)
879	albedo(ig,2,islope) = albedodat(ig)
880	emis(ig,islope) = emissiv
881	endif
882	else if ((co2_ice(ig,islope) > 1.e-3) .and. (d_co2ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
883	ave = 0.
884	do t = 1,timelen
885	ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
886	enddo
887	tsurf_avg(ig,islope) = ave/timelen
888	endif
889	enddo
890	enddo
891
892	do t = 1,timelen
893	tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_avg - Tsurfavg_before_saved
894	enddo
895	! for the start
896	do ig = 1,ngrid
897	do islope = 1,nslope
898	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_avg(ig,islope))
899	enddo
900	enddo
901
902	if (soil_pem) then
903
904	! II_d.2 Update soil temperature
905	write(,)"Updating soil temperature"
906	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
907	do islope = 1,nslope
908	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
909	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
910
911	do t = 1,timelen
912	Tsoil_locslope(:,1:nsoilmx) = tsoil_PEM(:,1:nsoilmx,islope) + tsoil_anom(:,:,islope,t)
913	Tsoil_locslope(:,nsoilmx + 1:) = tsoil_PEM(:,nsoilmx + 1:,islope)
914	do ig = 1,ngrid
915	do isoil = 1,nsoilmx_PEM
916	watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)mmol(igcm_h2o_vap)/(mugazr)
917	if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
918	enddo
919	enddo
920	enddo
921	enddo
922	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
923
924	write(,) "Update of soil temperature done"
925
926	deallocate(Tsoil_locslope)
927
928	! II_d.3 Update the ice table
929	if (icetable_equilibrium) then
930	write(,) "Compute ice table (equilibrium method)"
931	icetable_thickness_old = icetable_thickness
932	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
933	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
934	else if (icetable_dynamic) then
935	write(,) "Compute ice table (dynamic method)"
936	ice_porefilling_old = ice_porefilling
937	allocate(porefill(nsoilmx_PEM))
938	do ig = 1,ngrid
939	do islope = 1,nslope
940	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(ig),sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2),ice_porefilling(ig,:,islope),porefill,ssi_depth)
941	icetable_depth(ig,islope) = ssi_depth
942	ice_porefilling(ig,:,islope) = porefill
943	enddo
944	enddo
945	deallocate(porefill)
946	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
947	endif
948
949	! II_d.4 Update the soil thermal properties
950	call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
951
952	! II_d.5 Update the mass of the regolith adsorbed
953	if (adsorption_pem) then
954	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice, &
955	h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
956	h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
957
958	totmassco2_adsorbded = 0.
959	totmassh2o_adsorbded = 0.
960	do ig = 1,ngrid
961	do islope = 1,nslope
962	do l = 1,nsoilmx_PEM
963	if (l == 1) then
964	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
965	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
966	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
967	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
968	else
969	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
970	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
971	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
972	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
973	endif
974	enddo
975	enddo
976	enddo
977	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
978	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
979	endif
980	endif !soil_pem
981
982	!------------------------
983	! II Run
984	! II_e Outputs
985	!------------------------
986	call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/global_avg_press_new/))
987	do islope = 1,nslope
988	write(str2(1:2),'(i2.2)') islope
989	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
990	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
991	call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
992	call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
993	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
994	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
995	if (icetable_equilibrium) then
996	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
997	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
998	else if (icetable_dynamic) then
999	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
1000	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
1001	endif
1002
1003	if (soil_pem) then
1004	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
1005	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
1006	if (adsorption_pem) then
1007	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
1008	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
1009	endif
1010	endif
1011	enddo
1012
1013	!------------------------
1014	! II Run
1015	! II_f Update the tendencies
1016	!------------------------
1017	call recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,global_avg_press_PCM,global_avg_press_new)
1018
1019	!------------------------
1020	! II Run
1021	! II_g Checking the stopping criterion
1022	!------------------------
1023
1024	write(,) "Checking the stopping criteria..."
1025	call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
1026	call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
1027	i_myear_leg = i_myear_leg + dt
1028	i_myear = i_myear + dt
1029	if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
1030	if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
1031	call system_clock(c2)
1032	if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
1033	if (stopPEM > 0) then
1034	select case (stopPEM)
1035	case(1)
1036	write(,) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
1037	case(2)
1038	write(,) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
1039	case(3)
1040	write(,) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
1041	case(4)
1042	write(,) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
1043	case(5)
1044	write(,) "STOPPING because maximum number of iterations is reached (possibly due to orbital parameters):", stopPEM
1045	case(6)
1046	write(,) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
1047	case(7)
1048	write(,) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM
1049	case(8)
1050	write(,) "STOPPING because the layering algorithm met an hasty end:", stopPEM
1051	case default
1052	write(,) "STOPPING with unknown stopping criterion code:", stopPEM
1053	end select
1054	exit
1055	else
1056	write(,) "The PEM can continue!"
1057	write(,) "Number of iterations of the PEM: i_myear_leg =", i_myear_leg
1058	write(,) "Number of simulated Martian years: i_myear =", i_myear
1059	endif
1060
1061	global_avg_press_old = global_avg_press_new
1062
1063	enddo ! big time iteration loop of the pem
1064	!------------------------------ END RUN --------------------------------
1065
1066	!------------------------------- OUTPUT --------------------------------
1067	!------------------------
1068	! III Output
1069	! III_a Update surface value for the PCM start files
1070	!------------------------
1071	! III_a.1 Ice update (for startfi)
1072
1073	watercap = 0.
1074	perennial_co2ice = co2_ice
1075	do ig = 1,ngrid
1076	! H2O ice metamorphism
1077	if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_h2oice_threshold) then
1078	h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1079	qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
1080	endif
1081
1082	! Is H2O ice still considered as an infinite reservoir for the PCM?
1083	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
1084	! There is enough ice to be considered as an infinite reservoir
1085	watercaptag(ig) = .true.
1086	else
1087	! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
1088	watercaptag(ig) = .false.
1089	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1090	h2o_ice(ig,:) = 0.
1091	endif
1092
1093	! CO2 ice metamorphism
1094	if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_co2ice_threshold) then
1095	perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1096	qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
1097	endif
1098	enddo
1099
1100	! III_a.2 Tsoil update (for startfi)
1101	if (soil_pem) then
1102	call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
1103	tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen)
1104	#ifndef CPP_STD
1105	flux_geo = fluxgeo
1106	#endif
1107	endif
1108	deallocate(tsoil_anom)
1109
1110	! III_a.4 Pressure (for start)
1111	ps = ps*global_avg_press_new/global_avg_press_PCM
1112	ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
1113
1114	! III_a.5 Tracer (for start)
1115	allocate(zplev_new(ngrid,nlayer + 1))
1116
1117	do l = 1,nlayer + 1
1118	zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
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_PCM(ig,l) - zplev_PCM(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_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
1133	+ ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_PCM(ig,l) - zplev_PCM(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
1140	! Conserving the tracers mass for PCM start files
1141	do nnq = 1,nqtot
1142	do ig = 1,ngrid
1143	do l = 1,llm - 1
1144	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
1145	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1146	q(ig,l,nnq) = 1.
1147	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1148	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1149	endif
1150	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1151	enddo
1152	enddo
1153	enddo
1154
1155	if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
1156
1157	!------------------------
1158	! III Output
1159	! III_b Write "restart.nc" and "restartfi.nc"
1160	!------------------------
1161	! III_b.1 Write "restart.nc"
1162	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1163	pday = day_ini
1164	ztime_fin = time_phys
1165
1166	allocate(p(ip1jmp1,nlayer + 1))
1167	#ifndef CPP_1D
1168	! Correction on teta due to surface pressure changes
1169	do l = 1,nlayer
1170	do i = 1,ip1jmp1
1171	teta(i,l)= teta(i,l)(ps0(i)/ps(i))*kappa
1172	enddo
1173	enddo
1174	! Correction on atmospheric pressure
1175	call pression(ip1jmp1,ap,bp,ps,p)
1176	! Correction on the mass of atmosphere
1177	call massdair(p,masse)
1178	call dynredem0("restart.nc",day_ini,phis)
1179	call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps)
1180	write(,) "restart.nc has been written"
1181	#else
1182	call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1183	write(,) "restart1D.txt has been written"
1184	#endif
1185
1186	! III_b.2 Write the "restartfi.nc"
1187	#ifndef CPP_STD
1188	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1189	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1190	inertiedat,def_slope,subslope_dist)
1191	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1192	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1193	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1194	wstar,watercap,perennial_co2ice)
1195	#else
1196	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1197	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
1198	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1199	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1200	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1201	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
1202	tsea_ice,sea_ice)
1203	#endif
1204	write(,) "restartfi.nc has been written"
1205
1206	!------------------------
1207	! III Output
1208	! III_c Write the "restartpem.nc"
1209	!------------------------
1210	if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
1211	call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1212	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
1213	TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
1214	co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
1215	write(,) "restartpem.nc has been written"
1216
1217	call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
1218
1219	write(,) "The PEM has run for", i_myear_leg, "Martian years."
1220	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1221	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1222	write(,) "PEM: so far, so good!"
1223
1224	if (layering_algo) then
1225	do islope = 1,nslope
1226	do i = 1,ngrid
1227	call del_layering(stratif(i,islope))
1228	enddo
1229	enddo
1230	deallocate(new_str,new_lag1,new_lag2,current1,current2)
1231	endif
1232	deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
1233	deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
1234	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
1235	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
1236	deallocate(icetable_thickness_old,ice_porefilling_old)
1237	deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
1238	!----------------------------- END OUTPUT ------------------------------
1239
1240	END PROGRAM pem

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

Download in other formats:

Original Format