Context Navigation

pem.F90 @ 3983

Last change on this file since 3983 was 3983, checked in by jbclement, 7 days ago

PEM:

Removing completely the ice metamorphism computed by a threshold at the end of the PCM (which was commented).
Addition of a module "metamorphism" to compute the PCM frost at the PEM beginning and give it back to the PCM at the PEM end. The frost is considered as the ice given by the PCM "startfi.nc" which is above the yearly minimum. Thereby, metamorphism is performed through this operation.
Ice reservoirs representation in the PEM is modernized.

JBC

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

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

Download in other formats:

Original Format