Context Navigation

pem.F90 @ 3345

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

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

Download in other formats:

Original Format