Context Navigation

pem.F90 @ 3296

Last change on this file since 3296 was 3287, checked in by jbclement, 8 months ago

PEM:

Small correction to make the 3D PEM be able to compile.
Improvement of "launch_pem.sh": a file "kill_launch_pem.sh" is now automatically created which allows the user to kill the process of the launching script in case.

JBC

File size: 55.5 KB

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

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

Download in other formats:

Original Format