Context Navigation

pem.F90 @ 3961

Last change on this file since 3961 was 3961, checked in by jbclement, 8 days ago
PEM: Clearing all the tags and code related to the Generic PCM. JBC
File size: 62.1 KB

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

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

Download in other formats:

Original Format