Context Navigation

pemetat0.F90 @ 3498

Last change on this file since 3498 was 3498, checked in by jbclement, 2 weeks ago

PEM:

Correction of the variable name for the ice table depth in "pemetat0.F90". So it is now got as intended from the "startpem.nc" file;
Renaming of the tendencies in the PEM with the prefix 'd_' instead of 'tend_';
Modification of the PEM time step type from integer to real. As a consequence, all time variables are now of real type. This change adds the possibility to consider fractions of year as time step.

JBC

File size: 29.0 KB

Line
1	MODULE pemetat0_mod
2
3	implicit none
4
5	!=======================================================================
6	contains
7	!=======================================================================
8
9	SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness, &
10	ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice, &
11	global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, &
12	m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
13
14	use iostart_PEM, only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
15	use comsoil_h_PEM, only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
16	use comsoil_h, only: volcapa, inertiedat
17	use adsorption_mod, only: regolith_adsorption, adsorption_pem
18	use ice_table_mod, only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
19	use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
20	use soil_thermalproperties_mod, only: update_soil_thermalproperties
21	use tracer_mod, only: mmol, igcm_h2o_vap ! tracer names and molar masses
22	use abort_pem_mod, only: abort_pem
23	use compute_soiltemp_mod, only: ini_tsoil_pem, compute_tsoil_pem
24	use comslope_mod, only: def_slope_mean, subslope_dist
25	use layering_mod, only: layering, array2stratif, nb_str_max, layering_algo
26
27	#ifndef CPP_STD
28	use comcstfi_h, only: r, mugaz, pi
29	use surfdat_h, only: watercaptag, watercap, perennial_co2ice
30	#else
31	use comcstfi_mod, only: r, mugaz, pi
32	#endif
33
34	implicit none
35
36	include "callkeys.h"
37
38	character(*), intent(in) :: filename ! name of the startfi_PEM.nc
39	integer, intent(in) :: ngrid ! # of physical grid points
40	integer, intent(in) :: nsoil_PCM ! # of vertical grid points in the PCM
41	integer, intent(in) :: nsoil_PEM ! # of vertical grid points in the PEM
42	integer, intent(in) :: nslope ! # of sub-grid slopes
43	integer, intent(in) :: timelen ! # time samples
44	real, intent(in) :: timestep ! time step [s]
45	real, intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
46	real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr1 ! surface temperature at the first year of PCM call [K]
47	real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr2 ! surface temperature at the second year of PCM call [K]
48	real, dimension(ngrid,timelen), intent(in) :: q_co2 ! MMR tracer co2 [kg/kg]
49	real, dimension(ngrid,timelen), intent(in) :: q_h2o ! MMR tracer h2o [kg/kg]
50	real, dimension(ngrid,timelen), intent(in) :: ps_inst ! surface pressure [Pa]
51	real, dimension(ngrid,nslope), intent(in) :: d_h2oice ! tendencies on h2o ice
52	real, dimension(ngrid,nslope), intent(in) :: d_co2ice ! tendencies on co2 ice
53	real, dimension(ngrid,nslope), intent(in) :: watersurf_avg ! surface water ice density, yearly averaged (kg/m^3)
54	! outputs
55	real, dimension(ngrid), intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
56	real, dimension(ngrid), intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
57	real, dimension(ngrid,nslope), intent(out) :: h2o_ice ! h2o ice amount [kg/m^2]
58	real, dimension(ngrid,nslope), intent(out) :: co2_ice ! co2 ice amount [kg/m^2]
59	type(layering), dimension(ngrid,nslope), intent(inout) :: stratif ! stratification (layerings)
60	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! soil (mid-layer) thermal inertia in the PEM grid [SI]
61	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM ! soil (mid-layer) temperature [K]
62	real, dimension(ngrid,nslope), intent(inout) :: ice_table_depth ! Ice table depth [m]
63	real, dimension(ngrid,nslope), intent(inout) :: ice_table_thickness ! Ice table thickness [m]
64	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling ! Subsurface ice pore filling [m3/m3]
65	real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst ! instantaneous soil (mid-layer) temperature [K]
66	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
67	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
68	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg ! surface water ice density, yearly averaged (kg/m^3)
69	! local
70	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM ! soil temperature saved in the start [K]
71	real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM ! soil thermal inertia saved in the start [SI]
72	logical :: found, found2 ! check if variables are found in the start
73	integer :: iloop, ig, islope, it, isoil, k ! index for loops
74	real :: kcond ! Thermal conductivity, intermediate variable [SI]
75	real :: delta ! Depth of the interface regolith-breccia, breccia -bedrock [m]
76	character(2) :: num ! intermediate string to read PEM start sloped variables
77	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1 ! intermediate soil temperature during yr 1 [K]
78	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2 ! intermediate soil temperature during yr 2 [K]
79	logical :: startpem_file ! boolean to check if we read the startfile or not
80	real, dimension(:,:,:,:), allocatable :: stratif_array ! Array for stratification (layerings)
81
82	#ifdef CPP_STD
83	logical, dimension(ngrid) :: watercaptag
84	watercaptag = .false.
85	#endif
86
87	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
88	!!!
89	!!! Purpose: read start_pem. Need a specific iostart_PEM
90	!!!
91	!!! Order: 0. Previous year of the PEM run
92	!!! Ice initialization
93	!!! 1. Thermal Inertia
94	!!! 2. Soil Temperature
95	!!! 3. Ice table
96	!!! 4. Mass of CO2 & H2O adsorbed
97	!!!
98	!!! /!\ This order must be respected !
99	!!! Author: LL
100	!!!
101	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
102
103	!0.1 Check if the start_PEM exist.
104
105	inquire(file = filename,exist = startpem_file)
106
107	write(,)'Is there a "startpem" file?',startpem_file
108
109	!0.2 Set to default values
110	ice_table_depth = -1. ! by default, no ice table
111	ice_table_thickness = -1.
112	!1. Run
113	if (startpem_file) then
114	write(,)'There is a "startpem.nc"!'
115	! open pem initial state file:
116	call open_startphy(filename)
117
118	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
119	#ifndef CPP_STD
120	! h2o ice
121	h2o_ice = 0.
122	call get_field("h2o_ice",h2o_ice,found)
123	if (.not. found) then
124	write(,)'Pemetat0: failed loading <h2o_ice>'
125	write(,)'will reconstruct the values from watercaptag'
126	write(,)'with default value ''ini_huge_h2oice'''
127	do ig = 1,ngrid
128	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
129	enddo
130	else
131	! The variations of infinite reservoirs during the PCM years are taken into account
132	h2o_ice = h2o_ice + watercap
133	endif
134
135	! co2 ice
136	co2_ice = perennial_co2ice
137	#endif
138
139	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140	! Stratification (layerings)
141	nb_str_max = 1
142	if (layering_algo) then
143	found = inquire_dimension("nb_str_max")
144	if (.not. found) then
145	write(,) 'Pemetat0: failed loading <nb_str_max>!'
146	write(,) '''nb_str_max'' is set to 3!'
147	else
148	nb_str_max = int(inquire_dimension_length('nb_str_max'))
149	endif
150	allocate(stratif_array(ngrid,nslope,nb_str_max,6))
151	stratif_array = 0.
152	do islope = 1,nslope
153	write(num,'(i2.2)') islope
154	call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
155	call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
156	found = found .or. found2
157	call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
158	found = found .or. found2
159	call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
160	found = found .or. found2
161	call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
162	found = found .or. found2
163	call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
164	found = found .or. found2
165	if (.not. found) then
166	write(,) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
167	endif ! not found
168	enddo ! islope
169	call array2stratif(stratif_array,ngrid,nslope,stratif)
170	deallocate(stratif_array)
171	endif
172
173	if (soil_pem) then
174
175	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
176	!1. Thermal Inertia
177	! a. General case
178	do islope = 1,nslope
179	write(num,'(i2.2)') islope
180	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
181	if (.not. found) then
182	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
183	write(,)'will reconstruct the values of TI_PEM'
184
185	do ig = 1,ngrid
186	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
187	!!! transition
188	delta = depth_breccia
189	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
190	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
191	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
192	do iloop=index_breccia + 2,index_bedrock
193	TI_PEM(ig,iloop,islope) = TI_breccia
194	enddo
195	else ! we keep the high TI values
196	do iloop = index_breccia + 1,index_bedrock
197	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
198	enddo
199	endif ! TI PEM and breccia comparison
200	!! transition
201	delta = depth_bedrock
202	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
203	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
204	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
205	do iloop = index_bedrock + 2,nsoil_PEM
206	TI_PEM(ig,iloop,islope) = TI_bedrock
207	enddo
208	enddo
209	else ! found
210	do iloop = nsoil_PCM + 1,nsoil_PEM
211	TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope) ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
212	enddo
213	endif ! not found
214	enddo ! islope
215
216	write(,) 'PEMETAT0: THERMAL INERTIA done'
217
218	! b. Special case for inertiedat, inertiedat_PEM
219	call get_field("inertiedat_PEM",inertiedat_PEM,found)
220	if (.not. found) then
221	do iloop = 1,nsoil_PCM
222	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
223	enddo
224	!!! zone de transition
225	delta = depth_breccia
226	do ig = 1,ngrid
227	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
228	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
229	(((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
230	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
231
232	do iloop = index_breccia + 2,index_bedrock
233	inertiedat_PEM(ig,iloop) = TI_breccia
234	enddo
235	else
236	do iloop=index_breccia + 1,index_bedrock
237	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
238	enddo
239	endif ! comparison ti breccia
240	enddo ! ig
241
242	!!! zone de transition
243	delta = depth_bedrock
244	do ig = 1,ngrid
245	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
246	(((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
247	((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
248	enddo
249
250	do iloop = index_bedrock + 2, nsoil_PEM
251	do ig = 1,ngrid
252	inertiedat_PEM(ig,iloop) = TI_bedrock
253	enddo
254	enddo
255	endif ! not found
256
257	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
258	!2. Soil Temperature
259	do islope=1,nslope
260	write(num,fmt='(i2.2)') islope
261	call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
262	if (.not. found) then
263	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
264	write(,)'will reconstruct the values of Tsoil'
265	! do ig = 1,ngrid
266	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
267	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
268	! do iloop=index_breccia+2,index_bedrock
269	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
270	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
271	! enddo
272	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
273	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
274	!
275	! do iloop=index_bedrock+2,nsoil_PEM
276	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
277	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
278	! enddo
279	! enddo
280	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
281	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
282	else
283	! predictor corrector: restart from year 1 of the PCM and build the evolution of
284	! tsoil at depth
285	tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
286	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
287	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
288	tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
289	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
290
291	do iloop = nsoil_PCM+1,nsoil_PEM
292	tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
293	enddo
294	endif !found
295
296	do it = 1,timelen
297	do isoil = nsoil_PCM+1,nsoil_PEM
298	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
299	enddo
300	enddo
301	do isoil = nsoil_PCM+1,nsoil_PEM
302	do ig = 1,ngrid
303	watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)mmol(igcm_h2o_vap)/(mugazr)
304	enddo
305	enddo
306	enddo ! islope
307	write(,) 'PEMETAT0: SOIL TEMP done'
308
309	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310	!3. Ice Table
311	if (icetable_equilibrium) then
312	call get_field("ice_table_depth",ice_table_depth,found)
313	if (.not. found) then
314	write(,)'PEM settings: failed loading <ice_table_depth>'
315	write(,)'will reconstruct the values of the ice table given the current state'
316	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
317	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
318	do islope = 1,nslope
319	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
320	enddo
321	endif
322	write(,) 'PEMETAT0: ICE TABLE done'
323	else if (icetable_dynamic) then
324	call get_field("ice_table_depth",ice_table_depth,found)
325	if (.not. found) then
326	write(,)'PEM settings: failed loading <ice_table_depth>'
327	write(,)'will reconstruct the values of the ice table given the current state'
328	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
329	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
330	do islope = 1,nslope
331	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
332	enddo
333	endif
334	call get_field("ice_porefilling",ice_porefilling,found)
335	if (.not. found) then
336	write(,)'PEM settings: failed loading <ice_porefilling>'
337	ice_porefilling = 0.
338	endif
339	write(,) 'PEMETAT0: ICE TABLE done'
340	endif
341
342	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
343	!4. CO2 & H2O Adsorption
344	if (adsorption_pem) then
345	do islope = 1,nslope
346	write(num,fmt='(i2.2)') islope
347	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
348	if (.not. found) then
349	m_co2_regolith_phys = 0.
350	exit
351	endif
352	enddo
353	do islope=1,nslope
354	write(num,fmt='(i2.2)') islope
355	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
356	if (.not. found2) then
357	m_h2o_regolith_phys = 0.
358	exit
359	endif
360	enddo
361
362	call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
363	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
364
365	if (.not. found) deltam_co2_regolith_phys = 0.
366	if (.not. found2) deltam_h2o_regolith_phys = 0.
367
368	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
369	endif ! adsorption_pem
370	endif ! soil_pem
371
372	call close_startphy
373
374	else !No startpem, let's build all by hand
375
376	write(,)'There is no "startpem.nc"!'
377
378	! h2o ice
379	h2o_ice = 0.
380	write(,)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
381	do ig = 1,ngrid
382	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
383	enddo
384
385	! co2 ice
386	write(,)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
387	co2_ice = perennial_co2ice
388
389	! Stratification (layerings)
390	nb_str_max = 1
391	if (layering_algo) write(,)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
392
393	if (soil_pem) then
394
395	!a) Thermal inertia
396	do islope = 1,nslope
397	do ig = 1,ngrid
398	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
399	!!! transition
400	delta = depth_breccia
401	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
402	(((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
403	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
404	do iloop=index_breccia + 2,index_bedrock
405	TI_PEM(ig,iloop,islope) = TI_breccia
406	enddo
407	else ! we keep the high ti values
408	do iloop=index_breccia + 1,index_bedrock
409	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
410	enddo
411	endif
412	!! transition
413	delta = depth_bedrock
414	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
415	(((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
416	((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
417	do iloop = index_bedrock + 2,nsoil_PEM
418	TI_PEM(ig,iloop,islope) = TI_bedrock
419	enddo
420	enddo
421	enddo
422
423	do iloop = 1,nsoil_PCM
424	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
425	enddo
426	!!! zone de transition
427	delta = depth_breccia
428	do ig = 1,ngrid
429	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
430	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
431	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
432	((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
433	do iloop = index_breccia + 2,index_bedrock
434	inertiedat_PEM(ig,iloop) = TI_breccia
435	enddo
436	else
437	do iloop = index_breccia + 1,index_bedrock
438	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
439	enddo
440	endif
441	enddo
442
443	!!! zone de transition
444	delta = depth_bedrock
445	do ig = 1,ngrid
446	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
447	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
448	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
449	enddo
450
451	do iloop = index_bedrock + 2,nsoil_PEM
452	do ig = 1,ngrid
453	inertiedat_PEM(ig,iloop) = TI_bedrock
454	enddo
455	enddo
456
457	write(,) 'PEMETAT0: TI done'
458
459	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
460	!b) Soil temperature
461	do islope = 1,nslope
462	! do ig = 1,ngrid
463	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
464	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
465	!
466	! do iloop=index_breccia+2,index_bedrock
467	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
468	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
469	! enddo
470	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
471	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
472	!
473	! do iloop=index_bedrock+2,nsoil_PEM
474	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
475	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
476	! enddo
477	! enddo
478	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
479	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
480
481	! First raw initialization
482	do it = 1,timelen
483	do isoil = nsoil_PCM + 1,nsoil_PEM
484	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
485	enddo
486	enddo
487
488	do it = 1,timelen
489	do isoil = nsoil_PCM + 1,nsoil_PEM
490	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
491	enddo
492	enddo
493
494	do isoil = nsoil_PCM + 1,nsoil_PEM
495	do ig = 1,ngrid
496	watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)mmol(igcm_h2o_vap)/(mugazr)
497	enddo
498	enddo
499	enddo !islope
500	write(,) 'PEMETAT0: TSOIL done'
501
502	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
503	!c) Ice table
504	if (icetable_equilibrium) then
505	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
506	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
507	do islope = 1,nslope
508	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
509	enddo
510	write(,) 'PEMETAT0: Ice table done'
511	else if (icetable_dynamic) then
512	ice_porefilling = 0.
513	ice_table_depth = 0.
514	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
515	do islope = 1,nslope
516	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
517	enddo
518	write(,) 'PEMETAT0: Ice table done'
519	endif
520
521	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
522	!d) Regolith adsorbed
523	if (adsorption_pem) then
524	m_co2_regolith_phys = 0.
525	m_h2o_regolith_phys = 0.
526	call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
527	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
528	deltam_co2_regolith_phys = 0.
529	deltam_h2o_regolith_phys = 0.
530	endif
531
532	write(,) 'PEMETAT0: CO2 adsorption done'
533	endif !soil_pem
534	endif ! of if (startphy_file)
535
536	if (soil_pem) then ! Sanity check
537	do ig = 1,ngrid
538	do islope = 1,nslope
539	do iloop = 1,nsoil_PEM
540	if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
541	enddo
542	enddo
543	enddo
544	endif !soil_pem
545
546	END SUBROUTINE pemetat0
547
548	END MODULE pemetat0_mod

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

Download in other formats:

Original Format