Context Navigation

pemetat0.F90 @ 3991

Last change on this file since 3991 was 3991, checked in by jbclement, 5 days ago
PEM: Apply documentation template everywhere: standardized headers format with short description, separators between functions/subroutines, normalized code sections, aligned dependencies/arguments/variables declaration. JBC
File size: 28.3 KB

Line
1	MODULE pemetat0_mod
2	!-----------------------------------------------------------------------
3	! NAME
4	! pemetat0_mod
5	!
6	! DESCRIPTION
7	! Initialize PEM state from start files and PCM/XIOS data.
8	!
9	! AUTHORS & DATE
10	! L. Lange
11	! JB Clement, 2023-2025
12	!
13	! NOTES
14	!
15	!-----------------------------------------------------------------------
16
17	! DECLARATION
18	! -----------
19	implicit none
20
21	contains
22	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23
24	!=======================================================================
25	SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
26	tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,d_h2oice,d_co2ice,co2_ice,h2o_ice, &
27	watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,m_h2o_regolith_phys,deltam_h2o_regolith_phys,layerings_map)
28	!-----------------------------------------------------------------------
29	! NAME
30	! pemetat0
31	!
32	! DESCRIPTION
33	! Read PEM start file (if present) and initialize model state.
34	! Reconstructs when missing.
35	!
36	! AUTHORS & DATE
37	! L. Lange
38	! JB Clement, 2023-2025
39	!
40	! NOTES
41	! Order of operations (must be respected):
42	! 0) Previous-year context and ice fields
43	! 1) Thermal inertia
44	! 2) Soil temperature
45	! 3) Ice table (equilibrium or dynamic)
46	! 4) Regolith adsorption (CO2 & H2O)
47	!-----------------------------------------------------------------------
48
49	! DEPENDENCIES
50	! ------------
51	use iostart_pem, only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
52	use soil, only: do_soil, layer_PEM, fluxgeo, inertiedat_PEM, h2oice_huge_ini, depth_breccia, depth_bedrock, index_breccia, index_bedrock
53	use comsoil_h, only: inertiedat
54	use sorption, only: regolith_adsorption, do_sorption
55	use ice_table, only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
56	use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
57	use soil_thermal_inertia, only: update_soil_TI
58	use stop_pem, only: stop_clean
59	use soil_temp, only: ini_tsoil_pem, compute_tsoil
60	use layered_deposits, only: layering, array2map, nb_str_max, layering_algo, add_stratum, ini_layering
61	use callkeys_mod, only: startphy_file
62	use surf_ice, only: rho_co2ice, rho_h2oice
63	use phys_constants, only: pi, r, mugaz
64	use surfdat_h, only: watercaptag, perennial_co2ice, qsurf
65	use metamorphism, only: h2o_frost4PCM, co2_frost4PCM, iPCM_h2ofrost, iPCM_co2frost
66	use tracers, only: mmol, iPCM_qh2o
67
68	! DECLARATION
69	! -----------
70	implicit none
71
72	! ARGUMENTS
73	! ---------
74	character(*), intent(in) :: filename ! Name of the startfi_PEM.nc
75	integer, intent(in) :: ngrid ! # of physical grid points
76	integer, intent(in) :: nsoil_PCM ! # of vertical grid points in the PCM
77	integer, intent(in) :: nsoil_PEM ! # of vertical grid points in the PEM
78	integer, intent(in) :: nslope ! # of sub-grid slopes
79	integer, intent(in) :: timelen ! # time samples
80	real, intent(in) :: timestep ! Time step [s]
81	real, intent(in) :: ps_avg_global ! Mean average pressure on the planet [Pa]
82	real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr1 ! Surface temperature at the first year of PCM call [K]
83	real, dimension(ngrid,nslope), intent(in) :: tsurf_avg_yr2 ! Surface temperature at the second year of PCM call [K]
84	real, dimension(ngrid,timelen), intent(in) :: q_co2 ! MMR tracer co2 [kg/kg]
85	real, dimension(ngrid,timelen), intent(in) :: q_h2o ! MMR tracer h2o [kg/kg]
86	real, dimension(ngrid,timelen), intent(in) :: ps_timeseries ! Surface pressure [Pa]
87	real, dimension(ngrid,nslope), intent(in) :: d_h2oice ! Tendencies on h2o ice
88	real, dimension(ngrid,nslope), intent(in) :: d_co2ice ! Tendencies on co2 ice
89	real, dimension(ngrid,nslope), intent(in) :: watersurf_avg ! Surface water ice density, yearly averaged (kg/m^3)
90	real, dimension(ngrid), intent(out) :: deltam_co2_regolith_phys ! Mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
91	real, dimension(ngrid), intent(out) :: deltam_h2o_regolith_phys ! Mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
92	real, dimension(ngrid,nslope), intent(out) :: h2o_ice ! H2O ice amount [kg/m^2]
93	real, dimension(ngrid,nslope), intent(out) :: co2_ice ! CO2 ice amount [kg/m^2]
94	type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map ! Layerings
95	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Soil (mid-layer) thermal inertia in the PEM grid [SI]
96	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM ! Soil (mid-layer) temperature [K]
97	real, dimension(ngrid,nslope), intent(inout) :: icetable_depth ! Ice table depth [m]
98	real, dimension(ngrid,nslope), intent(inout) :: icetable_thickness ! Ice table thickness [m]
99	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling ! Subsurface ice pore filling [m3/m3]
100	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! Mass of co2 adsorbed [kg/m^2]
101	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! Mass of h2o adsorbed [kg/m^2]
102	real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg ! Surface water ice density, yearly averaged (kg/m^3)
103
104	! LOCAL VARIABLES
105	! ---------------
106	real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem ! Soil temperature saved in the start [K]
107	real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM ! Soil thermal inertia saved in the start [SI]
108	logical :: found, found2 ! Check if variables are found in the start
109	integer :: iloop, ig, islope ! Index for loops
110	real :: delta ! Depth of the interface regolith-breccia, breccia -bedrock [m]
111	character(2) :: num ! Intermediate string to read PEM start sloped variables
112	logical :: startpem_file ! Boolean to check if we read the startfile or not
113	real, dimension(:,:,:,:), allocatable :: layerings_array ! Array for layerings
114
115	! CODE
116	! ----
117	!0.1 Check if the start_PEM exist.
118
119	inquire(file = filename,exist = startpem_file)
120
121	write(,)'Is there a "startpem" file?',startpem_file
122
123	!0.2 Set to default values
124	icetable_depth = -1. ! by default, no ice table
125	icetable_thickness = -1.
126	ice_porefilling = 0.
127	!1. Run
128	if (startpem_file) then
129	write(,)'There is a "startpem.nc"!'
130	! open pem initial state file:
131	call open_startphy(filename)
132
133	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
134	! H2O ice
135	h2o_ice(:,:) = 0.
136	call get_field("h2o_ice",h2o_ice,found)
137	if (.not. found) then
138	write(,)'Pemetat0: failed loading <h2o_ice>'
139	write(,)'will reconstruct the values from watercaptag'
140	write(,)'with default value ''h2oice_huge_ini'''
141	do ig = 1,ngrid
142	if (watercaptag(ig)) h2o_ice(ig,:) = h2oice_huge_ini + qsurf(ig,iPCM_h2ofrost,:) - h2o_frost4PCM(ig,:)
143	enddo
144	endif
145
146	! CO2 ice
147	co2_ice(:,:) = 0.
148	call get_field("co2_ice",co2_ice,found)
149	if (.not. found) then
150	write(,)'Pemetat0: failed loading <co2_ice>'
151	write(,)'will reconstruct the values from perennial_co2ice'
152	co2_ice(:,:) = perennial_co2ice(:,:) + qsurf(:,iPCM_co2frost,:) - co2_frost4PCM(:,:)
153	endif
154
155	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
156	! Layerings
157	nb_str_max = 68
158	if (layering_algo) then
159	found = inquire_dimension("nb_str_max")
160	if (.not. found) then
161	write(,) 'Pemetat0: failed loading <nb_str_max>!'
162	write(,) '''nb_str_max'' is set to 68 (sub-surface layers)!'
163	else
164	nb_str_max = int(inquire_dimension_length('nb_str_max'))
165	endif
166	allocate(layerings_array(ngrid,nslope,nb_str_max,6))
167	layerings_array = 0.
168	do islope = 1,nslope
169	write(num,'(i2.2)') islope
170	call get_field('stratif_slope'//num//'_top_elevation',layerings_array(:,islope,:,1),found)
171	if (.not. found) then
172	write(,) 'Pemetat0: failed loading <stratif_slope'//num//'_top_elevation>!'
173	exit
174	endif
175	call get_field('stratif_slope'//num//'_h_co2ice',layerings_array(:,islope,:,2),found)
176	if (.not. found) then
177	write(,) 'Pemetat0: failed loading <stratif_slope'//num//'_h_co2ice>!'
178	exit
179	endif
180	call get_field('stratif_slope'//num//'_h_h2oice',layerings_array(:,islope,:,3),found)
181	if (.not. found) then
182	write(,) 'Pemetat0: failed loading <stratif_slope'//num//'_h_h2oice>!'
183	exit
184	endif
185	call get_field('stratif_slope'//num//'_h_dust',layerings_array(:,islope,:,4),found)
186	if (.not. found) then
187	write(,) 'Pemetat0: failed loading <stratif_slope'//num//'_h_dust>!'
188	exit
189	endif
190	call get_field('stratif_slope'//num//'_h_pore',layerings_array(:,islope,:,5),found)
191	if (.not. found) then
192	write(,) 'Pemetat0: failed loading <stratif_slope'//num//'_h_pore>!'
193	exit
194	endif
195	call get_field('stratif_slope'//num//'_poreice_volfrac',layerings_array(:,islope,:,6),found)
196	if (.not. found) then
197	write(,) 'Pemetat0: failed loading <stratif_slope'//num//'_poreice_volfrac>!'
198	exit
199	endif
200	enddo ! islope
201	if (.not. found) then
202	write(,) 'So the layerings are initialized with sub-surface strata.'
203	write(,) 'Ice is added with ''h2oice_huge_ini'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.'
204	do ig = 1,ngrid
205	if (watercaptag(ig)) then
206	do islope = 1,nslope
207	call ini_layering(layerings_map(ig,islope))
208	call add_stratum(layerings_map(ig,islope),1.05h2oice_huge_ini/rho_h2oice,0.,h2oice_huge_ini/rho_h2oice,0.05h2oice_huge_ini/rho_h2oice,0.,0.)
209	enddo
210	else
211	do islope = 1,nslope
212	call ini_layering(layerings_map(ig,islope))
213	if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),1.05perennial_co2ice(ig,islope)/rho_co2ice,perennial_co2ice(ig,islope)/rho_co2ice,0.,0.05perennial_co2ice(ig,islope)/rho_co2ice,0.,0.)
214	enddo
215	endif
216	enddo
217	else
218	call array2map(layerings_array,ngrid,nslope,layerings_map)
219	endif
220	deallocate(layerings_array)
221	endif
222
223	if (do_soil) then
224
225	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
226	!1. Thermal Inertia
227	! a. General case
228	do islope = 1,nslope
229	write(num,'(i2.2)') islope
230	call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
231	if (.not. found) then
232	write(,)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
233	write(,)'will reconstruct the values of TI_PEM'
234
235	do ig = 1,ngrid
236	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
237	!!! transition
238	delta = depth_breccia
239	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
240	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
241	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
242	do iloop = index_breccia + 2,index_bedrock
243	TI_PEM(ig,iloop,islope) = TI_breccia
244	enddo
245	else ! we keep the high TI values
246	do iloop = index_breccia + 1,index_bedrock
247	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
248	enddo
249	endif ! TI PEM and breccia comparison
250	!! transition
251	delta = depth_bedrock
252	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
253	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
254	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
255	do iloop = index_bedrock + 2,nsoil_PEM
256	TI_PEM(ig,iloop,islope) = TI_bedrock
257	enddo
258	enddo
259	else ! found
260	do iloop = nsoil_PCM + 1,nsoil_PEM
261	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.
262	enddo
263	endif ! not found
264	enddo ! islope
265
266	write(,) 'PEMETAT0: TI done'
267
268	! b. Special case for inertiedat, inertiedat_PEM
269	call get_field("inertiedat_PEM",inertiedat_PEM,found)
270	if (.not. found) then
271	do iloop = 1,nsoil_PCM
272	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
273	enddo
274	!!! zone de transition
275	delta = depth_breccia
276	do ig = 1,ngrid
277	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
278	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
279	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
280	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
281
282	do iloop = index_breccia + 2,index_bedrock
283	inertiedat_PEM(ig,iloop) = TI_breccia
284	enddo
285	else
286	do iloop = index_breccia + 1,index_bedrock
287	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
288	enddo
289	endif ! comparison ti breccia
290	enddo ! ig
291
292	!!! zone de transition
293	delta = depth_bedrock
294	do ig = 1,ngrid
295	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
296	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
297	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
298	enddo
299
300	do iloop = index_bedrock + 2,nsoil_PEM
301	do ig = 1,ngrid
302	inertiedat_PEM(ig,iloop) = TI_bedrock
303	enddo
304	enddo
305	endif ! not found
306
307	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
308	!2. Soil Temperature
309	do islope = 1,nslope
310	write(num,fmt='(i2.2)') islope
311	call get_field("tsoil_PEM_slope"//num,tsoil_startpem(:,:,islope),found)
312	if (.not. found) then
313	write(,)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
314	write(,)'will reconstruct the values of Tsoil'
315	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
316	call compute_tsoil(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
317	else
318	! predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startpem is adapted firstly
319	! for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth
320	call compute_tsoil(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
321	call compute_tsoil(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
322	call compute_tsoil(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
323
324	tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) = tsoil_startpem(:,nsoil_PCM + 1:nsoil_PEM,islope)
325	endif !found
326
327	watersoil_avg(:,nsoil_PCM + 1:nsoil_PEM,islope) = exp(beta_clap_h2o/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) + alpha_clap_h2o)/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)mmol(iPCM_qh2o)/(mugazr)
328	enddo ! islope
329	write(,) 'PEMETAT0: TSOIL done'
330
331	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
332	!3. Ice Table
333	if (icetable_equilibrium) then
334	call get_field("icetable_depth",icetable_depth,found)
335	if (.not. found) then
336	write(,)'PEM settings: failed loading <icetable_depth>'
337	write(,)'will reconstruct the values of the ice table given the current state'
338	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
339	call update_soil_TI(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
340	do islope = 1,nslope
341	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
342	enddo
343	endif
344	write(,) 'PEMETAT0: ICE TABLE done'
345	else if (icetable_dynamic) then
346	call get_field("ice_porefilling",ice_porefilling,found)
347	if (.not. found) write(,)'PEM settings: failed loading <ice_porefilling>'
348	call get_field("icetable_depth",icetable_depth,found)
349	if (.not. found) then
350	write(,)'PEM settings: failed loading <icetable_depth>'
351	write(,)'will reconstruct the values of the ice table given the current state'
352	call update_soil_TI(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
353	do islope = 1,nslope
354	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
355	enddo
356	endif
357	write(,) 'PEMETAT0: ICE TABLE done'
358	endif
359
360	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
361	!4. CO2 & H2O Adsorption
362	if (do_sorption) then
363	do islope = 1,nslope
364	write(num,fmt='(i2.2)') islope
365	call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
366	if (.not. found) then
367	m_co2_regolith_phys = 0.
368	exit
369	endif
370	enddo
371	do islope = 1,nslope
372	write(num,fmt='(i2.2)') islope
373	call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
374	if (.not. found2) then
375	m_h2o_regolith_phys = 0.
376	exit
377	endif
378	enddo
379
380	call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
381	m_h2o_regolith_phys,deltam_h2o_regolith_phys,m_co2_regolith_phys,deltam_co2_regolith_phys)
382
383	if (.not. found) deltam_co2_regolith_phys = 0.
384	if (.not. found2) deltam_h2o_regolith_phys = 0.
385
386	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
387	endif ! do_sorption
388	endif ! do_soil
389
390	call close_startphy
391
392	else ! No startpem, let's build all by hand
393
394	write(,)'There is no "startpem.nc"!'
395
396	! h2o ice
397	h2o_ice(:,:) = 0.
398	write(,)'So ''h2o_ice'' is initialized with default value ''h2oice_huge_ini'' where ''watercaptag'' is true.'
399	do ig = 1,ngrid
400	if (watercaptag(ig)) h2o_ice(ig,:) = h2oice_huge_ini + qsurf(ig,iPCM_h2ofrost,:) - h2o_frost4PCM(ig,:)
401	enddo
402
403	! co2 ice
404	write(,)'So ''co2_ice'' is initialized with ''perennial_co2ice'' found in the PCM.'
405	co2_ice(:,:) = perennial_co2ice(:,:) + qsurf(:,iPCM_co2frost,:) - co2_frost4PCM(:,:)
406
407	! Layerings
408	nb_str_max = 68
409	if (layering_algo) then
410	write(,)'So the layerings are initialized with sub-surface strata.'
411	write(,)'Ice is added with ''h2oice_huge_ini'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.'
412	do ig = 1,ngrid
413	if (watercaptag(ig)) then
414	do islope = 1,nslope
415	call ini_layering(layerings_map(ig,islope))
416	call add_stratum(layerings_map(ig,islope),1.05h2oice_huge_ini/rho_h2oice,0.,h2oice_huge_ini/rho_h2oice,0.05h2oice_huge_ini/rho_h2oice,0.,0.)
417	enddo
418	else
419	do islope = 1,nslope
420	call ini_layering(layerings_map(ig,islope))
421	if (perennial_co2ice(ig,islope) > 0.) call add_stratum(layerings_map(ig,islope),1.05perennial_co2ice(ig,islope)/rho_co2ice,perennial_co2ice(ig,islope)/rho_co2ice,0.,0.05perennial_co2ice(ig,islope)/rho_co2ice,0.,0.)
422	enddo
423	endif
424	enddo
425	endif
426
427	if (do_soil) then
428
429	!a) Thermal inertia
430	do islope = 1,nslope
431	do ig = 1,ngrid
432	if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
433	!!! transition
434	delta = depth_breccia
435	TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
436	(((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
437	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
438	do iloop = index_breccia + 2,index_bedrock
439	TI_PEM(ig,iloop,islope) = TI_breccia
440	enddo
441	else ! we keep the high ti values
442	do iloop = index_breccia + 1,index_bedrock
443	TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
444	enddo
445	endif
446	!! transition
447	delta = depth_bedrock
448	TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
449	(((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
450	((layer_PEM(index_bedrock + 1) - delta)/(TI_breccia**2))))
451	do iloop = index_bedrock + 2,nsoil_PEM
452	TI_PEM(ig,iloop,islope) = TI_bedrock
453	enddo
454	enddo
455	enddo
456
457	do iloop = 1,nsoil_PCM
458	inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
459	enddo
460	!!! zone de transition
461	delta = depth_breccia
462	do ig = 1,ngrid
463	if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
464	inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
465	(((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
466	((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
467	do iloop = index_breccia + 2,index_bedrock
468	inertiedat_PEM(ig,iloop) = TI_breccia
469	enddo
470	else
471	do iloop = index_breccia + 1,index_bedrock
472	inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
473	enddo
474	endif
475	enddo
476
477	!!! zone de transition
478	delta = depth_bedrock
479	do ig = 1,ngrid
480	inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
481	(((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
482	((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
483	enddo
484
485	do iloop = index_bedrock + 2,nsoil_PEM
486	do ig = 1,ngrid
487	inertiedat_PEM(ig,iloop) = TI_bedrock
488	enddo
489	enddo
490
491	write(,) 'PEMETAT0: TI done'
492
493	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
494	!b) Soil temperature
495	do islope = 1,nslope
496	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
497	call compute_tsoil(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
498
499	! First raw initialization
500	watersoil_avg(:,nsoil_PCM + 1:nsoil_PEM,islope) = exp(beta_clap_h2o/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope) + alpha_clap_h2o)/tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)*mmol(iPCM_qh2o)/(mugazr)
501	enddo !islope
502	write(,) 'PEMETAT0: TSOIL done'
503
504	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
505	!c) Ice table
506	if (icetable_equilibrium) then
507	call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
508	call update_soil_TI(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
509	do islope = 1,nslope
510	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
511	enddo
512	write(,) 'PEMETAT0: ICE TABLE done'
513	else if (icetable_dynamic) then
514	call update_soil_TI(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
515	do islope = 1,nslope
516	call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
517	enddo
518	write(,) 'PEMETAT0: ICE TABLE done'
519	endif
520
521	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
522	!d) Regolith adsorbed
523	if (do_sorption) then
524	m_co2_regolith_phys = 0.
525	m_h2o_regolith_phys = 0.
526	call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
527	m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
528	deltam_co2_regolith_phys = 0.
529	deltam_h2o_regolith_phys = 0.
530	endif
531
532	write(,) 'PEMETAT0: CO2 & H2O adsorption done'
533	endif !do_soil
534	endif ! of if (startphy_file)
535
536	if (do_soil) then ! Sanity check
537	do ig = 1,ngrid
538	do islope = 1,nslope
539	do iloop = 1,nsoil_PEM
540	if (isnan(tsoil_PEM(ig,iloop,islope))) call stop_clean("PEM - pemetat0","NaN detected in Tsoil",1)
541	enddo
542	enddo
543	enddo
544	endif !do_soil
545
546	END SUBROUTINE pemetat0
547	!========================================================================
548
549	END MODULE pemetat0_mod

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

Download in other formats:

Original Format