Home

Context Navigation

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

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

Rev	Line
[1496]	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	!
[1514]	23	!
[1496]	24	! Former revisions:
	25	! -----------------
	26	! $Id: land_surface_model.f90 1514 2014-12-19 09:14:55Z hoffmann $
[1514]	27	! Bugfix: REAL constants provided with KIND-attribute in call of
	28	! intrinsic function MAX and MIN
[1496]	29	!
[1501]	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	!
[1497]	35	! 1496 2014-12-02 17:25:50Z maronga
	36	! Initial revision
	37	!
[1496]	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.
[1500]	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
[1496]	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, &
[1500]	62	ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
[1496]	63
	64	USE control_parameters, &
	65	ONLY: dt_3d, humidity, intermediate_timestep_count, &
[1500]	66	intermediate_timestep_count_max, precipitation, pt_surface, &
	67	rho_surface, surface_pressure, timestep_scheme, tsc
[1496]	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)
[1500]	139	soil_moisture = 0.0_wp !: soil moisture content (m3/m3)
[1496]	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
[1500]	552
	553	IF ( humidity ) THEN
	554	LE = rho_l * l_v * qsws
	555	ELSE
	556	LE = 0.0_wp
	557	ENDIF
	558
[1496]	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)
[1500]	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
[1496]	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
[1500]	747	k = nzb_s_inner(j,i)
[1496]	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
[1500]	756
[1496]	757	!
[1500]	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)
[1496]	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 )
[1513]	812	f2 = MAX(f2,1.0E-20_wp)
[1496]	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)
[1513]	819	c_liq(j,i) = MIN(1.0_wp, m_liq(j,i)/m_liq_max)
[1496]	820
	821	q_s = 0.622_wp * e_s / surface_pressure
[1500]	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
[1496]	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
[1500]	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
[1496]	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
[1500]	866	IF ( humidity ) THEN
	867
[1496]	868	!
[1500]	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)
[1496]	873
	874	!
[1500]	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
[1496]	878
[1500]	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
[1496]	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
[1500]	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) )
[1496]	932
[1500]	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
[1496]	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
[1500]	951	!
	952	!-- Calculate the true surface resistance
[1496]	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	!
[1500]	961	!-- Calculate fluxes in the atmosphere
	962	shf(j,i) = H(j,i) / rho_cp
	963
	964	!
[1496]	965	!-- Calculate change in liquid water reservoir due to dew fall or
[1500]	966	!-- evaporation of liquid water
	967	IF ( humidity ) THEN
[1496]	968	!
[1500]	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
[1496]	974	ENDIF
[1500]	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
[1496]	987
	988	!
[1500]	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
[1496]	997
[1500]	998	tend = - LE_liq(j,i) * drho_l_lv
[1496]	999
[1500]	1000	m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend &
[1496]	1001	+ tsc(3) * tm_liq_m(j,i) )
	1002
	1003	!
[1500]	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)
[1496]	1007
	1008	!
[1500]	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)
[1496]	1012
	1013
	1014	!
[1500]	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
[1496]	1024	ENDIF
	1025
	1026	!
[1500]	1027	!-- Calculate moisture flux in the atmosphere
	1028	qsws(j,i) = LE(j,i) / rho_lv
[1496]	1029
[1500]	1030	ENDIF
	1031
[1496]	1032	ENDDO
[1500]	1033	ENDDO
[1496]	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
[1513]	1072	Ke = 1.0 + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
[1496]	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 |