Context Navigation

pemetat0.F90 @ 3210

Last change on this file since 3210 was 3206, checked in by jbclement, 22 months ago
PEM: Cleanings of unused variables/arguments and bad type conversions. JBC
File size: 25.1 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,nslope) :: tsoil_tmp_yr1 ! intermediate soil temperature during yr 1 [K]
76	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2 ! intermediate soil temperature during yr 2 [K]
77	logical :: startpem_file ! boolean to check if we read the startfile or not
78
79	#ifdef CPP_STD
80	logical, dimension(ngrid) :: watercaptag
81	watercaptag = .false.
82	#endif
83
84	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
85	!!!
86	!!! Purpose: read start_pem. Need a specific iostart_PEM
87	!!!
88	!!! Order: 0. Previous year of the PEM run
89	!!! Ice initialization
90	!!! 1. Thermal Inertia
91	!!! 2. Soil Temperature
92	!!! 3. Ice table
93	!!! 4. Mass of CO2 & H2O adsorbed
94	!!!
95	!!! /!\ This order must be respected !
96	!!! Author: LL
97	!!!
98	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
99
100	!0.1 Check if the start_PEM exist.
101
102	inquire(file = filename,exist = startpem_file)
103
104	write(,)'Is start PEM?',startpem_file
105
106	!0.2 Set to default values
107	ice_table = -1. ! by default, no ice table
108	ice_table_thickness = -1.
109	!1. Run
110	if (startpem_file) then
111	! open pem initial state file:
112	call open_startphy(filename)
113
114	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
115	#ifndef CPP_STD
116	! h2o ice
117	h2o_ice = 0.
118	call get_field("h2o_ice",h2o_ice,found)
119	if (.not. found) then
120	write(,)'Pemetat0: failed loading <h2o_ice>'
121	write(,)'will reconstruct the values from watercaptag'
122	write(,)'with default value ''ini_huge_h2oice'''
123	do ig = 1,ngrid
124	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
125	enddo
126	else
127	! The variations of infinite reservoirs during the PCM years are taken into account
128	h2o_ice = h2o_ice + watercap
129	endif
130
131	! co2 ice
132	co2_ice = perennial_co2ice
133	#endif
134
135	if (soil_pem) then
136
137	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
138	!1. Thermal Inertia
139	! a. General case
140	do islope = 1,nslope
141	write(num,'(i2.2)') islope
142	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
143	if (.not. found) then
144	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
145	write(,)'will reconstruct the values of TI_PEM'
146
147	do ig = 1,ngrid
148	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
149	!!! transition
150	delta = depth_breccia
151	TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
152	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
153	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
154	do iloop=index_breccia+2,index_bedrock
155	TI_PEM(ig,iloop,islope) = TI_breccia
156	enddo
157	else ! we keep the high TI values
158	do iloop = index_breccia + 1,index_bedrock
159	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
160	enddo
161	endif ! TI PEM and breccia comparison
162	!! transition
163	delta = depth_bedrock
164	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
165	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
166	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
167	do iloop = index_bedrock + 2,nsoil_PEM
168	TI_PEM(ig,iloop,islope) = TI_bedrock
169	enddo
170	enddo
171	else ! found
172	do iloop = nsoil_PCM + 1,nsoil_PEM
173	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.
174	enddo
175	endif ! not found
176	enddo ! islope
177
178	write(,) 'PEMETAT0: THERMAL INERTIA done'
179
180	! b. Special case for inertiedat, inertiedat_PEM
181	call get_field("inertiedat_PEM",inertiedat_PEM,found)
182	if (.not.found) then
183	do iloop = 1,nsoil_PCM
184	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
185	enddo
186	!!! zone de transition
187	delta = depth_breccia
188	do ig = 1,ngrid
189	if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
190	inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
191	(((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
192	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
193
194	do iloop = index_breccia+2,index_bedrock
195	inertiedat_PEM(ig,iloop) = TI_breccia
196	enddo
197	else
198	do iloop=index_breccia+1,index_bedrock
199	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
200	enddo
201	endif ! comparison ti breccia
202	enddo ! ig
203
204	!!! zone de transition
205	delta = depth_bedrock
206	do ig = 1,ngrid
207	inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
208	(((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
209	((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
210	enddo
211
212	do iloop = index_bedrock + 2, nsoil_PEM
213	do ig = 1,ngrid
214	inertiedat_PEM(ig,iloop) = TI_bedrock
215	enddo
216	enddo
217	endif ! not found
218
219	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
220	!2. Soil Temperature
221	do islope=1,nslope
222	write(num,fmt='(i2.2)') islope
223	call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
224	if (.not.found) then
225	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
226	write(,)'will reconstruct the values of Tsoil'
227	! do ig = 1,ngrid
228	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
229	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
230	! do iloop=index_breccia+2,index_bedrock
231	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
232	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
233	! enddo
234	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
235	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
236	!
237	! do iloop=index_bedrock+2,nsoil_PEM
238	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
239	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
240	! enddo
241	! enddo
242	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
243	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
244	else
245	! predictor corrector: restart from year 1 of the PCM and build the evolution of
246	! tsoil at depth
247	tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
248	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
249	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
250	tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
251	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
252
253	do iloop = nsoil_PCM+1,nsoil_PEM
254	tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
255	enddo
256	endif !found
257
258	do it = 1,timelen
259	do isoil = nsoil_PCM+1,nsoil_PEM
260	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
261	enddo
262	enddo
263	do isoil = nsoil_PCM+1,nsoil_PEM
264	do ig = 1,ngrid
265	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)
266	enddo
267	enddo
268	enddo ! islope
269	write(,) 'PEMETAT0: SOIL TEMP done'
270
271	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
272	!3. Ice Table
273	if(icetable_equilibrium) then
274	call get_field("ice_table",ice_table,found)
275	if (.not. found) then
276	write(,)'PEM settings: failed loading <ice_table>'
277	write(,)'will reconstruct the values of the ice table given the current state'
278	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
279	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
280	do islope = 1,nslope
281	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
282	enddo
283	endif
284	write(,) 'PEMETAT0: ICE TABLE done'
285	endif
286
287	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
288	!4. CO2 & H2O Adsorption
289	if (adsorption_pem) then
290	do islope = 1,nslope
291	write(num,fmt='(i2.2)') islope
292	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
293	if (.not. found) then
294	m_co2_regolith_phys = 0.
295	exit
296	endif
297	enddo
298	do islope=1,nslope
299	write(num,fmt='(i2.2)') islope
300	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
301	if (.not. found2) then
302	m_h2o_regolith_phys = 0.
303	exit
304	endif
305	enddo
306
307	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, &
308	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
309
310	if (.not. found) deltam_co2_regolith_phys = 0.
311	if (.not.found2) deltam_h2o_regolith_phys = 0.
312
313	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
314	endif ! adsorption_pem
315	endif ! soil_pem
316
317	call close_startphy
318
319	else !No startfi, let's build all by hand
320
321	! h2o ice
322	h2o_ice = 0.
323	write(,)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_huge_h2oice''.'
324	do ig = 1,ngrid
325	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
326	enddo
327
328	! co2 ice
329	co2_ice = perennial_co2ice
330
331	if (soil_pem) then
332
333	!a) Thermal inertia
334	do islope = 1,nslope
335	do ig = 1,ngrid
336	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
337	!!! transition
338	delta = depth_breccia
339	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
340	(((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
341	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
342	do iloop=index_breccia + 2,index_bedrock
343	TI_PEM(ig,iloop,islope) = TI_breccia
344	enddo
345	else ! we keep the high ti values
346	do iloop=index_breccia + 1,index_bedrock
347	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
348	enddo
349	endif
350	!! transition
351	delta = depth_bedrock
352	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
353	(((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
354	((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
355	do iloop = index_bedrock + 2,nsoil_PEM
356	TI_PEM(ig,iloop,islope) = TI_bedrock
357	enddo
358	enddo
359	enddo
360
361	do iloop = 1,nsoil_PCM
362	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
363	enddo
364	!!! zone de transition
365	delta = depth_breccia
366	do ig = 1,ngrid
367	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
368	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
369	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
370	((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
371	do iloop = index_breccia + 2,index_bedrock
372	inertiedat_PEM(ig,iloop) = TI_breccia
373	enddo
374	else
375	do iloop = index_breccia + 1,index_bedrock
376	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
377	enddo
378	endif
379	enddo
380
381	!!! zone de transition
382	delta = depth_bedrock
383	do ig = 1,ngrid
384	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
385	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
386	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
387	enddo
388
389	do iloop = index_bedrock + 2,nsoil_PEM
390	do ig = 1,ngrid
391	inertiedat_PEM(ig,iloop) = TI_bedrock
392	enddo
393	enddo
394
395	write(,) 'PEMETAT0: TI done'
396
397	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
398	!b) Soil temperature
399	do islope = 1,nslope
400	! do ig = 1,ngrid
401	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
402	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
403	!
404	! do iloop=index_breccia+2,index_bedrock
405	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
406	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
407	! enddo
408	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
409	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
410	!
411	! do iloop=index_bedrock+2,nsoil_PEM
412	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
413	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
414	! enddo
415	! enddo
416	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
417	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
418
419	! First raw initialization
420	do it = 1,timelen
421	do isoil = nsoil_PCM + 1,nsoil_PEM
422	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
423	enddo
424	enddo
425
426	do it = 1,timelen
427	do isoil = nsoil_PCM + 1,nsoil_PEM
428	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
429	enddo
430	enddo
431
432	do isoil = nsoil_PCM + 1,nsoil_PEM
433	do ig = 1,ngrid
434	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)
435	enddo
436	enddo
437	enddo !islope
438	write(,) 'PEMETAT0: TSOIL done'
439
440	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
441	!c) Ice table
442	if(icetable_equilibrium) then
443	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
444	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
445	do islope = 1,nslope
446	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
447	enddo
448	write(,) 'PEMETAT0: Ice table done'
449	endif
450
451	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452	!d) Regolith adsorbed
453	if (adsorption_pem) then
454	m_co2_regolith_phys = 0.
455	m_h2o_regolith_phys = 0.
456	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, &
457	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
458	deltam_co2_regolith_phys = 0.
459	deltam_h2o_regolith_phys = 0.
460	endif
461
462	write(,) 'PEMETAT0: CO2 adsorption done'
463	endif !soil_pem
464	endif ! of if (startphy_file)
465
466	if (soil_pem) then ! Sanity check
467	do ig = 1,ngrid
468	do islope = 1,nslope
469	do iloop = 1,nsoil_PEM
470	if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
471	enddo
472	enddo
473	enddo
474	endif !soil_pem
475
476	END SUBROUTINE pemetat0
477
478	END MODULE pemetat0_mod

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

Download in other formats:

Original Format