Context Navigation

pem.F90 @ 3840

Last change on this file since 3840 was 3840, checked in by jbclement, 18 hours ago

PEM:

Correction of a bug in the launching script.
Update of "visu_evol_layering.py", in particular to show value at cursor for 2D heatmaps.
Few cleanings.

JBC

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

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

Download in other formats:

Original Format