Context Navigation

pem.F90 @ 3995

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

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

Download in other formats:

Original Format