Context Navigation

pem.F90 @ 3985

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

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

Download in other formats:

Original Format