Context Navigation

pemetat0.F90 @ 3757

Last change on this file since 3757 was 3673, checked in by jbclement, 5 months ago

PEM:

Correction for the update of average pressure at the end.
Addition of pore-filled ice information for 'stratum' in the layering algorithm.
Deactivating the soil temperature shifting because of bad values given by interpolation.

JBC

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

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

Download in other formats:

Original Format