Context Navigation

pem.F90 @ 3069

Last change on this file since 3069 was 3069, checked in by jbclement, 2 years ago
Mars PCM: Related to commit r3066, correction of a bug to write/read a restart/start in 1D and more adaptations of the code. JBC
File size: 57.3 KB

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

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

Download in other formats:

Original Format