Context Navigation

pemetat0.F90 @ 3610

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

PEM:

Complete rework of "read_data_PCM_mod.F90" to reduce as much as possible memory consumption which caused problems and the model failure.
In case of very low surface shifting, checking if the index is inside the bounds of 'tsoil' and 'mlayer' and replace them by the value at surface if necessary.
Few cleanings.

JBC

File size: 25.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,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, isoil ! 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 3!'
142	else
143	nb_str_max = int(inquire_dimension_length('nb_str_max'))
144	endif
145	allocate(stratif_array(ngrid,nslope,nb_str_max,6))
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//'_air_volfrac',stratif_array(:,islope,:,6),found2)
159	found = found .or. found2
160	if (.not. found) then
161	write(,) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
162	endif ! not found
163	enddo ! islope
164	call array2stratif(stratif_array,ngrid,nslope,stratif)
165	deallocate(stratif_array)
166	endif
167
168	if (soil_pem) then
169
170	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
171	!1. Thermal Inertia
172	! a. General case
173	do islope = 1,nslope
174	write(num,'(i2.2)') islope
175	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
176	if (.not. found) then
177	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
178	write(,)'will reconstruct the values of TI_PEM'
179
180	do ig = 1,ngrid
181	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
182	!!! transition
183	delta = depth_breccia
184	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
185	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
186	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
187	do iloop=index_breccia + 2,index_bedrock
188	TI_PEM(ig,iloop,islope) = TI_breccia
189	enddo
190	else ! we keep the high TI values
191	do iloop = index_breccia + 1,index_bedrock
192	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
193	enddo
194	endif ! TI PEM and breccia comparison
195	!! transition
196	delta = depth_bedrock
197	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
198	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
199	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
200	do iloop = index_bedrock + 2,nsoil_PEM
201	TI_PEM(ig,iloop,islope) = TI_bedrock
202	enddo
203	enddo
204	else ! found
205	do iloop = nsoil_PCM + 1,nsoil_PEM
206	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.
207	enddo
208	endif ! not found
209	enddo ! islope
210
211	write(,) 'PEMETAT0: TI done'
212
213	! b. Special case for inertiedat, inertiedat_PEM
214	call get_field("inertiedat_PEM",inertiedat_PEM,found)
215	if (.not. found) then
216	do iloop = 1,nsoil_PCM
217	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
218	enddo
219	!!! zone de transition
220	delta = depth_breccia
221	do ig = 1,ngrid
222	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
223	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
224	(((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
225	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
226
227	do iloop = index_breccia + 2,index_bedrock
228	inertiedat_PEM(ig,iloop) = TI_breccia
229	enddo
230	else
231	do iloop=index_breccia + 1,index_bedrock
232	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
233	enddo
234	endif ! comparison ti breccia
235	enddo ! ig
236
237	!!! zone de transition
238	delta = depth_bedrock
239	do ig = 1,ngrid
240	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
241	(((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
242	((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
243	enddo
244
245	do iloop = index_bedrock + 2, nsoil_PEM
246	do ig = 1,ngrid
247	inertiedat_PEM(ig,iloop) = TI_bedrock
248	enddo
249	enddo
250	endif ! not found
251
252	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
253	!2. Soil Temperature
254	do islope = 1,nslope
255	write(num,fmt='(i2.2)') islope
256	call get_field("tsoil_PEM_slope"//num,tsoil_startpem(:,:,islope),found)
257	if (.not. found) then
258	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
259	write(,)'will reconstruct the values of Tsoil'
260	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
261	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
262	else
263	! predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startpem is adapted firstly
264	! for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth
265	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
266	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
267	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
268
269	do iloop = nsoil_PCM + 1,nsoil_PEM
270	tsoil_PEM(:,iloop,islope) = tsoil_startpem(:,iloop,islope)
271	enddo
272	endif !found
273
274	do isoil = nsoil_PCM + 1,nsoil_PEM
275	watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)mmol(igcm_h2o_vap)/(mugazr)
276	enddo
277	enddo ! islope
278	write(,) 'PEMETAT0: TSOIL done'
279
280	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
281	!3. Ice Table
282	if (icetable_equilibrium) then
283	call get_field("icetable_depth",icetable_depth,found)
284	if (.not. found) then
285	write(,)'PEM settings: failed loading <icetable_depth>'
286	write(,)'will reconstruct the values of the ice table given the current state'
287	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
288	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)
289	do islope = 1,nslope
290	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
291	enddo
292	endif
293	write(,) 'PEMETAT0: ICE TABLE done'
294	else if (icetable_dynamic) then
295	call get_field("ice_porefilling",ice_porefilling,found)
296	if (.not. found) write(,)'PEM settings: failed loading <ice_porefilling>'
297	call get_field("icetable_depth",icetable_depth,found)
298	if (.not. found) then
299	write(,)'PEM settings: failed loading <icetable_depth>'
300	write(,)'will reconstruct the values of the ice table given the current state'
301	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)
302	do islope = 1,nslope
303	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
304	enddo
305	endif
306	write(,) 'PEMETAT0: ICE TABLE done'
307	endif
308
309	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310	!4. CO2 & H2O Adsorption
311	if (adsorption_pem) then
312	do islope = 1,nslope
313	write(num,fmt='(i2.2)') islope
314	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
315	if (.not. found) then
316	m_co2_regolith_phys = 0.
317	exit
318	endif
319	enddo
320	do islope = 1,nslope
321	write(num,fmt='(i2.2)') islope
322	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
323	if (.not. found2) then
324	m_h2o_regolith_phys = 0.
325	exit
326	endif
327	enddo
328
329	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, &
330	m_h2o_regolith_phys,deltam_h2o_regolith_phys,m_co2_regolith_phys,deltam_co2_regolith_phys)
331
332	if (.not. found) deltam_co2_regolith_phys = 0.
333	if (.not. found2) deltam_h2o_regolith_phys = 0.
334
335	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
336	endif ! adsorption_pem
337	endif ! soil_pem
338
339	call close_startphy
340
341	else !No startpem, let's build all by hand
342
343	write(,)'There is no "startpem.nc"!'
344
345	! h2o ice
346	h2o_ice = 0.
347	write(,)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
348	do ig = 1,ngrid
349	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
350	enddo
351
352	! co2 ice
353	write(,)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
354	co2_ice = perennial_co2ice
355
356	! Stratification (layerings)
357	nb_str_max = 1
358	if (layering_algo) write(,)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
359
360	if (soil_pem) then
361
362	!a) Thermal inertia
363	do islope = 1,nslope
364	do ig = 1,ngrid
365	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
366	!!! transition
367	delta = depth_breccia
368	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
369	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
370	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
371	do iloop = index_breccia + 2,index_bedrock
372	TI_PEM(ig,iloop,islope) = TI_breccia
373	enddo
374	else ! we keep the high ti values
375	do iloop = index_breccia + 1,index_bedrock
376	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
377	enddo
378	endif
379	!! transition
380	delta = depth_bedrock
381	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
382	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
383	((layer_PEM(index_bedrock + 1) - delta)/(TI_breccia**2))))
384	do iloop = index_bedrock + 2,nsoil_PEM
385	TI_PEM(ig,iloop,islope) = TI_bedrock
386	enddo
387	enddo
388	enddo
389
390	do iloop = 1,nsoil_PCM
391	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
392	enddo
393	!!! zone de transition
394	delta = depth_breccia
395	do ig = 1,ngrid
396	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
397	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
398	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
399	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
400	do iloop = index_breccia + 2,index_bedrock
401	inertiedat_PEM(ig,iloop) = TI_breccia
402	enddo
403	else
404	do iloop = index_breccia + 1,index_bedrock
405	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
406	enddo
407	endif
408	enddo
409
410	!!! zone de transition
411	delta = depth_bedrock
412	do ig = 1,ngrid
413	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
414	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
415	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
416	enddo
417
418	do iloop = index_bedrock + 2,nsoil_PEM
419	do ig = 1,ngrid
420	inertiedat_PEM(ig,iloop) = TI_bedrock
421	enddo
422	enddo
423
424	write(,) 'PEMETAT0: TI done'
425
426	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
427	!b) Soil temperature
428	do islope = 1,nslope
429	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
430	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
431
432	! First raw initialization
433	do isoil = nsoil_PCM + 1,nsoil_PEM
434	do ig = 1,ngrid
435	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)
436	enddo
437	enddo
438	enddo !islope
439	write(,) 'PEMETAT0: TSOIL done'
440
441	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
442	!c) Ice table
443	if (icetable_equilibrium) then
444	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
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	else if (icetable_dynamic) then
451	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)
452	do islope = 1,nslope
453	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
454	enddo
455	write(,) 'PEMETAT0: ICE TABLE done'
456	endif
457
458	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
459	!d) Regolith adsorbed
460	if (adsorption_pem) then
461	m_co2_regolith_phys = 0.
462	m_h2o_regolith_phys = 0.
463	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, &
464	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
465	deltam_co2_regolith_phys = 0.
466	deltam_h2o_regolith_phys = 0.
467	endif
468
469	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
470	endif !soil_pem
471	endif ! of if (startphy_file)
472
473	if (soil_pem) then ! Sanity check
474	do ig = 1,ngrid
475	do islope = 1,nslope
476	do iloop = 1,nsoil_PEM
477	if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
478	enddo
479	enddo
480	enddo
481	endif !soil_pem
482
483	END SUBROUTINE pemetat0
484
485	END MODULE pemetat0_mod

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

Download in other formats:

Original Format