Context Navigation

pem.F90 @ 3989

Last change on this file since 3989 was 3989, checked in by jbclement, 6 days ago

PEM:
Massive structural refactor of the PEM codebase for improved readability, consistency and maintainability. The goal is to modernize, standardize and consolidate the code while removing legacy complexity. In detail, this change:

Performs large-scale cleanup removing unused code, obsolete routines, duplicated functionality and deprecated initialization logic;
Removes "*_mod" wrappers;
Replaces mixed naming conventions with clearer variable names, domain-based file/module names and purpose-based routine names;
Adds native reading/writing capabilities to the PEM removing the cumbersome dependency on Mars PCM subroutines;
Centralizes soil/ice/orbital/PEM configuration logic into dedicated modules;
Simplifies routines and updates calls/interfaces to match the new structure.

This refactor significantly clarifies the codebase and provides a cleaner foundation for forthcoming developments.
JBC

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

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

Download in other formats:

Original Format