Context Navigation

pem.F90 @ 3979

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

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

Download in other formats:

Original Format