Context Navigation

pem.F90 @ 3974

Last change on this file since 3974 was 3974, checked in by jbclement, 29 hours ago
Mars PCM: Simplification of the orbital parameters initialization (resolve ticket #109): removing redundant actions in subroutine 'iniorbit' in regard of what is done in "tabfi.F" and moving the subroutine 'lsp2solp' from "tabfi.F" to "planete_h.F90" + some cleanings. JBC
File size: 60.8 KB

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

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

Download in other formats:

Original Format