Context Navigation

pem.F90 @ 3516

Last change on this file since 3516 was 3516, checked in by jbclement, 7 days ago
PEM: Small corrections related to r3498 (time step from integer to real) for the launching script + making a proper initilization of 'year_bp_ini' in case of there is no orbital evolution. JBC
File size: 61.5 KB

Line
1	!------------------------
2	! I Initialization
3	! I_a Read the "run.def"
4	! I_b Read the "start.nc" and "startfi.nc"
5	! I_c Subslope parametrisation
6	! I_d Read the PCM data and convert them to the physical grid
7	! I_e Initialization of the PEM variable and soil
8	! I_f Compute tendencies
9	! I_g 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_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
226	real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
227	real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table
228	real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table
229
230	! Some variables for the PEM run
231	real, parameter :: year_step = 1 ! Timestep for the pem
232	real :: i_myear_leg ! Number of iteration
233	real :: n_myear_leg ! Maximum number of iterations before stopping
234	real :: i_myear ! Global number of Martian years of the chained simulations
235	real :: n_myear ! Maximum number of Martian years of the chained simulations
236	real :: timestep ! Timestep [s]
237	character(20) :: job_id ! Job id provided as argument passed on the command line when the program was invoked
238	logical :: timewall ! Flag to use the time limit stopping criterion in case of a PEM job
239	integer(kind=8) :: cr ! Number of clock ticks per second (count rate)
240	integer(kind=8) :: c1, c2 ! Counts of processor clock
241	character(100) :: chtimelimit ! Time limit for the PEM job outputted by the SLURM command
242	real :: timelimit ! Time limit for the PEM job in seconds
243	real, parameter :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default
244	integer :: cstat, days, hours, minutes, seconds
245	character(1) :: sep
246
247	#ifdef CPP_STD
248	real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
249	real :: albedo_h2o_frost ! albedo of h2o frost
250	real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
251	real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
252	real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
253	real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
254	real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
255	real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case
256	real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case
257	real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case
258	real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case
259	real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
260	logical, dimension(:), allocatable :: watercaptag ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
261	real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case
262	real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case
263	#endif
264
265	#ifdef CPP_1D
266	integer :: nsplit_phys
267	integer, parameter :: jjm_value = jjm - 1
268	integer :: day_step
269
270	! Dummy variables to use the subroutine 'init_testphys1d'
271	logical :: therestart1D, therestartfi
272	integer :: ndt, day0
273	real :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau
274	real, dimension(:), allocatable :: zqsat
275	real, dimension(:,:,:), allocatable :: dq, dqdyn
276	real, dimension(nlayer) :: play, w
277	real, dimension(nlayer + 1) :: plev
278	#else
279	integer, parameter :: jjm_value = jjm
280	real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
281	real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
282	#endif
283
284	! Loop variables
285	integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
286
287	! Elapsed time with system clock
288	call system_clock(count_rate = cr)
289	call system_clock(c1)
290	timewall = .true.
291	timelimit = 86400 ! 86400 seconds = 24 h by default
292	if (command_argument_count() > 0) then
293	! Read the job id passed as argument to the program
294	call get_command_argument(1,job_id)
295	! Execute the system command
296	call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
297	if (cstat /= 0) then
298	call execute_command_line('qstat -f '//trim(job_id)//' \| grep "Walltime" \| awk ''{print $3}'' > tmp_cmdout.txt', cmdstat = cstat)
299	if (cstat > 0) then
300	error stop 'pem: command execution failed!'
301	else if (cstat < 0) then
302	error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!'
303	endif
304	endif
305	! Read the output
306	open(1,file = 'tmp_cmdout.txt',status = 'old')
307	read(1,'(a)') chtimelimit
308	close(1)
309	chtimelimit = trim(chtimelimit)
310	call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat)
311	if (cstat > 0) then
312	error stop 'pem: command execution failed!'
313	else if (cstat < 0) then
314	error stop 'pem: command execution not supported!'
315	endif
316	if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS"
317	read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds
318	timelimit = days86400 + hours3600 + minutes*60 + seconds
319	else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS"
320	read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds
321	timelimit = hours3600 + minutes60 + seconds
322	else ! 'chtimelimit' format is "MM:SS"
323	read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds
324	timelimit = minutes*60 + seconds
325	endif
326	else
327	timewall = .false.
328	endif
329
330	! Parallel variables
331	#ifndef CPP_STD
332	is_sequential = .true.
333	is_parallel = .false.
334	is_mpi_root = .true.
335	is_omp_root = .true.
336	is_master = .true.
337	#endif
338
339	! Some constants
340	day_ini = 0 ! test
341	time_phys = 0. ! test
342	ngrid = ngridmx
343	A = (1/m_co2 - 1/m_noco2)
344	B = 1/m_noco2
345	year_day = 669
346	daysec = 88775.
347	timestep = year_daydaysecyear_step
348
349	!----------------------------- INITIALIZATION --------------------------
350	!------------------------
351	! I Initialization
352	! I_a Read the "run.def"
353	!------------------------
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))
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	call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, &
379	time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
380	play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
381	ps(2) = ps(1)
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_evol.nc
391	!------------------------
392	! I_b.1 Read "start.nc"
393	allocate(ps_start_PCM(ngrid))
394	#ifndef CPP_1D
395	call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
396
397	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
398
399	call iniconst !new
400	call inigeom
401
402	allocate(ap(nlayer + 1))
403	allocate(bp(nlayer + 1))
404	status = nf90_open(start_name,NF90_NOWRITE,ncid)
405	status = nf90_inq_varid(ncid,"ap",apvarid)
406	status = nf90_get_var(ncid,apvarid,ap)
407	status = nf90_inq_varid(ncid,"bp",bpvarid)
408	status = nf90_get_var(ncid,bpvarid,bp)
409	status = nf90_close(ncid)
410
411	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)
412	#else
413	ps_start_PCM(1) = ps(1)
414	#endif
415
416	! Save initial surface pressure
417	ps0 = ps
418
419	! In the PCM, these values are given to the physic by the dynamic.
420	! Here we simply read them in the "startfi.nc" file
421	status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
422
423	status = nf90_inq_varid(ncid,"longitude",lonvarid)
424	status = nf90_get_var(ncid,lonvarid,longitude)
425
426	status = nf90_inq_varid(ncid,"latitude",latvarid)
427	status = nf90_get_var(ncid,latvarid,latitude)
428
429	status = nf90_inq_varid(ncid,"area",areavarid)
430	status = nf90_get_var(ncid,areavarid,cell_area)
431
432	status = nf90_inq_varid(ncid,"soildepth",sdvarid)
433	status = nf90_get_var(ncid,sdvarid,mlayer)
434
435	status = nf90_close(ncid)
436
437	! I_b.2 Read the "startfi.nc"
438	! First we read the initial state (starfi.nc)
439	#ifndef CPP_STD
440	call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
441	tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, &
442	watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
443
444	! Remove unphysical values of surface tracer
445	where (qsurf < 0.) qsurf = 0.
446
447	call surfini(ngrid,nslope,qsurf)
448	#else
449	call phys_state_var_init(nq)
450	if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
451	call initracer(ngrid,nq)
452	call iniaerosol()
453	allocate(tsurf_read_generic(ngrid))
454	allocate(qsurf_read_generic(ngrid,nq))
455	allocate(tsoil_read_generic(ngrid,nsoilmx))
456	allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
457	allocate(emis_read_generic(ngrid))
458	allocate(tsurf(ngrid,1))
459	allocate(qsurf(ngrid,nq,1))
460	allocate(tsoil(ngrid,nsoilmx,1))
461	allocate(emis(ngrid,1))
462	allocate(watercap(ngrid,1))
463	allocate(watercaptag(ngrid))
464	allocate(albedo_read_generic(ngrid,2))
465	allocate(albedo(ngrid,2,1))
466	allocate(inertiesoil(ngrid,nsoilmx,1))
467	call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
468	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
469	qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, &
470	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
471	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
472
473	nslope = 1
474	call ini_comslope_h(ngrid,1)
475
476	qsurf(:,:,1) = qsurf_read_generic
477	tsurf(:,1) = tsurf_read_generic
478	tsoil(:,:,1) = tsoil_read_generic
479	emis(:,1) = emis_read_generic
480	watercap(:,1) = 0.
481	watercaptag(:) = .false.
482	albedo(:,1,1) = albedo_read_generic(:,1)
483	albedo(:,2,1) = albedo_read_generic(:,2)
484	inertiesoil(:,:,1) = inertiedat
485
486	if (nslope == 1) then
487	def_slope(1) = 0
488	def_slope(2) = 0
489	def_slope_mean = 0
490	subslope_dist(:,1) = 1.
491	endif
492
493	! Remove unphysical values of surface tracer
494	qsurf(:,:,1) = qsurf_read_generic
495	where (qsurf < 0.) qsurf = 0.
496	#endif
497
498	do nnq = 1,nqtot ! Why not using ini_tracer ?
499	if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
500	if (noms(nnq) == "h2o_vap") then
501	igcm_h2o_vap = nnq
502	mmol(igcm_h2o_vap) = 18.
503	endif
504	if (noms(nnq) == "co2") igcm_co2 = nnq
505	enddo
506	r = 8.314511*1000./mugaz
507
508	!------------------------
509	! I Initialization
510	! I_c Subslope parametrisation
511	!------------------------
512	! Define some slope statistics
513	iflat = 1
514	do islope = 2,nslope
515	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
516	enddo
517
518	write(,) 'Flat slope for islope = ',iflat
519	write(,) 'corresponding criterium = ',def_slope_mean(iflat)
520
521	allocate(flag_co2flow(ngrid,nslope))
522	allocate(flag_co2flow_mesh(ngrid))
523	allocate(flag_h2oflow(ngrid,nslope))
524	allocate(flag_h2oflow_mesh(ngrid))
525
526	flag_co2flow = 0
527	flag_co2flow_mesh = 0
528	flag_h2oflow = 0
529	flag_h2oflow_mesh = 0
530
531	!------------------------
532	! I Initialization
533	! I_d Read the PCM data and convert them to the physical grid
534	!------------------------
535	! 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
536	call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
537
538	allocate(tsoil_avg(ngrid,nsoilmx,nslope))
539	allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope))
540	allocate(vmr_co2_PCM(ngrid,timelen))
541	allocate(ps_timeseries(ngrid,timelen))
542	allocate(min_co2_ice(ngrid,nslope,2))
543	allocate(min_h2o_ice(ngrid,nslope,2))
544	allocate(tsurf_avg_yr1(ngrid,nslope))
545	allocate(tsurf_avg(ngrid,nslope))
546	allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
547	allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen))
548	allocate(q_co2_PEM_phys(ngrid,timelen))
549	allocate(q_h2o_PEM_phys(ngrid,timelen))
550	allocate(watersurf_density_avg(ngrid,nslope))
551	allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
552	allocate(Tsurfavg_before_saved(ngrid,nslope))
553	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
554	allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
555	allocate(delta_co2_adsorbded(ngrid))
556	allocate(co2ice_disappeared(ngrid,nslope))
557	allocate(icetable_thickness_prev_iter(ngrid,nslope))
558	allocate(delta_h2o_icetablesublim(ngrid))
559	allocate(delta_h2o_adsorbded(ngrid))
560	allocate(vmr_co2_PEM_phys(ngrid,timelen))
561
562	write(,) "Downloading data Y1..."
563	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), &
564	tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, &
565	watersurf_density_avg,watersoil_density_timeseries)
566	write(,) "Downloading data Y1 done!"
567
568	! 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
569	write(,) "Downloading data Y2..."
570	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), &
571	tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, &
572	watersurf_density_avg,watersoil_density_timeseries)
573	write(,) "Downloading data Y2 done!"
574
575	!------------------------
576	! I Initialization
577	! I_e Initialization of the PEM variables and soil
578	!------------------------
579	call end_comsoil_h_PEM
580	call ini_comsoil_h_PEM(ngrid,nslope)
581	call end_adsorption_h_PEM
582	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
583	call end_ice_table
584	call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
585
586	if (soil_pem) then
587	do t = 1,timelen
588	tsoil_anom(:,:,:,t) = tsoil_anom(:,:,:,t) - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - <Tsoil> to recompute properly tsoil in the restart
589	enddo
590	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
591	tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
592	tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_anom
593	watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
594	do l = nsoilmx + 1,nsoilmx_PEM
595	tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
596	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
597	enddo
598	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
599	endif !soil_pem
600	deallocate(tsoil_avg)
601
602	!------------------------
603	! I Initialization
604	! I_f Compute tendencies
605	!------------------------
606	allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope))
607	allocate(d_co2ice_ini(ngrid,nslope))
608
609	! Compute the tendencies of the evolution of ice over the years
610	call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
611	call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
612	d_co2ice_ini = d_co2ice
613	deallocate(min_co2_ice,min_h2o_ice)
614
615	!------------------------
616	! I Initialization
617	! I_g Save the initial situation
618	!------------------------
619	allocate(zplev_PCM(ngrid,nlayer + 1))
620	Total_surface = 0.
621	do ig = 1,ngrid
622	Total_surface = Total_surface + cell_area(ig)
623	zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
624	enddo
625	global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
626	global_avg_press_PCM = global_avg_press_old
627	global_avg_press_new = global_avg_press_old
628	write(,) "Total surface of the planet =", Total_surface
629	write(,) "Initial global average pressure =", global_avg_press_PCM
630
631	!------------------------
632	! I Initialization
633	! I_h Read the "startpem.nc"
634	!------------------------
635	allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
636
637	allocate(stratif(ngrid,nslope))
638	if (layering_algo) then
639	do islope = 1,nslope
640	do i = 1,ngrid
641	call ini_layering(stratif(i,islope))
642	enddo
643	enddo
644	endif
645
646	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
647	icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys, &
648	ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice, &
649	global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys, &
650	delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
651
652	! We save the places where h2o ice is sublimating
653	! We compute the surface of h2o ice sublimating
654	allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
655	co2ice_ini_surf = 0.
656	h2oice_ini_surf = 0.
657	ini_co2ice_sublim = .false.
658	ini_h2oice_sublim = .false.
659	is_co2ice_ini = .false.
660	co2ice_disappeared = .false.
661	do i = 1,ngrid
662	do islope = 1,nslope
663	if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
664	if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
665	ini_co2ice_sublim(i,islope) = .true.
666	co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
667	endif
668	if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
669	ini_h2oice_sublim(i,islope) = .true.
670	h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
671	endif
672	enddo
673	enddo
674	write(,) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf
675	write(,) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
676
677	delta_h2o_icetablesublim = 0.
678
679	if (adsorption_pem) then
680	totmassco2_adsorbded = 0.
681	totmassh2o_adsorbded = 0.
682	do ig = 1,ngrid
683	do islope = 1,nslope
684	do l = 1,nsoilmx_PEM - 1
685	if (l==1) then
686	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
687	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
688	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
689	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
690	else
691	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
692	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
693	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_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	endif
696	enddo
697	enddo
698	enddo
699	write(,) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
700	write(,) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
701	endif ! adsorption
702	deallocate(tsurf_avg_yr1)
703
704	!------------------------
705	! I Initialization
706	! I_i Compute orbit criterion
707	!------------------------
708	#ifndef CPP_STD
709	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
710	#else
711	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
712	#endif
713
714	n_myear_leg = Max_iter_pem
715	if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
716
717	!-------------------------- END INITIALIZATION -------------------------
718
719	!-------------------------------- RUN ----------------------------------
720	!------------------------
721	! II Run
722	! II_a Update pressure, ice and tracers
723	!------------------------
724	i_myear_leg = 0
725	stopPEM = 0
726	if (layering_algo) then
727	allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
728	new_str = .true.
729	new_lag1 = .true.
730	new_lag2 = .true.
731	do islope = 1,nslope
732	do ig = 1,ngrid
733	current1(ig,islope)%p => stratif(ig,islope)%top
734	current2(ig,islope)%p => stratif(ig,islope)%top
735	enddo
736	enddo
737	endif
738
739	do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
740	! II.a.1. Compute updated global pressure
741	write(,) "Recomputing the new pressure..."
742	do i = 1,ngrid
743	do islope = 1,nslope
744	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
745	enddo
746	enddo
747
748	if (adsorption_pem) then
749	do i = 1,ngrid
750	global_avg_press_new = global_avg_press_new - gcell_area(i)delta_co2_adsorbded(i)/Total_surface
751	enddo
752	endif
753	write(,) 'Global average pressure old time step',global_avg_press_old
754	write(,) 'Global average pressure new time step',global_avg_press_new
755
756	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
757	allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
758	write(,) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
759	do l = 1,nlayer + 1
760	do ig = 1,ngrid
761	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
762	enddo
763	enddo
764
765	! II.a.3. Surface pressure timeseries
766	write(,) "Recomputing the surface pressure timeserie adapted to the new pressure..."
767	do ig = 1,ngrid
768	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
769	enddo
770
771	! II.a.4. New pressure levels timeseries
772	allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
773	write(,) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
774	do l = 1,nlayer + 1
775	do ig = 1,ngrid
776	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
777	enddo
778	enddo
779
780	! II.a.5. Tracers timeseries
781	write(,) "Recomputing of tracer VMR timeseries for the new pressure..."
782
783	l = 1
784	do ig = 1,ngrid
785	do t = 1,timelen
786	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))/ &
787	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
788	if (q_h2o_PEM_phys(ig,t) < 0) then
789	q_h2o_PEM_phys(ig,t) = 1.e-30
790	else if (q_h2o_PEM_phys(ig,t) > 1) then
791	q_h2o_PEM_phys(ig,t) = 1.
792	endif
793	enddo
794	enddo
795
796	do ig = 1,ngrid
797	do t = 1,timelen
798	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))/ &
799	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
800	+ ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
801	- (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
802	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
803	if (q_co2_PEM_phys(ig,t) < 0) then
804	q_co2_PEM_phys(ig,t) = 1.e-30
805	else if (q_co2_PEM_phys(ig,t) > 1) then
806	q_co2_PEM_phys(ig,t) = 1.
807	endif
808	mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
809	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
810	enddo
811	enddo
812
813	deallocate(zplev_new_timeseries,zplev_old_timeseries)
814
815	!------------------------
816	! II Run
817	! II_b Evolution of ice
818	!------------------------
819	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
820	call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
821	if (layering_algo) then
822	do islope = 1,nslope
823	do ig = 1,ngrid
824	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)
825	enddo
826	enddo
827	endif
828
829	!------------------------
830	! II Run
831	! II_c Flow of glaciers
832	!------------------------
833	if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
834	global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
835	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)
836
837	!------------------------
838	! II Run
839	! II_d Update surface and soil temperatures
840	!------------------------
841	! II_d.1 Update Tsurf
842	write(,) "Updating the new Tsurf"
843	bool_sublim = .false.
844	Tsurfavg_before_saved = tsurf_avg
845	do ig = 1,ngrid
846	do islope = 1,nslope
847	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
848	co2ice_disappeared(ig,islope) = .true.
849	if (latitude_deg(ig) > 0) then
850	do ig_loop = ig,ngrid
851	do islope_loop = islope,iflat,-1
852	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
853	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
854	bool_sublim = .true.
855	exit
856	endif
857	enddo
858	if (bool_sublim) exit
859	enddo
860	else
861	do ig_loop = ig,1,-1
862	do islope_loop = islope,iflat
863	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
864	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
865	bool_sublim = .true.
866	exit
867	endif
868	enddo
869	if (bool_sublim) exit
870	enddo
871	endif
872	if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
873	albedo(ig,1,islope) = albedo_h2o_frost
874	albedo(ig,2,islope) = albedo_h2o_frost
875	else
876	albedo(ig,1,islope) = albedodat(ig)
877	albedo(ig,2,islope) = albedodat(ig)
878	emis(ig,islope) = emissiv
879	endif
880	else if ((co2_ice(ig,islope) > 1.e-3) .and. (d_co2ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
881	ave = 0.
882	do t = 1,timelen
883	ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
884	enddo
885	tsurf_avg(ig,islope) = ave/timelen
886	endif
887	enddo
888	enddo
889
890	do t = 1,timelen
891	tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_avg - Tsurfavg_before_saved
892	enddo
893	! for the start
894	do ig = 1,ngrid
895	do islope = 1,nslope
896	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_avg(ig,islope))
897	enddo
898	enddo
899
900	if (soil_pem) then
901
902	! II_d.2 Update soil temperature
903	write(,)"Updating soil temperature"
904	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
905	do islope = 1,nslope
906	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
907	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
908
909	do t = 1,timelen
910	Tsoil_locslope(:,1:nsoilmx) = tsoil_PEM(:,1:nsoilmx,islope) + tsoil_anom(:,:,islope,t)
911	Tsoil_locslope(:,nsoilmx + 1:) = tsoil_PEM(:,nsoilmx + 1:,islope)
912	do ig = 1,ngrid
913	do isoil = 1,nsoilmx_PEM
914	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)
915	if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
916	enddo
917	enddo
918	enddo
919	enddo
920	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
921
922	write(,) "Update of soil temperature done"
923
924	deallocate(Tsoil_locslope)
925
926	! II_d.3 Update the ice table
927	if (icetable_equilibrium) then
928	write(,) "Compute ice table (equilibrium method)"
929	icetable_thickness_prev_iter = icetable_thickness
930	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
931	call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
932	else if (icetable_dynamic) then
933	write(,) "Compute ice table (dynamic method)"
934	allocate(porefill(nsoilmx_PEM))
935	do ig = 1,ngrid
936	do islope = 1,nslope
937	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)
938	icetable_depth(ig,islope) = ssi_depth
939	ice_porefilling(ig,:,islope) = porefill
940	enddo
941	enddo
942	deallocate(porefill)
943	!call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
944	endif
945
946	! II_d.4 Update the soil thermal properties
947	call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
948
949	! II_d.5 Update the mass of the regolith adsorbed
950	if (adsorption_pem) then
951	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice, &
952	h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
953	h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
954
955	totmassco2_adsorbded = 0.
956	totmassh2o_adsorbded = 0.
957	do ig = 1,ngrid
958	do islope = 1,nslope
959	do l = 1,nsoilmx_PEM
960	if (l == 1) then
961	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
962	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
963	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
964	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
965	else
966	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
967	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
968	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
969	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
970	endif
971	enddo
972	enddo
973	enddo
974	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
975	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
976	endif
977	endif !soil_pem
978
979	!------------------------
980	! II Run
981	! II_e Outputs
982	!------------------------
983	call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/global_avg_press_new/))
984	do islope = 1,nslope
985	write(str2(1:2),'(i2.2)') islope
986	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
987	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
988	call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
989	call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
990	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
991	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
992	if (icetable_equilibrium) then
993	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
994	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
995	else if (icetable_dynamic) 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	endif
999
1000	if (soil_pem) then
1001	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
1002	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
1003	if (adsorption_pem) then
1004	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
1005	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
1006	endif
1007	endif
1008	enddo
1009
1010	!------------------------
1011	! II Run
1012	! II_f Update the tendencies
1013	!------------------------
1014	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)
1015
1016	!------------------------
1017	! II Run
1018	! II_g Checking the stopping criterion
1019	!------------------------
1020
1021	write(,) "Checking the stopping criteria..."
1022	call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
1023	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)
1024	i_myear_leg = i_myear_leg + dt
1025	i_myear = i_myear + dt
1026	if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
1027	if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
1028	call system_clock(c2)
1029	if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
1030	if (stopPEM > 0) then
1031	select case (stopPEM)
1032	case(1)
1033	write(,) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
1034	case(2)
1035	write(,) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
1036	case(3)
1037	write(,) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
1038	case(4)
1039	write(,) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
1040	case(5)
1041	write(,) "STOPPING because maximum number of iterations is reached (possibly due to orbital parameters):", stopPEM
1042	case(6)
1043	write(,) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
1044	case(7)
1045	write(,) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM
1046	case(8)
1047	write(,) "STOPPING because the layering algorithm met an hasty end:", stopPEM
1048	case default
1049	write(,) "STOPPING with unknown stopping criterion code:", stopPEM
1050	end select
1051	exit
1052	else
1053	write(,) "The PEM can continue!"
1054	write(,) "Number of iterations of the PEM: i_myear_leg =", i_myear_leg
1055	write(,) "Number of simulated Martian years: i_myear =", i_myear
1056	endif
1057
1058	global_avg_press_old = global_avg_press_new
1059
1060	enddo ! big time iteration loop of the pem
1061	!------------------------------ END RUN --------------------------------
1062
1063	!------------------------------- OUTPUT --------------------------------
1064	!------------------------
1065	! III Output
1066	! III_a Update surface value for the PCM start files
1067	!------------------------
1068	! III_a.1 Ice update (for startfi)
1069
1070	watercap = 0.
1071	perennial_co2ice = co2_ice
1072	do ig = 1,ngrid
1073	! H2O ice metamorphism
1074	if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_h2oice_threshold) then
1075	h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1076	qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
1077	endif
1078
1079	! Is H2O ice still considered as an infinite reservoir for the PCM?
1080	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
1081	! There is enough ice to be considered as an infinite reservoir
1082	watercaptag(ig) = .true.
1083	else
1084	! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
1085	watercaptag(ig) = .false.
1086	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1087	h2o_ice(ig,:) = 0.
1088	endif
1089
1090	! CO2 ice metamorphism
1091	if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_co2ice_threshold) then
1092	perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1093	qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
1094	endif
1095	enddo
1096
1097	! III_a.2 Tsoil update (for startfi)
1098	if (soil_pem) then
1099	call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
1100	tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen)
1101	#ifndef CPP_STD
1102	flux_geo = fluxgeo
1103	#endif
1104	endif
1105	deallocate(tsoil_anom)
1106
1107	! III_a.4 Pressure (for start)
1108	ps = ps*global_avg_press_new/global_avg_press_PCM
1109	ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
1110
1111	! III_a.5 Tracer (for start)
1112	allocate(zplev_new(ngrid,nlayer + 1))
1113
1114	do l = 1,nlayer + 1
1115	zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
1116	enddo
1117
1118	do nnq = 1,nqtot
1119	if (noms(nnq) /= "co2") then
1120	do l = 1,llm - 1
1121	do ig = 1,ngrid
1122	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))
1123	enddo
1124	q(:,llm,nnq) = q(:,llm - 1,nnq)
1125	enddo
1126	else
1127	do l = 1,llm - 1
1128	do ig = 1,ngrid
1129	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)) &
1130	+ ((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))
1131	enddo
1132	q(:,llm,nnq) = q(:,llm - 1,nnq)
1133	enddo
1134	endif
1135	enddo
1136
1137	! Conserving the tracers mass for PCM start files
1138	do nnq = 1,nqtot
1139	do ig = 1,ngrid
1140	do l = 1,llm - 1
1141	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
1142	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1143	q(ig,l,nnq) = 1.
1144	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1145	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1146	endif
1147	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1148	enddo
1149	enddo
1150	enddo
1151
1152	if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
1153
1154	!------------------------
1155	! III Output
1156	! III_b Write "restart.nc" and "restartfi.nc"
1157	!------------------------
1158	! III_b.1 Write "restart.nc"
1159	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1160	pday = day_ini
1161	ztime_fin = time_phys
1162
1163	allocate(p(ip1jmp1,nlayer + 1))
1164	#ifndef CPP_1D
1165	! Correction on teta due to surface pressure changes
1166	do l = 1,nlayer
1167	do i = 1,ip1jmp1
1168	teta(i,l)= teta(i,l)(ps0(i)/ps(i))*kappa
1169	enddo
1170	enddo
1171	! Correction on atmospheric pressure
1172	call pression(ip1jmp1,ap,bp,ps,p)
1173	! Correction on the mass of atmosphere
1174	call massdair(p,masse)
1175	call dynredem0("restart.nc",day_ini,phis)
1176	call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps)
1177	write(,) "restart.nc has been written"
1178	#else
1179	call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1180	write(,) "restart1D.txt has been written"
1181	#endif
1182
1183	! III_b.2 Write the "restartfi.nc"
1184	#ifndef CPP_STD
1185	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1186	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1187	inertiedat,def_slope,subslope_dist)
1188	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1189	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1190	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1191	wstar,watercap,perennial_co2ice)
1192	#else
1193	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1194	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
1195	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1196	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1197	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1198	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
1199	tsea_ice,sea_ice)
1200	#endif
1201	write(,) "restartfi.nc has been written"
1202
1203	!------------------------
1204	! III Output
1205	! III_c Write the "restartpem.nc"
1206	!------------------------
1207	if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
1208	call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1209	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
1210	TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
1211	co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
1212	write(,) "restartpem.nc has been written"
1213
1214	call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
1215
1216	write(,) "The PEM has run for", i_myear_leg, "Martian years."
1217	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1218	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1219	write(,) "PEM: so far, so good!"
1220
1221	if (layering_algo) then
1222	do islope = 1,nslope
1223	do i = 1,ngrid
1224	call del_layering(stratif(i,islope))
1225	enddo
1226	enddo
1227	deallocate(new_str,new_lag1,new_lag2,current1,current2)
1228	endif
1229	deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
1230	deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
1231	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
1232	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,icetable_thickness_prev_iter)
1233	deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
1234	!----------------------------- END OUTPUT ------------------------------
1235
1236	END PROGRAM pem

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

Download in other formats:

Original Format