Context Navigation

pemetat0.F90 @ 3563

Last change on this file since 3563 was 3563, checked in by jbclement, 5 days ago
PEM: Correction of Norbert Schorghofer's code due to missing initialization and bad shape array as subroutine argument + some cleanings. JBC
File size: 26.4 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, &
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) :: icetable_depth ! Ice table depth [m]
63	real, dimension(ngrid,nslope), intent(inout) :: icetable_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	icetable_depth = -1. ! by default, no ice table
111	icetable_thickness = -1.
112	ice_porefilling = 0.
113	!1. Run
114	if (startpem_file) then
115	write(,)'There is a "startpem.nc"!'
116	! open pem initial state file:
117	call open_startphy(filename)
118
119	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
120	#ifndef CPP_STD
121	! h2o ice
122	h2o_ice = 0.
123	call get_field("h2o_ice",h2o_ice,found)
124	if (.not. found) then
125	write(,)'Pemetat0: failed loading <h2o_ice>'
126	write(,)'will reconstruct the values from watercaptag'
127	write(,)'with default value ''ini_huge_h2oice'''
128	do ig = 1,ngrid
129	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
130	enddo
131	else
132	! The variations of infinite reservoirs during the PCM years are taken into account
133	h2o_ice = h2o_ice + watercap
134	endif
135
136	! co2 ice
137	co2_ice = perennial_co2ice
138	#endif
139
140	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141	! Stratification (layerings)
142	nb_str_max = 1
143	if (layering_algo) then
144	found = inquire_dimension("nb_str_max")
145	if (.not. found) then
146	write(,) 'Pemetat0: failed loading <nb_str_max>!'
147	write(,) '''nb_str_max'' is set to 3!'
148	else
149	nb_str_max = int(inquire_dimension_length('nb_str_max'))
150	endif
151	allocate(stratif_array(ngrid,nslope,nb_str_max,6))
152	stratif_array = 0.
153	do islope = 1,nslope
154	write(num,'(i2.2)') islope
155	call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
156	call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
157	found = found .or. found2
158	call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
159	found = found .or. found2
160	call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
161	found = found .or. found2
162	call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
163	found = found .or. found2
164	call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
165	found = found .or. found2
166	if (.not. found) then
167	write(,) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
168	endif ! not found
169	enddo ! islope
170	call array2stratif(stratif_array,ngrid,nslope,stratif)
171	deallocate(stratif_array)
172	endif
173
174	if (soil_pem) then
175
176	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
177	!1. Thermal Inertia
178	! a. General case
179	do islope = 1,nslope
180	write(num,'(i2.2)') islope
181	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
182	if (.not. found) then
183	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
184	write(,)'will reconstruct the values of TI_PEM'
185
186	do ig = 1,ngrid
187	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
188	!!! transition
189	delta = depth_breccia
190	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
191	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
192	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
193	do iloop=index_breccia + 2,index_bedrock
194	TI_PEM(ig,iloop,islope) = TI_breccia
195	enddo
196	else ! we keep the high TI values
197	do iloop = index_breccia + 1,index_bedrock
198	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
199	enddo
200	endif ! TI PEM and breccia comparison
201	!! transition
202	delta = depth_bedrock
203	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
204	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
205	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
206	do iloop = index_bedrock + 2,nsoil_PEM
207	TI_PEM(ig,iloop,islope) = TI_bedrock
208	enddo
209	enddo
210	else ! found
211	do iloop = nsoil_PCM + 1,nsoil_PEM
212	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.
213	enddo
214	endif ! not found
215	enddo ! islope
216
217	write(,) 'PEMETAT0: THERMAL INERTIA done'
218
219	! b. Special case for inertiedat, inertiedat_PEM
220	call get_field("inertiedat_PEM",inertiedat_PEM,found)
221	if (.not. found) then
222	do iloop = 1,nsoil_PCM
223	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
224	enddo
225	!!! zone de transition
226	delta = depth_breccia
227	do ig = 1,ngrid
228	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
229	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
230	(((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
231	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
232
233	do iloop = index_breccia + 2,index_bedrock
234	inertiedat_PEM(ig,iloop) = TI_breccia
235	enddo
236	else
237	do iloop=index_breccia + 1,index_bedrock
238	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
239	enddo
240	endif ! comparison ti breccia
241	enddo ! ig
242
243	!!! zone de transition
244	delta = depth_bedrock
245	do ig = 1,ngrid
246	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
247	(((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
248	((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
249	enddo
250
251	do iloop = index_bedrock + 2, nsoil_PEM
252	do ig = 1,ngrid
253	inertiedat_PEM(ig,iloop) = TI_bedrock
254	enddo
255	enddo
256	endif ! not found
257
258	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
259	!2. Soil Temperature
260	do islope=1,nslope
261	write(num,fmt='(i2.2)') islope
262	call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
263	if (.not. found) then
264	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
265	write(,)'will reconstruct the values of Tsoil'
266	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
267	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
268	else
269	! predictor corrector: restart from year 1 of the PCM and build the evolution of tsoil at depth
270	tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
271	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
272	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
273	tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
274	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
275
276	do iloop = nsoil_PCM + 1,nsoil_PEM
277	tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
278	enddo
279	endif !found
280
281	do it = 1,timelen
282	tsoil_inst(:,nsoil_PCM + 1:nsoil_PEM,islope,it) = tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)
283	enddo
284	do isoil = nsoil_PCM + 1,nsoil_PEM
285	watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)mmol(igcm_h2o_vap)/(mugazr)
286	enddo
287	enddo ! islope
288	write(,) 'PEMETAT0: SOIL TEMP done'
289
290	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
291	!3. Ice Table
292	if (icetable_equilibrium) then
293	call get_field("icetable_depth",icetable_depth,found)
294	if (.not. found) then
295	write(,)'PEM settings: failed loading <icetable_depth>'
296	write(,)'will reconstruct the values of the ice table given the current state'
297	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
298	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
299	do islope = 1,nslope
300	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
301	enddo
302	endif
303	write(,) 'PEMETAT0: ICE TABLE done'
304	else if (icetable_dynamic) then
305	call get_field("ice_porefilling",ice_porefilling,found)
306	if (.not. found) write(,)'PEM settings: failed loading <ice_porefilling>'
307	call get_field("icetable_depth",icetable_depth,found)
308	if (.not. found) then
309	write(,)'PEM settings: failed loading <icetable_depth>'
310	write(,)'will reconstruct the values of the ice table given the current state'
311	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
312	do islope = 1,nslope
313	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
314	enddo
315	endif
316	write(,) 'PEMETAT0: ICE TABLE done'
317	endif
318
319	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
320	!4. CO2 & H2O Adsorption
321	if (adsorption_pem) then
322	do islope = 1,nslope
323	write(num,fmt='(i2.2)') islope
324	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
325	if (.not. found) then
326	m_co2_regolith_phys = 0.
327	exit
328	endif
329	enddo
330	do islope=1,nslope
331	write(num,fmt='(i2.2)') islope
332	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
333	if (.not. found2) then
334	m_h2o_regolith_phys = 0.
335	exit
336	endif
337	enddo
338
339	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, &
340	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
341
342	if (.not. found) deltam_co2_regolith_phys = 0.
343	if (.not. found2) deltam_h2o_regolith_phys = 0.
344
345	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
346	endif ! adsorption_pem
347	endif ! soil_pem
348
349	call close_startphy
350
351	else !No startpem, let's build all by hand
352
353	write(,)'There is no "startpem.nc"!'
354
355	! h2o ice
356	h2o_ice = 0.
357	write(,)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
358	do ig = 1,ngrid
359	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
360	enddo
361
362	! co2 ice
363	write(,)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
364	co2_ice = perennial_co2ice
365
366	! Stratification (layerings)
367	nb_str_max = 1
368	if (layering_algo) write(,)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
369
370	if (soil_pem) then
371
372	!a) Thermal inertia
373	do islope = 1,nslope
374	do ig = 1,ngrid
375	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
376	!!! transition
377	delta = depth_breccia
378	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
379	(((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
380	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
381	do iloop=index_breccia + 2,index_bedrock
382	TI_PEM(ig,iloop,islope) = TI_breccia
383	enddo
384	else ! we keep the high ti values
385	do iloop=index_breccia + 1,index_bedrock
386	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
387	enddo
388	endif
389	!! transition
390	delta = depth_bedrock
391	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
392	(((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
393	((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
394	do iloop = index_bedrock + 2,nsoil_PEM
395	TI_PEM(ig,iloop,islope) = TI_bedrock
396	enddo
397	enddo
398	enddo
399
400	do iloop = 1,nsoil_PCM
401	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
402	enddo
403	!!! zone de transition
404	delta = depth_breccia
405	do ig = 1,ngrid
406	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
407	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
408	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
409	((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
410	do iloop = index_breccia + 2,index_bedrock
411	inertiedat_PEM(ig,iloop) = TI_breccia
412	enddo
413	else
414	do iloop = index_breccia + 1,index_bedrock
415	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
416	enddo
417	endif
418	enddo
419
420	!!! zone de transition
421	delta = depth_bedrock
422	do ig = 1,ngrid
423	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
424	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
425	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
426	enddo
427
428	do iloop = index_bedrock + 2,nsoil_PEM
429	do ig = 1,ngrid
430	inertiedat_PEM(ig,iloop) = TI_bedrock
431	enddo
432	enddo
433
434	write(,) 'PEMETAT0: TI done'
435
436	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
437	!b) Soil temperature
438	do islope = 1,nslope
439	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
440	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
441
442	! First raw initialization
443	do it = 1,timelen
444	do isoil = nsoil_PCM + 1,nsoil_PEM
445	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
446	enddo
447	enddo
448
449	do it = 1,timelen
450	do isoil = nsoil_PCM + 1,nsoil_PEM
451	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
452	enddo
453	enddo
454
455	do isoil = nsoil_PCM + 1,nsoil_PEM
456	do ig = 1,ngrid
457	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)
458	enddo
459	enddo
460	enddo !islope
461	write(,) 'PEMETAT0: TSOIL done'
462
463	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
464	!c) Ice table
465	if (icetable_equilibrium) then
466	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
467	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
468	do islope = 1,nslope
469	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
470	enddo
471	write(,) 'PEMETAT0: Ice table done'
472	else if (icetable_dynamic) then
473	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
474	do islope = 1,nslope
475	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
476	enddo
477	write(,) 'PEMETAT0: Ice table done'
478	endif
479
480	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481	!d) Regolith adsorbed
482	if (adsorption_pem) then
483	m_co2_regolith_phys = 0.
484	m_h2o_regolith_phys = 0.
485	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, &
486	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
487	deltam_co2_regolith_phys = 0.
488	deltam_h2o_regolith_phys = 0.
489	endif
490
491	write(,) 'PEMETAT0: CO2 adsorption done'
492	endif !soil_pem
493	endif ! of if (startphy_file)
494
495	if (soil_pem) then ! Sanity check
496	do ig = 1,ngrid
497	do islope = 1,nslope
498	do iloop = 1,nsoil_PEM
499	if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
500	enddo
501	enddo
502	enddo
503	endif !soil_pem
504
505	END SUBROUTINE pemetat0
506
507	END MODULE pemetat0_mod

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

Download in other formats:

Original Format