Home

Context Navigation

source: palm/trunk/SOURCE/land_surface_model.f90 @ 1520

Last change on this file since 1520 was 1514, checked in by heinze, 10 years ago
last commmit documented
Property svn:keywords set to `Id`
File size: 53.7 KB

Line
1	MODULE land_surface_model_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: land_surface_model.f90 1514 2014-12-19 09:14:55Z hoffmann $
27	! Bugfix: REAL constants provided with KIND-attribute in call of
28	! intrinsic function MAX and MIN
29	!
30	! 1500 2014-12-03 17:42:41Z maronga
31	! Corrected calculation of aerodynamic resistance (r_a).
32	! Precipitation is now added to liquid water reservoir using LE_liq.
33	! Added support for dry runs.
34	!
35	! 1496 2014-12-02 17:25:50Z maronga
36	! Initial revision
37	!
38	!
39	! Description:
40	! ------------
41	! Land surface model, consisting of a solver for the energy balance at the
42	! surface and a four layer soil scheme. The scheme is similar to the TESSEL
43	! scheme implemented in the ECMWF IFS model, with modifications according to
44	! H-TESSEL. The implementation is based on the formulation implemented in the
45	! DALES model.
46	!
47	! To do list:
48	! -----------
49	! - Add support for binary I/O support
50	! - Add support for lsm data output
51	! - Check for time step criterion
52	! - Check use with RK-2 and Euler time-stepping
53	! - Adaption for use with cloud physics (liquid water potential temperature)
54	! - Check reaction of plants at wilting point and at atmospheric saturation
55	! - Consider partial absorption of the net shortwave radiation by the skin layer
56	! - Allow for water surfaces, check performance for bare soils
57	!------------------------------------------------------------------------------!
58	USE arrays_3d, &
59	ONLY: pt, pt_p, q, q_p, qsws, rif, shf, ts, us, z0, z0h
60
61	USE cloud_parameters, &
62	ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
63
64	USE control_parameters, &
65	ONLY: dt_3d, humidity, intermediate_timestep_count, &
66	intermediate_timestep_count_max, precipitation, pt_surface, &
67	rho_surface, surface_pressure, timestep_scheme, tsc
68
69	USE indices, &
70	ONLY: nxlg, nxrg, nyng, nysg, nzb_s_inner
71
72	USE kinds
73
74	USE radiation_model_mod, &
75	ONLY: Rn, SW_in, sigma_SB
76
77
78	IMPLICIT NONE
79
80	!
81	!-- LSM model constants
82	INTEGER(iwp), PARAMETER :: soil_layers = 4 !: number of soil layers (fixed for now)
83
84	REAL(wp), PARAMETER :: &
85	b_CH = 6.04_wp, & ! Clapp & Hornberger exponent
86	lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil
87	lambda_h_sm = 3.44_wp, & ! heat conductivity of the soil matrix
88	lambda_h_water = 0.57_wp, & ! heat conductivity of water
89	psi_sat = -0.388_wp, & ! soil matrix potential at saturation
90	rhoC_soil = 2.19E6_wp, & ! volumetric heat capacity of soil
91	rhoC_water = 4.20E6_wp, & ! volumetric heat capacity of water
92	m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m)
93
94
95	!
96	!-- LSM variables
97	INTEGER(iwp) :: veg_type = 2, & !: vegetation type, 0: user-defined, 1-19: generic (see list)
98	soil_type = 3 !: soil type, 0: user-defined, 1-6: generic (see list)
99
100	LOGICAL :: conserve_water_content = .TRUE., & !: open or closed bottom surface for the soil model
101	land_surface = .FALSE. !: flag parameter indicating wheather the lsm is used
102
103	! value 9999999.9_wp -> generic available or user-defined value must be set
104	! otherwise -> no generic variable and user setting is optional
105	REAL(wp) :: alpha_VanGenuchten = 0.0_wp, & !: NAMELIST alpha_VG
106	canopy_resistance_coefficient = 0.0_wp, & !: NAMELIST gD
107	C_skin = 20000.0_wp, & !: Skin heat capacity
108	drho_l_lv, & !: (rho_l * l_v)**-1
109	exn, & !: value of the Exner function
110	e_s = 0.0_wp, & !: saturation water vapour pressure
111	field_capacity = 0.0_wp, & !: NAMELIST m_fc
112	f_shortwave_incoming = 9999999.9_wp, & !: NAMELIST f_SW_in
113	hydraulic_conductivity = 0.0_wp, & !: NAMELIST gamma_w_sat
114	Ke = 0.0_wp, & !: Kersten number
115	lambda_skin_stable = 9999999.9_wp, & !: NAMELIST lambda_skin_s
116	lambda_skin_unstable = 9999999.9_wp, & !: NAMELIST lambda_skin_u
117	leaf_area_index = 9999999.9_wp, & !: NAMELIST LAI
118	l_VanGenuchten = 0.0_WP, & !: NAMELIST l_VG
119	min_canopy_resistance = 110.0_wp, & !: NAMELIST r_s_min
120	m_total = 0.0_wp, & !: weighed total water content of the soil (m3/m3)
121	n_VanGenuchten = 0.0_WP, & !: NAMELIST n_VG
122	q_s = 0.0_wp, & !: saturation specific humidity
123	residual_moisture = 0.0_wp, & !: NAMELIST m_res
124	rho_cp, & !: rho_surface * cp
125	rho_lv, & !: rho * l_v
126	rd_d_rv, & !: r_d / r_v
127	saturation_moisture = 0.0_wp, & !: NAMELIST m_sat
128	vegetation_coverage = 9999999.9_wp, & !: NAMELIST c_veg
129	wilting_point = 0.0_wp !: NAMELIST m_wilt
130
131	REAL(wp), DIMENSION(0:soil_layers-1) :: &
132	ddz_soil, & !: 1/dz_soil
133	ddz_soil_stag, & !: 1/dz_soil_stag
134	dz_soil, & !: soil grid spacing (center-center)
135	dz_soil_stag, & !: soil grid spacing (edge-edge)
136	root_extr = 0.0_wp, & !: root extraction
137	root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers
138	soil_level = (/0.07_wp, 0.28_wp, 1.00_wp, 2.89_wp/), & !: soil layer depths (m)
139	soil_moisture = 0.0_wp !: soil moisture content (m3/m3)
140
141	REAL(wp), DIMENSION(0:soil_layers) :: &
142	soil_temperature = 9999999.9_wp !: soil temperature (K)
143
144	#if defined( __nopointer )
145	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0, & !: skin temperature (K)
146	T_0_p, & !: progn. skin temperature (K)
147	m_liq, & !: liquid water reservoir (m)
148	m_liq_p !: progn. liquid water reservoir (m)
149	#else
150	REAL(wp), DIMENSION(:,:), POINTER :: T_0, &
151	T_0_p, &
152	m_liq, &
153	m_liq_p
154
155	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0_1, T_0_2, &
156	m_liq_1, m_liq_2
157	#endif
158
159	!
160	!-- Temporal tendencies for time stepping
161	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tT_0_m, & !: skin temperature tendency (K)
162	tm_liq_m !: liquid water reservoir tendency (m)
163
164	!
165	!-- Energy balance variables
166	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
167	alpha_VG, & !: coef. of Van Genuchten
168	c_liq, & !: liquid water coverage (of vegetated area)
169	c_veg, & !: vegetation coverage
170	f_SW_in, & !: ?
171	G, & !: surface soil heat flux
172	H, & !: surface flux of sensible heat
173	gamma_w_sat, & !: hydraulic conductivity at saturation
174	gD, & !: coefficient for dependence of r_canopy on water vapour pressure deficit
175	LAI, & !: leaf area index
176	LE, & !: surface flux of latent heat (total)
177	LE_veg, & !: surface flux of latent heat (vegetation portion)
178	LE_soil, & !: surface flux of latent heat (soil portion)
179	LE_liq, & !: surface flux of latent heat (liquid water portion)
180	lambda_h_sat, & !: heat conductivity for dry soil
181	lambda_skin_s, & !: coupling between skin and soil (depends on vegetation type)
182	lambda_skin_u, & !: coupling between skin and soil (depends on vegetation type)
183	l_VG, & !: coef. of Van Genuchten
184	m_fc, & !: soil moisture at field capacity (m3/m3)
185	m_res, & !: residual soil moisture
186	m_sat, & !: saturation soil moisture (m3/m3)
187	m_wilt, & !: soil moisture at permanent wilting point (m3/m3)
188	n_VG, & !: coef. Van Genuchten
189	r_a, & !: aerodynamic resistance
190	r_canopy, & !: canopy resistance
191	r_soil, & !: soil resitance
192	r_soil_min, & !: minimum soil resistance
193	r_s, & !: total surface resistance (combination of r_soil and r_canopy)
194	r_s_min !: minimum canopy resistance
195
196	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
197	lambda_h, & !: heat conductivity of soil (?)
198	lambda_w, & !: hydraulic diffusivity of soil (?)
199	gamma_w, & !: hydraulic conductivity of soil (?)
200	rhoC_total !: volumetric heat capacity of the actual soil matrix (?)
201
202	#if defined( __nopointer )
203	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
204	T_soil, & !: Soil temperature (K)
205	T_soil_p, & !: Prog. soil temperature (K)
206	m_soil, & !: Soil moisture (m3/m3)
207	m_soil_p !: Prog. soil moisture (m3/m3)
208	#else
209	REAL(wp), DIMENSION(:,:,:), POINTER :: &
210	T_soil, T_soil_p, &
211	m_soil, m_soil_p
212
213	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
214	T_soil_1, T_soil_2, &
215	m_soil_1, m_soil_2
216
217
218	#endif
219
220
221	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
222	tT_soil_m, & !: T_soil storage array
223	tm_soil_m, & !: m_soil storage array
224	root_fr !: root fraction (sum=1)
225
226	!
227	!-- Land surface parameters according to the following classes (veg_type)
228	!-- (0 user defined)
229	!-- 1 crops, mixed farming
230	!-- 2 short grass
231	!-- 3 evergreen needleleaf trees
232	!-- 4 deciduous needleleaf trees
233	!-- 5 evergreen broadleaf trees
234	!-- 6 deciduous broadleaf trees
235	!-- 7 tall grass
236	!-- 8 desert
237	!-- 9 tundra
238	!-- 10 irrigated crops
239	!-- 11 semidesert
240	!-- 12 ice caps and glaciers
241	!-- 13 bogs and marshes
242	!-- 14 inland water
243	!-- 15 ocean
244	!-- 16 evergreen shrubs
245	!-- 17 deciduous shrubs
246	!-- 18 mixed forest/woodland
247	!-- 19 interrupted forest
248
249	!
250	!-- Land surface parameters I r_s_min, LAI, c_veg, gD
251	REAL(wp), DIMENSION(0:3,1:19) :: veg_pars = RESHAPE( (/ &
252	180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 1
253	110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & ! 2
254	500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 3
255	500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 4
256	175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 5
257	240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & ! 6
258	100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & ! 7
259	250.0_wp, 0.50_wp, 0.00_wp, 0.00_wp, & ! 8
260	80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & ! 9
261	180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10
262	150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, & ! 11
263	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
264	240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, & ! 13
265	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
266	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
267	225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, & ! 16
268	225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17
269	250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18
270	175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp & ! 19
271	/), (/ 4, 19 /) )
272
273	!
274	!-- Land surface parameters II z0, z0h
275	REAL(wp), DIMENSION(0:1,1:19) :: roughness_par = RESHAPE( (/ &
276	0.25_wp, 0.25E-2_wp, & ! 1
277	0.20_wp, 0.20E-2_wp, & ! 2
278	2.00_wp, 2.00_wp, & ! 3
279	2.00_wp, 2.00_wp, & ! 4
280	2.00_wp, 2.00_wp, & ! 5
281	2.00_wp, 2.00_wp, & ! 6
282	0.47_wp, 0.47E-2_wp, & ! 7
283	0.013_wp, 0.013E-2_wp, & ! 8
284	0.034_wp, 0.034E-2_wp, & ! 9
285	0.5_wp, 0.50E-2_wp, & ! 10
286	0.17_wp, 0.17E-2_wp, & ! 11
287	1.3E-3_wp, 1.3E-4_wp, & ! 12
288	0.83_wp, 0.83E-2_wp, & ! 13
289	0.00_wp, 0.00E-2_wp, & ! 14
290	0.00_wp, 0.00E-2_wp, & ! 15
291	0.10_wp, 0.10E-2_wp, & ! 16
292	0.25_wp, 0.25E-2_wp, & ! 17
293	2.00_wp, 2.00E-2_wp, & ! 18
294	1.10_wp, 1.10E-2_wp & ! 19
295	/), (/ 2, 19 /) )
296
297	!
298	!-- Land surface parameters III lambda_skin_s, lambda_skin_u, f_SW_in
299	REAL(wp), DIMENSION(0:2,1:19) :: skin_pars = RESHAPE( (/ &
300	10.0_wp, 10.0_wp, 0.05_wp, & ! 1
301	10.0_wp, 10.0_wp, 0.05_wp, & ! 2
302	20.0_wp, 15.0_wp, 0.03_wp, & ! 3
303	20.0_wp, 15.0_wp, 0.03_wp, & ! 4
304	20.0_wp, 15.0_wp, 0.03_wp, & ! 5
305	20.0_wp, 15.0_wp, 0.03_wp, & ! 6
306	10.0_wp, 10.0_wp, 0.05_wp, & ! 7
307	15.0_wp, 15.0_wp, 0.00_wp, & ! 8
308	10.0_wp, 10.0_wp, 0.05_wp, & ! 9
309	10.0_wp, 10.0_wp, 0.05_wp, & ! 10
310	10.0_wp, 10.0_wp, 0.05_wp, & ! 11
311	58.0_wp, 58.0_wp, 0.00_wp, & ! 12
312	10.0_wp, 10.0_wp, 0.05_wp, & ! 13
313	1.0E20_wp, 1.0E20_wp, 0.00_wp, & ! 14
314	1.0E20_wp, 1.0E20_wp, 0.00_wp, & ! 15
315	10.0_wp, 10.0_wp, 0.05_wp, & ! 16
316	10.0_wp, 10.0_wp, 0.05_wp, & ! 17
317	20.0_wp, 15.0_wp, 0.03_wp, & ! 18
318	20.0_wp, 15.0_wp, 0.03_wp & ! 19
319	/), (/ 3, 19 /) )
320
321	!
322	!-- Root distribution (sum = 1) level 1, level 2, level 3, level 4,
323	REAL(wp), DIMENSION(0:3,1:19) :: root_distribution = RESHAPE( (/ &
324	0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 1
325	0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 2
326	0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & ! 3
327	0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & ! 4
328	0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & ! 5
329	0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & ! 6
330	0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & ! 7
331	1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 8
332	0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 9
333	0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 10
334	0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 11
335	0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
336	0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 13
337	0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
338	0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
339	0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16
340	0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17
341	0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18
342	0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 19
343	/), (/ 4, 19 /) )
344
345	!
346	!-- Soil parameters according to the following porosity classes (soil_type)
347	!-- (0 user defined)
348	!-- 1 coarse
349	!-- 2 medium
350	!-- 3 medium-fine
351	!-- 4 fine
352	!-- 5 very fine
353	!-- 6 organic
354	!
355	!-- Soil parameters I alpha_VG, l_VG, n_VG, gamma_w_sat
356	REAL(wp), DIMENSION(0:3,1:6) :: soil_pars = RESHAPE( (/ &
357	3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, & ! 1
358	3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, & ! 2
359	0.83_wp, -0.588_wp, 1.25_wp, 0.26E-6_wp, & ! 3
360	3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, & ! 4
361	2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, & ! 5
362	1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp & ! 6
363	/), (/ 4, 6 /) )
364
365	!
366	!-- Soil parameters II m_sat, m_fc, m_wilt, m_res
367	REAL(wp), DIMENSION(0:3,1:6) :: m_soil_pars = RESHAPE( (/ &
368	0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
369	0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
370	0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
371	0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
372	0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
373	0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6
374	/), (/ 4, 6 /) )
375
376
377	SAVE
378
379
380	PRIVATE
381
382
383	PUBLIC alpha_VanGenuchten, C_skin, canopy_resistance_coefficient, &
384	conserve_water_content, field_capacity, f_shortwave_incoming, &
385	hydraulic_conductivity, init_lsm, lambda_skin_stable, &
386	lambda_skin_unstable, land_surface, leaf_area_index, &
387	lsm_energy_balance, lsm_soil_model, l_VanGenuchten, &
388	min_canopy_resistance, n_VanGenuchten, residual_moisture, &
389	root_fraction, saturation_moisture, soil_level, soil_moisture, &
390	soil_temperature, soil_type, vegetation_coverage, veg_type, &
391	wilting_point
392
393	#if defined( __nopointer )
394	PUBLIC m_liq, m_liq_p, m_soil, m_soil_p, T_0, T_0_p, T_soil, T_soil_p
395	#else
396	PUBLIC m_liq, m_liq_1, m_liq_2, m_liq_p, m_soil, m_soil_1, m_soil_2, &
397	m_soil_p, T_0, T_0_1, T_0_2, T_0_p, T_soil, T_soil_1, T_soil_2, &
398	T_soil_p
399	#endif
400
401
402	INTERFACE init_lsm
403	MODULE PROCEDURE init_lsm
404	END INTERFACE init_lsm
405
406	INTERFACE lsm_energy_balance
407	MODULE PROCEDURE lsm_energy_balance
408	END INTERFACE lsm_energy_balance
409
410	INTERFACE lsm_soil_model
411	MODULE PROCEDURE lsm_soil_model
412	END INTERFACE lsm_soil_model
413
414
415	CONTAINS
416
417
418	!------------------------------------------------------------------------------!
419	! Description:
420	! ------------
421	!-- Initialization of the land surface model
422	!------------------------------------------------------------------------------!
423	SUBROUTINE init_lsm
424
425
426	IMPLICIT NONE
427
428	INTEGER(iwp) :: i !: running index
429	INTEGER(iwp) :: j !: running index
430	INTEGER(iwp) :: k !: running index
431
432
433	!
434	!-- Calculate frequently used parameters
435	rho_cp = cp * rho_surface
436	rd_d_rv = r_d / r_v
437	rho_lv = rho_surface * l_v
438	drho_l_lv = 1.0 / (rho_l * l_v)
439
440	!
441	!-- Allocate skin and soil temperature / humidity
442	#if defined( __nopointer )
443	ALLOCATE ( T_0(nysg:nyng,nxlg:nxrg) )
444	ALLOCATE ( T_0_p(nysg:nyng,nxlg:nxrg) )
445	#else
446	ALLOCATE ( T_0_1(nysg:nyng,nxlg:nxrg) )
447	ALLOCATE ( T_0_2(nysg:nyng,nxlg:nxrg) )
448	#endif
449
450	ALLOCATE ( tT_0_m(nysg:nyng,nxlg:nxrg) )
451
452	#if defined( __nopointer )
453	ALLOCATE ( T_soil(0:soil_layers,nysg:nyng,nxlg:nxrg) )
454	ALLOCATE ( T_soil_p(0:soil_layers,nysg:nyng,nxlg:nxrg) )
455	#else
456	ALLOCATE ( T_soil_1(0:soil_layers,nysg:nyng,nxlg:nxrg) )
457	ALLOCATE ( T_soil_2(0:soil_layers,nysg:nyng,nxlg:nxrg) )
458	#endif
459
460	ALLOCATE ( tT_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
461
462	#if defined( __nopointer )
463	ALLOCATE ( m_liq(nysg:nyng,nxlg:nxrg) )
464	ALLOCATE ( m_liq_p(nysg:nyng,nxlg:nxrg) )
465	#else
466	ALLOCATE ( m_liq_1(nysg:nyng,nxlg:nxrg) )
467	ALLOCATE ( m_liq_2(nysg:nyng,nxlg:nxrg) )
468	#endif
469
470	ALLOCATE ( tm_liq_m(nysg:nyng,nxlg:nxrg) )
471
472	#if defined( __nopointer )
473	ALLOCATE ( m_soil(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
474	ALLOCATE ( m_soil_p(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
475	#else
476	ALLOCATE ( m_soil_1(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
477	ALLOCATE ( m_soil_2(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
478	#endif
479
480	ALLOCATE ( tm_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
481
482
483	#if ! defined( __nopointer )
484	!
485	!-- Initial assignment of the pointers
486	T_soil => T_soil_1; T_soil_p => T_soil_2
487	T_0 => T_0_1; T_0_p => T_0_2
488	m_soil => m_soil_1; m_soil_p => m_soil_2
489	m_liq => m_liq_1; m_liq_p => m_liq_2
490	#endif
491
492	T_0 = 0.0_wp
493	T_0_p = 0.0_wp
494	tT_0_m = 0.0_wp
495
496	T_soil = 0.0_wp
497	T_soil_p = 0.0_wp
498	tT_soil_m = 0.0_wp
499
500	m_liq = 0.0_wp
501	m_liq_p = 0.0_wp
502	tm_liq_m = 0.0_wp
503
504	m_soil = 0.0_wp
505	m_soil_p = 0.0_wp
506	tm_soil_m = 0.0_wp
507
508	!
509	!-- Allocate 2D vegetation model arrays
510	ALLOCATE ( alpha_VG(nysg:nyng,nxlg:nxrg) )
511	ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
512	ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
513	ALLOCATE ( f_SW_in(nysg:nyng,nxlg:nxrg) )
514	ALLOCATE ( G(nysg:nyng,nxlg:nxrg) )
515	ALLOCATE ( H(nysg:nyng,nxlg:nxrg) )
516	ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
517	ALLOCATE ( gD(nysg:nyng,nxlg:nxrg) )
518	ALLOCATE ( LAI(nysg:nyng,nxlg:nxrg) )
519	ALLOCATE ( LE(nysg:nyng,nxlg:nxrg) )
520	ALLOCATE ( LE_veg(nysg:nyng,nxlg:nxrg) )
521	ALLOCATE ( LE_soil(nysg:nyng,nxlg:nxrg) )
522	ALLOCATE ( LE_liq(nysg:nyng,nxlg:nxrg) )
523	ALLOCATE ( l_VG(nysg:nyng,nxlg:nxrg) )
524	ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
525	ALLOCATE ( lambda_skin_u(nysg:nyng,nxlg:nxrg) )
526	ALLOCATE ( lambda_skin_s(nysg:nyng,nxlg:nxrg) )
527	ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
528	ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
529	ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
530	ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
531	ALLOCATE ( n_VG(nysg:nyng,nxlg:nxrg) )
532	ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
533	ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
534	ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
535	ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
536	ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
537	ALLOCATE ( r_s_min(nysg:nyng,nxlg:nxrg) )
538
539	!
540	!-- Set initial and default values
541	c_liq = 0.0_wp
542	c_veg = 0.0_wp
543	f_SW_in = 0.05_wp
544	gD = 0.0_wp
545	LAI = 0.0_wp
546	lambda_skin_u = 10.0_wp
547	lambda_skin_s = 10.0_wp
548
549
550	G = 0.0_wp
551	H = rho_cp * shf
552
553	IF ( humidity ) THEN
554	LE = rho_l * l_v * qsws
555	ELSE
556	LE = 0.0_wp
557	ENDIF
558
559	LE_veg = 0.0_wp
560	LE_soil = LE
561	LE_liq = 0.0_wp
562
563	r_a = 50.0_wp
564	r_canopy = 0.0_wp
565	r_soil = 0.0_wp
566	r_soil_min = 50.0_wp
567	r_s = 110.0_wp
568	r_s_min = min_canopy_resistance
569
570	!
571	!-- Allocate 3D soil model arrays
572	ALLOCATE ( root_fr(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
573	ALLOCATE ( lambda_h(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
574	ALLOCATE ( rhoC_total(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
575
576	lambda_h = 0.0_wp
577	!
578	!-- If required, allocate humidity-related variables for the soil model
579	IF ( humidity ) THEN
580	ALLOCATE ( lambda_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
581	ALLOCATE ( gamma_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
582
583	lambda_w = 0.0_wp
584	ENDIF
585
586	!
587	!-- Calculate grid spacings. Temperature and moisture are defined at
588	!-- the center of the soil layers, whereas gradients/fluxes are defined
589	!-- at the edges (_stag)
590	dz_soil_stag(0) = soil_level(0)
591
592	DO k = 1, soil_layers-1
593	dz_soil_stag(k) = soil_level(k) - soil_level(k-1)
594	ENDDO
595
596	DO k = 0, soil_layers-2
597	dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k))
598	ENDDO
599	dz_soil(soil_layers-1) = dz_soil_stag(soil_layers-1)
600
601	ddz_soil = 1.0 / dz_soil
602	ddz_soil_stag = 1.0 / dz_soil_stag
603	!
604	!-- Initialize soil
605	IF ( soil_type .NE. 0 ) THEN
606	alpha_VG = soil_pars(0,soil_type)
607	l_VG = soil_pars(1,soil_type)
608	n_VG = soil_pars(2,soil_type)
609	gamma_w_sat = soil_pars(3,soil_type)
610	m_sat = m_soil_pars(0,soil_type)
611	m_fc = m_soil_pars(1,soil_type)
612	m_wilt = m_soil_pars(2,soil_type)
613	m_res = m_soil_pars(3,soil_type)
614	ELSE
615	alpha_VG = alpha_VanGenuchten
616	l_VG = l_VanGenuchten
617	n_VG = n_VanGenuchten
618	gamma_w_sat = hydraulic_conductivity
619	m_sat = saturation_moisture
620	m_fc = field_capacity
621	m_wilt = wilting_point
622	m_res = residual_moisture
623	ENDIF
624
625	!
626	!-- Map user settings of T and q for each soil layer
627	!-- (make sure that the soil moisture does not drop below the permanent
628	!-- wilting point) -> problems with devision by zero)
629	DO k = 0, soil_layers-1
630	T_soil(k,:,:) = soil_temperature(k)
631	m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:))
632	ENDDO
633	T_soil(soil_layers,:,:) = soil_temperature(soil_layers)
634
635
636	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
637	T_0 = pt_surface * exn
638
639	T_soil_p = T_soil
640	m_soil_p = m_soil
641
642	!
643	!-- Calculate saturation soil heat conductivity
644	lambda_h_sat(:,:) = lambda_h_sm ** (1.0_wp - m_sat(:,:)) * &
645	lambda_h_water ** m_sat(:,:)
646
647	!
648	!-- Initialize vegetation
649	IF ( veg_type .NE. 0 ) THEN
650
651	r_s_min = veg_pars(0,veg_type)
652	LAI = veg_pars(1,veg_type)
653	c_veg = veg_pars(2,veg_type)
654	gD = veg_pars(3,veg_type)
655	lambda_skin_s = skin_pars(0,veg_type)
656	lambda_skin_u = skin_pars(1,veg_type)
657	f_SW_in = skin_pars(2,veg_type)
658	z0 = roughness_par(0,veg_type)
659	z0h = roughness_par(1,veg_type)
660
661
662	DO k = 0, soil_layers-1
663	root_fr(k,:,:) = root_distribution(k,veg_type)
664	ENDDO
665
666	ELSE
667
668	DO k = 0, soil_layers-1
669	root_fr(k,:,:) = root_fraction(k)
670	ENDDO
671
672	ENDIF
673
674	!
675	!-- Possibly do user-defined actions (e.g. define heterogeneous land surface)
676	CALL user_init_land_surface
677
678	!
679	!-- Set artifical values for ts and us so that r_a has its initial value for
680	!-- the first time step
681	DO i = nxlg, nxrg
682	DO j = nysg, nyng
683	k = nzb_s_inner(j,i)
684	!
685	!-- Assure that r_a cannot be zero at model start
686	IF ( pt(k+1,j,i) == pt(k,j,i) ) pt(k+1,j,i) = pt(k+1,j,i) + 1.0E-10_wp
687
688	us(j,i) = 0.1_wp
689	ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
690	shf(j,i) = - us(j,i) * ts(j,i)
691	ENDDO
692	ENDDO
693
694	!
695	!-- Calculate humidity at the surface
696	IF ( humidity ) THEN
697	CALL calc_q0
698	ENDIF
699
700	RETURN
701
702	END SUBROUTINE init_lsm
703
704
705
706	!------------------------------------------------------------------------------!
707	! Description:
708	! ------------
709	!
710	!------------------------------------------------------------------------------!
711	SUBROUTINE lsm_energy_balance
712
713
714	IMPLICIT NONE
715
716	INTEGER(iwp) :: i !: running index
717	INTEGER(iwp) :: j !: running index
718	INTEGER(iwp) :: k, ks !: running index
719
720	REAL(wp) :: f1, & !: resistance correction term 1
721	f2, & !: resistance correction term 2
722	f3, & !: resistance correction term 3
723	m_min, & !: minimum soil moisture
724	T_1, & !: actual temperature at first grid point
725	e, & !: water vapour pressure
726	e_s, & !: water vapour saturation pressure
727	e_s_dT, & !: derivate of e_s with respect to T
728	tend, & !: tendency
729	dq_s_dT, & !: derivate of q_s with respect to T
730	coef_1, & !: coef. for prognostic equation
731	coef_2, & !: coef. for prognostic equation
732	f_LE, & !: factor for LE
733	f_LE_veg, & !: factor for LE_veg
734	f_LE_soil, & !: factor for LE_soil
735	f_LE_liq, & !: factor for LE_liq
736	f_H, & !: factor for H
737	lambda_skin, & !: Current value of lambda_skin
738	m_liq_max !: maxmimum value of the liquid water reservoir
739
740	!
741	!-- Calculate the exner function for the current time step
742	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
743
744
745	DO i = nxlg, nxrg
746	DO j = nysg, nyng
747	k = nzb_s_inner(j,i)
748
749	!
750	!-- Set lambda_skin according to stratification
751	IF ( rif(j,i) >= 0.0_wp ) THEN
752	lambda_skin = lambda_skin_s(j,i)
753	ELSE
754	lambda_skin = lambda_skin_u(j,i)
755	ENDIF
756
757	!
758	!-- First step: calculate aerodyamic resistance. As pt, us, ts
759	!-- are not available for the prognostic time step, data from the last
760	!-- time step is used here. Note that this formulation is the
761	!-- equivalent to the ECMWF formulation using drag coefficients
762	r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20)
763
764	!
765	!-- Second step: calculate canopy resistance r_canopy
766	!-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
767
768	!-- f1: correction for incoming shortwave radiation
769	f1 = MIN(1.0_wp, ( 0.004_wp * SW_in(j,i) + 0.05_wp ) / &
770	(0.81_wp * (0.004_wp * SW_in(j,i) + 1.0_wp) ) )
771
772	!
773	!-- f2: correction for soil moisture f2=0 for very dry soil
774	m_total = 0.0_wp
775	DO ks = 0, soil_layers-1
776	m_total = m_total + root_fr(ks,j,i) * m_soil(ks,j,i)
777	ENDDO
778
779	IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LE. m_fc(j,i) ) THEN
780	f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
781	ELSE
782	f2 = 1.0E-20_wp
783	ENDIF
784
785	!
786	!-- Calculate water vapour pressure at saturation
787	!-- (T_0 should be replaced by liquid water temp?!)
788	e_s = 0.01 * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) - 273.16_wp )&
789	/ ( T_0(j,i) - 35.86_wp ) )
790
791	!
792	!-- f3: correction for vapour pressure deficit
793	IF ( gD(j,i) .NE. 0.0_wp ) THEN
794	!
795	!-- Calculate vapour pressure
796	e = q_p(k+1,j,i) * surface_pressure / 0.622
797	f3 = EXP ( -gD(j,i) * (e_s - e) )
798	ELSE
799	f3 = 1.0_wp
800	ENDIF
801
802	!
803	!-- To do: check for very dry soil -> r_canopy goes to infinity
804	r_canopy(j,i) = r_s_min(j,i) / (LAI(j,i) * f1 * f2 * f3 + 1.0E-20)
805
806	!
807	!-- Third step: calculate bare soil resistance r_soil
808	m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * &
809	m_res(j,i)
810
811	f2 = ( m_soil(0,j,i) - m_min ) / ( m_fc(j,i) - m_min )
812	f2 = MAX(f2,1.0E-20_wp)
813
814	r_soil(j,i) = r_soil_min(j,i) / f2
815
816	!
817	!-- Calculate fraction of liquid water reservoir
818	m_liq_max = m_max_depth * LAI(j,i)
819	c_liq(j,i) = MIN(1.0_wp, m_liq(j,i)/m_liq_max)
820
821	q_s = 0.622_wp * e_s / surface_pressure
822
823	!
824	!-- In case of dew fall, set resistances to zero.
825	!-- To do: what does that physically reasoning behind this?
826	IF ( humidity ) THEN
827	IF ( q_s .LE. q_p(k+1,j,i) ) THEN
828	r_canopy(j,i) = 0.0_wp
829	r_soil(j,i) = 0.0_wp
830	ENDIF
831	ENDIF
832
833
834	!
835	!-- Calculate coefficients for the total evapotranspiration
836	f_LE_veg = rho_lv * c_veg(j,i) * (1.0 - c_liq(j,i)) / (r_a(j,i) &
837	+ r_canopy(j,i))
838	f_LE_soil = rho_lv * (1.0 - c_veg(j,i)) / (r_a(j,i) + r_soil(j,i))
839	f_LE_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
840
841
842	!
843	!-- If soil moisture is below wilting point, plants do no longer
844	!-- transpirate.
845	IF ( m_soil(k,j,i) .LT. m_wilt(j,i) ) THEN
846	f_LE_veg = 0.0
847	ENDIF
848
849	f_H = rho_cp / r_a(j,i)
850	f_LE = f_LE_veg + f_LE_soil + f_LE_liq
851
852	!
853	!-- Calculate derivative of q_s for Taylor series expansion
854	e_s_dT = e_s * ( 17.269_wp / (T_0(j,i) - 35.86_wp) - &
855	17.269_wp*(T_0(j,i) - 273.16_wp) / (T_0(j,i) &
856	- 35.86_wp)**2 )
857
858	dq_s_dT = 0.622_wp * e_s_dT / surface_pressure
859
860	T_1 = pt_p(k+1,j,i) * exn
861
862	!
863	!-- Add LW up so that it can be removed in prognostic equation
864	Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4
865
866	IF ( humidity ) THEN
867
868	!
869	!-- Numerator of the prognostic equation
870	coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H &
871	/ exn * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT &
872	* T_0(j,i) ) + lambda_skin * T_soil(0,j,i)
873
874	!
875	!-- Denominator of the prognostic equation
876	coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT &
877	+ lambda_skin + f_H / exn
878
879	ELSE
880
881	!
882	!-- Numerator of the prognostic equation
883	coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H / &
884	exn * T_1 + lambda_skin * T_soil(0,j,i)
885
886	!
887	!-- Denominator of the prognostic equation
888	coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 &
889	+ lambda_skin + f_H / exn
890
891	ENDIF
892
893	tend = 0.0_wp
894
895	!
896	!-- Implicit solution when the skin layer has no heat capacity,
897	!-- otherwise use RK3 scheme.
898	T_0_p(j,i) = ( coef_1 * dt_3d * tsc(2) + C_skin * T_0(j,i) ) / &
899	( C_skin + coef_2 * dt_3d * tsc(2) )
900
901	!
902	!-- Add RK3 term
903	T_0_p(j,i) = T_0_p(j,i) + dt_3d * tsc(3) * tT_soil_m(0,j,i)
904
905	!
906	!-- Calculate true tendency
907	tend = (T_0_p(j,i) - T_0(j,i) - tsc(3) * tT_0_m(j,i)) / (dt_3d &
908	* tsc(2))
909
910	!
911	!-- Calculate T_0 tendencies for the next Runge-Kutta step
912	IF ( timestep_scheme(1:5) == 'runge' ) THEN
913	IF ( intermediate_timestep_count == 1 ) THEN
914	tT_0_m(j,i) = tend
915	ELSEIF ( intermediate_timestep_count < &
916	intermediate_timestep_count_max ) THEN
917	tT_0_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tT_0_m(j,i)
918	ENDIF
919	ENDIF
920
921	pt_p(k,j,i) = T_0_p(j,i) / exn
922	!
923	!-- Calculate fluxes
924	Rn(j,i) = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i)**4 &
925	- 4.0_wp * sigma_SB * T_0(j,i)*3 T_0_p(j,i)
926	G(j,i) = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i))
927	H(j,i) = - f_H * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
928
929	IF ( humidity ) THEN
930	LE(j,i) = - f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT &
931	* T_0(j,i) - dq_s_dT * T_0_p(j,i) )
932
933	LE_veg(j,i) = - f_LE_veg * ( q_p(k+1,j,i) - q_s + dq_s_dT &
934	* T_0(j,i) - dq_s_dT * T_0_p(j,i) )
935	LE_soil(j,i) = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT &
936	* T_0(j,i) - dq_s_dT * T_0_p(j,i) )
937	LE_liq(j,i) = - f_LE_liq * ( q_p(k+1,j,i) - q_s + dq_s_dT &
938	* T_0(j,i) - dq_s_dT * T_0_p(j,i) )
939	ENDIF
940
941	! IF ( i == 1 .AND. j == 1 ) THEN
942	! PRINT*, "Rn", Rn(j,i)
943	! PRINT*, "H", H(j,i)
944	! PRINT*, "LE", LE(j,i)
945	! PRINT*, "LE_liq", LE_liq(j,i)
946	! PRINT*, "LE_veg", LE_veg(j,i)
947	! PRINT*, "LE_soil", LE_soil(j,i)
948	! PRINT*, "G", G(j,i)
949	! ENDIF
950
951	!
952	!-- Calculate the true surface resistance
953	IF ( LE(j,i) .EQ. 0.0 ) THEN
954	r_s(j,i) = 1.0E10
955	ELSE
956	r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dT * T_0(j,i)&
957	- dq_s_dT * T_0_p(j,i) ) / LE(j,i) - r_a(j,i)
958	ENDIF
959
960	!
961	!-- Calculate fluxes in the atmosphere
962	shf(j,i) = H(j,i) / rho_cp
963
964	!
965	!-- Calculate change in liquid water reservoir due to dew fall or
966	!-- evaporation of liquid water
967	IF ( humidity ) THEN
968	!
969	!-- If precipitation is activated, add rain water to LE_liq.
970	!-- precipitation_rate is given in mm.
971	IF ( precipitation ) THEN
972	LE_liq(j,i) = LE_liq(j,i) + precipitation_rate(j,i) &
973	* 0.001_wp * rho_l * l_v
974	ENDIF
975	!
976	!-- If the air is saturated, check the reservoir water level
977	IF ( q_s .LE. q_p(k+1,j,i)) THEN
978	!
979	!-- Check if reservoir is full (avoid values > m_liq_max)
980	!-- In that case, LE_liq goes to LE_soil. In this case
981	!-- LE_veg is zero anyway (because c_liq = 1), so that tend is
982	!-- zero and no further check is needed
983	IF ( m_liq(j,i) .EQ. m_liq_max ) THEN
984	LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i)
985	LE_liq(j,i) = 0.0_wp
986	ENDIF
987
988	!
989	!-- In case LE_veg becomes negative (unphysical behavior), let
990	!-- the water enter the liquid water reservoir as dew on the
991	!-- plant
992	IF ( LE_veg(j,i) .LT. 0.0_wp ) THEN
993	LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)
994	LE_veg(j,i) = 0.0_wp
995	ENDIF
996	ENDIF
997
998	tend = - LE_liq(j,i) * drho_l_lv
999
1000	m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend &
1001	+ tsc(3) * tm_liq_m(j,i) )
1002
1003	!
1004	!-- Check if reservoir is overfull -> reduce to maximum
1005	!-- (conservation of water is violated here)
1006	m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max)
1007
1008	!
1009	!-- Check if reservoir is empty (avoid values < 0.0)
1010	!-- (conservation of water is violated here)
1011	m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp)
1012
1013
1014	!
1015	!-- Calculate m_liq tendencies for the next Runge-Kutta step
1016	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1017	IF ( intermediate_timestep_count == 1 ) THEN
1018	tm_liq_m(j,i) = tend
1019	ELSEIF ( intermediate_timestep_count < &
1020	intermediate_timestep_count_max ) THEN
1021	tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp &
1022	* tm_liq_m(j,i)
1023	ENDIF
1024	ENDIF
1025
1026	!
1027	!-- Calculate moisture flux in the atmosphere
1028	qsws(j,i) = LE(j,i) / rho_lv
1029
1030	ENDIF
1031
1032	ENDDO
1033	ENDDO
1034
1035
1036
1037	END SUBROUTINE lsm_energy_balance
1038
1039
1040	!------------------------------------------------------------------------------!
1041	! Description:
1042	! ------------
1043	!
1044	!------------------------------------------------------------------------------!
1045	SUBROUTINE lsm_soil_model
1046
1047
1048	IMPLICIT NONE
1049
1050	INTEGER(iwp) :: i !: running index
1051	INTEGER(iwp) :: j !: running index
1052	INTEGER(iwp) :: k !: running index
1053
1054	REAL(wp) :: h_VG !: Van Genuchten coef. h
1055
1056	REAL(wp), DIMENSION(0:soil_layers-1) :: gamma_temp, & !: temp. gamma
1057	lambda_temp, & !: temp. lambda
1058	tend !: tendency
1059
1060	DO i = nxlg, nxrg
1061	DO j = nysg, nyng
1062	DO k = 0, soil_layers-1
1063	!
1064	!-- Calculate volumetric heat capacity of the soil, taking into
1065	!-- account water content
1066	rhoC_total(k,j,i) = (rhoC_soil * (1.0 - m_sat(j,i)) &
1067	+ rhoC_water * m_soil(k,j,i))
1068
1069	!
1070	!-- Calculate soil heat conductivity at the center of the soil
1071	!-- layers
1072	Ke = 1.0 + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
1073	lambda_temp(k) = Ke * (lambda_h_sat(j,i) + lambda_h_dry) + &
1074	lambda_h_dry
1075
1076	ENDDO
1077
1078	!
1079	!-- Calculate soil heat conductivity (lambda_h) at the _stag level
1080	!-- using linear interpolation
1081	DO k = 0, soil_layers-2
1082
1083	lambda_h(k,j,i) = lambda_temp(k) + &
1084	( lambda_temp(k+1) - lambda_temp(k) ) &
1085	* 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
1086
1087	ENDDO
1088	lambda_h(soil_layers-1,j,i) = lambda_temp(soil_layers-1)
1089
1090	!
1091	!-- Prognostic equation for soil temperature T_soil
1092	tend(:) = 0.0_wp
1093	tend(0) = (1.0/rhoC_total(0,j,i)) * &
1094	( lambda_h(0,j,i) * ( T_soil(1,j,i) - T_soil(0,j,i) ) &
1095	* ddz_soil(0) + G(j,i) ) * ddz_soil_stag(0)
1096
1097	DO k = 1, soil_layers-1
1098	tend(k) = (1.0/rhoC_total(k,j,i)) &
1099	* ( lambda_h(k,j,i) &
1100	* ( T_soil(k+1,j,i) - T_soil(k,j,i) ) &
1101	* ddz_soil(k) &
1102	- lambda_h(k-1,j,i) &
1103	* ( T_soil(k,j,i) - T_soil(k-1,j,i) ) &
1104	* ddz_soil(k-1) &
1105	) * ddz_soil_stag(k)
1106	ENDDO
1107
1108	T_soil_p(0:soil_layers-1,j,i) = T_soil(0:soil_layers-1,j,i) &
1109	+ dt_3d * ( tsc(2) &
1110	* tend(:) + tsc(3) &
1111	* tT_soil_m(:,j,i) )
1112
1113	!
1114	!-- Calculate T_soil tendencies for the next Runge-Kutta step
1115	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1116	IF ( intermediate_timestep_count == 1 ) THEN
1117	DO k = 0, soil_layers-1
1118	tT_soil_m(k,j,i) = tend(k)
1119	ENDDO
1120	ELSEIF ( intermediate_timestep_count < &
1121	intermediate_timestep_count_max ) THEN
1122	DO k = 0, soil_layers-1
1123	tT_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp &
1124	* tT_soil_m(k,j,i)
1125	ENDDO
1126	ENDIF
1127	ENDIF
1128
1129
1130	DO k = 0, soil_layers-1
1131	!
1132	!-- Calculate soil diffusivity at the center of the soil layers
1133	lambda_temp(k) = (- b_CH * gamma_w_sat(j,i) * psi_sat &
1134	/ m_sat(j,i) ) * ( MAX(m_soil(k,j,i), &
1135	m_wilt(j,i)) / m_sat(j,i) )**(b_CH + 2.0_wp)
1136
1137	!
1138	!-- Calculate the hydraulic conductivity after Van Genuchten (1980)
1139	h_VG = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - &
1140	MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_VG(j,i) &
1141	/ (n_VG(j,i)-1.0_wp)) - 1.0_wp &
1142	)**(1.0_wp/n_VG(j,i)) / alpha_VG(j,i)
1143
1144	gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + &
1145	(alpha_VG(j,i)h_VG)n_VG(j,i))*(1.0_wp &
1146	-1.0_wp/n_VG(j,i)) - (alpha_VG(j,i)*h_VG &
1147	)(n_VG(j,i)-1.0_wp))2 ) &
1148	/ ( (1.0_wp + (alpha_VG(j,i)h_VG)*n_VG(j,i) &
1149	)*((1.0_wp - 1.0_wp/n_VG(j,i))(l_VG(j,i) &
1150	+ 2.0)) )
1151
1152	ENDDO
1153
1154
1155	IF ( humidity ) THEN
1156	!
1157	!-- Calculate soil diffusivity (lambda_w) at the _stag level
1158	!-- using linear interpolation
1159	DO k = 0, soil_layers-2
1160
1161	lambda_w(k,j,i) = lambda_temp(k) + &
1162	( lambda_temp(k+1) - lambda_temp(k) ) &
1163	* 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
1164	gamma_w(k,j,i) = gamma_temp(k) + &
1165	( gamma_temp(k+1) - gamma_temp(k) ) &
1166	* 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
1167
1168	ENDDO
1169
1170	!
1171	!
1172	!-- In case of a closed bottom (= water content is conserved), set
1173	!-- hydraulic conductivity to zero to that no water will be lost
1174	!-- in the bottom layer.
1175	IF ( conserve_water_content ) THEN
1176	gamma_w(soil_layers-1,j,i) = 0.0_wp
1177	ELSE
1178	gamma_w(soil_layers-1,j,i) = lambda_temp(soil_layers-1)
1179	ENDIF
1180
1181	!-- The root extraction (= root_extr * LE_veg / (rho_l * l_v))
1182	!-- ensures the mass conservation for water. The transpiration of
1183	!-- plants equals the cumulative withdrawals by the roots in the
1184	!-- soil. The scheme takes into account the availability of water
1185	!-- in the soil layers as well as the root fraction in the
1186	!-- respective layer
1187
1188	!
1189	!-- Calculate the root extraction (ECMWF 7.69, with some
1190	!-- modifications)
1191	m_total = 0.0_wp
1192	DO k = 0, soil_layers-1
1193	m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) * &
1194	dz_soil_stag(k)
1195
1196	ENDDO
1197
1198	!
1199	!-- For conservation of mass, the sum of root_extr must be 1
1200	DO k = 0, soil_layers-1
1201	root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) &
1202	* dz_soil_stag(k) / m_total
1203	ENDDO
1204
1205
1206	!
1207	!-- Prognostic equation for soil water content m_soil
1208	tend(:) = 0.0_wp
1209	tend(0) = ( lambda_w(0,j,i) * ( m_soil(1,j,i) - m_soil(0,j,i) )&
1210	* ddz_soil(0) - gamma_w(0,j,i) - ( root_extr(0) &
1211	* LE_veg(j,i) + LE_soil(j,i) ) * drho_l_lv &
1212	) * ddz_soil_stag(0)
1213
1214	DO k = 1, soil_layers-2
1215	tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i) &
1216	- m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&
1217	- lambda_w(k-1,j,i) * (m_soil(k,j,i) - &
1218	m_soil(k-1,j,i)) * ddz_soil(k-1) &
1219	+ gamma_w(k-1,j,i) - (root_extr(k) * LE_veg(j,i)&
1220	* drho_l_lv) &
1221	) * ddz_soil_stag(k)
1222
1223	ENDDO
1224	tend(soil_layers-1) = ( - gamma_w(soil_layers-1,j,i) &
1225	- lambda_w(soil_layers-2,j,i) &
1226	* (m_soil(soil_layers-1,j,i) &
1227	- m_soil(soil_layers-2,j,i)) &
1228	* ddz_soil(soil_layers-2) &
1229	+ gamma_w(soil_layers-2,j,i) - ( &
1230	root_extr(soil_layers-1) &
1231	* LE_veg(j,i) * drho_l_lv ) &
1232	) * ddz_soil_stag(soil_layers-1)
1233
1234	m_soil_p(0:soil_layers-1,j,i) = m_soil(0:soil_layers-1,j,i) &
1235	+ dt_3d * ( tsc(2) * tend(:) &
1236	+ tsc(3) * tm_soil_m(:,j,i) )
1237
1238	!
1239	!-- Account for dry soils (find a better solution here!)
1240	m_soil_p(:,j,i) = MAX(m_soil_p(:,j,i),0.0_wp)
1241
1242	!
1243	!-- Calculate m_soil tendencies for the next Runge-Kutta step
1244	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1245	IF ( intermediate_timestep_count == 1 ) THEN
1246	DO k = 0, soil_layers-1
1247	tm_soil_m(k,j,i) = tend(k)
1248	ENDDO
1249	ELSEIF ( intermediate_timestep_count < &
1250	intermediate_timestep_count_max ) THEN
1251	DO k = 0, soil_layers-1
1252	tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp &
1253	* tm_soil_m(k,j,i)
1254	ENDDO
1255	ENDIF
1256	ENDIF
1257
1258	ENDIF
1259
1260	ENDDO
1261	ENDDO
1262
1263	!
1264	!-- Calculate surface specific humidity
1265	IF ( humidity ) THEN
1266	CALL calc_q0
1267	ENDIF
1268
1269
1270	END SUBROUTINE lsm_soil_model
1271
1272
1273	!------------------------------------------------------------------------------!
1274	! Description:
1275	! ------------
1276	!
1277	!------------------------------------------------------------------------------!
1278	SUBROUTINE calc_q0
1279
1280	IMPLICIT NONE
1281
1282	INTEGER :: i !: running index
1283	INTEGER :: j !: running index
1284	INTEGER :: k !: running index
1285	REAL(wp) :: resistance !: aerodynamic and soil resistance term
1286
1287	DO i = nxlg, nxrg
1288	DO j = nysg, nyng
1289	k = nzb_s_inner(j,i)
1290	!
1291	!-- Temporary solution as long as T_0 is prescribed
1292
1293	pt_p(k,j,i) = T_0(j,i) / exn
1294	!
1295	!-- Calculate water vapour pressure at saturation
1296	e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) - &
1297	273.16_wp ) / ( T_0(j,i) - &
1298	35.86_wp ) )
1299
1300	!
1301	!-- Calculate specific humidity at saturation
1302	q_s = 0.622_wp * e_s / surface_pressure
1303
1304
1305	resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
1306
1307	!
1308	!-- Calculate specific humidity at surface
1309	q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance) &
1310	* q_p(k+1,j,i)
1311
1312	ENDDO
1313	ENDDO
1314
1315	END SUBROUTINE calc_q0
1316
1317
1318	END MODULE land_surface_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |