Context Navigation

pem.F90 @ 3305

Last change on this file since 3305 was 3297, checked in by jbclement, 21 months ago

PEM:
Integration of the module "layering_mod.F90" with the rest of the PEM:

The linked list data structure representative of layered deposits is converted into an array which can be outputed in the "restratpem.nc" files. This array has dimensions (ngrid,nslope,nb_str_max,6) where 'nb_str_max' is the maximum number of 'stratum' through the layerings and '6' is the number of properties of 'stratum';
this structure can also be read from "startpem.nc" files to initialize PEM runs;
The layering algorithm is now used in the main PEM loop to make the layerings evolve.

JBC

File size: 57.1 KB

Line
1	!------------------------
2	! I Initialization
3	! I_a Read the "run.def"
4	! I_b Read the "start_evol.nc" and "startfi_evol.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_evol.nc" and "restartfi_evol.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
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_evol.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_evol.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_evol.nc',exist = therestartfi)
316	if (.not. therestartfi) then
317	write(,) 'There is no "startfi_evol.nc" file!'
318	error stop 'Initialization cannot be done for the 1D PEM.'
319	endif
320
321	call init_testphys1d('start1D_evol.txt','startfi_evol.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_evol.nc" and starfi_evol.nc
333	!------------------------
334	! I_b.1 Read "start_evol.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_evol.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_evol.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	do islope = 1,nslope
592	do i = 1,ngrid
593	call ini_layering(stratif(i,islope))
594	enddo
595	enddo
596
597	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
598	porefillingice_thickness,tsurf_avg_yr1,tsurf_ave,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries, &
599	tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM, &
600	watersurf_density_ave,watersoil_density_PEM_ave,co2_adsorbded_phys,delta_co2_adsorbded, &
601	h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
602
603	allocate(co2_ice_ini(ngrid,nslope))
604	co2_ice_ini = co2_ice
605
606	delta_h2o_icetablesublim = 0.
607
608	if (adsorption_pem) then
609	totmassco2_adsorbded = 0.
610	totmassh2o_adsorbded = 0.
611	do ig = 1,ngrid
612	do islope = 1,nslope
613	do l = 1,nsoilmx_PEM - 1
614	if (l==1) then
615	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
616	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
617	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
618	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
619	else
620	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
621	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
622	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_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	endif
625	enddo
626	enddo
627	enddo
628	write(,) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
629	write(,) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
630	endif ! adsorption
631	deallocate(tsurf_avg_yr1)
632
633	!------------------------
634	! I Initialization
635	! I_i Compute orbit criterion
636	!------------------------
637	#ifndef CPP_STD
638	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
639	#else
640	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
641	#endif
642
643	if (evol_orbit_pem) then
644	call orbit_param_criterion(i_myear,year_iter_max)
645	else
646	year_iter_max = Max_iter_pem
647	endif
648	!-------------------------- END INITIALIZATION -------------------------
649
650	!-------------------------------- RUN ----------------------------------
651	!------------------------
652	! II Run
653	! II_a Update pressure, ice and tracers
654	!------------------------
655	year_iter = 0
656	stopPEM = 0
657	allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
658	new_str = .true.
659	new_lag1 = .true.
660	new_lag2 = .true.
661	do islope = 1,nslope
662	do ig = 1,ngrid
663	current1(ig,islope)%p => stratif(ig,islope)%top
664	current2(ig,islope)%p => stratif(ig,islope)%top
665	enddo
666	enddo
667
668	do while (year_iter < 10 .and. i_myear < n_myear)
669	! II.a.1. Compute updated global pressure
670	write(,) "Recomputing the new pressure..."
671	do i = 1,ngrid
672	do islope = 1,nslope
673	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
674	enddo
675	enddo
676
677	if (adsorption_pem) then
678	do i = 1,ngrid
679	global_avg_press_new = global_avg_press_new - gcell_area(i)delta_co2_adsorbded(i)/Total_surface
680	enddo
681	endif
682	write(,) 'Global average pressure old time step',global_avg_press_old
683	write(,) 'Global average pressure new time step',global_avg_press_new
684
685	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
686	allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
687	write(,) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
688	do l = 1,nlayer + 1
689	do ig = 1,ngrid
690	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
691	enddo
692	enddo
693
694	! II.a.3. Surface pressure timeseries
695	write(,) "Recomputing the surface pressure timeserie adapted to the new pressure..."
696	do ig = 1,ngrid
697	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
698	enddo
699
700	! II.a.4. New pressure levels timeseries
701	allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
702	write(,) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
703	do l = 1,nlayer + 1
704	do ig = 1,ngrid
705	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
706	enddo
707	enddo
708
709	! II.a.5. Tracers timeseries
710	write(,) "Recomputing of tracer VMR timeseries for the new pressure..."
711
712	l = 1
713	do ig = 1,ngrid
714	do t = 1,timelen
715	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))/ &
716	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
717	if (q_h2o_PEM_phys(ig,t) < 0) then
718	q_h2o_PEM_phys(ig,t) = 1.e-30
719	else if (q_h2o_PEM_phys(ig,t) > 1) then
720	q_h2o_PEM_phys(ig,t) = 1.
721	endif
722	enddo
723	enddo
724
725	do ig = 1,ngrid
726	do t = 1,timelen
727	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))/ &
728	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
729	+ ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
730	- (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
731	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
732	if (q_co2_PEM_phys(ig,t) < 0) then
733	q_co2_PEM_phys(ig,t) = 1.e-30
734	else if (q_co2_PEM_phys(ig,t) > 1) then
735	q_co2_PEM_phys(ig,t) = 1.
736	endif
737	mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
738	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
739	enddo
740	enddo
741
742	deallocate(zplev_new_timeseries,zplev_old_timeseries)
743
744	!------------------------
745	! II Run
746	! II_b Evolution of ice
747	!------------------------
748	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
749	call evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
750	do islope = 1,nslope
751	do ig = 1,ngrid
752	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)
753	enddo
754	enddo
755
756	!------------------------
757	! II Run
758	! II_c Flow of glaciers
759	!------------------------
760	if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
761	global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
762	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)
763
764	!------------------------
765	! II Run
766	! II_d Update surface and soil temperatures
767	!------------------------
768	! II_d.1 Update Tsurf
769	write(,) "Updating the new Tsurf"
770	bool_sublim = .false.
771	Tsurfavg_before_saved = tsurf_ave
772	do ig = 1,ngrid
773	do islope = 1,nslope
774	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
775	if (latitude_deg(ig) > 0) then
776	do ig_loop = ig,ngrid
777	do islope_loop = islope,iflat,-1
778	if (co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
779	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
780	bool_sublim = .true.
781	exit
782	endif
783	enddo
784	if (bool_sublim) exit
785	enddo
786	else
787	do ig_loop = ig,1,-1
788	do islope_loop = islope,iflat
789	if(co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
790	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
791	bool_sublim = .true.
792	exit
793	endif
794	enddo
795	if (bool_sublim) exit
796	enddo
797	endif
798	co2_ice_ini(ig,islope) = 0
799	if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
800	albedo(ig,1,islope) = albedo_h2o_frost
801	albedo(ig,2,islope) = albedo_h2o_frost
802	else
803	albedo(ig,1,islope) = albedodat(ig)
804	albedo(ig,2,islope) = albedodat(ig)
805	emis(ig,islope) = emissiv
806	endif
807	else if ((co2_ice(ig,islope) > 1.e-3) .and. (tend_co2_ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
808	ave = 0.
809	do t = 1,timelen
810	if (co2_ice_PCM(ig,islope,t) > 1.e-3) then
811	ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
812	else
813	ave = ave + tsurf_PCM_timeseries(ig,islope,t)
814	endif
815	enddo
816	tsurf_ave(ig,islope) = ave/timelen
817	! set the surface albedo to be the ice albedo
818	if (latitude_deg(ig) > 0) then
819	icap = 1
820	else
821	icap = 2
822	endif
823	albedo(ig,1,islope) = albedice(icap)
824	albedo(ig,2,islope) = albedice(icap)
825	emis(ig,islope) = emisice(icap)
826	endif
827	enddo
828	enddo
829
830	do t = 1,timelen
831	tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_ave - Tsurfavg_before_saved
832	enddo
833	! for the start
834	do ig = 1,ngrid
835	do islope = 1,nslope
836	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_ave(ig,islope))
837	enddo
838	enddo
839
840	if (soil_pem) then
841
842	! II_d.2 Update soil temperature
843	allocate(TI_locslope(ngrid,nsoilmx_PEM))
844	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
845	allocate(Tsurf_locslope(ngrid))
846	write(,)"Updating soil temperature"
847
848	! Soil averaged
849	do islope = 1,nslope
850	TI_locslope = TI_PEM(:,:,islope)
851	do t = 1,timelen
852	Tsoil_locslope = tsoil_phys_PEM_timeseries(:,:,islope,t)
853	Tsurf_locslope = tsurf_PCM_timeseries(:,islope,t)
854	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
855	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
856	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope
857	do ig = 1,ngrid
858	do isoil = 1,nsoilmx_PEM
859	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)
860	if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1)
861	enddo
862	enddo
863	enddo
864	enddo
865	tsoil_PEM = sum(tsoil_phys_PEM_timeseries,4)/timelen
866	watersoil_density_PEM_ave = sum(watersoil_density_PEM_timeseries,4)/timelen
867
868	write(,) "Update of soil temperature done"
869
870	deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
871
872	! II_d.3 Update the ice table
873	if (icetable_equilibrium) then
874	write(,) "Compute ice table"
875	porefillingice_thickness_prev_iter = porefillingice_thickness
876	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
877	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
878	endif
879	! II_d.4 Update the soil thermal properties
880	call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
881
882	! II_d.5 Update the mass of the regolith adsorbed
883	if (adsorption_pem) then
884	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice, &
885	h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
886	h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
887
888	totmassco2_adsorbded = 0.
889	totmassh2o_adsorbded = 0.
890	do ig = 1,ngrid
891	do islope =1, nslope
892	do l = 1,nsoilmx_PEM
893	if (l==1) then
894	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
895	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
896	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l)) &
897	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
898	else
899	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
900	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
901	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
902	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
903	endif
904	enddo
905	enddo
906	enddo
907	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
908	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
909	endif
910	endif !soil_pem
911
912	!------------------------
913	! II Run
914	! II_e Outputs
915	!------------------------
916	call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_avg_press_new/))
917	do islope = 1,nslope
918	write(str2(1:2),'(i2.2)') islope
919	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
920	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
921	call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
922	call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
923	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
924	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
925	if(icetable_equilibrium) then
926	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
927	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
928	endif
929	if(soil_pem) then
930	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
931	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
932	if (adsorption_pem) then
933	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
934	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
935	endif
936	endif
937	enddo
938
939	!------------------------
940	! II Run
941	! II_f Update the tendencies
942	!------------------------
943	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, &
944	global_avg_press_PCM,global_avg_press_new)
945
946	!------------------------
947	! II Run
948	! II_g Checking the stopping criterion
949	!------------------------
950	call stopping_crit_h2o_ice(cell_area,h2o_surf_ini,h2o_ice,stopPEM,ngrid)
951	call stopping_crit_co2(cell_area,co2_surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
952
953	year_iter = year_iter + dt_pem
954	i_myear = i_myear + dt_pem
955
956	write(,) "Checking the stopping criteria..."
957	if (year_iter >= year_iter_max) stopPEM = 5
958	if (i_myear >= n_myear) stopPEM = 6
959	if (stopPEM > 0) then
960	select case (stopPEM)
961	case(1)
962	write(,) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
963	case(2)
964	write(,) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
965	case(3)
966	write(,) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
967	case(4)
968	write(,) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
969	case(5)
970	write(,) "STOPPING because maximum number of iterations due to orbital parameters is reached:", stopPEM
971	case(6)
972	write(,) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
973	case default
974	write(,) "STOPPING with unknown stopping criterion code:", stopPEM
975	end select
976	exit
977	else
978	write(,) "The PEM can continue!"
979	write(,) "Number of iterations of the PEM: year_iter =", year_iter
980	write(,) "Number of simulated Martian years: i_myear =", i_myear
981	endif
982
983	global_avg_press_old = global_avg_press_new
984
985	enddo ! big time iteration loop of the pem
986	!------------------------------ END RUN --------------------------------
987
988	!------------------------------- OUTPUT --------------------------------
989	!------------------------
990	! III Output
991	! III_a Update surface value for the PCM start files
992	!------------------------
993	! III_a.1 Ice update (for startfi)
994
995	watercap = 0.
996	perennial_co2ice = co2_ice
997	do ig = 1,ngrid
998	! H2O ice metamorphism
999	if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_h2oice_threshold) then
1000	h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
1001	qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
1002	endif
1003
1004	! Is H2O ice still considered as an infinite reservoir for the PCM?
1005	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
1006	! There is enough ice to be considered as an infinite reservoir
1007	watercaptag(ig) = .true.
1008	else
1009	! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
1010	watercaptag(ig) = .false.
1011	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1012	h2o_ice(ig,:) = 0.
1013	endif
1014
1015	! CO2 ice metamorphism
1016	if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_co2ice_threshold) then
1017	perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
1018	qsurf(ig,igcm_co2,:) = metam_co2ice_threshold/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
1019	endif
1020	enddo
1021
1022	! III_a.2 Tsoil update (for startfi)
1023	if (soil_pem) then
1024	call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
1025	tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom
1026	deallocate(tsoil_anom)
1027	#ifndef CPP_STD
1028	flux_geo = fluxgeo
1029	#endif
1030	endif
1031
1032	! III_a.4 Pressure (for start)
1033	ps = ps*global_avg_press_new/global_avg_press_PCM
1034	ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
1035
1036	! III_a.5 Tracer (for start)
1037	allocate(zplev_new(ngrid,nlayer + 1))
1038
1039	do l = 1,nlayer + 1
1040	zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
1041	enddo
1042
1043	do nnq = 1,nqtot
1044	if (noms(nnq) /= "co2") then
1045	do l = 1,llm - 1
1046	do ig = 1,ngrid
1047	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))
1048	enddo
1049	q(:,llm,nnq) = q(:,llm - 1,nnq)
1050	enddo
1051	else
1052	do l = 1,llm - 1
1053	do ig = 1,ngrid
1054	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)) &
1055	+ ((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))
1056	enddo
1057	q(:,llm,nnq) = q(:,llm - 1,nnq)
1058	enddo
1059	endif
1060	enddo
1061
1062	! Conserving the tracers mass for PCM start files
1063	do nnq = 1,nqtot
1064	do ig = 1,ngrid
1065	do l = 1,llm - 1
1066	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
1067	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1068	q(ig,l,nnq) = 1.
1069	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1070	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1071	endif
1072	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1073	enddo
1074	enddo
1075	enddo
1076
1077	if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
1078
1079	!------------------------
1080	! III Output
1081	! III_b Write "restart_evol.nc" and "restartfi_evol.nc"
1082	!------------------------
1083	! III_b.1 Write "restart_evol.nc"
1084	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1085	pday = day_ini
1086	ztime_fin = time_phys
1087
1088	allocate(p(ip1jmp1,nlayer + 1))
1089	#ifndef CPP_1D
1090	call pression (ip1jmp1,ap,bp,ps,p)
1091	call massdair(p,masse)
1092	call dynredem0("restart_evol.nc",day_ini,phis)
1093	call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
1094	write(,) "restart_evol.nc has been written"
1095	#else
1096	call writerestart1D('restart1D_evol.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1097	write(,) "restart1D_evol.txt has been written"
1098	#endif
1099
1100	! III_b.2 Write the "restartfi_evol.nc"
1101	#ifndef CPP_STD
1102	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
1103	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1104	inertiedat,def_slope,subslope_dist)
1105	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1106	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1107	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1108	wstar,watercap,perennial_co2ice)
1109	#else
1110	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
1111	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
1112	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1113	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1114	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1115	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
1116	tsea_ice,sea_ice)
1117	#endif
1118	write(,) "restartfi_evol.nc has been written"
1119
1120	!------------------------
1121	! III Output
1122	! III_c Write the "restartpem.nc"
1123	!------------------------
1124	nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
1125	call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
1126	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
1127	TI_PEM,porefillingice_depth,porefillingice_thickness, &
1128	co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
1129	write(,) "restartpem.nc has been written"
1130
1131	call info_PEM(year_iter,stopPEM,i_myear,n_myear)
1132
1133	write(,) "The PEM has run for", year_iter, "Martian years."
1134	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1135	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1136	write(,) "LL & RV & JBC: so far, so good!"
1137
1138	do islope = 1,nslope
1139	do i = 1,ngrid
1140	call del_layering(stratif(i,islope))
1141	enddo
1142	enddo
1143	deallocate(stratif,new_str,new_lag1,new_lag2,current1,current2)
1144	deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
1145	deallocate(co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfavg_before_saved)
1146	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
1147	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
1148	deallocate(co2_ice_ini)
1149	!----------------------------- END OUTPUT ------------------------------
1150
1151	END PROGRAM pem

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

Download in other formats:

Original Format