Context Navigation

pem.F90 @ 3093

Last change on this file since 3093 was 3088, checked in by jbclement, 2 years ago

PEM:

Addition of the bash script "concat_diagpem.sh" in deftank to concatenate all the "diagpem.nc" outputs into one file;
Addition of "writediagpem_mod.F90" inspired by the subroutine "writediagfi.F" in phymars/;
Renaming of PEM outputs ('startfi_PEM' -> 'startpem' and 'diagfi' -> 'diagpem');
Gathering of all "diagpem.nc" outputs at the same place at the end of the main PEM loop;
Update of "launch_pem.sh", "README" and "modify_startfi_orbit.sh" in the deftank.

JBC

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