Context Navigation

pem.F90 @ 3130

Last change on this file since 3130 was 3130, checked in by jbclement, 19 months ago

PEM:
The perennial co2 ice is now taken into account with co2 frost (qsurf) to compute the tendency and to make the update + Rework of how co2 frost is converted to perennial co2 ice at the end of the PEM run + Correction of the value of 'threshold_co2_frost2perennial' to correspond to 10 m + Perennial co2 ice is now handled outside 'paleoclimate' in "phyetat0_mod.F90" of the Mars PCM + Some cleanings.

/!\ Commit for the PEM management of co2 ice before a rework of ice management in the PEM!
JBC

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

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

Download in other formats:

Original Format