Context Navigation

pemetat0.F90 @ 3297

Last change on this file since 3297 was 3297, checked in by jbclement, 8 months ago

PEM:
Integration of the module "layering_mod.F90" with the rest of the PEM:

The linked list data structure representative of layered deposits is converted into an array which can be outputed in the "restratpem.nc" files. This array has dimensions (ngrid,nslope,nb_str_max,6) where 'nb_str_max' is the maximum number of 'stratum' through the layerings and '6' is the number of properties of 'stratum';
this structure can also be read from "startpem.nc" files to initialize PEM runs;
The layering algorithm is now used in the main PEM loop to make the layerings evolve.

JBC

File size: 27.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,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
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
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) :: tend_h2o_ice ! tendencies on h2o ice
52	real, dimension(ngrid,nslope), intent(in) :: tend_co2_ice ! 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) :: ice_table ! Ice table depth [m]
63	real, dimension(ngrid,nslope), intent(inout) :: ice_table_thickness ! Ice table thickness [m]
64	real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst ! instantaneous soil (mid-layer) temperature [K]
65	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
66	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
67	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg ! surface water ice density, yearly averaged (kg/m^3)
68	! local
69	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM ! soil temperature saved in the start [K]
70	real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM ! soil thermal inertia saved in the start [SI]
71	logical :: found, found2 ! check if variables are found in the start
72	integer :: iloop, ig, islope, it, isoil, k ! index for loops
73	real :: kcond ! Thermal conductivity, intermediate variable [SI]
74	real :: delta ! Depth of the interface regolith-breccia, breccia -bedrock [m]
75	character(2) :: num ! intermediate string to read PEM start sloped variables
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	logical :: startpem_file ! boolean to check if we read the startfile or not
79	real, dimension(:,:,:,:), allocatable :: stratif_array ! Array for stratification (layerings)
80	real :: nb_str_max
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	write(,)'There is a "startpem.nc"!'
115	! open pem initial state file:
116	call open_startphy(filename)
117
118	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
119	! Stratification (layerings)
120	found = inquire_dimension("nb_str_max")
121	if (.not. found) then
122	write(,) 'Pemetat0: failed loading <nb_str_max>!'
123	write(,) '''nb_str_max'' is set to 10!'
124	nb_str_max = 10.
125	else
126	nb_str_max = inquire_dimension_length('nb_str_max')
127	endif
128	allocate(stratif_array(ngrid,nslope,int(nb_str_max),6))
129	stratif_array = 0.
130	do islope = 1,nslope
131	write(num,'(i2.2)') islope
132	call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
133	call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
134	found = found .or. found2
135	call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
136	found = found .or. found2
137	call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
138	found = found .or. found2
139	call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
140	found = found .or. found2
141	call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
142	found = found .or. found2
143	if (.not. found) then
144	write(,) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
145	endif ! not found
146	enddo ! islope
147	call array2stratif(stratif_array,ngrid,nslope,stratif)
148	deallocate(stratif_array)
149
150	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151	#ifndef CPP_STD
152	! h2o ice
153	h2o_ice = 0.
154	call get_field("h2o_ice",h2o_ice,found)
155	if (.not. found) then
156	write(,)'Pemetat0: failed loading <h2o_ice>'
157	write(,)'will reconstruct the values from watercaptag'
158	write(,)'with default value ''ini_huge_h2oice'''
159	do ig = 1,ngrid
160	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
161	enddo
162	else
163	! The variations of infinite reservoirs during the PCM years are taken into account
164	h2o_ice = h2o_ice + watercap
165	endif
166
167	! co2 ice
168	co2_ice = perennial_co2ice
169	#endif
170
171	if (soil_pem) then
172
173	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174	!1. Thermal Inertia
175	! a. General case
176	do islope = 1,nslope
177	write(num,'(i2.2)') islope
178	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
179	if (.not. found) then
180	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
181	write(,)'will reconstruct the values of TI_PEM'
182
183	do ig = 1,ngrid
184	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
185	!!! transition
186	delta = depth_breccia
187	TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
188	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
189	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
190	do iloop=index_breccia+2,index_bedrock
191	TI_PEM(ig,iloop,islope) = TI_breccia
192	enddo
193	else ! we keep the high TI values
194	do iloop = index_breccia + 1,index_bedrock
195	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
196	enddo
197	endif ! TI PEM and breccia comparison
198	!! transition
199	delta = depth_bedrock
200	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
201	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
202	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
203	do iloop = index_bedrock + 2,nsoil_PEM
204	TI_PEM(ig,iloop,islope) = TI_bedrock
205	enddo
206	enddo
207	else ! found
208	do iloop = nsoil_PCM + 1,nsoil_PEM
209	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.
210	enddo
211	endif ! not found
212	enddo ! islope
213
214	write(,) 'PEMETAT0: THERMAL INERTIA done'
215
216	! b. Special case for inertiedat, inertiedat_PEM
217	call get_field("inertiedat_PEM",inertiedat_PEM,found)
218	if (.not.found) then
219	do iloop = 1,nsoil_PCM
220	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
221	enddo
222	!!! zone de transition
223	delta = depth_breccia
224	do ig = 1,ngrid
225	if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
226	inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
227	(((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
228	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
229
230	do iloop = index_breccia+2,index_bedrock
231	inertiedat_PEM(ig,iloop) = TI_breccia
232	enddo
233	else
234	do iloop=index_breccia+1,index_bedrock
235	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
236	enddo
237	endif ! comparison ti breccia
238	enddo ! ig
239
240	!!! zone de transition
241	delta = depth_bedrock
242	do ig = 1,ngrid
243	inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
244	(((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
245	((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
246	enddo
247
248	do iloop = index_bedrock + 2, nsoil_PEM
249	do ig = 1,ngrid
250	inertiedat_PEM(ig,iloop) = TI_bedrock
251	enddo
252	enddo
253	endif ! not found
254
255	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
256	!2. Soil Temperature
257	do islope=1,nslope
258	write(num,fmt='(i2.2)') islope
259	call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
260	if (.not.found) then
261	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
262	write(,)'will reconstruct the values of Tsoil'
263	! do ig = 1,ngrid
264	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
265	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
266	! do iloop=index_breccia+2,index_bedrock
267	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
268	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
269	! enddo
270	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
271	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
272	!
273	! do iloop=index_bedrock+2,nsoil_PEM
274	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
275	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
276	! enddo
277	! enddo
278	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
279	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
280	else
281	! predictor corrector: restart from year 1 of the PCM and build the evolution of
282	! tsoil at depth
283	tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
284	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
285	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
286	tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
287	call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
288
289	do iloop = nsoil_PCM+1,nsoil_PEM
290	tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
291	enddo
292	endif !found
293
294	do it = 1,timelen
295	do isoil = nsoil_PCM+1,nsoil_PEM
296	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
297	enddo
298	enddo
299	do isoil = nsoil_PCM+1,nsoil_PEM
300	do ig = 1,ngrid
301	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)
302	enddo
303	enddo
304	enddo ! islope
305	write(,) 'PEMETAT0: SOIL TEMP done'
306
307	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
308	!3. Ice Table
309	if(icetable_equilibrium) then
310	call get_field("ice_table",ice_table,found)
311	if (.not. found) then
312	write(,)'PEM settings: failed loading <ice_table>'
313	write(,)'will reconstruct the values of the ice table given the current state'
314	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
315	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
316	do islope = 1,nslope
317	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
318	enddo
319	endif
320	write(,) 'PEMETAT0: ICE TABLE done'
321	endif
322
323	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
324	!4. CO2 & H2O Adsorption
325	if (adsorption_pem) then
326	do islope = 1,nslope
327	write(num,fmt='(i2.2)') islope
328	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
329	if (.not. found) then
330	m_co2_regolith_phys = 0.
331	exit
332	endif
333	enddo
334	do islope=1,nslope
335	write(num,fmt='(i2.2)') islope
336	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
337	if (.not. found2) then
338	m_h2o_regolith_phys = 0.
339	exit
340	endif
341	enddo
342
343	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, &
344	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
345
346	if (.not. found) deltam_co2_regolith_phys = 0.
347	if (.not.found2) deltam_h2o_regolith_phys = 0.
348
349	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
350	endif ! adsorption_pem
351	endif ! soil_pem
352
353	call close_startphy
354
355	else !No startfi, let's build all by hand
356
357	write(,)'There is no "startpem.nc"!'
358
359	! Stratification (layerings)
360	write(,)'So the stratification (layerings) is initialized with nothing more than the 3 sub-surface strata.'
361	nb_str_max = 3
362
363	! h2o ice
364	h2o_ice = 0.
365	write(,)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
366	do ig = 1,ngrid
367	if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)cos(pidef_slope_mean(:)*180.)
368	enddo
369
370	! co2 ice
371	write(,)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
372	co2_ice = perennial_co2ice
373
374	if (soil_pem) then
375
376	!a) Thermal inertia
377	do islope = 1,nslope
378	do ig = 1,ngrid
379	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
380	!!! transition
381	delta = depth_breccia
382	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
383	(((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
384	((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
385	do iloop=index_breccia + 2,index_bedrock
386	TI_PEM(ig,iloop,islope) = TI_breccia
387	enddo
388	else ! we keep the high ti values
389	do iloop=index_breccia + 1,index_bedrock
390	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
391	enddo
392	endif
393	!! transition
394	delta = depth_bedrock
395	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
396	(((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
397	((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
398	do iloop = index_bedrock + 2,nsoil_PEM
399	TI_PEM(ig,iloop,islope) = TI_bedrock
400	enddo
401	enddo
402	enddo
403
404	do iloop = 1,nsoil_PCM
405	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
406	enddo
407	!!! zone de transition
408	delta = depth_breccia
409	do ig = 1,ngrid
410	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
411	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
412	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
413	((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
414	do iloop = index_breccia + 2,index_bedrock
415	inertiedat_PEM(ig,iloop) = TI_breccia
416	enddo
417	else
418	do iloop = index_breccia + 1,index_bedrock
419	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
420	enddo
421	endif
422	enddo
423
424	!!! zone de transition
425	delta = depth_bedrock
426	do ig = 1,ngrid
427	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
428	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
429	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
430	enddo
431
432	do iloop = index_bedrock + 2,nsoil_PEM
433	do ig = 1,ngrid
434	inertiedat_PEM(ig,iloop) = TI_bedrock
435	enddo
436	enddo
437
438	write(,) 'PEMETAT0: TI done'
439
440	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
441	!b) Soil temperature
442	do islope = 1,nslope
443	! do ig = 1,ngrid
444	! kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
445	! tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
446	!
447	! do iloop=index_breccia+2,index_bedrock
448	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
449	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
450	! enddo
451	! kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
452	! tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
453	!
454	! do iloop=index_bedrock+2,nsoil_PEM
455	! kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
456	! tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
457	! enddo
458	! enddo
459	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
460	call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
461
462	! First raw initialization
463	do it = 1,timelen
464	do isoil = nsoil_PCM + 1,nsoil_PEM
465	tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
466	enddo
467	enddo
468
469	do it = 1,timelen
470	do isoil = nsoil_PCM + 1,nsoil_PEM
471	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
472	enddo
473	enddo
474
475	do isoil = nsoil_PCM + 1,nsoil_PEM
476	do ig = 1,ngrid
477	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)
478	enddo
479	enddo
480	enddo !islope
481	write(,) 'PEMETAT0: TSOIL done'
482
483	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
484	!c) Ice table
485	if(icetable_equilibrium) then
486	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
487	call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
488	do islope = 1,nslope
489	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
490	enddo
491	write(,) 'PEMETAT0: Ice table done'
492	endif
493
494	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
495	!d) Regolith adsorbed
496	if (adsorption_pem) then
497	m_co2_regolith_phys = 0.
498	m_h2o_regolith_phys = 0.
499	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, &
500	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
501	deltam_co2_regolith_phys = 0.
502	deltam_h2o_regolith_phys = 0.
503	endif
504
505	write(,) 'PEMETAT0: CO2 adsorption done'
506	endif !soil_pem
507	endif ! of if (startphy_file)
508
509	if (soil_pem) then ! Sanity check
510	do ig = 1,ngrid
511	do islope = 1,nslope
512	do iloop = 1,nsoil_PEM
513	if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
514	enddo
515	enddo
516	enddo
517	endif !soil_pem
518
519	END SUBROUTINE pemetat0
520
521	END MODULE pemetat0_mod

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

Download in other formats:

Original Format