Context Navigation

pem.F90 @ 3071

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

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

Download in other formats:

Original Format