Context Navigation

pem.F90 @ 3967

Last change on this file since 3967 was 3967, checked in by jbclement, 2 days ago
PEM: New management of stopping criteria: giving explicit names and collecting everything in the same place. JBC
File size: 60.9 KB

Line
1	!------------------------
2	! I Initialization
3	! I_a Read the "run.def"
4	! I_b Read the "start.nc" and "startfi.nc"
5	! I_c Subslope parametrisation
6	! I_d Read the PCM data and convert them to the physical grid
7	! I_e Initialization of the PEM variable and soil
8	! I_f Compute tendencies
9	! I_g Compute global surface pressure
10	! I_h Read the "startpem.nc"
11	! I_i Compute orbit criterion
12
13	! II Run
14	! II_a Update pressure, ice and tracers
15	! II_b Evolution of ice
16	! II_c Flow of glaciers
17	! II_d Update surface and soil temperatures
18	! II_e Outputs
19	! II_f Update the tendencies
20	! II_g Checking the stopping criterion
21
22	! III Output
23	! III_a Update surface values for the PCM start files
24	! III_b Write the "restart.nc" and "restartfi.nc"
25	! III_c Write the "restartpem.nc"
26	!------------------------
27
28	PROGRAM pem
29
30	use phyetat0_mod, only: phyetat0
31	use phyredem, only: physdem0, physdem1
32	use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
33	use turb_mod, only: q2, wstar
34	use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h
35	use logic_mod, only: iflag_phys
36	use mod_const_mpi, only: COMM_LMDZ
37	use infotrac
38	use geometry_mod, only: latitude_deg
39	use conf_pem_mod, only: conf_pem
40	use pemredem, only: pemdem0, pemdem1
41	use glaciers_mod, only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
42	metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
43	use stopping_crit_mod, only: stopping_crit_h2o_ice, stopping_crit_co2, 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 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	type(stopFlags) :: stopCrit ! Stopping criteria
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	if (layering_algo) then
615	allocate(h2o_ice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
616	new_str = .true.
617	new_lag = .true.
618	do islope = 1,nslope
619	do ig = 1,ngrid
620	current(ig,islope)%p => layerings_map(ig,islope)%top
621	enddo
622	enddo
623	endif
624
625	do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
626	! II.a.1. Compute updated global pressure
627	write(,'(a,f10.2)') ' *** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
628	write(,) "> Updating the surface pressure"
629	ps_avg_global_old = ps_avg_global_new
630	do i = 1,ngrid
631	do islope = 1,nslope
632	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
633	enddo
634	enddo
635	if (adsorption_pem) then
636	do i = 1,ngrid
637	ps_avg_global_new = ps_avg_global_new - gcell_area(i)delta_co2_adsorbed(i)/total_surface
638	enddo
639	endif
640	ps_avg = ps_avg*ps_avg_global_new/ps_avg_global_old
641	write(,) 'Global average pressure old time step:',ps_avg_global_old
642	write(,) 'Global average pressure new time step:',ps_avg_global_new
643
644	! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
645	write(,) "> Updating the surface pressure timeseries for the new pressure"
646	allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
647	do l = 1,nlayer + 1
648	do ig = 1,ngrid
649	zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
650	enddo
651	enddo
652	ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old
653	write(,) "> Updating the pressure levels timeseries for the new pressure"
654	allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
655	do l = 1,nlayer + 1
656	do ig = 1,ngrid
657	zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
658	enddo
659	enddo
660
661	! II.a.3. Tracers timeseries
662	write(,) "> Updating the tracer VMR timeseries for the new pressure"
663	allocate(vmr_co2_PEM_phys(ngrid,timelen))
664	l = 1
665	do ig = 1,ngrid
666	do t = 1,timelen
667	! H2O
668	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))/ &
669	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
670	if (q_h2o_PEM_phys(ig,t) < 0) then
671	q_h2o_PEM_phys(ig,t) = 1.e-30
672	else if (q_h2o_PEM_phys(ig,t) > 1) then
673	q_h2o_PEM_phys(ig,t) = 1.
674	endif
675	! CO2
676	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))/ &
677	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
678	+ ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
679	- (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/ &
680	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
681	if (q_co2_PEM_phys(ig,t) < 0) then
682	q_co2_PEM_phys(ig,t) = 1.e-30
683	else if (q_co2_PEM_phys(ig,t) > 1) then
684	q_co2_PEM_phys(ig,t) = 1.
685	endif
686	mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
687	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
688	enddo
689	enddo
690	deallocate(zplev_timeseries_new,zplev_timeseries_old)
691
692	!------------------------
693	! II Run
694	! II_b Evolution of ice
695	!------------------------
696	allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
697	if (layering_algo) then
698	h2o_ice_depth_old = h2o_ice_depth
699
700	allocate(d_h2oice_new(ngrid,nslope))
701	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)
702	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)
703
704	do islope = 1,nslope
705	do ig = 1,ngrid
706	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)
707	!call print_layering(layerings_map(ig,islope))
708	co2_ice(ig,islope) = 0.
709	h2o_ice(ig,islope) = 0.
710	h2o_ice_depth(ig,islope) = 0.
711	if (is_co2ice_str(layerings_map(ig,islope)%top)) then
712	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice*rho_co2ice
713	else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
714	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice*rho_h2oice
715	else
716	call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
717	endif
718	enddo
719	enddo
720	deallocate(d_h2oice_new)
721	else
722	zlag = 0.
723	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopCrit)
724	call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
725	endif
726
727	!------------------------
728	! II Run
729	! II_c Flow of glaciers
730	!------------------------
731	allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
732	if (co2ice_flow .and. nslope > 1) then
733	call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
734	ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
735	if (layering_algo) then
736	do islope = 1,nslope
737	do ig = 1,ngrid
738	layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)/rho_co2ice
739	enddo
740	enddo
741	endif
742	endif
743	if (h2oice_flow .and. nslope > 1) then
744	call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
745	if (layering_algo) then
746	do islope = 1,nslope
747	do ig = 1,ngrid
748	layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)/rho_h2oice
749	enddo
750	enddo
751	endif
752	endif
753
754	!------------------------
755	! II Run
756	! II_d Update surface and soil temperatures
757	!------------------------
758	! II_d.1 Update Tsurf
759	call update_tsurf_nearest_baresoil(ngrid,nslope,iim,jjm_value,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
760
761	if (soil_pem) then
762	! II_d.2 Shifting soil temperature to surface
763	call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
764
765	! II_d.3 Update soil temperature
766	write(,)"> Updating soil temperature profile"
767	allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen))
768	tsoil_PEM_timeseries_old = tsoil_PEM_timeseries
769	do islope = 1,nslope
770	tsoil_avg_old = tsoil_PEM(:,:,islope)
771	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
772	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
773
774	do t = 1,timelen
775	do ig = 1,ngrid
776	do isoil = 1,nsoilmx_PEM
777	! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
778	tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
779	! Update of watersoil density
780	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)
781	if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
782	enddo
783	enddo
784	enddo
785	enddo
786	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
787	deallocate(tsoil_avg_old)
788
789	! II_d.4 Update the ice table
790	allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope))
791	if (icetable_equilibrium) then
792	write(,) "> Updating ice table (equilibrium method)"
793	icetable_thickness_old = icetable_thickness
794	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
795	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
796	else if (icetable_dynamic) then
797	write(,) "> Updating ice table (dynamic method)"
798	ice_porefilling_old = ice_porefilling
799	icetable_depth_old = icetable_depth
800	allocate(porefill(nsoilmx_PEM))
801	do ig = 1,ngrid
802	do islope = 1,nslope
803	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)
804	icetable_depth(ig,islope) = ssi_depth
805	ice_porefilling(ig,:,islope) = porefill
806	enddo
807	enddo
808	deallocate(porefill)
809	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
810	endif
811	deallocate(icetable_thickness_old,ice_porefilling_old)
812
813	! II_d.5 Update the soil thermal properties
814	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)
815
816	! II_d.6 Update the mass of the regolith adsorbed
817	totmass_adsco2 = 0.
818	totmass_adsh2o = 0.
819	if (adsorption_pem) then
820	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
821	tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
822	h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
823	do ig = 1,ngrid
824	do islope = 1,nslope
825	do l = 1,nsoilmx_PEM
826	if (l == 1) then
827	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
828	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
829	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
830	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
831	else
832	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
833	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
834	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
835	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
836	endif
837	enddo
838	enddo
839	enddo
840	write(,) "Total mass of CO2 in the regolith =", totmass_adsco2
841	write(,) "Total mass of H2O in the regolith =", totmass_adsh2o
842	endif
843	endif !soil_pem
844	deallocate(zshift_surf,zlag)
845
846	!------------------------
847	! II Run
848	! II_e Outputs
849	!------------------------
850	call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
851	do islope = 1,nslope
852	write(str2(1:2),'(i2.2)') islope
853	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
854	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
855	call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
856	call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
857	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope)))
858	call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope)))
859	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
860	if (icetable_equilibrium) then
861	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
862	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope))
863	else if (icetable_dynamic) then
864	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
865	endif
866
867	if (soil_pem) then
868	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
869	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
870	if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
871	if (adsorption_pem) then
872	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope))
873	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope))
874	endif
875	endif
876	enddo
877	deallocate(flag_co2flow,flag_h2oflow)
878
879	! Checking mass balance for CO2
880	if (abs(CO2cond_ps - 1.) < 1.e-10) then
881	totmass_co2ice = 0.
882	totmass_atmco2 = 0.
883	do ig = 1,ngrid
884	totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
885	do islope = 1,nslope
886	totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)cell_area(ig)subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
887	enddo
888	enddo
889	totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,1.e-10)
890	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, ' %'
891	if (abs((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini) > 0.01) then
892	write(,) ' /!\ Warning: mass balance is not conserved!'
893	totmass_ini = max(totmass_atmco2_ini,1.e-10)
894	write(,) ' Atmospheric CO2 mass balance = ', 100.*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini, ' %'
895	totmass_ini = max(totmass_co2ice_ini,1.e-10)
896	write(,) ' CO2 ice mass balance = ', 100.*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini, ' %'
897	totmass_ini = max(totmass_adsco2_ini,1.e-10)
898	write(,) ' Adsorbed CO2 mass balance = ', 100.*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini, ' %'
899	endif
900	endif
901
902	!------------------------
903	! II Run
904	! II_f Update the tendencies
905	!------------------------
906	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)
907	write(,) "> Updating the H2O sub-surface ice tendency due to lag layer"
908	if (layering_algo) then
909	do ig = 1,ngrid
910	do islope = 1,nslope
911	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))
912	enddo
913	enddo
914	!~ else
915	!~ do ig = 1,ngrid
916	!~ do islope = 1,nslope
917	!~ 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))
918	!~ enddo
919	!~ enddo
920	endif
921	if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old)
922	deallocate(vmr_co2_PEM_phys)
923
924	!------------------------
925	! II Run
926	! II_g Checking the stopping criterion
927	!------------------------
928	write(,) "> Checking the stopping criteria"
929	call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopCrit,ngrid)
930	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)
931	i_myear_leg = i_myear_leg + dt
932	i_myear = i_myear + dt
933	if (stopCrit%stop_code() == 0 .and. i_myear_leg >= n_myear_leg) stopCrit%max_iter_reached = .true.
934	if (stopCrit%stop_code() == 0 .and. i_myear >= n_myear) stopCrit%max_years_reached = .true.
935	call system_clock(c2)
936	if (stopCrit%stop_code() == 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopCrit%time_limit_reached = .true.
937	if (stopCrit%is_any_set()) then
938	write(,) stopCrit%stop_message()
939	exit
940	else
941	write(,'(a,f10.2,a)') ' *** The chained simulation has run for ',i_myear,' Martian years.'
942	write(,) '**** The PEM can continue!'
943	write(,) '****'
944	endif
945	enddo ! big time iteration loop of the pem
946	deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed)
947	deallocate(watersoil_density_PEM_avg,watersurf_density_avg)
948	deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries)
949	deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim)
950	deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
951	deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
952	if (layering_algo) deallocate(h2o_ice_depth_old,new_str,new_lag,current)
953	!------------------------------ END RUN --------------------------------
954
955	!------------------------------- OUTPUT --------------------------------
956	!------------------------
957	! III Output
958	! III_a Update surface values for the PCM start files
959	!------------------------
960	write(,)
961	write(,) '******* PEM finalization *******'
962	! III_a.1 Ice update for start file
963	write(,) '> Reconstructing perennial ice and frost for the PCM'
964	watercap = 0.
965	perennial_co2ice = co2_ice
966	do ig = 1,ngrid
967	! H2O ice metamorphism
968	!if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_h2oice_threshold) then
969	! h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
970	! qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
971	!endif
972
973	! Is H2O ice still considered as an infinite reservoir for the PCM?
974	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
975	! There is enough ice to be considered as an infinite reservoir
976	watercaptag(ig) = .true.
977	else
978	! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
979	watercaptag(ig) = .false.
980	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
981	h2o_ice(ig,:) = 0.
982	endif
983
984	! CO2 ice metamorphism
985	!if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_co2ice_threshold) then
986	! perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
987	! qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
988	!endif
989	enddo
990
991	! III.a.3. Tsurf update for start file
992	write(,) '> Reconstructing the surface temperature for the PCM'
993	tsurf = tsurf_avg + tsurf_dev
994	deallocate(tsurf_dev)
995
996	! III_a.4 Tsoil update for start file
997	if (soil_pem) then
998	write(,) '> Reconstructing the soil temperature profile for the PCM'
999	inertiesoil = TI_PEM(:,:nsoilmx,:)
1000	! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value
1001	do isoil = 1,nsoilmx
1002	tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:)
1003	enddo
1004	tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev
1005	flux_geo = fluxgeo
1006	endif
1007	deallocate(tsurf_avg,tsoil_dev)
1008
1009	! III_a.5 Pressure update for start file
1010	write(,) '> Reconstructing the pressure for the PCM'
1011	allocate(ps_start(ngrid))
1012	! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop
1013	ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini
1014	deallocate(ps_avg,ps_dev)
1015
1016	! III_a.6 Tracers update for start file
1017	write(,) '> Reconstructing the tracer VMR for the PCM'
1018	allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
1019	do l = 1,nlayer + 1
1020	zplev_start0(:,l) = ap(l) + bp(l)*ps_start0
1021	zplev_new(:,l) = ap(l) + bp(l)*ps_start
1022	enddo
1023
1024	do nnq = 1,nqtot
1025	if (noms(nnq) /= "co2") then
1026	do l = 1,llm - 1
1027	do ig = 1,ngrid
1028	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))
1029	enddo
1030	q(:,llm,nnq) = q(:,llm - 1,nnq)
1031	enddo
1032	else
1033	do l = 1,llm - 1
1034	do ig = 1,ngrid
1035	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)) &
1036	+ ((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))
1037	enddo
1038	q(:,llm,nnq) = q(:,llm - 1,nnq)
1039	enddo
1040	endif
1041	enddo
1042	deallocate(zplev_start0)
1043
1044	! Conserving the tracers mass for start file
1045	do nnq = 1,nqtot
1046	do ig = 1,ngrid
1047	do l = 1,llm - 1
1048	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
1049	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1050	q(ig,l,nnq) = 1.
1051	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1052	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1053	endif
1054	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1055	enddo
1056	enddo
1057	enddo
1058	deallocate(zplev_new)
1059
1060	! III_a.7 Albedo update for start file
1061	write(,) '> Reconstructing the albedo for the PCM'
1062	do ig = 1,ngrid
1063	if (latitude(ig) < 0.) then
1064	icap = 2 ! Southern hemisphere
1065	else
1066	icap = 1 ! Northern hemisphere
1067	endif
1068	do islope = 1,nslope
1069	! Bare ground
1070	albedo(ig,:,islope) = albedodat(ig)
1071	emis(ig,islope) = emissiv
1072
1073	! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant
1074	! H2O ice
1075	if (h2o_ice(ig,islope) > 0.) then
1076	albedo(ig,:,islope) = albedo_h2o_cap
1077	emis(ig,islope) = 1.
1078	endif
1079	! CO2 ice
1080	if (co2_ice(ig,islope) > 0.) then
1081	albedo(ig,:,islope) = albedo_perennialco2(icap)
1082	emis(ig,islope) = emisice(icap)
1083	endif
1084	! H2O frost
1085	if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
1086	albedo(ig,:,islope) = albedo_h2o_frost
1087	emis(ig,islope) = 1.
1088	endif
1089	! CO2 frost
1090	if (qsurf(ig,igcm_co2,islope) > 0.) then
1091	albedo(ig,:,islope) = albedice(icap)
1092	emis(ig,islope) = emisice(icap)
1093	endif
1094	enddo
1095	enddo
1096
1097	! III_a.8 Orbital parameters update for start file
1098	write(,) '> Setting the new orbital parameters'
1099	if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
1100
1101	!------------------------
1102	! III Output
1103	! III_b Write "restart.nc" and "restartfi.nc"
1104	!------------------------
1105	! III_b.1 Write "restart.nc"
1106	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1107	pday = day_ini
1108	ztime_fin = time_phys
1109	#ifndef CPP_1D
1110	write(,) '> Writing "restart.nc"'
1111	! Correction on teta due to surface pressure changes
1112	allocate(pdyn(ip1jmp1))
1113	call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn)
1114	do i = 1,ip1jmp1
1115	teta(i,:) = teta(i,:)pdyn(i)*rcp
1116	enddo
1117	! Correction on atmospheric pressure
1118	allocate(p(ip1jmp1,nlayer + 1))
1119	call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start,pdyn)
1120	call pression(ip1jmp1,ap,bp,pdyn,p)
1121	! Correction on the mass of atmosphere
1122	call massdair(p,masse)
1123	call dynredem0("restart.nc",day_ini,phis)
1124	call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,pdyn)
1125	deallocate(ap,bp,p,pdyn)
1126	#else
1127	write(,) '> Writing "restart1D.txt"'
1128	call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1129	#endif
1130	deallocate(ps_start0,ps_start)
1131
1132	! III_b.2 Write the "restartfi.nc"
1133	write(,) '> Writing "restartfi.nc"'
1134	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1135	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1136	inertiedat,def_slope,subslope_dist)
1137	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1138	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1139	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1140	wstar,watercap,perennial_co2ice)
1141
1142	!------------------------
1143	! III Output
1144	! III_c Write the "restartpem.nc"
1145	!------------------------
1146	write(,) '> Writing "restartpem.nc"'
1147	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)
1148	call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1149	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
1150	co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
1151
1152	call info_PEM(i_myear_leg,stopCrit%stop_code(),i_myear,n_myear)
1153
1154	write(,)
1155	write(,) '**** PEM final information *****'
1156	write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
1157	write(,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myearconvert_years, " Earth years."
1158	write(,'(a,f16.4,a)') " + The reached date is now ", (year_bp_ini + i_myear)convert_years, " Earth years."
1159	write(,) "+ PEM: so far, so good!"
1160	write(,) '************************************'
1161
1162	if (layering_algo) then
1163	do islope = 1,nslope
1164	do i = 1,ngrid
1165	call del_layering(layerings_map(i,islope))
1166	enddo
1167	enddo
1168	endif
1169	deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
1170	deallocate(co2_ice,h2o_ice,layerings_map)
1171	!----------------------------- END OUTPUT ------------------------------
1172
1173	END PROGRAM pem

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

Download in other formats:

Original Format