Context Navigation

pemetat0.F90 @ 3296

Last change on this file since 3296 was 3214, checked in by jbclement, 10 months ago

PEM:

It is now possible to set the number of initial PCM calls independently of the number of "inter-PEM" PCM calls. It is useful to get a stable situation for the start of the simulation.
Correction of a bug: 'reshape_XIOS_output' treated the first two PCM runs instead of the last two. The numbering has been adapted in "launch_pem.sh" to get it right.
PEM outputs ("diagpem.nc" and "diagsoil_pem.nc") has been improved: 'Time' now evolves according to 'dt_pem' (usually year by year) and 'ecritpem' is a full variable of "time_evol_mod.F90".

JBC

File size: 25.2 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	write(,)'There is a "startpem.nc"!'
112	! open pem initial state file:
113	call open_startphy(filename)
114
115	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116	#ifndef CPP_STD
117	! h2o ice
118	h2o_ice = 0.
119	call get_field("h2o_ice",h2o_ice,found)
120	if (.not. found) then
121	write(,)'Pemetat0: failed loading <h2o_ice>'
122	write(,)'will reconstruct the values from watercaptag'
123	write(,)'with default value ''ini_huge_h2oice'''
124	do ig = 1,ngrid
125	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
126	enddo
127	else
128	! The variations of infinite reservoirs during the PCM years are taken into account
129	h2o_ice = h2o_ice + watercap
130	endif
131
132	! co2 ice
133	co2_ice = perennial_co2ice
134	#endif
135
136	if (soil_pem) then
137
138	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139	!1. Thermal Inertia
140	! a. General case
141	do islope = 1,nslope
142	write(num,'(i2.2)') islope
143	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
144	if (.not. found) then
145	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
146	write(,)'will reconstruct the values of TI_PEM'
147
148	do ig = 1,ngrid
149	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
150	!!! transition
151	delta = depth_breccia
152	TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
153	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
154	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
155	do iloop=index_breccia+2,index_bedrock
156	TI_PEM(ig,iloop,islope) = TI_breccia
157	enddo
158	else ! we keep the high TI values
159	do iloop = index_breccia + 1,index_bedrock
160	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
161	enddo
162	endif ! TI PEM and breccia comparison
163	!! transition
164	delta = depth_bedrock
165	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
166	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
167	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
168	do iloop = index_bedrock + 2,nsoil_PEM
169	TI_PEM(ig,iloop,islope) = TI_bedrock
170	enddo
171	enddo
172	else ! found
173	do iloop = nsoil_PCM + 1,nsoil_PEM
174	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.
175	enddo
176	endif ! not found
177	enddo ! islope
178
179	write(,) 'PEMETAT0: THERMAL INERTIA done'
180
181	! b. Special case for inertiedat, inertiedat_PEM
182	call get_field("inertiedat_PEM",inertiedat_PEM,found)
183	if (.not.found) then
184	do iloop = 1,nsoil_PCM
185	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
186	enddo
187	!!! zone de transition
188	delta = depth_breccia
189	do ig = 1,ngrid
190	if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
191	inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
192	(((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
193	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
194
195	do iloop = index_breccia+2,index_bedrock
196	inertiedat_PEM(ig,iloop) = TI_breccia
197	enddo
198	else
199	do iloop=index_breccia+1,index_bedrock
200	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
201	enddo
202	endif ! comparison ti breccia
203	enddo ! ig
204
205	!!! zone de transition
206	delta = depth_bedrock
207	do ig = 1,ngrid
208	inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
209	(((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
210	((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
211	enddo
212
213	do iloop = index_bedrock + 2, nsoil_PEM
214	do ig = 1,ngrid
215	inertiedat_PEM(ig,iloop) = TI_bedrock
216	enddo
217	enddo
218	endif ! not found
219
220	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
221	!2. Soil Temperature
222	do islope=1,nslope
223	write(num,fmt='(i2.2)') islope
224	call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
225	if (.not.found) then
226	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
227	write(,)'will reconstruct the values of Tsoil'
228	! do ig = 1,ngrid
229	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
230	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
231	! do iloop=index_breccia+2,index_bedrock
232	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
233	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
234	! enddo
235	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
236	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
237	!
238	! do iloop=index_bedrock+2,nsoil_PEM
239	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
240	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
241	! enddo
242	! enddo
243	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
244	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
245	else
246	! predictor corrector: restart from year 1 of the PCM and build the evolution of
247	! tsoil at depth
248	tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
249	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
250	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
251	tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
252	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
253
254	do iloop = nsoil_PCM+1,nsoil_PEM
255	tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
256	enddo
257	endif !found
258
259	do it = 1,timelen
260	do isoil = nsoil_PCM+1,nsoil_PEM
261	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
262	enddo
263	enddo
264	do isoil = nsoil_PCM+1,nsoil_PEM
265	do ig = 1,ngrid
266	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)
267	enddo
268	enddo
269	enddo ! islope
270	write(,) 'PEMETAT0: SOIL TEMP done'
271
272	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
273	!3. Ice Table
274	if(icetable_equilibrium) then
275	call get_field("ice_table",ice_table,found)
276	if (.not. found) then
277	write(,)'PEM settings: failed loading <ice_table>'
278	write(,)'will reconstruct the values of the ice table given the current state'
279	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
280	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
281	do islope = 1,nslope
282	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
283	enddo
284	endif
285	write(,) 'PEMETAT0: ICE TABLE done'
286	endif
287
288	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
289	!4. CO2 & H2O Adsorption
290	if (adsorption_pem) then
291	do islope = 1,nslope
292	write(num,fmt='(i2.2)') islope
293	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
294	if (.not. found) then
295	m_co2_regolith_phys = 0.
296	exit
297	endif
298	enddo
299	do islope=1,nslope
300	write(num,fmt='(i2.2)') islope
301	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
302	if (.not. found2) then
303	m_h2o_regolith_phys = 0.
304	exit
305	endif
306	enddo
307
308	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, &
309	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
310
311	if (.not. found) deltam_co2_regolith_phys = 0.
312	if (.not.found2) deltam_h2o_regolith_phys = 0.
313
314	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
315	endif ! adsorption_pem
316	endif ! soil_pem
317
318	call close_startphy
319
320	else !No startfi, let's build all by hand
321
322	! h2o ice
323	h2o_ice = 0.
324	write(,)'There is no "startpem.nc"!'
325	write(,)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
326	do ig = 1,ngrid
327	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
328	enddo
329
330	! co2 ice
331	write(,)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
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