Context Navigation

pem.F90 @ 3151

Last change on this file since 3151 was 3149, checked in by jbclement, 2 years ago

PEM:

Simplification of the algorithm managing the stopping criteria;
Complete rework of the ice management in the PEM (H2O & CO2);

Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').
Renaming of variables/subroutines for clarity;
Some cleanings throughout the code;
Small updates in files of the deftank.

JBC

File size: 52.6 KB

Line
1	!------------------------
2	! I Initialization
3	! I_a Read run.def
4	! I_b Read of start_evol.nc and starfi_evol.nc
5	! I_c Subslope parametrisation
6	! I_d Read PCM data and convert to the physical grid
7	! I_e Initialization of the PEM variable and soil
8	! I_f Compute tendencies
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 restart_evol.nc and restartfi_evol.nc
25	! III_c Write restartpem.nc
26	!------------------------
27
28	PROGRAM pem
29
30	use phyetat0_mod, only: phyetat0
31	use phyredem, only: physdem0, physdem1
32	use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
33	use turb_mod, only: q2, wstar
34	use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
35	use logic_mod, only: iflag_phys
36	use mod_const_mpi, only: COMM_LMDZ
37	use infotrac
38	use geometry_mod, only: latitude_deg
39	use conf_pem_mod, only: conf_pem
40	use pemredem, only: pemdem0, pemdem1
41	use glaciers_mod, only: flow_co2glaciers, flow_h2oglaciers, co2_ice_flow, h2o_ice_flow
42	use stopping_crit_mod, only: stopping_crit_h2o_ice, stopping_crit_co2
43	use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
44	m_noco2, inf_h2oice_threshold
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, 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: writediagpem
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 water 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 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 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 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 PCM data and convert 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 water ice is sublimating
553	! We compute the surface of water 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 consuption)
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 (co2_ice_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 (h2o_ice_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	write(,) "Compute ice table"
841	porefillingice_thickness_prev_iter = porefillingice_thickness
842	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
843	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
844
845	! II_d.4 Update the soil thermal properties
846	call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
847
848	! II_d.5 Update the mass of the regolith adsorbed
849	if (adsorption_pem) then
850	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice, &
851	h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
852	h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
853
854	totmassco2_adsorbded = 0.
855	totmassh2o_adsorbded = 0.
856	do ig = 1,ngrid
857	do islope =1, nslope
858	do l = 1,nsoilmx_PEM - 1
859	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l + 1) - layer_PEM(l)) &
860	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
861	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l + 1) - layer_PEM(l)) &
862	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
863	enddo
864	enddo
865	enddo
866	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
867	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
868	endif
869	endif !soil_pem
870
871	!------------------------
872	! II Run
873	! II_e Outputs
874	!------------------------
875	call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_avg_press_new/))
876	do islope = 1,nslope
877	write(str2(1:2),'(i2.2)') islope
878	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
879	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
880	call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
881	call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
882	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
883	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
884	enddo
885
886	!------------------------
887	! II Run
888	! II_f Update the tendencies
889	!------------------------
890	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, &
891	global_avg_press_PCM,global_avg_press_new)
892
893	!------------------------
894	! II Run
895	! II_g Checking the stopping criterion
896	!------------------------
897	call stopping_crit_h2o_ice(cell_area,h2o_surf_ini,h2o_ice,stopPEM,ngrid)
898	call stopping_crit_co2(cell_area,co2_surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
899
900	year_iter = year_iter + dt_pem
901	i_myear = i_myear + dt_pem
902
903	write(,) "Checking the stopping criteria..."
904	if (year_iter >= year_iter_max) stopPEM = 5
905	if (i_myear >= n_myear) stopPEM = 6
906	if (stopPEM > 0) then
907	select case (stopPEM)
908	case(1)
909	write(,) "STOPPING because surface of water ice sublimating is too low:", stopPEM, "See message above."
910	case(2)
911	write(,) "STOPPING because tendencies on water ice = 0:", stopPEM, "See message above."
912	case(3)
913	write(,) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
914	case(4)
915	write(,) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
916	case(5)
917	write(,) "STOPPING because maximum number of iterations due to orbital parameters is reached:", stopPEM
918	case(6)
919	write(,) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
920	case default
921	write(,) "STOPPING with unknown stopping criterion code:", stopPEM
922	end select
923	exit
924	else
925	write(,) "The PEM can continue!"
926	write(,) "Number of iterations of the PEM: year_iter =", year_iter
927	write(,) "Number of simulated Martian years: i_myear =", i_myear
928	endif
929
930	global_avg_press_old = global_avg_press_new
931
932	enddo ! big time iteration loop of the pem
933	!------------------------------ END RUN --------------------------------
934
935	!------------------------------- OUTPUT --------------------------------
936	!------------------------
937	! III Output
938	! III_a Update surface value for the PCM start files
939	!------------------------
940	! III_a.1 Ice update (for startfi)
941
942	! H2O ice
943	watercap = 0.
944	do ig = 1,ngrid
945	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
946	! If there is enough ice to be considered as an infinite reservoir
947	watercaptag(ig) = .true.
948	else
949	! If there too little ice to be considered as an infinite reservoir
950	! Ice is transferred to the frost
951	watercaptag(ig) = .false.
952	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
953	h2o_ice(ig,:) = 0.
954	endif
955	enddo
956
957	! III_a.2 Tsoil update (for startfi)
958	if (soil_pem) then
959	call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
960	tsoil = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
961	endif
962
963	! III_a.4 Pressure (for start)
964	ps = ps*global_avg_press_new/global_avg_press_PCM
965	ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
966
967	! III_a.5 Tracer (for start)
968	allocate(zplev_new(ngrid,nlayer + 1))
969
970	do l = 1,nlayer + 1
971	zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
972	enddo
973
974	do nnq = 1,nqtot
975	if (noms(nnq) /= "co2") then
976	do l = 1,llm - 1
977	do ig = 1,ngrid
978	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))
979	enddo
980	q(:,llm,nnq) = q(:,llm - 1,nnq)
981	enddo
982	else
983	do l = 1,llm - 1
984	do ig = 1,ngrid
985	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)) &
986	+ ((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))
987	enddo
988	q(:,llm,nnq) = q(:,llm - 1,nnq)
989	enddo
990	endif
991	enddo
992
993	! Conserving the tracers mass for PCM start files
994	do nnq = 1,nqtot
995	do ig = 1,ngrid
996	do l = 1,llm - 1
997	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
998	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
999	q(ig,l,nnq) = 1.
1000	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
1001	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
1002	endif
1003	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1004	enddo
1005	enddo
1006	enddo
1007
1008	if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
1009
1010	!------------------------
1011	! III Output
1012	! III_b Write restart_evol.nc and restartfi_evol.nc
1013	!------------------------
1014	! III_b.1 Write restart_evol.nc
1015	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
1016	pday = day_ini
1017	ztime_fin = time_phys
1018
1019	allocate(p(ip1jmp1,nlayer + 1))
1020	#ifndef CPP_1D
1021	call pression (ip1jmp1,ap,bp,ps,p)
1022	call massdair(p,masse)
1023	call dynredem0("restart_evol.nc",day_ini,phis)
1024	call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
1025	write(,) "restart_evol.nc has been written"
1026	#else
1027	call writerestart1D('restart1D_evol.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1028	write(,) "restart1D_evol.txt has been written"
1029	#endif
1030
1031	! III_b.2 Write restartfi_evol.nc
1032	#ifndef CPP_STD
1033	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
1034	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
1035	inertiedat,def_slope,subslope_dist)
1036	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1037	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
1038	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
1039	wstar,watercap,perennial_co2ice)
1040	#else
1041	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
1042	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
1043	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1044	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
1045	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1046	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
1047	tsea_ice,sea_ice)
1048	#endif
1049	write(,) "restartfi_evol.nc has been written"
1050
1051	!------------------------
1052	! III Output
1053	! III_c Write restartpem.nc
1054	!------------------------
1055	call pemdem0("restartpem.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
1056	float(day_ini),0.,nslope,def_slope,subslope_dist)
1057	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
1058	TI_PEM, porefillingice_depth,porefillingice_thickness, &
1059	co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice)
1060	write(,) "restartpem.nc has been written"
1061
1062	call info_PEM(year_iter,stopPEM,i_myear,n_myear)
1063
1064	write(,) "The PEM has run for", year_iter, "Martian years."
1065	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1066	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1067	write(,) "LL & RV & JBC: so far, so good!"
1068
1069	deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
1070	deallocate(co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfavg_before_saved)
1071	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
1072	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
1073	deallocate(co2_ice_ini)
1074	!----------------------------- END OUTPUT ------------------------------
1075
1076	END PROGRAM pem

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

Download in other formats:

Original Format