Context Navigation

pem.F90 @ 3171

Last change on this file since 3171 was 3171, checked in by llange, 11 months ago

PEM
Add the possibily to output the soil fields during a PEM run in a dedicated diagsoilpem

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

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

Download in other formats:

Original Format