Context Navigation

pemetat0.F90 @ 3543

Last change on this file since 3543 was 3537, checked in by jbclement, 2 weeks ago
PEM: Small fixes to initialize and output ice table-related variables and in the launching script. JBC
File size: 28.8 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	! do ig = 1,ngrid
267	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
268	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
269	! do iloop=index_breccia+2,index_bedrock
270	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
271	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
272	! enddo
273	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
274	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
275	!
276	! do iloop=index_bedrock+2,nsoil_PEM
277	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
278	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
279	! enddo
280	! enddo
281	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
282	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
283	else
284	! predictor corrector: restart from year 1 of the PCM and build the evolution of tsoil at depth
285	tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
286	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
287	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
288	tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
289	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
290
291	do iloop = nsoil_PCM + 1,nsoil_PEM
292	tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
293	enddo
294	endif !found
295
296	do it = 1,timelen
297	tsoil_inst(:,nsoil_PCM + 1:nsoil_PEM,islope,it) = tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)
298	enddo
299	do isoil = nsoil_PCM + 1,nsoil_PEM
300	watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)mmol(igcm_h2o_vap)/(mugazr)
301	enddo
302	enddo ! islope
303	write(,) 'PEMETAT0: SOIL TEMP done'
304
305	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
306	!3. Ice Table
307	if (icetable_equilibrium) then
308	call get_field("icetable_depth",icetable_depth,found)
309	if (.not. found) then
310	write(,)'PEM settings: failed loading <icetable_depth>'
311	write(,)'will reconstruct the values of the ice table given the current state'
312	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
313	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)
314	do islope = 1,nslope
315	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
316	enddo
317	endif
318	write(,) 'PEMETAT0: ICE TABLE done'
319	else if (icetable_dynamic) then
320	call get_field("ice_porefilling",ice_porefilling,found)
321	if (.not. found) write(,)'PEM settings: failed loading <ice_porefilling>'
322	call get_field("icetable_depth",icetable_depth,found)
323	if (.not. found) then
324	write(,)'PEM settings: failed loading <icetable_depth>'
325	write(,)'will reconstruct the values of the ice table given the current state'
326	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)
327	do islope = 1,nslope
328	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
329	enddo
330	endif
331	write(,) 'PEMETAT0: ICE TABLE done'
332	endif
333
334	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
335	!4. CO2 & H2O Adsorption
336	if (adsorption_pem) then
337	do islope = 1,nslope
338	write(num,fmt='(i2.2)') islope
339	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
340	if (.not. found) then
341	m_co2_regolith_phys = 0.
342	exit
343	endif
344	enddo
345	do islope=1,nslope
346	write(num,fmt='(i2.2)') islope
347	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
348	if (.not. found2) then
349	m_h2o_regolith_phys = 0.
350	exit
351	endif
352	enddo
353
354	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, &
355	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
356
357	if (.not. found) deltam_co2_regolith_phys = 0.
358	if (.not. found2) deltam_h2o_regolith_phys = 0.
359
360	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
361	endif ! adsorption_pem
362	endif ! soil_pem
363
364	call close_startphy
365
366	else !No startpem, let's build all by hand
367
368	write(,)'There is no "startpem.nc"!'
369
370	! h2o ice
371	h2o_ice = 0.
372	write(,)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
373	do ig = 1,ngrid
374	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
375	enddo
376
377	! co2 ice
378	write(,)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
379	co2_ice = perennial_co2ice
380
381	! Stratification (layerings)
382	nb_str_max = 1
383	if (layering_algo) write(,)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
384
385	if (soil_pem) then
386
387	!a) Thermal inertia
388	do islope = 1,nslope
389	do ig = 1,ngrid
390	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
391	!!! transition
392	delta = depth_breccia
393	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
394	(((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
395	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
396	do iloop=index_breccia + 2,index_bedrock
397	TI_PEM(ig,iloop,islope) = TI_breccia
398	enddo
399	else ! we keep the high ti values
400	do iloop=index_breccia + 1,index_bedrock
401	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
402	enddo
403	endif
404	!! transition
405	delta = depth_bedrock
406	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
407	(((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
408	((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
409	do iloop = index_bedrock + 2,nsoil_PEM
410	TI_PEM(ig,iloop,islope) = TI_bedrock
411	enddo
412	enddo
413	enddo
414
415	do iloop = 1,nsoil_PCM
416	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
417	enddo
418	!!! zone de transition
419	delta = depth_breccia
420	do ig = 1,ngrid
421	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
422	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
423	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
424	((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
425	do iloop = index_breccia + 2,index_bedrock
426	inertiedat_PEM(ig,iloop) = TI_breccia
427	enddo
428	else
429	do iloop = index_breccia + 1,index_bedrock
430	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
431	enddo
432	endif
433	enddo
434
435	!!! zone de transition
436	delta = depth_bedrock
437	do ig = 1,ngrid
438	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
439	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
440	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
441	enddo
442
443	do iloop = index_bedrock + 2,nsoil_PEM
444	do ig = 1,ngrid
445	inertiedat_PEM(ig,iloop) = TI_bedrock
446	enddo
447	enddo
448
449	write(,) 'PEMETAT0: TI done'
450
451	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452	!b) Soil temperature
453	do islope = 1,nslope
454	! do ig = 1,ngrid
455	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
456	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
457	!
458	! do iloop=index_breccia+2,index_bedrock
459	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
460	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
461	! enddo
462	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
463	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
464	!
465	! do iloop=index_bedrock+2,nsoil_PEM
466	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
467	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
468	! enddo
469	! enddo
470	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
471	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
472
473	! First raw initialization
474	do it = 1,timelen
475	do isoil = nsoil_PCM + 1,nsoil_PEM
476	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
477	enddo
478	enddo
479
480	do it = 1,timelen
481	do isoil = nsoil_PCM + 1,nsoil_PEM
482	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
483	enddo
484	enddo
485
486	do isoil = nsoil_PCM + 1,nsoil_PEM
487	do ig = 1,ngrid
488	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)
489	enddo
490	enddo
491	enddo !islope
492	write(,) 'PEMETAT0: TSOIL done'
493
494	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
495	!c) Ice table
496	if (icetable_equilibrium) then
497	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
498	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)
499	do islope = 1,nslope
500	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
501	enddo
502	write(,) 'PEMETAT0: Ice table done'
503	else if (icetable_dynamic) then
504	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)
505	do islope = 1,nslope
506	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
507	enddo
508	write(,) 'PEMETAT0: Ice table done'
509	endif
510
511	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
512	!d) Regolith adsorbed
513	if (adsorption_pem) then
514	m_co2_regolith_phys = 0.
515	m_h2o_regolith_phys = 0.
516	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, &
517	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
518	deltam_co2_regolith_phys = 0.
519	deltam_h2o_regolith_phys = 0.
520	endif
521
522	write(,) 'PEMETAT0: CO2 adsorption done'
523	endif !soil_pem
524	endif ! of if (startphy_file)
525
526	if (soil_pem) then ! Sanity check
527	do ig = 1,ngrid
528	do islope = 1,nslope
529	do iloop = 1,nsoil_PEM
530	if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
531	enddo
532	enddo
533	enddo
534	endif !soil_pem
535
536	END SUBROUTINE pemetat0
537
538	END MODULE pemetat0_mod

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

Download in other formats:

Original Format