Context Navigation

pemetat0.F90 @ 3321

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

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

Download in other formats:

Original Format