Context Navigation

pem.F90 @ 3116

Last change on this file since 3116 was 3114, checked in by jbclement, 21 months ago
PEM: Following r3113, addition of 'nqsoil' and 'qsoil' in the arguments of the subroutines 'phyetat0' and 'physdem1' to be able to compile. JBC
File size: 57.4 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 PCM 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 PCM
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_PEM_mod, only: info_PEM
66	use interpolate_TIPEM_TIGCM_mod, only: interpolate_TIPEM_TIGCM
67	use nb_time_step_PCM_mod, only: nb_time_step_PCM
68	use pemetat0_mod, only: pemetat0
69	use read_data_PCM_mod, only: read_data_PCM
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, nqsoil, qsoil
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, iniorbit
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 PCM
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 PCM
129	real :: ptimestep ! Same as PCM
130	real :: ztime_fin ! Same as PCM
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 PCM [kg/m^2]
169	real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
170	real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
171	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 PCM [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 PCM [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 PCM [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 PCM [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 PCM
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 PCM [K]
210	real, dimension(:,:), allocatable :: TI_locslope ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI]
211	real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: intermediate when computing Tsoil [K]
212	real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
213	real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
214	real, dimension(:,:), allocatable :: watersurf_density_ave ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
215	real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
216	real, dimension(:,:,:), allocatable :: watersoil_density_PEM_ave ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
217	real, dimension(:,:), allocatable :: 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 PCM, 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,nqsoil,day_ini,time_phys,tsurf, &
370	tsoil,albedo,emis,q2,qsurf,qsoil,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(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
386	allocate(emis_read_generic(ngrid))
387	allocate(tsurf(ngrid,1))
388	allocate(qsurf(ngrid,nq,1))
389	allocate(tsoil(ngrid,nsoilmx,1))
390	allocate(emis(ngrid,1))
391	allocate(watercap(ngrid,1))
392	allocate(watercaptag(ngrid))
393	allocate(albedo_read_generic(ngrid,2))
394	allocate(albedo(ngrid,2,1))
395	allocate(inertiesoil(ngrid,nsoilmx,1))
396	call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
397	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
398	qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, &
399	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
400	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
401
402	nslope = 1
403	call ini_comslope_h(ngrid,1)
404
405	qsurf(:,:,1) = qsurf_read_generic(:,:)
406	tsurf(:,1) = tsurf_read_generic(:)
407	tsoil(:,:,1) = tsoil_read_generic(:,:)
408	emis(:,1) = emis_read_generic(:)
409	watercap(:,1) = 0.
410	watercaptag(:) = .false.
411	albedo(:,1,1) = albedo_read_generic(:,1)
412	albedo(:,2,1) = albedo_read_generic(:,2)
413	inertiesoil(:,:,1) = inertiedat(:,:)
414
415	if (nslope == 1) then
416	def_slope(1) = 0
417	def_slope(2) = 0
418	def_slope_mean = 0
419	subslope_dist(:,1) = 1.
420	endif
421
422	! Remove unphysical values of surface tracer
423	qsurf(:,:,1) = qsurf_read_generic(:,:)
424	where (qsurf < 0.) qsurf = 0.
425	#endif
426
427	do nnq = 1,nqtot ! Why not using ini_tracer ?
428	if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
429	if (noms(nnq) == "h2o_vap") then
430	igcm_h2o_vap = nnq
431	mmol(igcm_h2o_vap)=18.
432	endif
433	if (noms(nnq) == "co2") igcm_co2 = nnq
434	enddo
435	r = 8.314511*1000./mugaz
436
437	!------------------------
438	! I Initialization
439	! I_c Subslope parametrisation
440	!------------------------
441	! Define some slope statistics
442	iflat = 1
443	do islope = 2,nslope
444	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
445	enddo
446
447	write(,) 'Flat slope for islope = ',iflat
448	write(,) 'corresponding criterium = ',def_slope_mean(iflat)
449
450	allocate(flag_co2flow(ngrid,nslope))
451	allocate(flag_co2flow_mesh(ngrid))
452	allocate(flag_h2oflow(ngrid,nslope))
453	allocate(flag_h2oflow_mesh(ngrid))
454
455	flag_co2flow(:,:) = 0
456	flag_co2flow_mesh(:) = 0
457	flag_h2oflow(:,:) = 0
458	flag_h2oflow_mesh(:) = 0
459
460	!------------------------
461	! I Initialization
462	! I_d READ PCM data and convert to the physical grid
463	!------------------------
464	! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value
465	call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
466
467	allocate(tsoil_ave(ngrid,nsoilmx,nslope))
468	allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
469	allocate(vmr_co2_gcm(ngrid,timelen))
470	allocate(ps_timeseries(ngrid,timelen))
471	allocate(min_co2_ice_1(ngrid,nslope))
472	allocate(min_h2o_ice_1(ngrid,nslope))
473	allocate(min_co2_ice_2(ngrid,nslope))
474	allocate(min_h2o_ice_2(ngrid,nslope))
475	allocate(tsurf_ave_yr1(ngrid,nslope))
476	allocate(tsurf_ave(ngrid,nslope))
477	allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
478	allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
479	allocate(q_co2_PEM_phys(ngrid,timelen))
480	allocate(q_h2o_PEM_phys(ngrid,timelen))
481	allocate(co2_ice_GCM(ngrid,nslope,timelen))
482	allocate(watersurf_density_ave(ngrid,nslope))
483	allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
484	allocate(Tsurfave_before_saved(ngrid,nslope))
485	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
486	allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
487	allocate(delta_co2_adsorbded(ngrid))
488	allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
489	allocate(delta_h2o_icetablesublim(ngrid))
490	allocate(delta_h2o_adsorbded(ngrid))
491	allocate(vmr_co2_pem_phys(ngrid,timelen))
492
493	write(,) "Downloading data Y1..."
494	call read_data_PCM("data_PCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
495	tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
496	co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
497	write(,) "Downloading data Y1 done"
498
499	! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value
500	write(,) "Downloading data Y2"
501	call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
502	tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
503	co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
504	write(,) "Downloading data Y2 done"
505
506	!------------------------
507	! I Initialization
508	! I_e Initialization of the PEM variables and soil
509	!------------------------
510	call end_comsoil_h_PEM
511	call ini_comsoil_h_PEM(ngrid,nslope)
512	call end_adsorption_h_PEM
513	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
514	call end_ice_table_porefilling
515	call ini_ice_table_porefilling(ngrid,nslope)
516
517	if (soil_pem) then
518	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
519	do l = 1,nsoilmx
520	tsoil_PEM(:,l,:) = tsoil_ave(:,l,:)
521	tsoil_phys_PEM_timeseries(:,l,:,:) = tsoil_GCM_timeseries(:,l,:,:)
522	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,l,:,:)
523	enddo
524	do l = nsoilmx + 1,nsoilmx_PEM
525	tsoil_PEM(:,l,:) = tsoil_ave(:,nsoilmx,:)
526	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
527	enddo
528	watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
529	endif !soil_pem
530	deallocate(tsoil_ave,tsoil_GCM_timeseries)
531
532	!------------------------
533	! I Initialization
534	! I_f Compute tendencies & Save initial situation
535	!------------------------
536	allocate(tendencies_co2_ice(ngrid,nslope))
537	allocate(tendencies_co2_ice_ini(ngrid,nslope))
538	allocate(tendencies_h2o_ice(ngrid,nslope))
539
540	! Compute the tendencies of the evolution of ice over the years
541	call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice)
542	tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:)
543	call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice)
544
545	deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2)
546
547	!------------------------
548	! I Initialization
549	! I_g Save initial PCM situation
550	!------------------------
551	allocate(initial_co2_ice_sublim(ngrid,nslope))
552	allocate(initial_co2_ice(ngrid,nslope))
553	allocate(initial_h2o_ice(ngrid,nslope))
554
555	! We save the places where water ice is sublimating
556	! We compute the surface of water ice sublimating
557	ini_surf_co2 = 0.
558	ini_surf_h2o = 0.
559	Total_surface = 0.
560	do i = 1,ngrid
561	Total_surface = Total_surface + cell_area(i)
562	do islope = 1,nslope
563	if (tendencies_co2_ice(i,islope) < 0) then
564	initial_co2_ice_sublim(i,islope) = 1.
565	ini_surf_co2 = ini_surf_co2 + cell_area(i)*subslope_dist(i,islope)
566	else
567	initial_co2_ice_sublim(i,islope) = 0.
568	endif
569	if (qsurf(i,igcm_co2,islope) > 0) then
570	initial_co2_ice(i,islope) = 1.
571	else
572	initial_co2_ice(i,islope) = 0.
573	endif
574	if (tendencies_h2o_ice(i,islope) < 0) then
575	initial_h2o_ice(i,islope) = 1.
576	ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
577	else
578	initial_h2o_ice(i,islope) = 0.
579	endif
580	enddo
581	enddo
582
583	write(,) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
584	write(,) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
585	write(,) "Total surface of the planet=", Total_surface
586	allocate(zplev_gcm(ngrid,nlayer + 1))
587
588	do ig = 1,ngrid
589	zplev_gcm(ig,:) = ap(:) + bp(:)*ps_start_GCM(ig)
590	enddo
591
592	global_ave_press_old = 0.
593	do i = 1,ngrid
594	global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
595	enddo
596
597	global_ave_press_GCM = global_ave_press_old
598	global_ave_press_new = global_ave_press_old
599	write(,) "Initial global average pressure=", global_ave_press_GCM
600
601	!------------------------
602	! I Initialization
603	! I_h Read the startpem.nc
604	!------------------------
605	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
606	porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries, &
607	tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:), &
608	qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave, &
609	co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
610
611	delta_h2o_icetablesublim(:) = 0.
612
613	do ig = 1,ngrid
614	do islope = 1,nslope
615	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
616	enddo
617	enddo
618
619	if (adsorption_pem) then
620	totmassco2_adsorbded = 0.
621	totmassh2o_adsorbded = 0.
622	do ig = 1,ngrid
623	do islope = 1,nslope
624	do l = 1,nsoilmx_PEM - 1
625	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
626	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
627	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
628	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
629	enddo
630	enddo
631	enddo
632
633	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
634	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
635	endif ! adsorption
636	deallocate(tsurf_ave_yr1)
637
638	!------------------------
639	! I Initialization
640	! I_i Compute orbit criterion
641	!------------------------
642	#ifndef CPP_STD
643	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
644	#else
645	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
646	#endif
647
648	if (evol_orbit_pem) then
649	call orbit_param_criterion(i_myear,year_iter_max)
650	else
651	year_iter_max = Max_iter_pem
652	endif
653	!-------------------------- END INITIALIZATION -------------------------
654
655	!-------------------------------- RUN ----------------------------------
656	!------------------------
657	! II Run
658	! II_a Update pressure, ice and tracers
659	!------------------------
660	year_iter = 0
661
662	do while (year_iter < year_iter_max .and. i_myear < n_myear)
663	! II.a.1. Compute updated global pressure
664	write(,) "Recomputing the new pressure..."
665	do i = 1,ngrid
666	do islope = 1,nslope
667	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
668	enddo
669	enddo
670
671	if (adsorption_pem) then
672	do i = 1,ngrid
673	global_ave_press_new = global_ave_press_new - gcell_area(i)delta_co2_adsorbded(i)/Total_surface
674	enddo
675	endif
676	write(,) 'Global average pressure old time step',global_ave_press_old
677	write(,) 'Global average pressure new time step',global_ave_press_new
678
679	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
680	allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
681	write(,) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
682	do l = 1,nlayer + 1
683	do ig = 1,ngrid
684	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
685	enddo
686	enddo
687
688	! II.a.3. Surface pressure timeseries
689	write(,) "Recomputing the surface pressure timeserie adapted to the new pressure..."
690	do ig = 1,ngrid
691	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
692	enddo
693
694	! II.a.4. New pressure levels timeseries
695	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
696	write(,) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
697	do l = 1,nlayer + 1
698	do ig = 1,ngrid
699	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
700	enddo
701	enddo
702
703	! II.a.5. Tracers timeseries
704	write(,) "Recomputing of tracer VMR timeseries for the new pressure..."
705
706	l = 1
707	do ig = 1,ngrid
708	do t = 1,timelen
709	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))/ &
710	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
711	if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30
712	if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1.
713	enddo
714	enddo
715
716	do ig = 1,ngrid
717	do t = 1,timelen
718	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))/ &
719	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
720	+ ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
721	- (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
722	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
723	if (q_co2_PEM_phys(ig,t) < 0) then
724	q_co2_PEM_phys(ig,t) = 1.e-30
725	elseif (q_co2_PEM_phys(ig,t) > 1) then
726	q_co2_PEM_phys(ig,t) = 1.
727	endif
728	mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
729	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
730	enddo
731	enddo
732
733	deallocate(zplev_new_timeseries,zplev_old_timeseries)
734
735	!------------------------
736	! II Run
737	! II_b Evolution of the ice
738	!------------------------
739	write(,) "Evolution of h2o ice"
740	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)
741
742	write(,) "Evolution of co2 ice"
743	call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
744
745	!------------------------
746	! II Run
747	! II_c CO2 & H2O glaciers flows
748	!------------------------
749	write(,) "CO2 glacier flows"
750	if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
751	global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
752
753	write(,) "H2O glacier flows"
754	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)
755
756	!------------------------
757	! II Run
758	! II_d Update surface and soil temperatures
759	!------------------------
760	! II_d.1 Update Tsurf
761	write(,) "Updating the new Tsurf"
762	bool_sublim = .false.
763	Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
764	do ig = 1,ngrid
765	do islope = 1,nslope
766	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
767	if (latitude_deg(ig) > 0) then
768	do ig_loop = ig,ngrid
769	do islope_loop = islope,iflat,-1
770	if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
771	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
772	bool_sublim = .true.
773	exit
774	endif
775	enddo
776	if (bool_sublim) exit
777	enddo
778	else
779	do ig_loop = ig,1,-1
780	do islope_loop = islope,iflat
781	if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
782	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
783	bool_sublim = .true.
784	exit
785	endif
786	enddo
787	if (bool_sublim) exit
788	enddo
789	endif
790	initial_co2_ice(ig,islope) = 0
791	if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
792	albedo(ig,1,islope) = albedo_h2o_frost
793	albedo(ig,2,islope) = albedo_h2o_frost
794	else
795	albedo(ig,1,islope) = albedodat(ig)
796	albedo(ig,2,islope) = albedodat(ig)
797	emis(ig,islope) = emissiv
798	endif
799	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
800	ave = 0.
801	do t = 1,timelen
802	if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
803	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
804	else
805	ave = ave + tsurf_GCM_timeseries(ig,islope,t)
806	endif
807	enddo
808	tsurf_ave(ig,islope) = ave/timelen
809	! set the surface albedo to be the ice albedo
810	if (latitude_deg(ig) > 0) then
811	icap = 1
812	else
813	icap = 2
814	endif
815	albedo(ig,1,islope) = albedice(icap)
816	albedo(ig,2,islope) = albedice(icap)
817	emis(ig,islope) = emisice(icap)
818	endif
819	enddo
820	enddo
821
822	do t = 1,timelen
823	tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:))
824	enddo
825	! for the start
826	do ig = 1,ngrid
827	do islope = 1,nslope
828	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
829	enddo
830	enddo
831
832	if (soil_pem) then
833
834	! II_d.2 Update soil temperature
835	allocate(TI_locslope(ngrid,nsoilmx_PEM))
836	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
837	allocate(Tsurf_locslope(ngrid))
838	write(,)"Updating soil temperature"
839
840	! Soil averaged
841	do islope = 1,nslope
842	TI_locslope(:,:) = TI_PEM(:,:,islope)
843	do t = 1,timelen
844	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
845	Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
846	call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
847	call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
848	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
849	do ig = 1,ngrid
850	do isoil = 1,nsoilmx_PEM
851	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)
852	if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1)
853	enddo
854	enddo
855	enddo
856	enddo
857	tsoil_PEM(:,:,:) = sum(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
858	watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
859
860	write(,) "Update of soil temperature done"
861
862	deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
863	write(,) "Compute ice table"
864
865	! II_d.3 Update the ice table
866	porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:)
867	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
868
869	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
870
871	write(,) "Update soil propreties"
872
873	! II_d.4 Update the soil thermal properties
874	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)
875
876	! II_d.5 Update the mass of the regolith adsorbded
877	if (adsorption_pem) then
878	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
879	qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, &
880	q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
881
882	totmassco2_adsorbded = 0.
883	totmassh2o_adsorbded = 0.
884	do ig = 1,ngrid
885	do islope =1, nslope
886	do l = 1,nsoilmx_PEM - 1
887	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
888	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
889	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
890	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
891	enddo
892	enddo
893	enddo
894	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
895	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
896	endif
897	endif !soil_pem
898
899	!------------------------
900	! II Run
901	! II_e Outputs
902	!------------------------
903	call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_ave_press_new/))
904	do islope = 1,nslope
905	write(str2(1:2),'(i2.2)') islope
906	call writediagpem(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
907	call writediagpem(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
908	call writediagpem(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
909	call writediagpem(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
910	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
911	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
912	enddo
913
914	!------------------------
915	! II Run
916	! II_f Update the tendencies
917	!------------------------
918	write(,) "Adaptation of the new co2 tendencies to the current pressure"
919	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, &
920	global_ave_press_GCM,global_ave_press_new)
921
922	!------------------------
923	! II Run
924	! II_g Checking the stopping criterion
925	!------------------------
926	call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
927
928	call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
929	initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
930
931	year_iter = year_iter + dt_pem
932	i_myear = i_myear + dt_pem
933
934	write(,) "Checking all the stopping criterion."
935	if (STOPPING_water) then
936	write(,) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
937	criterion_stop = 1
938	endif
939	if (STOPPING_1_water) then
940	write(,) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
941	criterion_stop = 1
942	endif
943	if (STOPPING_co2) then
944	write(,) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
945	criterion_stop = 2
946	endif
947	if (STOPPING_pressure) then
948	write(,) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
949	criterion_stop = 3
950	endif
951	if (year_iter >= year_iter_max) then
952	write(,) "STOPPING because maximum number of iterations reached"
953	criterion_stop = 4
954	endif
955	if (i_myear >= n_myear) then
956	write(,) "STOPPING because maximum number of Martian years to be simulated reached"
957	criterion_stop = 5
958	endif
959
960	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then
961	exit
962	else
963	write(,) "We continue!"
964	write(,) "Number of iterations of the PEM: year_iter =", year_iter
965	write(,) "Number of simulated Martian years: i_myear =", i_myear
966	endif
967
968	global_ave_press_old = global_ave_press_new
969
970	enddo ! big time iteration loop of the pem
971	!------------------------------ END RUN --------------------------------
972
973	!------------------------------- OUTPUT --------------------------------
974	!------------------------
975	! III Output
976	! III_a Update surface value for the PCM start files
977	!------------------------
978	! III_a.1 Ice update (for startfi)
979
980	! H2O ice
981	do ig = 1,ngrid
982	if (watercaptag(ig)) then
983	watercap_sum = 0.
984	do islope = 1,nslope
985	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
986	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
987	else
988	! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
989	watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.))
990	qsurf(ig,igcm_h2o_ice,islope)=0.
991	endif
992	watercap_sum = watercap_sum+watercap(ig,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
993	watercap(ig,islope) = 0.
994	enddo
995	water_reservoir(ig) = water_reservoir(ig) + watercap_sum
996	endif
997	enddo
998
999	do ig = 1,ngrid
1000	water_sum = 0.
1001	do islope = 1,nslope
1002	water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
1003	enddo
1004	if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
1005	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
1006	watercaptag(ig) = .true.
1007	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
1008	do islope = 1,nslope
1009	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.cos(pidef_slope_mean(islope)/180.)
1010	enddo
1011	endif
1012	else ! let's check that the infinite source of water disapear
1013	if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then
1014	watercaptag(ig) = .false.
1015	do islope = 1,nslope
1016	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
1017	enddo
1018	water_reservoir(ig) = 0.
1019	endif
1020	endif
1021	enddo
1022
1023	! CO2 ice
1024	do ig = 1,ngrid
1025	do islope = 1,nslope
1026	if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then
1027	perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
1028	qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope)
1029	albedo(ig,1,islope) = albedo_perenialco2
1030	albedo(ig,2,islope) = albedo_perenialco2
1031	endif
1032	enddo
1033	enddo
1034
1035	! III_a.2 Tsoil update (for startfi)
1036	if (soil_pem) then
1037	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
1038	tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
1039	endif
1040
1041	! III_a.4 Pressure (for start)
1042	ps(:) = ps(:)*global_ave_press_new/global_ave_press_GCM
1043	ps_start_GCM(:) = ps_start_GCM(:)*global_ave_press_new/global_ave_press_GCM
1044
1045	! III_a.5 Tracer (for start)
1046	allocate(zplev_new(ngrid,nlayer + 1))
1047
1048	do l = 1,nlayer + 1
1049	zplev_new(:,l) = ap(l) + bp(l)*ps_start_GCM(:)
1050	enddo
1051
1052	do nnq = 1,nqtot
1053	if (noms(nnq) /= "co2") then
1054	do l = 1,llm - 1
1055	do ig = 1,ngrid
1056	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))
1057	enddo
1058	q(:,llm,nnq) = q(:,llm - 1,nnq)
1059	enddo
1060	else
1061	do l = 1,llm - 1
1062	do ig = 1,ngrid
1063	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)) &
1064	+ ((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))
1065	enddo
1066	q(:,llm,nnq) = q(:,llm - 1,nnq)
1067	enddo
1068	endif
1069	enddo
1070
1071	! Conserving the tracers mass for PCM start files
1072	do nnq = 1,nqtot
1073	do ig = 1,ngrid
1074	do l = 1,llm - 1
1075	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
1076	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1077	q(ig,l,nnq) = 1.
1078	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1079	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1080	endif
1081	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1082	enddo
1083	enddo
1084	enddo
1085
1086	if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
1087
1088	!------------------------
1089	! III Output
1090	! III_b Write restart_evol.nc and restartfi_evol.nc
1091	!------------------------
1092	! III_b.1 Write restart_evol.nc
1093	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1094	pday = day_ini
1095	ztime_fin = time_phys
1096
1097	allocate(p(ip1jmp1,nlayer + 1))
1098	#ifndef CPP_1D
1099	call pression (ip1jmp1,ap,bp,ps,p)
1100	call massdair(p,masse)
1101	call dynredem0("restart_evol.nc",day_ini,phis)
1102	call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
1103	write(,) "restart_evol.nc has been written"
1104	#else
1105	call writerestart1D('restart1D_evol.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1106	write(,) "restart1D_evol.txt has been written"
1107	#endif
1108
1109	! III_b.2 Write restartfi_evol.nc
1110	#ifndef CPP_STD
1111	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
1112	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1113	inertiedat,def_slope,subslope_dist)
1114
1115	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1116	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1117	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1118	wstar,watercap,perenial_co2ice)
1119	#else
1120	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
1121	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
1122	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1123
1124	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1125	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1126	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
1127	tsea_ice,sea_ice)
1128	#endif
1129	write(,) "restartfi_evol.nc has been written"
1130
1131	!------------------------
1132	! III Output
1133	! III_c Write restartpem.nc
1134	!------------------------
1135	call pemdem0("restartpem.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
1136	float(day_ini),0.,nslope,def_slope,subslope_dist)
1137	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
1138	TI_PEM, porefillingice_depth,porefillingice_thickness, &
1139	co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
1140	write(,) "restartpem.nc has been written"
1141
1142	call info_PEM(year_iter,criterion_stop,i_myear,n_myear)
1143
1144	write(,) "The PEM has run for", year_iter, "Martian years."
1145	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1146	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1147	write(,) "LL & RV & JBC: so far, so good!"
1148
1149	deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
1150	deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
1151	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
1152	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
1153	!----------------------------- END OUTPUT ------------------------------
1154
1155	END PROGRAM pem

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

Download in other formats:

Original Format