Context Navigation

ini_soil_mod.F90 @ 2849

Last change on this file since 2849 was 2842, checked in by romain.vande, 3 years ago
Mars PEM: PEM is now adapted to run with XIOS diurnal averages (when they will work properly) RV
File size: 11.8 KB

Line
1	MODULE ini_soil_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine ini_icetable(timelen,ngrid,nsoil_PEM, &
8	therm_i, timestep,tsurf_ave,tsoil_ave,tsurf_inst, tsoil_inst,q_co2,q_h2o,ps,ice_table)
9
10	#ifndef CPP_STD
11
12	use vertical_layers_mod, only: ap,bp
13	use comsoil_h_PEM, only: fluxgeo,layer_PEM,inertiedat_PEM
14	use comsoil_h,only: volcapa, nsoilmx
15	#endif
16
17	implicit none
18
19	!-----------------------------------------------------------------------
20	! Author: LL
21	! Purpose: Compute soil temperature using an implict 1st order scheme
22	!
23	! Note: depths of layers and mid-layers, soil thermal inertia and
24	! heat capacity are commons in comsoil_PEM.h
25	! A convergence loop is added until equilibrium
26	!-----------------------------------------------------------------------
27
28	#include "dimensions.h"
29
30	!-----------------------------------------------------------------------
31	! arguments
32	! ---------
33	! inputs:
34	integer,intent(in) :: timelen ! Time length in for time-series data
35	integer,intent(in) :: ngrid ! number of (horizontal) grid-points
36	integer,intent(in) :: nsoil_PEM ! number of soil layers in the PEM
37
38	real,intent(in) :: timestep ! time step
39	real,intent(in) :: tsurf_ave(ngrid) ! surface temperature
40
41	real,intent(in) :: q_co2(ngrid,timelen) ! MMR tracer co2 [kg/kg]
42	real,intent(in) :: q_h2o(ngrid,timelen) ! MMR tracer h2o [kg/kg]
43	real,intent(in) :: ps(ngrid,timelen) ! surface pressure [Pa]
44	real,intent(in) :: tsurf_inst(ngrid,timelen) ! soil (mid-layer) temperature
45
46	! outputs:
47	real,intent(inout) :: tsoil_ave(ngrid,nsoil_PEM) ! soil (mid-layer) temperature.
48	real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,timelen) ! soil (mid-layer) temperature
49	real,intent(out) :: ice_table(ngrid) ! ice table [m]
50	real,intent(inout) :: therm_i(ngrid,nsoil_PEM) ! thermal inertia
51
52	! local variables:
53	integer ig,isoil,it,k,iref
54	REAL :: error_depth = 2.
55	REAL :: tsoil_saved(nsoil_PEM)
56	integer :: countmax = 20
57	integer :: countloop
58	REAL :: tol_error = 0.1
59	REAL :: ice_inertia = 2120.
60	REAL :: alph_PEM(nsoil_PEM-1)
61	REAL :: beta_PEM(nsoil_PEM-1)
62	real :: rhoc(nsoil_PEM)
63	real :: volcapa_ice = 1.43e7
64
65	real :: k_soil(nsoil_PEM)
66	real :: d1(nsoil_PEM),d2(nsoil_PEM),d3(nsoil_PEM),re(nsoil_PEM)
67	real :: Tcol_saved(nsoil_PEM)
68	real :: tsoil_prev(nsoil_PEM)
69	real :: tsurf_prev
70	real :: icedepth_prev
71
72	real :: m_h2o = 18.01528E-3
73	real :: m_co2 = 44.01E-3
74	real :: m_noco2 = 33.37E-3
75	real :: A,B,z1,z2
76	real :: alpha = -6143.7
77	real :: beta = 29.9074
78
79	real,allocatable :: mass_mean(:) ! mean mass above the surface
80	real,allocatable :: zplev(:) ! pressure above the surface
81	real,allocatable :: pvapor(:) ! partial pressure above the surface
82
83	real,allocatable :: rhovapor(:)
84	real :: rhovapor_avg ! mean vapor_density at the surface yearly averaged
85
86	real :: psv_surf
87	real,allocatable :: rho_soil(:) ! water vapor in the soil
88	real,allocatable :: rho_soil_avg(:) ! water vapor in the soil yearly averaged
89
90	real,allocatable :: diff_rho(:) ! difference of vapor content
91
92	#ifndef CPP_STD
93
94	A =(1/m_co2 - 1/m_noco2)
95	B=1/m_noco2
96
97	ice_table(:) = 1.
98	do ig = 1,ngrid
99	error_depth = 1.
100	countloop = 0
101
102	do while(( error_depth.gt.tol_error).and.(countloop.lt.countmax).and.(ice_table(ig).gt.-1e-20))
103
104	countloop = countloop +1
105	Tcol_saved(:) = tsoil_ave(ig,:)
106
107	!2. Compute ice table
108
109	! 2.1 Compute water density at the surface, yearly averaged
110	allocate(mass_mean(timelen))
111	! 1.1 Compute the partial pressure of vapor
112	! a. the molecular mass into the column
113	mass_mean(:) = 1/(A*q_co2(ig,:) +B)
114	! b. pressure level
115	allocate(zplev(timelen))
116	do it = 1,timelen
117	zplev(it) = ap(1) + bp(1)*ps(ig,it)
118	enddo
119
120	! c. Vapor pressure
121	allocate(pvapor(timelen))
122	pvapor(:) = mass_mean(:)/m_h2oq_h2o(ig,:)zplev(:)
123	deallocate(zplev)
124	deallocate(mass_mean)
125
126	! d! Check if there is frost at the surface and then compute the density
127	! at the surface
128	allocate(rhovapor(timelen))
129
130	do it = 1,timelen
131	psv_surf = exp(alpha/tsurf_inst(ig,it) +beta)
132	rhovapor(it) = min(psv_surf,pvapor(it))/tsurf_inst(ig,it)
133	enddo
134	deallocate(pvapor)
135	rhovapor_avg =SUM(rhovapor(:),1)/timelen
136	deallocate(rhovapor)
137
138	! 2.2 Compute water density at the soil layer, yearly averaged
139
140	allocate(rho_soil(timelen))
141	allocate(rho_soil_avg(nsoil_PEM))
142
143	do isoil = 1,nsoil_PEM
144	do it = 1,timelen
145	rho_soil(it) = exp(alpha/tsoil_inst(ig,isoil,it) +beta)/tsoil_inst(ig,isoil,it)
146	enddo
147	rho_soil_avg(isoil) = SUM(rho_soil(:),1)/timelen
148	enddo
149	deallocate(rho_soil)
150
151	!2.3 Final: compute ice table
152	icedepth_prev = ice_table(ig)
153	ice_table(ig) = -1
154	allocate(diff_rho(nsoil_PEM))
155	do isoil = 1,nsoil_PEM
156	diff_rho(isoil) = rhovapor_avg - rho_soil_avg(isoil)
157	enddo
158	deallocate(rho_soil_avg)
159
160	if(diff_rho(1) > 0) then
161	ice_table(ig) = 0.
162	else
163	do isoil = 1,nsoil_PEM -1
164	if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
165	z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
166	z2 = -layer_PEM(isoil+1)*z1 + diff_rho(isoil+1)
167	ice_table(ig) = -z2/z1
168	exit
169	endif
170	enddo
171	endif
172
173	deallocate(diff_rho)
174
175	!3. Update Soil Thermal Inertia
176
177	if (ice_table(ig).gt. 0.) then
178	if (ice_table(ig).lt. 1e-10) then
179	do isoil = 1,nsoil_PEM
180	therm_i(ig,isoil)=ice_inertia
181	enddo
182	else
183	! 4.1 find the index of the mixed layer
184	iref=0 ! initialize iref
185	do k=1,nsoil_PEM ! loop on layers
186	if (ice_table(ig).ge.layer_PEM(k)) then
187	iref=k ! pure regolith layer up to here
188	else
189	! correct iref was obtained in previous cycle
190	exit
191	endif
192	enddo
193
194	! 4.2 Build the new ti
195	do isoil=1,iref
196	therm_i(ig,isoil) =inertiedat_PEM(ig,isoil)
197	enddo
198
199	if (iref.lt.nsoil_PEM) then
200	if (iref.ne.0) then
201	! mixed layer
202	therm_i(ig,iref+1)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
203	(((ice_table(ig)-layer_PEM(iref))/(inertiedat_PEM(ig,iref)**2))+ &
204	((layer_PEM(iref+1)-ice_table(ig))/(ice_inertia**2))))
205
206	else ! first layer is already a mixed layer
207	! (ie: take layer(iref=0)=0)
208	therm_i(ig,1)=sqrt((layer_PEM(1))/ &
209	(((ice_table(ig))/(inertiedat_PEM(ig,1)**2))+ &
210	((layer_PEM(1)-ice_table(ig))/(ice_inertia**2))))
211	endif ! of if (iref.ne.0)
212	! lower layers of pure ice
213	do isoil=iref+2,nsoil_PEM
214	therm_i(ig,isoil)=ice_inertia
215	enddo
216	endif ! of if (iref.lt.(nsoilmx))
217	endif ! permanent glaciers
218
219	call soil_pem_1D(nsoil_PEM,.true.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM)
220	call soil_pem_1D(nsoil_PEM,.false.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM)
221
222	do it = 1,timelen
223	tsoil_inst(ig,:,it) = tsoil_inst(ig,:,it) - (Tcol_saved(:) - tsoil_ave(ig,:))
224	enddo
225
226	error_depth = abs(icedepth_prev - ice_table(ig))
227	endif
228
229	enddo
230
231	error_depth = 1.
232	countloop = 0
233	enddo
234
235	#endif
236
237	END SUBROUTINE ini_icetable
238	subroutine soil_pem_1D(nsoil,firstcall, &
239	therm_i, &
240	timestep,tsurf,tsoil,alph,beta)
241	#ifndef CPP_STD
242	use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, &
243	mu_PEM,fluxgeo
244	use comsoil_h,only: volcapa
245	#endif
246	implicit none
247
248	!-----------------------------------------------------------------------
249	! Author: LL
250	! Purpose: Compute soil temperature using an implict 1st order scheme
251	!
252	! Note: depths of layers and mid-layers, soil thermal inertia and
253	! heat capacity are commons in comsoil_PEM.h
254	! A convergence loop is added until equilibrium
255	!-----------------------------------------------------------------------
256
257	#include "dimensions.h"
258	!#include "dimphys.h"
259
260	!#include"comsoil.h"
261
262
263	!-----------------------------------------------------------------------
264	! arguments
265	! ---------
266	! inputs:
267	integer,intent(in) :: nsoil ! number of soil layers
268	logical,intent(in) :: firstcall ! identifier for initialization call
269	real,intent(in) :: therm_i(nsoil) ! thermal inertia
270	real,intent(in) :: timestep ! time step
271	real,intent(in) :: tsurf ! surface temperature
272
273	! outputs:
274	real,intent(inout) :: tsoil(nsoil) ! soil (mid-layer) temperature
275	real,intent(inout) :: alph(nsoil-1)
276	real,intent(inout) :: beta(nsoil-1)
277
278	! local variables:
279	integer ik
280	real :: thermdiff_PEM(nsoil-1)
281	real :: mthermdiff_PEM(0:nsoil-1)
282	real :: coefd_PEM(nsoil-1)
283	real :: coefq_PEM(0:nsoil-1)
284
285	#ifndef CPP_STD
286
287	! 0. Initialisations and preprocessing step
288	if (firstcall) then
289
290	! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
291	do ik=0,nsoil-1
292	mthermdiff_PEM(ik)=therm_i(ik+1)*therm_i(ik+1)/volcapa
293
294	enddo
295
296	! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
297	do ik=1,nsoil-1
298	thermdiff_PEM(ik)=((layer_PEM(ik)-mlayer_PEM(ik-1))*mthermdiff_PEM(ik) &
299	+(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ik-1)) &
300	/(mlayer_PEM(ik)-mlayer_PEM(ik-1))
301	enddo
302
303	! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alpha_k and capcal
304	! mu_PEM
305	mu_PEM=mlayer_PEM(0)/(mlayer_PEM(1)-mlayer_PEM(0))
306
307	! q_{1/2}
308	coefq_PEM(0)=volcapa*layer_PEM(1)/timestep
309	! q_{k+1/2}
310	do ik=1,nsoil-1
311	coefq_PEM(ik)=volcapa*(layer_PEM(ik+1)-layer_PEM(ik)) &
312	/timestep
313	enddo
314
315	! d_k
316	do ik=1,nsoil-1
317	coefd_PEM(ik)=thermdiff_PEM(ik)/(mlayer_PEM(ik)-mlayer_PEM(ik-1))
318	enddo
319
320	! alph_PEM_{N-1}
321	alph(nsoil-1)=coefd_PEM(nsoil-1)/ &
322	(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1))
323	! alph_PEM_k
324	do ik=nsoil-2,1,-1
325	alph(ik)=coefd_PEM(ik)/(coefq_PEM(ik)+coefd_PEM(ik+1)* &
326	(1.-alph(ik+1))+coefd_PEM(ik))
327	enddo
328
329	endif ! of if (firstcall)
330
331	IF (.not.firstcall) THEN
332	! 2. Compute soil temperatures
333	! First layer:
334	tsoil(1)=(tsurf+mu_PEMbeta(1) &
335	thermdiff_PEM(1)/mthermdiff_PEM(0))/ &
336	(1.+mu_PEM(1.0-alph(1))&
337	thermdiff_PEM(1)/mthermdiff_PEM(0))
338
339	! Other layers:
340	do ik=1,nsoil-1
341	tsoil(ik+1)=alph(ik)*tsoil(ik)+beta(ik)
342	enddo
343
344	ENDIF
345
346	! 2. Compute beta_PEM coefficients (preprocessing for next time step)
347	! Bottom layer, beta_PEM_{N-1}
348
349	beta(nsoil-1)=coefq_PEM(nsoil-1)*tsoil(nsoil) &
350	/(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1)) &
351	+ fluxgeo/(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1))
352	! Other layers
353	do ik=nsoil-2,1,-1
354	beta(ik)=(coefq_PEM(ik)*tsoil(ik+1)+ &
355	coefd_PEM(ik+1)*beta(ik+1))/ &
356	(coefq_PEM(ik)+coefd_PEM(ik+1)*(1.0-alph(ik+1)) &
357	+coefd_PEM(ik))
358	enddo
359	#endif
360	end
361
362	END MODULE ini_soil_mod

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

Download in other formats:

Original Format