Context Navigation

pemetat0.F90 @ 3188

Last change on this file since 3188 was 3178, checked in by llange, 11 months ago
PEM Cleaning of the several subroutine regarding soil temperatures: they are now gathered in an unique Tsoil module. (cosmetic commit) LL
File size: 25.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,ice_table,ice_table_thickness, &
10	tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,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)
13
14	use iostart_PEM, only: open_startphy, close_startphy, get_field, get_var
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
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
26	#ifndef CPP_STD
27	use comcstfi_h, only: r, mugaz, pi
28	use surfdat_h, only: watercaptag, watercap, perennial_co2ice
29	#else
30	use comcstfi_mod, only: r, mugaz, pi
31	#endif
32
33	implicit none
34
35	include "callkeys.h"
36
37	character(len=*), intent(in) :: filename ! name of the startfi_PEM.nc
38	integer, intent(in) :: ngrid ! # of physical grid points
39	integer, intent(in) :: nsoil_PCM ! # of vertical grid points in the PCM
40	integer, intent(in) :: nsoil_PEM ! # of vertical grid points in the PEM
41	integer, intent(in) :: nslope ! # of sub-grid slopes
42	integer, intent(in) :: timelen ! # time samples
43	real, intent(in) :: timestep ! time step [s]
44	real, intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
45	real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr1 ! surface temperature at the first year of PCM call [K]
46	real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr2 ! surface temperature at the second year of PCM call [K]
47	real, dimension(ngrid,timelen), intent(in) :: q_co2 ! MMR tracer co2 [kg/kg]
48	real, dimension(ngrid,timelen), intent(in) :: q_h2o ! MMR tracer h2o [kg/kg]
49	real, dimension(ngrid,timelen), intent(in) :: ps_inst ! surface pressure [Pa]
50	real, dimension(ngrid,nslope), intent(in) :: tend_h2o_ice ! tendencies on h2o ice
51	real, dimension(ngrid,nslope), intent(in) :: tend_co2_ice ! tendencies on co2 ice
52	real, dimension(ngrid,nslope), intent(in) :: watersurf_avg ! surface water ice density, yearly averaged (kg/m^3)
53	! outputs
54	real, dimension(ngrid), intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
55	real, dimension(ngrid), intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
56	real, dimension(ngrid,nslope), intent(out) :: h2o_ice ! h2o ice amount [kg/m^2]
57	real, dimension(ngrid,nslope), intent(out) :: co2_ice ! co2 ice amount [kg/m^2]
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) :: ice_table ! Ice table depth [m]
61	real, dimension(ngrid,nslope), intent(inout) :: ice_table_thickness ! Ice table thickness [m]
62	real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst ! instantaneous soil (mid-layer) temperature [K]
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 ! check if variables are found in the start
70	logical :: found2 ! check if variables are found in the start
71	integer :: iloop, ig, islope, it, isoil ! index for loops
72	real :: kcond ! Thermal conductivity, intermediate variable [SI]
73	real :: delta ! Depth of the interface regolith-breccia, breccia -bedrock [m]
74	character(2) :: num ! intermediate string to read PEM start sloped variables
75	real, dimension(ngrid,nsoil_PEM) :: tsoil_saved ! saved soil temperature [K]
76	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1 ! intermediate soil temperature during yr 1 [K]
77	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2 ! intermediate soil temperature during yr 2 [K]
78	real, dimension(ngrid,nsoil_PEM - 1) :: alph_tmp ! Intermediate for tsoil computation []
79	real, dimension(ngrid,nsoil_PEM - 1) :: beta_tmp ! Intermediate for tsoil computatio []
80	logical :: startpem_file ! boolean to check if we read the startfile or not
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 start PEM?',startpem_file
108
109	!0.2 Set to default values
110	ice_table = -1. ! by default, no ice table
111	ice_table_thickness = -1.
112	!1. Run
113	if (startpem_file) then
114	! open pem initial state file:
115	call open_startphy(filename)
116
117	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
118	#ifndef CPP_STD
119	! h2o ice
120	h2o_ice = 0.
121	call get_field("h2o_ice",h2o_ice,found)
122	if (.not. found) then
123	write(,)'Pemetat0: failed loading <h2o_ice>'
124	write(,)'will reconstruct the values from watercaptag'
125	write(,)'with default value ''ini_huge_h2oice'''
126	do ig = 1,ngrid
127	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
128	enddo
129	else
130	! The variations of infinite reservoirs during the PCM years are taken into account
131	h2o_ice = h2o_ice + watercap
132	endif
133
134	! co2 ice
135	co2_ice = perennial_co2ice
136	#endif
137
138	if (soil_pem) then
139
140	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141	!1. Thermal Inertia
142	! a. General case
143	do islope = 1,nslope
144	write(num,'(i2.2)') islope
145	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
146	if (.not. found) then
147	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
148	write(,)'will reconstruct the values of TI_PEM'
149
150	do ig = 1,ngrid
151	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
152	!!! transition
153	delta = depth_breccia
154	TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
155	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
156	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
157	do iloop=index_breccia+2,index_bedrock
158	TI_PEM(ig,iloop,islope) = TI_breccia
159	enddo
160	else ! we keep the high TI values
161	do iloop = index_breccia + 1,index_bedrock
162	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
163	enddo
164	endif ! TI PEM and breccia comparison
165	!! transition
166	delta = depth_bedrock
167	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
168	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
169	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
170	do iloop = index_bedrock + 2,nsoil_PEM
171	TI_PEM(ig,iloop,islope) = TI_bedrock
172	enddo
173	enddo
174	else ! found
175	do iloop = nsoil_PCM + 1,nsoil_PEM
176	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.
177	enddo
178	endif ! not found
179	enddo ! islope
180
181	write(,) 'PEMETAT0: THERMAL INERTIA done'
182
183	! b. Special case for inertiedat, inertiedat_PEM
184	call get_field("inertiedat_PEM",inertiedat_PEM,found)
185	if (.not.found) then
186	do iloop = 1,nsoil_PCM
187	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
188	enddo
189	!!! zone de transition
190	delta = depth_breccia
191	do ig = 1,ngrid
192	if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
193	inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
194	(((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
195	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
196
197	do iloop = index_breccia+2,index_bedrock
198	inertiedat_PEM(ig,iloop) = TI_breccia
199	enddo
200	else
201	do iloop=index_breccia+1,index_bedrock
202	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
203	enddo
204	endif ! comparison ti breccia
205	enddo ! ig
206
207	!!! zone de transition
208	delta = depth_bedrock
209	do ig = 1,ngrid
210	inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
211	(((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
212	((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
213	enddo
214
215	do iloop = index_bedrock + 2, nsoil_PEM
216	do ig = 1,ngrid
217	inertiedat_PEM(ig,iloop) = TI_bedrock
218	enddo
219	enddo
220	endif ! not found
221
222	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
223	!2. Soil Temperature
224	do islope=1,nslope
225	write(num,fmt='(i2.2)') islope
226	call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
227	if (.not.found) then
228	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
229	write(,)'will reconstruct the values of Tsoil'
230	! do ig = 1,ngrid
231	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
232	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
233	! do iloop=index_breccia+2,index_bedrock
234	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
235	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
236	! enddo
237	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
238	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
239	!
240	! do iloop=index_bedrock+2,nsoil_PEM
241	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
242	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
243	! enddo
244	! enddo
245	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
246	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
247	else
248	! predictor corrector: restart from year 1 of the PCM and build the evolution of
249	! tsoil at depth
250	tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
251	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
252	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
253	tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
254	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
255
256	do iloop = nsoil_PCM+1,nsoil_PEM
257	tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
258	enddo
259	endif !found
260
261	do it = 1,timelen
262	do isoil = nsoil_PCM+1,nsoil_PEM
263	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
264	enddo
265	enddo
266	do isoil = nsoil_PCM+1,nsoil_PEM
267	do ig = 1,ngrid
268	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)
269	enddo
270	enddo
271	enddo ! islope
272	write(,) 'PEMETAT0: SOIL TEMP done'
273
274	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275	!3. Ice Table
276	if(icetable_equilibrium) then
277	call get_field("ice_table",ice_table,found)
278	if (.not. found) then
279	write(,)'PEM settings: failed loading <ice_table>'
280	write(,)'will reconstruct the values of the ice table given the current state'
281	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
282	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
283	do islope = 1,nslope
284	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
285	enddo
286	endif
287	write(,) 'PEMETAT0: ICE TABLE done'
288	endif
289
290	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
291	!4. CO2 & H2O Adsorption
292	if (adsorption_pem) then
293	do islope = 1,nslope
294	write(num,fmt='(i2.2)') islope
295	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
296	if (.not. found) then
297	m_co2_regolith_phys = 0.
298	exit
299	endif
300	enddo
301	do islope=1,nslope
302	write(num,fmt='(i2.2)') islope
303	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
304	if (.not. found2) then
305	m_h2o_regolith_phys = 0.
306	exit
307	endif
308	enddo
309
310	call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
311	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
312
313	if (.not. found) deltam_co2_regolith_phys = 0.
314	if (.not.found2) deltam_h2o_regolith_phys = 0.
315
316	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
317	endif ! adsorption_pem
318	endif ! soil_pem
319
320	call close_startphy
321
322	else !No startfi, let's build all by hand
323
324	! h2o ice
325	h2o_ice = 0.
326	write(,)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_huge_h2oice''.'
327	do ig = 1,ngrid
328	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
329	enddo
330
331	! co2 ice
332	co2_ice = perennial_co2ice
333
334	if (soil_pem) then
335
336	!a) Thermal inertia
337	do islope = 1,nslope
338	do ig = 1,ngrid
339	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
340	!!! transition
341	delta = depth_breccia
342	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
343	(((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
344	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
345	do iloop=index_breccia + 2,index_bedrock
346	TI_PEM(ig,iloop,islope) = TI_breccia
347	enddo
348	else ! we keep the high ti values
349	do iloop=index_breccia + 1,index_bedrock
350	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
351	enddo
352	endif
353	!! transition
354	delta = depth_bedrock
355	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
356	(((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
357	((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
358	do iloop = index_bedrock + 2,nsoil_PEM
359	TI_PEM(ig,iloop,islope) = TI_bedrock
360	enddo
361	enddo
362	enddo
363
364	do iloop = 1,nsoil_PCM
365	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
366	enddo
367	!!! zone de transition
368	delta = depth_breccia
369	do ig = 1,ngrid
370	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
371	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
372	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
373	((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
374	do iloop = index_breccia + 2,index_bedrock
375	inertiedat_PEM(ig,iloop) = TI_breccia
376	enddo
377	else
378	do iloop = index_breccia + 1,index_bedrock
379	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
380	enddo
381	endif
382	enddo
383
384	!!! zone de transition
385	delta = depth_bedrock
386	do ig = 1,ngrid
387	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
388	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
389	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
390	enddo
391
392	do iloop = index_bedrock + 2,nsoil_PEM
393	do ig = 1,ngrid
394	inertiedat_PEM(ig,iloop) = TI_bedrock
395	enddo
396	enddo
397
398	write(,) 'PEMETAT0: TI done'
399
400	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
401	!b) Soil temperature
402	do islope = 1,nslope
403	! do ig = 1,ngrid
404	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
405	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
406	!
407	! do iloop=index_breccia+2,index_bedrock
408	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
409	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
410	! enddo
411	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
412	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
413	!
414	! do iloop=index_bedrock+2,nsoil_PEM
415	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
416	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
417	! enddo
418	! enddo
419	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
420	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
421
422	! First raw initialization
423	do it = 1,timelen
424	do isoil = nsoil_PCM + 1,nsoil_PEM
425	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
426	enddo
427	enddo
428
429	do it = 1,timelen
430	do isoil = nsoil_PCM + 1,nsoil_PEM
431	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
432	enddo
433	enddo
434
435	do isoil = nsoil_PCM + 1,nsoil_PEM
436	do ig = 1,ngrid
437	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)
438	enddo
439	enddo
440	enddo !islope
441	write(,) 'PEMETAT0: TSOIL done'
442
443	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
444	!c) Ice table
445	if(icetable_equilibrium) then
446	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
447	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
448	do islope = 1,nslope
449	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
450	enddo
451	write(,) 'PEMETAT0: Ice table done'
452	endif
453
454	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
455	!d) Regolith adsorbed
456	if (adsorption_pem) then
457	m_co2_regolith_phys = 0.
458	m_h2o_regolith_phys = 0.
459	call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
460	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
461	deltam_co2_regolith_phys = 0.
462	deltam_h2o_regolith_phys = 0.
463	endif
464
465	write(,) 'PEMETAT0: CO2 adsorption done'
466	endif !soil_pem
467	endif ! of if (startphy_file)
468
469	if (soil_pem) then ! Sanity check
470	do ig = 1,ngrid
471	do islope = 1,nslope
472	do iloop = 1,nsoil_PEM
473	if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
474	enddo
475	enddo
476	enddo
477	endif !soil_pem
478
479	END SUBROUTINE pemetat0
480
481	END MODULE pemetat0_mod

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

Download in other formats:

Original Format