Context Navigation

pem.F90 @ 3321

Last change on this file since 3321 was 3319, checked in by jbclement, 7 months ago
PEM: Small correction about the dimension of an array. JBC
File size: 57.1 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 tendencies
9	! I_g Save initial PCM situation
10	! I_h Read the "startpem.nc"
11	! I_i Compute orbit criterion
12
13	! II Run
14	! II_a Update pressure, ice and tracers
15	! II_b Evolution of ice
16	! II_c Flow of glaciers
17	! II_d Update surface and soil temperatures
18	! II_e Outputs
19	! II_f Update the tendencies
20	! II_g Checking the stopping criterion
21
22	! III Output
23	! III_a Update surface value for the PCM start files
24	! III_b Write the "restart.nc" and "restartfi.nc"
25	! III_c Write the "restartpem.nc"
26	!------------------------
27
28	PROGRAM pem
29
30	use phyetat0_mod, only: phyetat0
31	use phyredem, only: physdem0, physdem1
32	use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
33	use turb_mod, only: q2, wstar
34	use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h
35	use logic_mod, only: iflag_phys
36	use mod_const_mpi, only: COMM_LMDZ
37	use infotrac
38	use geometry_mod, only: latitude_deg
39	use conf_pem_mod, only: conf_pem
40	use pemredem, only: pemdem0, pemdem1
41	use glaciers_mod, only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
42	metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice
43	use stopping_crit_mod, only: stopping_crit_h2o_ice, stopping_crit_co2
44	use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
45	use evol_ice_mod, only: evol_co2_ice, evol_h2o_ice
46	use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
47	TI_PEM, & ! soil thermal inertia
48	tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
49	fluxgeo ! Geothermal flux for the PEM and PCM
50	use adsorption_mod, only: regolith_adsorption, adsorption_pem, & ! Bool to check if adsorption, main subroutine
51	ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
52	co2_adsorbded_phys, h2o_adsorbded_phys ! Mass of co2 and h2O adsorbded
53	use time_evol_mod, only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
54	use orbit_param_criterion_mod, only: orbit_param_criterion
55	use recomp_orb_param_mod, only: recomp_orb_param
56	use ice_table_mod, only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
57	ini_ice_table_porefilling, icetable_equilibrium, computeice_table_equilibrium,compute_massh2o_exchange_ssi
58	use soil_thermalproperties_mod, only: update_soil_thermalproperties
59	use time_phylmdz_mod, only: daysec, dtphys
60	use abort_pem_mod, only: abort_pem
61	use soil_settings_PEM_mod, only: soil_settings_PEM
62	use compute_tend_mod, only: compute_tend
63	use info_PEM_mod, only: info_PEM
64	use interpol_TI_PEM2PCM_mod, only: interpol_TI_PEM2PCM
65	use nb_time_step_PCM_mod, only: nb_time_step_PCM
66	use pemetat0_mod, only: pemetat0
67	use read_data_PCM_mod, only: read_data_PCM
68	use recomp_tend_co2_slope_mod, only: recomp_tend_co2_slope
69	use compute_soiltemp_mod, only: compute_tsoil_pem
70	use writediagpem_mod, only: writediagpem, writediagsoilpem
71	use co2condens_mod, only: CO2cond_ps
72	use layering_mod, only: tend_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
73
74	#ifndef CPP_STD
75	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil
76	use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
77	albedodat, zmea, zstd, zsig, zgam, zthe, &
78	albedo_h2o_frost,frost_albedo_threshold, &
79	emissiv, watercaptag, perennial_co2ice, emisice, albedice
80	use dimradmars_mod, only: totcloudfrac, albedo
81	use dust_param_mod, only: tauscaling
82	use tracer_mod, only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
83	use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
84	use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
85	use comcstfi_h, only: pi, rad, g, cpp, mugaz, r
86	use surfini_mod, only: surfini
87	#else
88	use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
89	use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, &
90	PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
91	ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
92	use aerosol_mod, only: iniaerosol
93	use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit
94	use comcstfi_mod, only: pi, rad, g, cpp, mugaz, r
95	#endif
96
97	#ifndef CPP_1D
98	use iniphysiq_mod, only: iniphysiq
99	use control_mod, only: iphysiq, day_step, nsplit_phys
100	#else
101	use time_phylmdz_mod, only: iphysiq, day_step
102	use regular_lonlat_mod, only: init_regular_lonlat
103	use physics_distribution_mod, only: init_physics_distribution
104	use mod_grid_phy_lmdz, only: regular_lonlat
105	use init_testphys1d_mod, only: init_testphys1d
106	use comvert_mod, only: ap, bp
107	use writerestart1D_mod, only: writerestart1D
108	#endif
109
110	implicit none
111
112	include "dimensions.h"
113	include "paramet.h"
114	include "comgeom.h"
115	include "iniprint.h"
116	include "callkeys.h"
117
118	integer ngridmx
119	parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
120
121	! Same variable names as in the PCM
122	integer, parameter :: nlayer = llm ! Number of vertical layer
123	integer :: ngrid ! Number of physical grid points
124	integer :: nq ! Number of tracer
125	integer :: day_ini ! First day of the simulation
126	real :: pday ! Physical day
127	real :: time_phys ! Same as in PCM
128	real :: ptimestep ! Same as in PCM
129	real :: ztime_fin ! Same as in PCM
130
131	! Variables to read start.nc
132	character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM
133
134	! Dynamic variables
135	real, dimension(ip1jm,llm) :: vcov ! vents covariants
136	real, dimension(ip1jmp1,llm) :: ucov ! vents covariants
137	real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle
138	real, dimension(:,:,:), allocatable :: q ! champs advectes
139	real, dimension(ip1jmp1) :: ps ! pression au sol
140	real, dimension(:), allocatable :: ps_start_PCM ! (ngrid) surface pressure
141	real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure
142	real, dimension(ip1jmp1,llm) :: masse ! masse d'air
143	real, dimension(ip1jmp1) :: phis ! geopotentiel au sol
144	real :: time_0
145
146	! Variables to read starfi.nc
147	character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM
148	character(2) :: str2
149	integer :: ncid, status ! Variable for handling opening of files
150	integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files
151	integer :: apvarid, bpvarid ! Variable ID for Netcdf files
152
153	! Variables to read starfi.nc and write restartfi.nc
154	real, dimension(:), allocatable :: longitude ! Longitude read in startfi_name and written in restartfi
155	real, dimension(:), allocatable :: latitude ! Latitude read in startfi_name and written in restartfi
156	real, dimension(:), allocatable :: cell_area ! Cell_area read in startfi_name and written in restartfi
157	real :: Total_surface ! Total surface of the planet
158
159	! Variables for h2o_ice evolution
160	real, dimension(:,:), allocatable :: h2o_ice ! h2o ice in the PEM
161	real, dimension(:,:,:), allocatable :: min_h2o_ice ! Minima of h2o ice at each point for the PCM years [kg/m^2]
162	real :: h2o_surf_ini ! Initial surface of sublimating h2o ice
163	real :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
164	real :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
165	real :: global_avg_press_new ! constant: Global average pressure of current time step
166	real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
167	real, dimension(:,:), allocatable :: zplev_PCM ! same but retrieved from the PCM [kg/m^2]
168	real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
169	real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
170	integer :: stopPEM ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
171	real, save :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg]
172	real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
173	integer :: timelen ! # time samples
174	real :: ave ! intermediate varibale to compute average
175	real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
176	real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1
177
178	! Variables for co2_ice evolution
179	real, dimension(:,:), allocatable :: co2_ice ! co2 ice in the PEM
180	real, dimension(:,:), allocatable :: co2_ice_ini ! Initial amount of co2 ice in the PEM
181	real, dimension(:,:,:), allocatable :: min_co2_ice ! Minimum of co2 ice at each point for the first year [kg/m^2]
182	real :: co2_surf_ini ! Initial surface of sublimating co2 ice
183	real, dimension(:,:), allocatable :: vmr_co2_PCM ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
184	real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM
185	real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
186
187	! Variables for stratification (layering) evolution
188	type(layering), dimension(:,:), allocatable :: stratif ! Layering (linked list of "stratum") for each grid point and slope
189	type(ptrarray), dimension(:,:), allocatable :: current1, current2 ! Current active stratum in the layering
190	logical, dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm
191
192	! Variables for slopes
193	real, dimension(:,:,:), allocatable :: co2_ice_PCM ! Physics x NSLOPE x Times field: co2 ice given by the PCM [kg/m^2]
194	real, dimension(:,:), allocatable :: tend_co2_ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
195	real, dimension(:,:), allocatable :: tend_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM
196	real, dimension(:,:), allocatable :: tend_h2o_ice ! physical point x slope field: Tendency of evolution of perennial h2o ice
197	real, dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow
198	real, dimension(:), allocatable :: flag_co2flow_mesh ! (ngrid) : Flag where there is a CO2 glacier flow
199	real, dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow
200	real, dimension(:), allocatable :: flag_h2oflow_mesh ! (ngrid) : Flag where there is a H2O glacier flow
201
202	! Variables for surface and soil
203	real, dimension(:,:), allocatable :: tsurf_ave ! Physic x SLOPE field: Averaged Surface Temperature [K]
204	real, dimension(:,:,:), allocatable :: tsoil_ave ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
205	real, dimension(:,:,:), allocatable :: tsoil_anom ! Amplitude between instataneous and yearly average soil temperature [K]
206	real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE XTULES field: Surface Temperature in timeseries [K]
207	real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SLOPE XTULES field: Non averaged Soil Temperature [K]
208	real, dimension(:,:,:,:), allocatable :: tsoil_PCM_timeseries ! IG x SLOPE XTULES field: Non averaged Soil Temperature [K]
209	real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
210	real, dimension(:,:), allocatable :: TI_locslope ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI]
211	real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: intermediate when computing Tsoil [K]
212	real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
213	real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
214	real, dimension(:,:), allocatable :: watersurf_density_ave ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
215	real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
216	real, dimension(:,:,:), allocatable :: watersoil_density_PEM_ave ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
217	real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K]
218	real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
219	real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
220	real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
221	real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
222	logical :: bool_sublim ! logical to check if there is sublimation or not
223	real, dimension(:,:), allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
224	real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
225
226	! Some variables for the PEM run
227	real, parameter :: year_step = 1 ! timestep for the pem
228	integer :: year_iter ! number of iteration
229	integer :: year_iter_max ! maximum number of iterations before stopping
230	integer :: i_myear ! Global number of Martian years of the chained simulations
231	integer :: n_myear ! Maximum number of Martian years of the chained simulations
232	real :: timestep ! timestep [s]
233
234	#ifdef CPP_STD
235	real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
236	real :: albedo_h2o_frost ! albedo of h2o frost
237	real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
238	real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
239	real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
240	real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
241	real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic
242	real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case
243	real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case
244	real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case
245	real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case
246	real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
247	logical, dimension(:), allocatable :: watercaptag ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
248	real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case
249	real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case
250	#endif
251
252	#ifdef CPP_1D
253	integer :: nsplit_phys
254	integer, parameter :: jjm_value = jjm - 1
255
256	! Dummy variables to use the subroutine 'init_testphys1d'
257	logical :: therestart1D, therestartfi
258	integer :: ndt, day0
259	real :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau
260	real, dimension(:), allocatable :: zqsat
261	real, dimension(:,:,:), allocatable :: dq, dqdyn
262	real, dimension(nlayer) :: play, w
263	real, dimension(nlayer + 1) :: plev
264	#else
265	integer, parameter :: jjm_value = jjm
266	real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
267	real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
268	#endif
269
270	! Loop variables
271	integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
272
273	! Parallel variables
274	#ifndef CPP_STD
275	is_sequential = .true.
276	is_parallel = .false.
277	is_mpi_root = .true.
278	is_omp_root = .true.
279	is_master = .true.
280	#endif
281
282	! Some constants
283	day_ini = 0 ! test
284	time_phys = 0. ! test
285	ngrid = ngridmx
286	A = (1/m_co2 - 1/m_noco2)
287	B = 1/m_noco2
288	year_day = 669
289	daysec = 88775.
290	timestep = year_day*daysec/year_step
291
292	!----------------------------- INITIALIZATION --------------------------
293	!------------------------
294	! I Initialization
295	! I_a Read the "run.def"
296	!------------------------
297	#ifndef CPP_1D
298	dtphys = 0
299	call conf_gcm(99,.true.)
300	call infotrac_init
301	nq = nqtot
302	allocate(q(ip1jmp1,llm,nqtot))
303	allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
304	#else
305	allocate(q(1,llm,nqtot))
306	allocate(longitude(1),latitude(1),cell_area(1))
307
308	therestart1D = .false. ! Default value
309	inquire(file = 'start1D_evol.txt',exist = therestart1D)
310	if (.not. therestart1D) then
311	write(,) 'There is no "start1D_evol.txt" file!'
312	error stop 'Initialization cannot be done for the 1D PEM.'
313	endif
314	therestartfi = .false. ! Default value
315	inquire(file = 'startfi.nc',exist = therestartfi)
316	if (.not. therestartfi) then
317	write(,) 'There is no "startfi.nc" file!'
318	error stop 'Initialization cannot be done for the 1D PEM.'
319	endif
320
321	call init_testphys1d('start1D_evol.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, &
322	time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
323	play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
324	ps(2) = ps(1)
325	nsplit_phys = 1
326	#endif
327
328	call conf_pem(i_myear,n_myear)
329
330	!------------------------
331	! I Initialization
332	! I_b Read of the "start.nc" and starfi_evol.nc
333	!------------------------
334	! I_b.1 Read "start.nc"
335	allocate(ps_start_PCM(ngrid))
336	#ifndef CPP_1D
337	call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
338
339	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
340
341	call iniconst !new
342	call inigeom
343
344	allocate(ap(nlayer + 1))
345	allocate(bp(nlayer + 1))
346	status = nf90_open(start_name,NF90_NOWRITE,ncid)
347	status = nf90_inq_varid(ncid,"ap",apvarid)
348	status = nf90_get_var(ncid,apvarid,ap)
349	status = nf90_inq_varid(ncid,"bp",bpvarid)
350	status = nf90_get_var(ncid,bpvarid,bp)
351	status = nf90_close(ncid)
352
353	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)
354	#else
355	ps_start_PCM(1) = ps(1)
356	#endif
357
358	! In the PCM, these values are given to the physic by the dynamic.
359	! Here we simply read them in the "startfi.nc" file
360	status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
361
362	status = nf90_inq_varid(ncid,"longitude",lonvarid)
363	status = nf90_get_var(ncid,lonvarid,longitude)
364
365	status = nf90_inq_varid(ncid,"latitude",latvarid)
366	status = nf90_get_var(ncid,latvarid,latitude)
367
368	status = nf90_inq_varid(ncid,"area",areavarid)
369	status = nf90_get_var(ncid,areavarid,cell_area)
370
371	status = nf90_inq_varid(ncid,"soildepth",sdvarid)
372	status = nf90_get_var(ncid,sdvarid,mlayer)
373
374	status = nf90_close(ncid)
375
376	! I_b.2 Read the "startfi.nc"
377	! First we read the initial state (starfi.nc)
378	#ifndef CPP_STD
379	call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
380	tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, &
381	watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
382
383	! Remove unphysical values of surface tracer
384	where (qsurf < 0.) qsurf = 0.
385
386	call surfini(ngrid,nslope,qsurf)
387	#else
388	call phys_state_var_init(nq)
389	if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
390	call initracer(ngrid,nq)
391	call iniaerosol()
392	allocate(tsurf_read_generic(ngrid))
393	allocate(qsurf_read_generic(ngrid,nq))
394	allocate(tsoil_read_generic(ngrid,nsoilmx))
395	allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
396	allocate(emis_read_generic(ngrid))
397	allocate(tsurf(ngrid,1))
398	allocate(qsurf(ngrid,nq,1))
399	allocate(tsoil(ngrid,nsoilmx,1))
400	allocate(emis(ngrid,1))
401	allocate(watercap(ngrid,1))
402	allocate(watercaptag(ngrid))
403	allocate(albedo_read_generic(ngrid,2))
404	allocate(albedo(ngrid,2,1))
405	allocate(inertiesoil(ngrid,nsoilmx,1))
406	call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
407	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
408	qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, &
409	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
410	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
411
412	nslope = 1
413	call ini_comslope_h(ngrid,1)
414
415	qsurf(:,:,1) = qsurf_read_generic
416	tsurf(:,1) = tsurf_read_generic
417	tsoil(:,:,1) = tsoil_read_generic
418	emis(:,1) = emis_read_generic
419	watercap(:,1) = 0.
420	watercaptag(:) = .false.
421	albedo(:,1,1) = albedo_read_generic(:,1)
422	albedo(:,2,1) = albedo_read_generic(:,2)
423	inertiesoil(:,:,1) = inertiedat
424
425	if (nslope == 1) then
426	def_slope(1) = 0
427	def_slope(2) = 0
428	def_slope_mean = 0
429	subslope_dist(:,1) = 1.
430	endif
431
432	! Remove unphysical values of surface tracer
433	qsurf(:,:,1) = qsurf_read_generic
434	where (qsurf < 0.) qsurf = 0.
435	#endif
436
437	do nnq = 1,nqtot ! Why not using ini_tracer ?
438	if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
439	if (noms(nnq) == "h2o_vap") then
440	igcm_h2o_vap = nnq
441	mmol(igcm_h2o_vap) = 18.
442	endif
443	if (noms(nnq) == "co2") igcm_co2 = nnq
444	enddo
445	r = 8.314511*1000./mugaz
446
447	!------------------------
448	! I Initialization
449	! I_c Subslope parametrisation
450	!------------------------
451	! Define some slope statistics
452	iflat = 1
453	do islope = 2,nslope
454	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
455	enddo
456
457	write(,) 'Flat slope for islope = ',iflat
458	write(,) 'corresponding criterium = ',def_slope_mean(iflat)
459
460	allocate(flag_co2flow(ngrid,nslope))
461	allocate(flag_co2flow_mesh(ngrid))
462	allocate(flag_h2oflow(ngrid,nslope))
463	allocate(flag_h2oflow_mesh(ngrid))
464
465	flag_co2flow = 0
466	flag_co2flow_mesh = 0
467	flag_h2oflow = 0
468	flag_h2oflow_mesh = 0
469
470	!------------------------
471	! I Initialization
472	! I_d Read the PCM data and convert them to the physical grid
473	!------------------------
474	! 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
475	call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
476
477	allocate(tsoil_ave(ngrid,nsoilmx,nslope))
478	allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
479	allocate(vmr_co2_PCM(ngrid,timelen))
480	allocate(ps_timeseries(ngrid,timelen))
481	allocate(min_co2_ice(ngrid,nslope,2))
482	allocate(min_h2o_ice(ngrid,nslope,2))
483	allocate(tsurf_avg_yr1(ngrid,nslope))
484	allocate(tsurf_ave(ngrid,nslope))
485	allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
486	allocate(tsoil_PCM_timeseries(ngrid,nsoilmx,nslope,timelen))
487	allocate(q_co2_PEM_phys(ngrid,timelen))
488	allocate(q_h2o_PEM_phys(ngrid,timelen))
489	allocate(co2_ice_PCM(ngrid,nslope,timelen))
490	allocate(watersurf_density_ave(ngrid,nslope))
491	allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
492	allocate(Tsurfavg_before_saved(ngrid,nslope))
493	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
494	allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
495	allocate(delta_co2_adsorbded(ngrid))
496	allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
497	allocate(delta_h2o_icetablesublim(ngrid))
498	allocate(delta_h2o_adsorbded(ngrid))
499	allocate(vmr_co2_PEM_phys(ngrid,timelen))
500
501	write(,) "Downloading data Y1..."
502	call read_data_PCM("data_PCM_Y1.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), &
503	tsurf_avg_yr1,tsoil_ave,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
504	co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries)
505	write(,) "Downloading data Y1 done!"
506
507	! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value
508	write(,) "Downloading data Y2..."
509	call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,2),min_h2o_ice(:,:,2), &
510	tsurf_ave,tsoil_ave,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
511	co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries)
512	write(,) "Downloading data Y2 done!"
513
514	!------------------------
515	! I Initialization
516	! I_e Initialization of the PEM variables and soil
517	!------------------------
518	call end_comsoil_h_PEM
519	call ini_comsoil_h_PEM(ngrid,nslope)
520	call end_adsorption_h_PEM
521	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
522	call end_ice_table_porefilling
523	call ini_ice_table_porefilling(ngrid,nslope)
524
525	if (soil_pem) then
526	allocate(tsoil_anom(ngrid,nsoilmx,nslope))
527	tsoil_anom = tsoil - tsoil_ave ! compute anomaly between Tsoil(t) in the startfi - <Tsoil> to recompute properly tsoil in the restart
528	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
529	tsoil_PEM(:,1:nsoilmx,:) = tsoil_ave
530	tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_PCM_timeseries
531	watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
532	do l = nsoilmx + 1,nsoilmx_PEM
533	tsoil_PEM(:,l,:) = tsoil_ave(:,nsoilmx,:)
534	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
535	enddo
536	watersoil_density_PEM_ave = sum(watersoil_density_PEM_timeseries,4)/timelen
537	endif !soil_pem
538	deallocate(tsoil_ave,tsoil_PCM_timeseries)
539
540	!------------------------
541	! I Initialization
542	! I_f Compute tendencies
543	!------------------------
544	allocate(tend_co2_ice(ngrid,nslope),tend_h2o_ice(ngrid,nslope))
545	allocate(tend_co2_ice_ini(ngrid,nslope))
546
547	! Compute the tendencies of the evolution of ice over the years
548	call compute_tend(ngrid,nslope,min_co2_ice,tend_co2_ice)
549	call compute_tend(ngrid,nslope,min_h2o_ice,tend_h2o_ice)
550	tend_co2_ice_ini = tend_co2_ice
551
552	!------------------------
553	! I Initialization
554	! I_g Save initial PCM situation
555	!------------------------
556	! We save the places where h2o ice is sublimating
557	! We compute the surface of h2o ice sublimating
558	co2_surf_ini = 0.
559	h2o_surf_ini = 0.
560	Total_surface = 0.
561	do i = 1,ngrid
562	Total_surface = Total_surface + cell_area(i)
563	do islope = 1,nslope
564	if (tend_co2_ice(i,islope) < 0.) co2_surf_ini = co2_surf_ini + cell_area(i)*subslope_dist(i,islope)
565	if (tend_h2o_ice(i,islope) < 0.) h2o_surf_ini = h2o_surf_ini + cell_area(i)*subslope_dist(i,islope)
566	enddo
567	enddo
568
569	write(,) "Total initial surface of co2 ice sublimating (slope) =", co2_surf_ini
570	write(,) "Total initial surface of h2o ice sublimating (slope) =", h2o_surf_ini
571	write(,) "Total surface of the planet =", Total_surface
572	allocate(zplev_PCM(ngrid,nlayer + 1))
573
574	do ig = 1,ngrid
575	zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
576	enddo
577
578	global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
579	global_avg_press_PCM = global_avg_press_old
580	global_avg_press_new = global_avg_press_old
581	write(,) "Initial global average pressure =", global_avg_press_PCM
582
583	!------------------------
584	! I Initialization
585	! I_h Read the "startpem.nc"
586	!------------------------
587	allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
588	deallocate(min_co2_ice,min_h2o_ice)
589
590	allocate(stratif(ngrid,nslope))
591	if (layering_algo) then
592	do islope = 1,nslope
593	do i = 1,ngrid
594	call ini_layering(stratif(i,islope))
595	enddo
596	enddo
597	endif
598
599	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
600	porefillingice_thickness,tsurf_avg_yr1,tsurf_ave,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries, &
601	tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM, &
602	watersurf_density_ave,watersoil_density_PEM_ave,co2_adsorbded_phys,delta_co2_adsorbded, &
603	h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
604
605	allocate(co2_ice_ini(ngrid,nslope))
606	co2_ice_ini = co2_ice
607
608	delta_h2o_icetablesublim = 0.
609
610	if (adsorption_pem) then
611	totmassco2_adsorbded = 0.
612	totmassh2o_adsorbded = 0.
613	do ig = 1,ngrid
614	do islope = 1,nslope
615	do l = 1,nsoilmx_PEM - 1
616	if (l==1) then
617	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
618	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
619	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
620	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
621	else
622	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
623	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
624	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
625	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
626	endif
627	enddo
628	enddo
629	enddo
630	write(,) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
631	write(,) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
632	endif ! adsorption
633	deallocate(tsurf_avg_yr1)
634
635	!------------------------
636	! I Initialization
637	! I_i Compute orbit criterion
638	!------------------------
639	#ifndef CPP_STD
640	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
641	#else
642	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
643	#endif
644
645	if (evol_orbit_pem) then
646	call orbit_param_criterion(i_myear,year_iter_max)
647	else
648	year_iter_max = Max_iter_pem
649	endif
650	!-------------------------- END INITIALIZATION -------------------------
651
652	!-------------------------------- RUN ----------------------------------
653	!------------------------
654	! II Run
655	! II_a Update pressure, ice and tracers
656	!------------------------
657	year_iter = 0
658	stopPEM = 0
659	if (layering_algo) then
660	allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
661	new_str = .true.
662	new_lag1 = .true.
663	new_lag2 = .true.
664	do islope = 1,nslope
665	do ig = 1,ngrid
666	current1(ig,islope)%p => stratif(ig,islope)%top
667	current2(ig,islope)%p => stratif(ig,islope)%top
668	enddo
669	enddo
670	endif
671
672	do while (year_iter < 10 .and. i_myear < n_myear)
673	! II.a.1. Compute updated global pressure
674	write(,) "Recomputing the new pressure..."
675	do i = 1,ngrid
676	do islope = 1,nslope
677	global_avg_press_new = global_avg_press_new - CO2cond_psgcell_area(i)tend_co2_ice(i,islope)subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
678	enddo
679	enddo
680
681	if (adsorption_pem) then
682	do i = 1,ngrid
683	global_avg_press_new = global_avg_press_new - gcell_area(i)delta_co2_adsorbded(i)/Total_surface
684	enddo
685	endif
686	write(,) 'Global average pressure old time step',global_avg_press_old
687	write(,) 'Global average pressure new time step',global_avg_press_new
688
689	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
690	allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
691	write(,) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
692	do l = 1,nlayer + 1
693	do ig = 1,ngrid
694	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
695	enddo
696	enddo
697
698	! II.a.3. Surface pressure timeseries
699	write(,) "Recomputing the surface pressure timeserie adapted to the new pressure..."
700	do ig = 1,ngrid
701	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
702	enddo
703
704	! II.a.4. New pressure levels timeseries
705	allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
706	write(,) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
707	do l = 1,nlayer + 1
708	do ig = 1,ngrid
709	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
710	enddo
711	enddo
712
713	! II.a.5. Tracers timeseries
714	write(,) "Recomputing of tracer VMR timeseries for the new pressure..."
715
716	l = 1
717	do ig = 1,ngrid
718	do t = 1,timelen
719	q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
720	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
721	if (q_h2o_PEM_phys(ig,t) < 0) then
722	q_h2o_PEM_phys(ig,t) = 1.e-30
723	else if (q_h2o_PEM_phys(ig,t) > 1) then
724	q_h2o_PEM_phys(ig,t) = 1.
725	endif
726	enddo
727	enddo
728
729	do ig = 1,ngrid
730	do t = 1,timelen
731	q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
732	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
733	+ ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
734	- (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
735	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
736	if (q_co2_PEM_phys(ig,t) < 0) then
737	q_co2_PEM_phys(ig,t) = 1.e-30
738	else if (q_co2_PEM_phys(ig,t) > 1) then
739	q_co2_PEM_phys(ig,t) = 1.
740	endif
741	mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
742	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
743	enddo
744	enddo
745
746	deallocate(zplev_new_timeseries,zplev_old_timeseries)
747
748	!------------------------
749	! II Run
750	! II_b Evolution of ice
751	!------------------------
752	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
753	call evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
754	if (layering_algo) then
755	do islope = 1,nslope
756	do ig = 1,ngrid
757	call make_layering(stratif(ig,islope),tend_co2_ice(ig,islope),tend_h2o_ice(ig,islope),tend_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
758	enddo
759	enddo
760	endif
761
762	!------------------------
763	! II Run
764	! II_c Flow of glaciers
765	!------------------------
766	if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
767	global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
768	if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
769
770	!------------------------
771	! II Run
772	! II_d Update surface and soil temperatures
773	!------------------------
774	! II_d.1 Update Tsurf
775	write(,) "Updating the new Tsurf"
776	bool_sublim = .false.
777	Tsurfavg_before_saved = tsurf_ave
778	do ig = 1,ngrid
779	do islope = 1,nslope
780	if (co2_ice_ini(ig,islope) > 0.5 .and. co2_ice(ig,islope) < 1.e-10) then ! co2 ice disappeared, look for closest point without co2ice
781	if (latitude_deg(ig) > 0) then
782	do ig_loop = ig,ngrid
783	do islope_loop = islope,iflat,-1
784	if (co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
785	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
786	bool_sublim = .true.
787	exit
788	endif
789	enddo
790	if (bool_sublim) exit
791	enddo
792	else
793	do ig_loop = ig,1,-1
794	do islope_loop = islope,iflat
795	if(co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
796	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
797	bool_sublim = .true.
798	exit
799	endif
800	enddo
801	if (bool_sublim) exit
802	enddo
803	endif
804	co2_ice_ini(ig,islope) = 0
805	if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
806	albedo(ig,1,islope) = albedo_h2o_frost
807	albedo(ig,2,islope) = albedo_h2o_frost
808	else
809	albedo(ig,1,islope) = albedodat(ig)
810	albedo(ig,2,islope) = albedodat(ig)
811	emis(ig,islope) = emissiv
812	endif
813	else if ((co2_ice(ig,islope) > 1.e-3) .and. (tend_co2_ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
814	ave = 0.
815	do t = 1,timelen
816	if (co2_ice_PCM(ig,islope,t) > 1.e-3) then
817	ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
818	else
819	ave = ave + tsurf_PCM_timeseries(ig,islope,t)
820	endif
821	enddo
822	tsurf_ave(ig,islope) = ave/timelen
823	! set the surface albedo to be the ice albedo
824	if (latitude_deg(ig) > 0) then
825	icap = 1
826	else
827	icap = 2
828	endif
829	albedo(ig,1,islope) = albedice(icap)
830	albedo(ig,2,islope) = albedice(icap)
831	emis(ig,islope) = emisice(icap)
832	endif
833	enddo
834	enddo
835
836	do t = 1,timelen
837	tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_ave - Tsurfavg_before_saved
838	enddo
839	! for the start
840	do ig = 1,ngrid
841	do islope = 1,nslope
842	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_ave(ig,islope))
843	enddo
844	enddo
845
846	if (soil_pem) then
847
848	! II_d.2 Update soil temperature
849	allocate(TI_locslope(ngrid,nsoilmx_PEM))
850	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
851	allocate(Tsurf_locslope(ngrid))
852	write(,)"Updating soil temperature"
853
854	! Soil averaged
855	do islope = 1,nslope
856	TI_locslope = TI_PEM(:,:,islope)
857	do t = 1,timelen
858	Tsoil_locslope = tsoil_phys_PEM_timeseries(:,:,islope,t)
859	Tsurf_locslope = tsurf_PCM_timeseries(:,islope,t)
860	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
861	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
862	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope
863	do ig = 1,ngrid
864	do isoil = 1,nsoilmx_PEM
865	watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)mmol(igcm_h2o_vap)/(mugazr)
866	if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1)
867	enddo
868	enddo
869	enddo
870	enddo
871	tsoil_PEM = sum(tsoil_phys_PEM_timeseries,4)/timelen
872	watersoil_density_PEM_ave = sum(watersoil_density_PEM_timeseries,4)/timelen
873
874	write(,) "Update of soil temperature done"
875
876	deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
877
878	! II_d.3 Update the ice table
879	if (icetable_equilibrium) then
880	write(,) "Compute ice table"
881	porefillingice_thickness_prev_iter = porefillingice_thickness
882	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
883	call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_ave, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
884	endif
885	! II_d.4 Update the soil thermal properties
886	call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
887
888	! II_d.5 Update the mass of the regolith adsorbed
889	if (adsorption_pem) then
890	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice, &
891	h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
892	h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
893
894	totmassco2_adsorbded = 0.
895	totmassh2o_adsorbded = 0.
896	do ig = 1,ngrid
897	do islope =1, nslope
898	do l = 1,nsoilmx_PEM
899	if (l==1) then
900	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
901	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
902	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
903	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
904	else
905	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
906	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
907	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
908	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
909	endif
910	enddo
911	enddo
912	enddo
913	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
914	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
915	endif
916	endif !soil_pem
917
918	!------------------------
919	! II Run
920	! II_e Outputs
921	!------------------------
922	call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_avg_press_new/))
923	do islope = 1,nslope
924	write(str2(1:2),'(i2.2)') islope
925	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
926	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
927	call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
928	call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
929	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
930	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
931	if(icetable_equilibrium) then
932	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
933	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
934	endif
935	if(soil_pem) then
936	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
937	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
938	if (adsorption_pem) then
939	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
940	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
941	endif
942	endif
943	enddo
944
945	!------------------------
946	! II Run
947	! II_f Update the tendencies
948	!------------------------
949	call recomp_tend_co2_slope(ngrid,nslope,timelen,tend_co2_ice,tend_co2_ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries, &
950	global_avg_press_PCM,global_avg_press_new)
951
952	!------------------------
953	! II Run
954	! II_g Checking the stopping criterion
955	!------------------------
956	call stopping_crit_h2o_ice(cell_area,h2o_surf_ini,h2o_ice,stopPEM,ngrid)
957	call stopping_crit_co2(cell_area,co2_surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
958
959	year_iter = year_iter + dt_pem
960	i_myear = i_myear + dt_pem
961
962	write(,) "Checking the stopping criteria..."
963	if (year_iter >= year_iter_max) stopPEM = 5
964	if (i_myear >= n_myear) stopPEM = 6
965	if (stopPEM > 0) then
966	select case (stopPEM)
967	case(1)
968	write(,) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
969	case(2)
970	write(,) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
971	case(3)
972	write(,) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
973	case(4)
974	write(,) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
975	case(5)
976	write(,) "STOPPING because maximum number of iterations due to orbital parameters is reached:", stopPEM
977	case(6)
978	write(,) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
979	case default
980	write(,) "STOPPING with unknown stopping criterion code:", stopPEM
981	end select
982	exit
983	else
984	write(,) "The PEM can continue!"
985	write(,) "Number of iterations of the PEM: year_iter =", year_iter
986	write(,) "Number of simulated Martian years: i_myear =", i_myear
987	endif
988
989	global_avg_press_old = global_avg_press_new
990
991	enddo ! big time iteration loop of the pem
992	!------------------------------ END RUN --------------------------------
993
994	!------------------------------- OUTPUT --------------------------------
995	!------------------------
996	! III Output
997	! III_a Update surface value for the PCM start files
998	!------------------------
999	! III_a.1 Ice update (for startfi)
1000
1001	watercap = 0.
1002	perennial_co2ice = co2_ice
1003	do ig = 1,ngrid
1004	! H2O ice metamorphism
1005	if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_h2oice_threshold) then
1006	h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1007	qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
1008	endif
1009
1010	! Is H2O ice still considered as an infinite reservoir for the PCM?
1011	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
1012	! There is enough ice to be considered as an infinite reservoir
1013	watercaptag(ig) = .true.
1014	else
1015	! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
1016	watercaptag(ig) = .false.
1017	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1018	h2o_ice(ig,:) = 0.
1019	endif
1020
1021	! CO2 ice metamorphism
1022	if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_co2ice_threshold) then
1023	perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1024	qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
1025	endif
1026	enddo
1027
1028	! III_a.2 Tsoil update (for startfi)
1029	if (soil_pem) then
1030	call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
1031	tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom
1032	deallocate(tsoil_anom)
1033	#ifndef CPP_STD
1034	flux_geo = fluxgeo
1035	#endif
1036	endif
1037
1038	! III_a.4 Pressure (for start)
1039	ps = ps*global_avg_press_new/global_avg_press_PCM
1040	ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
1041
1042	! III_a.5 Tracer (for start)
1043	allocate(zplev_new(ngrid,nlayer + 1))
1044
1045	do l = 1,nlayer + 1
1046	zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
1047	enddo
1048
1049	do nnq = 1,nqtot
1050	if (noms(nnq) /= "co2") then
1051	do l = 1,llm - 1
1052	do ig = 1,ngrid
1053	q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
1054	enddo
1055	q(:,llm,nnq) = q(:,llm - 1,nnq)
1056	enddo
1057	else
1058	do l = 1,llm - 1
1059	do ig = 1,ngrid
1060	q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
1061	+ ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_PCM(ig,l) - zplev_PCM(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
1062	enddo
1063	q(:,llm,nnq) = q(:,llm - 1,nnq)
1064	enddo
1065	endif
1066	enddo
1067
1068	! Conserving the tracers mass for PCM start files
1069	do nnq = 1,nqtot
1070	do ig = 1,ngrid
1071	do l = 1,llm - 1
1072	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
1073	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1074	q(ig,l,nnq) = 1.
1075	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1076	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1077	endif
1078	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1079	enddo
1080	enddo
1081	enddo
1082
1083	if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
1084
1085	!------------------------
1086	! III Output
1087	! III_b Write "restart.nc" and "restartfi.nc"
1088	!------------------------
1089	! III_b.1 Write "restart.nc"
1090	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1091	pday = day_ini
1092	ztime_fin = time_phys
1093
1094	allocate(p(ip1jmp1,nlayer + 1))
1095	#ifndef CPP_1D
1096	call pression (ip1jmp1,ap,bp,ps,p)
1097	call massdair(p,masse)
1098	call dynredem0("restart.nc",day_ini,phis)
1099	call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps)
1100	write(,) "restart.nc has been written"
1101	#else
1102	call writerestart1D('restart1D_evol.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1103	write(,) "restart1D_evol.txt has been written"
1104	#endif
1105
1106	! III_b.2 Write the "restartfi.nc"
1107	#ifndef CPP_STD
1108	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1109	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1110	inertiedat,def_slope,subslope_dist)
1111	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1112	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1113	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1114	wstar,watercap,perennial_co2ice)
1115	#else
1116	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
1117	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
1118	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1119	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1120	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1121	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
1122	tsea_ice,sea_ice)
1123	#endif
1124	write(,) "restartfi.nc has been written"
1125
1126	!------------------------
1127	! III Output
1128	! III_c Write the "restartpem.nc"
1129	!------------------------
1130	if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
1131	call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1132	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
1133	TI_PEM,porefillingice_depth,porefillingice_thickness, &
1134	co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
1135	write(,) "restartpem.nc has been written"
1136
1137	call info_PEM(year_iter,stopPEM,i_myear,n_myear)
1138
1139	write(,) "The PEM has run for", year_iter, "Martian years."
1140	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1141	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1142	write(,) "LL & RV & JBC: so far, so good!"
1143
1144	if (layering_algo) then
1145	do islope = 1,nslope
1146	do i = 1,ngrid
1147	call del_layering(stratif(i,islope))
1148	enddo
1149	enddo
1150	deallocate(new_str,new_lag1,new_lag2,current1,current2)
1151	endif
1152	deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
1153	deallocate(co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfavg_before_saved)
1154	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
1155	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
1156	deallocate(co2_ice_ini,stratif)
1157	!----------------------------- END OUTPUT ------------------------------
1158
1159	END PROGRAM pem

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

Download in other formats:

Original Format