Context Navigation

ini_soil_mod.F90 @ 2794

Last change on this file since 2794 was 2794, checked in by llange, 2 years ago

MARS PEM:

Add a PEMETAT0 that read "startfi_pem.nc"
Add the soil in the model: soil temperature, thermal properties, ice table
Add a routine that compute CO2 + H2O adsorption
Minor corrections in PEM.F90

File size: 11.7 KB

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

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

Download in other formats:

Original Format