Home

Context Navigation

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

Last change on this file since 1590 was 1590, checked in by maronga, 9 years ago
bugfix in definition of character arrays within the LSM and radiation code
Property svn:keywords set to `Id`
File size: 71.4 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	! -----------------
[1590]	22	! Bugfix: definition of character strings requires same length for all elements
[1552]	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: land_surface_model.f90 1590 2015-05-08 13:56:27Z maronga $
	27	! Bugfix: REAL constants provided with KIND-attribute in call of
[1586]	28	! 1585 2015-04-30 07:05:52Z maronga
	29	! Modifications for RRTMG. Changed tables to PARAMETER type.
	30	!
[1572]	31	! 1571 2015-03-12 16:12:49Z maronga
	32	! Removed upper-case variable names. Corrected distribution of precipitation to
	33	! the liquid water reservoir and the bare soil fractions.
	34	!
[1556]	35	! 1555 2015-03-04 17:44:27Z maronga
	36	! Added output of r_a and r_s
	37	!
[1554]	38	! 1553 2015-03-03 17:33:54Z maronga
	39	! Improved better treatment of roughness lengths. Added default soil temperature
	40	! profile
	41	!
[1552]	42	! 1551 2015-03-03 14:18:16Z maronga
[1551]	43	! Flux calculation is now done in prandtl_fluxes. Added support for data output.
	44	! Vertical indices have been replaced. Restart runs are now possible. Some
	45	! variables have beem renamed. Bugfix in the prognostic equation for the surface
	46	! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
	47	! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
	48	! the hydraulic conductivity. Bugfix for root fraction and extraction
	49	! calculation
[1496]	50	!
[1514]	51	! intrinsic function MAX and MIN
[1496]	52	!
[1501]	53	! 1500 2014-12-03 17:42:41Z maronga
	54	! Corrected calculation of aerodynamic resistance (r_a).
	55	! Precipitation is now added to liquid water reservoir using LE_liq.
	56	! Added support for dry runs.
	57	!
[1497]	58	! 1496 2014-12-02 17:25:50Z maronga
	59	! Initial revision
	60	!
[1496]	61	!
	62	! Description:
	63	! ------------
	64	! Land surface model, consisting of a solver for the energy balance at the
	65	! surface and a four layer soil scheme. The scheme is similar to the TESSEL
	66	! scheme implemented in the ECMWF IFS model, with modifications according to
	67	! H-TESSEL. The implementation is based on the formulation implemented in the
[1551]	68	! DALES and UCLA-LES models.
[1500]	69	!
	70	! To do list:
	71	! -----------
[1551]	72	! - Check dewfall parametrization for fog simulations
	73	! - Consider partial absorption of the net shortwave radiation by the surface layer
[1500]	74	! - Allow for water surfaces, check performance for bare soils
[1551]	75	! - Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3))
	76	! - Implement surface runoff model (required when performing long-term LES with
	77	! considerable precipitation
[1571]	78	!
[1551]	79	! Notes:
	80	! ------
	81	! - No time step criterion is required as long as the soil layers do not become
	82	! too thin
[1496]	83	!------------------------------------------------------------------------------!
	84	USE arrays_3d, &
[1585]	85	ONLY: pt, pt_p, q_p, qsws, rif, shf, ts, us, z0, z0h
[1496]	86
	87	USE cloud_parameters, &
[1551]	88	ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
[1496]	89
	90	USE control_parameters, &
[1551]	91	ONLY: cloud_physics, dt_3d, humidity, intermediate_timestep_count, &
	92	initializing_actions, intermediate_timestep_count_max, &
	93	max_masks, precipitation, pt_surface, rho_surface, &
	94	roughness_length, surface_pressure, timestep_scheme, tsc, &
[1585]	95	z0h_factor, time_since_reference_point
[1496]	96
	97	USE indices, &
[1585]	98	ONLY: nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner
[1496]	99
	100	USE kinds
	101
[1551]	102	USE netcdf_control, &
	103	ONLY: dots_label, dots_num, dots_unit
	104
[1496]	105	USE radiation_model_mod, &
[1585]	106	ONLY: radiation_scheme, rad_net, rad_sw_in, sigma_sb
[1496]	107
[1551]	108	USE statistics, &
	109	ONLY: hom, statistic_regions
[1496]	110
	111	IMPLICIT NONE
	112
	113	!
	114	!-- LSM model constants
[1551]	115	INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !: bottom of the soil model (to be switched)
	116	nzt_soil = 3, & !: top of the soil model (to be switched)
	117	nzs = 4 !: number of soil layers (fixed for now)
[1496]	118
[1551]	119	INTEGER(iwp) :: dots_soil = 0 !: starting index for timeseries output
	120
	121	INTEGER(iwp), DIMENSION(0:1) :: id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz, &
	122	id_dim_zs_3d, id_var_zs_xy, &
	123	id_var_zs_xz, id_var_zs_yz, id_var_zs_3d
	124
	125	INTEGER(iwp), DIMENSION(1:max_masks,0:1) :: id_dim_zs_mask, id_var_zs_mask
	126
[1496]	127	REAL(wp), PARAMETER :: &
[1551]	128	b_ch = 6.04_wp, & ! Clapp & Hornberger exponent
	129	lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil
[1496]	130	lambda_h_sm = 3.44_wp, & ! heat conductivity of the soil matrix
	131	lambda_h_water = 0.57_wp, & ! heat conductivity of water
	132	psi_sat = -0.388_wp, & ! soil matrix potential at saturation
[1551]	133	rho_c_soil = 2.19E6_wp, & ! volumetric heat capacity of soil
	134	rho_c_water = 4.20E6_wp, & ! volumetric heat capacity of water
[1496]	135	m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m)
	136
	137
	138	!
	139	!-- LSM variables
	140	INTEGER(iwp) :: veg_type = 2, & !: vegetation type, 0: user-defined, 1-19: generic (see list)
	141	soil_type = 3 !: soil type, 0: user-defined, 1-6: generic (see list)
	142
	143	LOGICAL :: conserve_water_content = .TRUE., & !: open or closed bottom surface for the soil model
[1551]	144	dewfall = .TRUE., & !: allow/inhibit dewfall
[1496]	145	land_surface = .FALSE. !: flag parameter indicating wheather the lsm is used
	146
	147	! value 9999999.9_wp -> generic available or user-defined value must be set
	148	! otherwise -> no generic variable and user setting is optional
[1551]	149	REAL(wp) :: alpha_vangenuchten = 9999999.9_wp, & !: NAMELIST alpha_vg
	150	canopy_resistance_coefficient = 9999999.9_wp, & !: NAMELIST g_d
	151	c_surface = 20000.0_wp, & !: Surface (skin) heat capacity
[1496]	152	drho_l_lv, & !: (rho_l * l_v)**-1
	153	exn, & !: value of the Exner function
	154	e_s = 0.0_wp, & !: saturation water vapour pressure
[1551]	155	field_capacity = 9999999.9_wp, & !: NAMELIST m_fc
	156	f_shortwave_incoming = 9999999.9_wp, & !: NAMELIST f_sw_in
	157	hydraulic_conductivity = 9999999.9_wp, & !: NAMELIST gamma_w_sat
	158	ke = 0.0_wp, & !: Kersten number
	159	lambda_surface_stable = 9999999.9_wp, & !: NAMELIST lambda_surface_s
	160	lambda_surface_unstable = 9999999.9_wp, & !: NAMELIST lambda_surface_u
	161	leaf_area_index = 9999999.9_wp, & !: NAMELIST lai
	162	l_vangenuchten = 9999999.9_wp, & !: NAMELIST l_vg
	163	min_canopy_resistance = 9999999.9_wp, & !: NAMELIST r_canopy_min
	164	min_soil_resistance = 50.0_wp, & !: NAMELIST r_soil_min
	165	m_total = 0.0_wp, & !: weighted total water content of the soil (m3/m3)
	166	n_vangenuchten = 9999999.9_wp, & !: NAMELIST n_vg
[1496]	167	q_s = 0.0_wp, & !: saturation specific humidity
[1551]	168	residual_moisture = 9999999.9_wp, & !: NAMELIST m_res
[1496]	169	rho_cp, & !: rho_surface * cp
	170	rho_lv, & !: rho * l_v
	171	rd_d_rv, & !: r_d / r_v
[1551]	172	saturation_moisture = 9999999.9_wp, & !: NAMELIST m_sat
[1496]	173	vegetation_coverage = 9999999.9_wp, & !: NAMELIST c_veg
[1551]	174	wilting_point = 9999999.9_wp, & !: NAMELIST m_wilt
	175	z0_eb = 9999999.9_wp, & !: NAMELIST z0 (lsm_par)
	176	z0h_eb = 9999999.9_wp !: NAMELIST z0h (lsm_par)
[1496]	177
[1551]	178	REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
[1496]	179	ddz_soil, & !: 1/dz_soil
	180	ddz_soil_stag, & !: 1/dz_soil_stag
	181	dz_soil, & !: soil grid spacing (center-center)
	182	dz_soil_stag, & !: soil grid spacing (edge-edge)
	183	root_extr = 0.0_wp, & !: root extraction
[1551]	184	root_fraction = (/9999999.9_wp, 9999999.9_wp, &
	185	9999999.9_wp, 9999999.9_wp /), & !: distribution of root surface area to the individual soil layers
	186	zs = (/0.07_wp, 0.28_wp, 1.00_wp, 2.89_wp/), & !: soil layer depths (m)
[1500]	187	soil_moisture = 0.0_wp !: soil moisture content (m3/m3)
[1496]	188
[1551]	189	REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) :: &
[1553]	190	soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp, 283.0_wp, &
	191	283.0_wp /) !: soil temperature (K)
[1496]	192
	193	#if defined( __nopointer )
[1551]	194	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface, & !: surface temperature (K)
	195	t_surface_p, & !: progn. surface temperature (K)
	196	m_liq_eb, & !: liquid water reservoir (m)
	197	m_liq_eb_av, & !: liquid water reservoir (m)
	198	m_liq_eb_p !: progn. liquid water reservoir (m)
[1496]	199	#else
[1551]	200	REAL(wp), DIMENSION(:,:), POINTER :: t_surface, &
	201	t_surface_p, &
	202	m_liq_eb, &
	203	m_liq_eb_p
[1496]	204
[1551]	205	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
	206	m_liq_eb_av, &
	207	m_liq_eb_1, m_liq_eb_2
[1496]	208	#endif
	209
	210	!
	211	!-- Temporal tendencies for time stepping
[1551]	212	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m, & !: surface temperature tendency (K)
	213	tm_liq_eb_m !: liquid water reservoir tendency (m)
[1496]	214
	215	!
	216	!-- Energy balance variables
	217	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
[1551]	218	alpha_vg, & !: coef. of Van Genuchten
	219	c_liq, & !: liquid water coverage (of vegetated area)
	220	c_liq_av, & !: average of c_liq
	221	c_soil_av, & !: average of c_soil
	222	c_veg, & !: vegetation coverage
	223	c_veg_av, & !: average of c_veg
	224	f_sw_in, & !: fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
	225	ghf_eb, & !: ground heat flux
	226	ghf_eb_av, & !: average of ghf_eb
	227	gamma_w_sat, & !: hydraulic conductivity at saturation
	228	g_d, & !: coefficient for dependence of r_canopy on water vapour pressure deficit
	229	lai, & !: leaf area index
	230	lai_av, & !: average of lai
	231	lambda_h_sat, & !: heat conductivity for dry soil
	232	lambda_surface_s, & !: coupling between surface and soil (depends on vegetation type)
	233	lambda_surface_u, & !: coupling between surface and soil (depends on vegetation type)
	234	l_vg, & !: coef. of Van Genuchten
	235	m_fc, & !: soil moisture at field capacity (m3/m3)
	236	m_res, & !: residual soil moisture
	237	m_sat, & !: saturation soil moisture (m3/m3)
	238	m_wilt, & !: soil moisture at permanent wilting point (m3/m3)
	239	n_vg, & !: coef. Van Genuchten
	240	qsws_eb, & !: surface flux of latent heat (total)
	241	qsws_eb_av, & !: average of qsws_eb
	242	qsws_liq_eb, & !: surface flux of latent heat (liquid water portion)
	243	qsws_liq_eb_av, & !: average of qsws_liq_eb
	244	qsws_soil_eb, & !: surface flux of latent heat (soil portion)
	245	qsws_soil_eb_av, & !: average of qsws_soil_eb
	246	qsws_veg_eb, & !: surface flux of latent heat (vegetation portion)
	247	qsws_veg_eb_av, & !: average of qsws_veg_eb
[1585]	248	rad_net_l, & !: local copy of rad_net (net radiation at surface)
[1551]	249	r_a, & !: aerodynamic resistance
[1555]	250	r_a_av, & !: avergae of r_a
[1551]	251	r_canopy, & !: canopy resistance
	252	r_soil, & !: soil resitance
	253	r_soil_min, & !: minimum soil resistance
[1555]	254	r_s, & !: total surface resistance (combination of r_soil and r_canopy)
	255	r_s_av, & !: avergae of r_s
[1551]	256	r_canopy_min, & !: minimum canopy (stomatal) resistance
	257	shf_eb, & !: surface flux of sensible heat
	258	shf_eb_av !: average of shf_eb
[1496]	259
	260	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
	261	lambda_h, & !: heat conductivity of soil (?)
	262	lambda_w, & !: hydraulic diffusivity of soil (?)
	263	gamma_w, & !: hydraulic conductivity of soil (?)
[1551]	264	rho_c_total !: volumetric heat capacity of the actual soil matrix (?)
[1496]	265
	266	#if defined( __nopointer )
	267	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
[1551]	268	t_soil, & !: Soil temperature (K)
	269	t_soil_av, & !: Average of t_soil
	270	t_soil_p, & !: Prog. soil temperature (K)
[1496]	271	m_soil, & !: Soil moisture (m3/m3)
[1551]	272	m_soil_av, & !: Average of m_soil
[1496]	273	m_soil_p !: Prog. soil moisture (m3/m3)
	274	#else
	275	REAL(wp), DIMENSION(:,:,:), POINTER :: &
[1551]	276	t_soil, t_soil_p, &
[1496]	277	m_soil, m_soil_p
	278
	279	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
[1551]	280	t_soil_av, t_soil_1, t_soil_2, &
	281	m_soil_av, m_soil_1, m_soil_2
[1496]	282	#endif
	283
	284
	285	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
[1551]	286	tt_soil_m, & !: t_soil storage array
[1496]	287	tm_soil_m, & !: m_soil storage array
	288	root_fr !: root fraction (sum=1)
	289
[1551]	290
[1496]	291	!
[1551]	292	!-- Predefined Land surface classes (veg_type)
[1585]	293	CHARACTER(26), DIMENSION(0:19), PARAMETER :: veg_type_name = (/ &
[1590]	294	'user defined ', & ! 0
	295	'crops, mixed farming ', & ! 1
	296	'short grass ', & ! 2
[1585]	297	'evergreen needleleaf trees', & ! 3
	298	'deciduous needleleaf trees', & ! 4
[1590]	299	'evergreen broadleaf trees ', & ! 5
	300	'deciduous broadleaf trees ', & ! 6
	301	'tall grass ', & ! 7
	302	'desert ', & ! 8
	303	'tundra ', & ! 9
	304	'irrigated crops ', & ! 10
	305	'semidesert ', & ! 11
	306	'ice caps and glaciers ', & ! 12
	307	'bogs and marshes ', & ! 13
	308	'inland water ', & ! 14
	309	'ocean ', & ! 15
	310	'evergreen shrubs ', & ! 16
	311	'deciduous shrubs ', & ! 17
	312	'mixed forest/woodland ', & ! 18
	313	'interrupted forest ' & ! 19
[1585]	314	/)
[1496]	315
	316	!
[1551]	317	!-- Soil model classes (soil_type)
[1585]	318	CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
	319	'user defined', & ! 0
[1590]	320	'coarse ', & ! 1
	321	'medium ', & ! 2
	322	'medium-fine ', & ! 3
	323	'fine ', & ! 4
	324	'very fine ', & ! 5
	325	'organic ', & ! 6
	326	'loamy (CH) ' & ! 7
[1585]	327	/)
[1551]	328	!
	329	!-- Land surface parameters according to the respective classes (veg_type)
	330
	331	!
	332	!-- Land surface parameters I
	333	!-- r_canopy_min, lai, c_veg, g_d
[1585]	334	REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: veg_pars = RESHAPE( (/ &
	335	180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 1
	336	110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & ! 2
	337	500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 3
	338	500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 4
	339	175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 5
	340	240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & ! 6
	341	100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & ! 7
	342	250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, & ! 8
	343	80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & ! 9
	344	180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10
	345	150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, & ! 11
	346	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
	347	240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, & ! 13
	348	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
	349	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
	350	225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, & ! 16
	351	225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17
	352	250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18
	353	175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp & ! 19
[1496]	354	/), (/ 4, 19 /) )
	355
	356	!
	357	!-- Land surface parameters II z0, z0h
[1585]	358	REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ &
	359	0.25_wp, 0.25E-2_wp, & ! 1
	360	0.20_wp, 0.20E-2_wp, & ! 2
	361	2.00_wp, 2.00_wp, & ! 3
	362	2.00_wp, 2.00_wp, & ! 4
	363	2.00_wp, 2.00_wp, & ! 5
	364	2.00_wp, 2.00_wp, & ! 6
	365	0.47_wp, 0.47E-2_wp, & ! 7
	366	0.013_wp, 0.013E-2_wp, & ! 8
	367	0.034_wp, 0.034E-2_wp, & ! 9
	368	0.5_wp, 0.50E-2_wp, & ! 10
	369	0.17_wp, 0.17E-2_wp, & ! 11
	370	1.3E-3_wp, 1.3E-4_wp, & ! 12
	371	0.83_wp, 0.83E-2_wp, & ! 13
	372	0.00_wp, 0.00E-2_wp, & ! 14
	373	0.00_wp, 0.00E-2_wp, & ! 15
	374	0.10_wp, 0.10E-2_wp, & ! 16
	375	0.25_wp, 0.25E-2_wp, & ! 17
	376	2.00_wp, 2.00E-2_wp, & ! 18
	377	1.10_wp, 1.10E-2_wp & ! 19
[1496]	378	/), (/ 2, 19 /) )
	379
	380	!
[1551]	381	!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
[1585]	382	REAL(wp), DIMENSION(0:2,1:19), PARAMETER :: surface_pars = RESHAPE( (/ &
	383	10.0_wp, 10.0_wp, 0.05_wp, & ! 1
	384	10.0_wp, 10.0_wp, 0.05_wp, & ! 2
	385	20.0_wp, 15.0_wp, 0.03_wp, & ! 3
	386	20.0_wp, 15.0_wp, 0.03_wp, & ! 4
	387	20.0_wp, 15.0_wp, 0.03_wp, & ! 5
	388	20.0_wp, 15.0_wp, 0.03_wp, & ! 6
	389	10.0_wp, 10.0_wp, 0.05_wp, & ! 7
	390	15.0_wp, 15.0_wp, 0.00_wp, & ! 8
	391	10.0_wp, 10.0_wp, 0.05_wp, & ! 9
	392	10.0_wp, 10.0_wp, 0.05_wp, & ! 10
	393	10.0_wp, 10.0_wp, 0.05_wp, & ! 11
	394	58.0_wp, 58.0_wp, 0.00_wp, & ! 12
	395	10.0_wp, 10.0_wp, 0.05_wp, & ! 13
	396	1.0E20_wp, 1.0E20_wp, 0.00_wp, & ! 14
	397	1.0E20_wp, 1.0E20_wp, 0.00_wp, & ! 15
	398	10.0_wp, 10.0_wp, 0.05_wp, & ! 16
	399	10.0_wp, 10.0_wp, 0.05_wp, & ! 17
	400	20.0_wp, 15.0_wp, 0.03_wp, & ! 18
	401	20.0_wp, 15.0_wp, 0.03_wp & ! 19
[1496]	402	/), (/ 3, 19 /) )
	403
	404	!
	405	!-- Root distribution (sum = 1) level 1, level 2, level 3, level 4,
[1585]	406	REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: root_distribution = RESHAPE( (/ &
	407	0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 1
	408	0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 2
	409	0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & ! 3
	410	0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & ! 4
	411	0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & ! 5
	412	0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & ! 6
	413	0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & ! 7
	414	1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 8
	415	0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 9
	416	0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 10
	417	0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 11
	418	0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
	419	0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 13
	420	0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
	421	0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
	422	0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16
	423	0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17
	424	0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18
	425	0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 19
[1496]	426	/), (/ 4, 19 /) )
	427
	428	!
	429	!-- Soil parameters according to the following porosity classes (soil_type)
[1551]	430
[1496]	431	!
[1551]	432	!-- Soil parameters I alpha_vg, l_vg, n_vg, gamma_w_sat
[1585]	433	REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/ &
[1496]	434	3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, & ! 1
	435	3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, & ! 2
	436	0.83_wp, -0.588_wp, 1.25_wp, 0.26E-6_wp, & ! 3
	437	3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, & ! 4
	438	2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, & ! 5
[1551]	439	1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp, & ! 6
	440	0.00_wp, 0.00_wp, 0.00_wp, 0.57E-6_wp & ! 7
	441	/), (/ 4, 7 /) )
[1496]	442
	443	!
	444	!-- Soil parameters II m_sat, m_fc, m_wilt, m_res
[1585]	445	REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
[1496]	446	0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
	447	0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
	448	0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
	449	0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
	450	0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
[1551]	451	0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
	452	0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp & ! 7
	453	/), (/ 4, 7 /) )
[1496]	454
	455
	456	SAVE
	457
	458
	459	PRIVATE
	460
	461
[1551]	462	!
	463	!-- Public parameters, constants and initial values
	464	PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient, &
	465	conserve_water_content, dewfall, field_capacity, &
	466	f_shortwave_incoming, hydraulic_conductivity, init_lsm, &
	467	init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable, &
	468	land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model, &
	469	lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance, &
	470	min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp, &
	471	rho_lv, root_fraction, saturation_moisture, soil_moisture, &
	472	soil_temperature, soil_type, soil_type_name, vegetation_coverage, &
	473	veg_type, veg_type_name, wilting_point, z0_eb, z0h_eb
[1496]	474
[1551]	475	!
	476	!-- Public grid and NetCDF variables
	477	PUBLIC dots_soil, id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz, &
	478	id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy, id_var_zs_xz, &
	479	id_var_zs_yz, id_var_zs_3d, id_var_zs_mask, nzb_soil, nzs, nzt_soil,&
	480	zs
[1496]	481
	482
[1551]	483	!
	484	!-- Public 2D output variables
	485	PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av, &
	486	lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av, &
	487	qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av, &
[1555]	488	r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
[1551]	489
	490
	491	!
	492	!-- Public prognostic variables
	493	PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
	494
[1496]	495	INTERFACE init_lsm
	496	MODULE PROCEDURE init_lsm
	497	END INTERFACE init_lsm
	498
	499	INTERFACE lsm_energy_balance
	500	MODULE PROCEDURE lsm_energy_balance
	501	END INTERFACE lsm_energy_balance
	502
	503	INTERFACE lsm_soil_model
	504	MODULE PROCEDURE lsm_soil_model
	505	END INTERFACE lsm_soil_model
	506
[1551]	507	INTERFACE lsm_swap_timelevel
	508	MODULE PROCEDURE lsm_swap_timelevel
	509	END INTERFACE lsm_swap_timelevel
[1496]	510
	511	CONTAINS
	512
	513
	514	!------------------------------------------------------------------------------!
	515	! Description:
	516	! ------------
[1551]	517	!-- Allocate land surface model arrays and define pointers
[1496]	518	!------------------------------------------------------------------------------!
[1551]	519	SUBROUTINE init_lsm_arrays
[1496]	520
	521
	522	IMPLICIT NONE
	523
	524	!
[1551]	525	!-- Allocate surface and soil temperature / humidity
[1496]	526	#if defined( __nopointer )
[1551]	527	ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
	528	ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
	529	ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
	530	ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
	531	ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
	532	ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
	533	ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
	534	ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
[1496]	535	#else
[1551]	536	ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
	537	ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
	538	ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
	539	ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
	540	ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
	541	ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
	542	ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
	543	ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
[1496]	544	#endif
	545
	546	!
[1551]	547	!-- Allocate intermediate timestep arrays
	548	ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
	549	ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
	550	ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
	551	ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
[1496]	552
	553	!
	554	!-- Allocate 2D vegetation model arrays
[1551]	555	ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
[1496]	556	ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
	557	ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
[1551]	558	ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
	559	ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
[1496]	560	ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
[1551]	561	ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
	562	ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
	563	ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
[1496]	564	ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
[1551]	565	ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
	566	ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
[1496]	567	ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
	568	ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
	569	ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
	570	ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
[1551]	571	ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
	572	ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
	573	ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
	574	ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
	575	ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
[1585]	576	ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
[1496]	577	ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
	578	ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
	579	ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
	580	ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
	581	ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
[1551]	582	ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
	583	ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
[1496]	584
[1551]	585	#if ! defined( __nopointer )
[1496]	586	!
[1585]	587	!-- Initial assignment of the pointers
[1551]	588	t_soil => t_soil_1; t_soil_p => t_soil_2
	589	t_surface => t_surface_1; t_surface_p => t_surface_2
	590	m_soil => m_soil_1; m_soil_p => m_soil_2
	591	m_liq_eb => m_liq_eb_1; m_liq_eb_p => m_liq_eb_2
	592	#endif
[1496]	593
	594
[1551]	595	END SUBROUTINE init_lsm_arrays
[1500]	596
[1551]	597	!------------------------------------------------------------------------------!
	598	! Description:
	599	! ------------
	600	!-- Initialization of the land surface model
	601	!------------------------------------------------------------------------------!
	602	SUBROUTINE init_lsm
	603
	604
	605	IMPLICIT NONE
	606
	607	INTEGER(iwp) :: i !: running index
	608	INTEGER(iwp) :: j !: running index
	609	INTEGER(iwp) :: k !: running index
	610
	611
	612	!
	613	!-- Calculate Exner function
	614	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
	615
	616
	617	!
	618	!-- If no cloud physics is used, rho_surface has not been calculated before
	619	IF ( .NOT. cloud_physics ) THEN
	620	rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
	621	ENDIF
	622
	623	!
	624	!-- Calculate frequently used parameters
	625	rho_cp = cp * rho_surface
	626	rd_d_rv = r_d / r_v
	627	rho_lv = rho_surface * l_v
	628	drho_l_lv = 1.0_wp / (rho_l * l_v)
	629
	630	!
	631	!-- Set inital values for prognostic quantities
	632	tt_surface_m = 0.0_wp
	633	tt_soil_m = 0.0_wp
	634	tm_liq_eb_m = 0.0_wp
	635	c_liq = 0.0_wp
	636
	637	ghf_eb = 0.0_wp
	638	shf_eb = rho_cp * shf
	639
[1500]	640	IF ( humidity ) THEN
[1551]	641	qsws_eb = rho_l * l_v * qsws
[1500]	642	ELSE
[1551]	643	qsws_eb = 0.0_wp
[1500]	644	ENDIF
	645
[1551]	646	qsws_liq_eb = 0.0_wp
	647	qsws_soil_eb = qsws_eb
	648	qsws_veg_eb = 0.0_wp
[1496]	649
	650	r_a = 50.0_wp
[1551]	651	r_s = 50.0_wp
[1496]	652	r_canopy = 0.0_wp
	653	r_soil = 0.0_wp
	654
	655	!
	656	!-- Allocate 3D soil model arrays
[1551]	657	ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
	658	ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
	659	ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
[1496]	660
	661	lambda_h = 0.0_wp
	662	!
	663	!-- If required, allocate humidity-related variables for the soil model
	664	IF ( humidity ) THEN
[1551]	665	ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
	666	ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
[1496]	667
	668	lambda_w = 0.0_wp
	669	ENDIF
	670
	671	!
	672	!-- Calculate grid spacings. Temperature and moisture are defined at
	673	!-- the center of the soil layers, whereas gradients/fluxes are defined
	674	!-- at the edges (_stag)
[1551]	675	dz_soil_stag(nzb_soil) = zs(nzb_soil)
[1496]	676
[1551]	677	DO k = nzb_soil+1, nzt_soil
	678	dz_soil_stag(k) = zs(k) - zs(k-1)
[1496]	679	ENDDO
	680
[1551]	681	DO k = nzb_soil, nzt_soil-1
[1496]	682	dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k))
	683	ENDDO
[1551]	684	dz_soil(nzt_soil) = dz_soil_stag(nzt_soil)
[1496]	685
[1551]	686	ddz_soil = 1.0_wp / dz_soil
	687	ddz_soil_stag = 1.0_wp / dz_soil_stag
	688
[1496]	689	!
[1551]	690	!-- Initialize standard soil types. It is possible to overwrite each
	691	!-- parameter by setting the respecticy NAMELIST variable to a
	692	!-- value /= 9999999.9.
	693	IF ( soil_type /= 0 ) THEN
	694
	695	IF ( alpha_vangenuchten == 9999999.9_wp ) THEN
	696	alpha_vangenuchten = soil_pars(0,soil_type)
	697	ENDIF
	698
	699	IF ( l_vangenuchten == 9999999.9_wp ) THEN
	700	l_vangenuchten = soil_pars(1,soil_type)
	701	ENDIF
	702
	703	IF ( n_vangenuchten == 9999999.9_wp ) THEN
	704	n_vangenuchten = soil_pars(2,soil_type)
	705	ENDIF
	706
	707	IF ( hydraulic_conductivity == 9999999.9_wp ) THEN
	708	hydraulic_conductivity = soil_pars(3,soil_type)
	709	ENDIF
	710
	711	IF ( saturation_moisture == 9999999.9_wp ) THEN
	712	saturation_moisture = m_soil_pars(0,soil_type)
	713	ENDIF
	714
	715	IF ( field_capacity == 9999999.9_wp ) THEN
	716	field_capacity = m_soil_pars(1,soil_type)
	717	ENDIF
	718
	719	IF ( wilting_point == 9999999.9_wp ) THEN
	720	wilting_point = m_soil_pars(2,soil_type)
	721	ENDIF
	722
	723	IF ( residual_moisture == 9999999.9_wp ) THEN
	724	residual_moisture = m_soil_pars(3,soil_type)
	725	ENDIF
	726
[1496]	727	ENDIF
	728
[1551]	729	alpha_vg = alpha_vangenuchten
	730	l_vg = l_vangenuchten
	731	n_vg = n_vangenuchten
	732	gamma_w_sat = hydraulic_conductivity
	733	m_sat = saturation_moisture
	734	m_fc = field_capacity
	735	m_wilt = wilting_point
	736	m_res = residual_moisture
	737	r_soil_min = min_soil_resistance
	738
[1496]	739	!
[1551]	740	!-- Initial run actions
	741	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
[1496]	742
[1551]	743	t_soil = 0.0_wp
	744	m_liq_eb = 0.0_wp
	745	m_soil = 0.0_wp
[1496]	746
[1551]	747	!
	748	!-- Map user settings of T and q for each soil layer
	749	!-- (make sure that the soil moisture does not drop below the permanent
	750	!-- wilting point) -> problems with devision by zero)
	751	DO k = nzb_soil, nzt_soil
	752	t_soil(k,:,:) = soil_temperature(k)
	753	m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:))
	754	soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
	755	ENDDO
	756	t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
[1496]	757
[1551]	758	!
	759	!-- Calculate surface temperature
	760	t_surface = pt_surface * exn
[1496]	761
	762	!
[1551]	763	!-- Set artifical values for ts and us so that r_a has its initial value for
	764	!-- the first time step
	765	DO i = nxlg, nxrg
	766	DO j = nysg, nyng
	767	k = nzb_s_inner(j,i)
	768
	769	!
	770	!-- Assure that r_a cannot be zero at model start
	771	IF ( pt(k+1,j,i) == pt(k,j,i) ) pt(k+1,j,i) = pt(k+1,j,i) &
	772	+ 1.0E-10_wp
	773
	774	us(j,i) = 0.1_wp
	775	ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
	776	shf(j,i) = - us(j,i) * ts(j,i)
	777	ENDDO
	778	ENDDO
	779
	780	!
	781	!-- Actions for restart runs
	782	ELSE
	783
	784	DO i = nxlg, nxrg
	785	DO j = nysg, nyng
	786	k = nzb_s_inner(j,i)
	787	t_surface(j,i) = pt(k,j,i) * exn
	788	ENDDO
	789	ENDDO
	790
	791	ENDIF
	792
	793	!
[1496]	794	!-- Calculate saturation soil heat conductivity
	795	lambda_h_sat(:,:) = lambda_h_sm ** (1.0_wp - m_sat(:,:)) * &
	796	lambda_h_water ** m_sat(:,:)
	797
	798
	799
	800
[1551]	801	DO k = nzb_soil, nzt_soil
	802	root_fr(k,:,:) = root_fraction(k)
	803	ENDDO
[1496]	804
[1551]	805	IF ( veg_type /= 0 ) THEN
	806	IF ( min_canopy_resistance == 9999999.9_wp ) THEN
	807	min_canopy_resistance = veg_pars(0,veg_type)
	808	ENDIF
	809	IF ( leaf_area_index == 9999999.9_wp ) THEN
	810	leaf_area_index = veg_pars(1,veg_type)
	811	ENDIF
	812	IF ( vegetation_coverage == 9999999.9_wp ) THEN
	813	vegetation_coverage = veg_pars(2,veg_type)
	814	ENDIF
	815	IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN
	816	canopy_resistance_coefficient= veg_pars(3,veg_type)
	817	ENDIF
	818	IF ( lambda_surface_stable == 9999999.9_wp ) THEN
	819	lambda_surface_stable = surface_pars(0,veg_type)
	820	ENDIF
	821	IF ( lambda_surface_unstable == 9999999.9_wp ) THEN
	822	lambda_surface_unstable = surface_pars(1,veg_type)
	823	ENDIF
	824	IF ( f_shortwave_incoming == 9999999.9_wp ) THEN
	825	f_shortwave_incoming = surface_pars(2,veg_type)
	826	ENDIF
	827	IF ( z0_eb == 9999999.9_wp ) THEN
	828	roughness_length = roughness_par(0,veg_type)
	829	z0_eb = roughness_par(0,veg_type)
	830	ENDIF
	831	IF ( z0h_eb == 9999999.9_wp ) THEN
	832	z0h_eb = roughness_par(1,veg_type)
	833	ENDIF
	834	z0h_factor = z0h_eb / z0_eb
[1496]	835
[1551]	836	IF ( ANY( root_fraction == 9999999.9_wp ) ) THEN
	837	DO k = nzb_soil, nzt_soil
	838	root_fr(k,:,:) = root_distribution(k,veg_type)
	839	root_fraction(k) = root_distribution(k,veg_type)
	840	ENDDO
	841	ENDIF
[1496]	842
[1553]	843	ELSE
	844
	845	IF ( z0_eb == 9999999.9_wp ) THEN
	846	z0_eb = roughness_length
	847	ENDIF
	848	IF ( z0h_eb == 9999999.9_wp ) THEN
	849	z0h_eb = z0_eb * z0h_factor
	850	ENDIF
	851
[1496]	852	ENDIF
	853
	854	!
[1551]	855	!-- Initialize vegetation
	856	r_canopy_min = min_canopy_resistance
	857	lai = leaf_area_index
	858	c_veg = vegetation_coverage
	859	g_d = canopy_resistance_coefficient
	860	lambda_surface_s = lambda_surface_stable
	861	lambda_surface_u = lambda_surface_unstable
	862	f_sw_in = f_shortwave_incoming
	863	z0 = z0_eb
	864	z0h = z0h_eb
	865
	866	!
[1496]	867	!-- Possibly do user-defined actions (e.g. define heterogeneous land surface)
	868	CALL user_init_land_surface
	869
[1500]	870
[1551]	871	t_soil_p = t_soil
	872	m_soil_p = m_soil
	873	m_liq_eb_p = m_liq_eb
[1496]	874
[1551]	875	!-- Store initial profiles of t_soil and m_soil (assuming they are
	876	!-- horizontally homogeneous on this PE)
	877	hom(nzb_soil:nzt_soil,1,90,:) = SPREAD( t_soil(nzb_soil:nzt_soil, &
	878	nysg,nxlg), 2, &
	879	statistic_regions+1 )
	880	hom(nzb_soil:nzt_soil,1,92,:) = SPREAD( m_soil(nzb_soil:nzt_soil, &
	881	nysg,nxlg), 2, &
	882	statistic_regions+1 )
	883
[1496]	884	!
	885	!-- Calculate humidity at the surface
	886	IF ( humidity ) THEN
[1551]	887	CALL calc_q_surface
[1496]	888	ENDIF
	889
[1551]	890
	891
	892	!
	893	!-- Add timeseries for land surface model
[1585]	894
[1551]	895	dots_label(dots_num+1) = "ghf_eb"
	896	dots_label(dots_num+2) = "shf_eb"
	897	dots_label(dots_num+3) = "qsws_eb"
	898	dots_label(dots_num+4) = "qsws_liq_eb"
	899	dots_label(dots_num+5) = "qsws_soil_eb"
	900	dots_label(dots_num+6) = "qsws_veg_eb"
	901	dots_unit(dots_num+1:dots_num+6) = "W/m2"
	902
[1555]	903	dots_label(dots_num+7) = "r_a"
	904	dots_label(dots_num+8) = "r_s"
	905	dots_unit(dots_num+7:dots_num+8) = "s/m"
	906
[1551]	907	dots_soil = dots_num + 1
[1555]	908	dots_num = dots_num + 8
[1551]	909
[1496]	910	RETURN
	911
	912	END SUBROUTINE init_lsm
	913
	914
	915
	916	!------------------------------------------------------------------------------!
	917	! Description:
	918	! ------------
[1551]	919	! Solver for the energy balance at the surface.
[1496]	920	!
[1551]	921	! Note: surface fluxes are calculated in the land surface model, but these are
	922	! not used in the atmospheric code. The fluxes are calculated afterwards in
	923	! prandtl_fluxes using the surface values of temperature and humidity as
	924	! provided by the land surface model. In this way, the fluxes in the land
	925	! surface model are not equal to the ones calculated in prandtl_fluxes
[1496]	926	!------------------------------------------------------------------------------!
	927	SUBROUTINE lsm_energy_balance
	928
	929
	930	IMPLICIT NONE
	931
	932	INTEGER(iwp) :: i !: running index
	933	INTEGER(iwp) :: j !: running index
	934	INTEGER(iwp) :: k, ks !: running index
	935
	936	REAL(wp) :: f1, & !: resistance correction term 1
	937	f2, & !: resistance correction term 2
	938	f3, & !: resistance correction term 3
	939	m_min, & !: minimum soil moisture
	940	T_1, & !: actual temperature at first grid point
	941	e, & !: water vapour pressure
	942	e_s, & !: water vapour saturation pressure
[1571]	943	e_s_dt, & !: derivate of e_s with respect to T
[1496]	944	tend, & !: tendency
[1571]	945	dq_s_dt, & !: derivate of q_s with respect to T
[1496]	946	coef_1, & !: coef. for prognostic equation
	947	coef_2, & !: coef. for prognostic equation
[1551]	948	f_qsws, & !: factor for qsws_eb
	949	f_qsws_veg, & !: factor for qsws_veg_eb
	950	f_qsws_soil, & !: factor for qsws_soil_eb
	951	f_qsws_liq, & !: factor for qsws_liq_eb
	952	f_shf, & !: factor for shf_eb
	953	lambda_surface, & !: Current value of lambda_surface
	954	m_liq_eb_max !: maxmimum value of the liq. water reservoir
[1496]	955
[1551]	956
[1496]	957	!
	958	!-- Calculate the exner function for the current time step
	959	exn = ( surface_pressure / 1000.0_wp )**0.286_wp
	960
[1585]	961
[1496]	962
	963	DO i = nxlg, nxrg
	964	DO j = nysg, nyng
[1500]	965	k = nzb_s_inner(j,i)
[1496]	966
	967	!
[1551]	968	!-- Set lambda_surface according to stratification
[1496]	969	IF ( rif(j,i) >= 0.0_wp ) THEN
[1551]	970	lambda_surface = lambda_surface_s(j,i)
[1496]	971	ELSE
[1551]	972	lambda_surface = lambda_surface_u(j,i)
[1496]	973	ENDIF
[1500]	974
[1496]	975	!
[1500]	976	!-- First step: calculate aerodyamic resistance. As pt, us, ts
	977	!-- are not available for the prognostic time step, data from the last
	978	!-- time step is used here. Note that this formulation is the
	979	!-- equivalent to the ECMWF formulation using drag coefficients
[1551]	980	r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + &
	981	1.0E-20_wp)
[1496]	982
	983	!
	984	!-- Second step: calculate canopy resistance r_canopy
	985	!-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
	986
[1551]	987	!-- f1: correction for incoming shortwave radiation (stomata close at
	988	!-- night)
[1585]	989	IF ( radiation_scheme /= 'constant' ) THEN
	990	f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /&
	991	(0.81_wp * (0.004_wp * rad_sw_in(k,j,i) &
	992	+ 1.0_wp)) )
[1551]	993	ELSE
	994	f1 = 1.0_wp
	995	ENDIF
[1496]	996
	997	!
[1551]	998	!-- f2: correction for soil moisture availability to plants (the
	999	!-- integrated soil moisture must thus be considered here)
	1000	!-- f2 = 0 for very dry soils
[1496]	1001	m_total = 0.0_wp
[1551]	1002	DO ks = nzb_soil, nzt_soil
	1003	m_total = m_total + root_fr(ks,j,i) &
	1004	* MAX(m_soil(ks,j,i),m_wilt(j,i))
[1496]	1005	ENDDO
	1006
[1551]	1007	IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT. m_fc(j,i) ) THEN
[1496]	1008	f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
[1551]	1009	ELSEIF ( m_total .GE. m_fc(j,i) ) THEN
	1010	f2 = 1.0_wp
[1496]	1011	ELSE
	1012	f2 = 1.0E-20_wp
	1013	ENDIF
	1014
	1015	!
	1016	!-- Calculate water vapour pressure at saturation
[1551]	1017	e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) &
	1018	- 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
[1496]	1019
	1020	!
	1021	!-- f3: correction for vapour pressure deficit
[1551]	1022	IF ( g_d(j,i) /= 0.0_wp ) THEN
[1496]	1023	!
	1024	!-- Calculate vapour pressure
[1551]	1025	e = q_p(k+1,j,i) * surface_pressure / 0.622_wp
	1026	f3 = EXP ( -g_d(j,i) * (e_s - e) )
[1496]	1027	ELSE
	1028	f3 = 1.0_wp
	1029	ENDIF
	1030
	1031	!
[1551]	1032	!-- Calculate canopy resistance. In case that c_veg is 0 (bare soils),
	1033	!-- this calculation is obsolete, as r_canopy is not used below.
[1496]	1034	!-- To do: check for very dry soil -> r_canopy goes to infinity
[1551]	1035	r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3 &
	1036	+ 1.0E-20_wp)
[1496]	1037
	1038	!
[1551]	1039	!-- Third step: calculate bare soil resistance r_soil. The Clapp &
	1040	!-- Hornberger parametrization does not consider c_veg.
	1041	IF ( soil_type /= 7 ) THEN
	1042	m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * &
	1043	m_res(j,i)
	1044	ELSE
	1045	m_min = m_wilt(j,i)
	1046	ENDIF
[1496]	1047
[1551]	1048	f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
[1513]	1049	f2 = MAX(f2,1.0E-20_wp)
[1551]	1050	f2 = MIN(f2,1.0_wp)
[1496]	1051
	1052	r_soil(j,i) = r_soil_min(j,i) / f2
	1053
	1054	!
	1055	!-- Calculate fraction of liquid water reservoir
[1551]	1056	m_liq_eb_max = m_max_depth * lai(j,i)
	1057	c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
	1058
[1496]	1059
[1551]	1060	!
	1061	!-- Calculate saturation specific humidity
[1496]	1062	q_s = 0.622_wp * e_s / surface_pressure
[1500]	1063
	1064	!
[1551]	1065	!-- In case of dewfall, set evapotranspiration to zero
	1066	!-- All super-saturated water is then removed from the air
	1067	IF ( humidity .AND. dewfall .AND. q_s .LE. q_p(k+1,j,i) ) THEN
	1068	r_canopy(j,i) = 0.0_wp
	1069	r_soil(j,i) = 0.0_wp
[1496]	1070	ENDIF
	1071
	1072	!
	1073	!-- Calculate coefficients for the total evapotranspiration
[1551]	1074	f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ &
	1075	(r_a(j,i) + r_canopy(j,i))
[1496]	1076
[1551]	1077	f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + &
	1078	r_soil(j,i))
	1079	f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
[1496]	1080
[1551]	1081
[1500]	1082	!
	1083	!-- If soil moisture is below wilting point, plants do no longer
	1084	!-- transpirate.
[1551]	1085	! IF ( m_soil(k,j,i) .LT. m_wilt(j,i) ) THEN
	1086	! f_qsws_veg = 0.0_wp
	1087	! ENDIF
[1496]	1088
	1089	!
[1551]	1090	!-- If dewfall is deactivated, vegetation, soil and liquid water
	1091	!-- reservoir are not allowed to take up water from the super-saturated
	1092	!-- air.
	1093	IF ( humidity ) THEN
	1094	IF ( q_s .LE. q_p(k+1,j,i) ) THEN
	1095	IF ( .NOT. dewfall ) THEN
	1096	f_qsws_veg = 0.0_wp
	1097	f_qsws_soil = 0.0_wp
	1098	f_qsws_liq = 0.0_wp
	1099	ENDIF
	1100	ENDIF
	1101	ENDIF
	1102
	1103	f_shf = rho_cp / r_a(j,i)
	1104	f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
	1105
	1106	!
[1496]	1107	!-- Calculate derivative of q_s for Taylor series expansion
[1571]	1108	e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) - &
[1551]	1109	17.269_wp*(t_surface(j,i) - 273.16_wp) &
	1110	/ (t_surface(j,i) - 35.86_wp)**2 )
[1496]	1111
[1571]	1112	dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
[1496]	1113
	1114	T_1 = pt_p(k+1,j,i) * exn
	1115
	1116	!
	1117	!-- Add LW up so that it can be removed in prognostic equation
[1585]	1118	rad_net_l(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
[1496]	1119
[1500]	1120	IF ( humidity ) THEN
	1121
[1496]	1122	!
[1500]	1123	!-- Numerator of the prognostic equation
[1585]	1124	coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
[1551]	1125	+ f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s &
[1571]	1126	+ dq_s_dt * t_surface(j,i) ) + lambda_surface &
[1551]	1127	* t_soil(nzb_soil,j,i)
[1496]	1128
	1129	!
[1500]	1130	!-- Denominator of the prognostic equation
[1571]	1131	coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws &
	1132	* dq_s_dt + lambda_surface + f_shf / exn
[1496]	1133
[1500]	1134	ELSE
	1135
	1136	!
	1137	!-- Numerator of the prognostic equation
[1585]	1138	coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
[1551]	1139	+ f_shf / exn * T_1 + lambda_surface &
	1140	* t_soil(nzb_soil,j,i)
[1500]	1141
	1142	!
	1143	!-- Denominator of the prognostic equation
[1571]	1144	coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 &
[1551]	1145	+ lambda_surface + f_shf / exn
[1500]	1146
	1147	ENDIF
	1148
[1496]	1149	tend = 0.0_wp
	1150
	1151	!
[1551]	1152	!-- Implicit solution when the surface layer has no heat capacity,
[1496]	1153	!-- otherwise use RK3 scheme.
[1551]	1154	t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface * &
	1155	t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
	1156	* tsc(2) )
[1496]	1157	!
	1158	!-- Add RK3 term
[1551]	1159	t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) &
	1160	* tt_surface_m(j,i)
[1496]	1161
	1162	!
	1163	!-- Calculate true tendency
[1551]	1164	tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) &
	1165	* tt_surface_m(j,i)) / (dt_3d * tsc(2))
[1496]	1166	!
[1551]	1167	!-- Calculate t_surface tendencies for the next Runge-Kutta step
[1496]	1168	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1169	IF ( intermediate_timestep_count == 1 ) THEN
[1551]	1170	tt_surface_m(j,i) = tend
[1496]	1171	ELSEIF ( intermediate_timestep_count < &
	1172	intermediate_timestep_count_max ) THEN
[1551]	1173	tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp &
	1174	* tt_surface_m(j,i)
[1496]	1175	ENDIF
	1176	ENDIF
	1177
[1551]	1178	pt_p(k,j,i) = t_surface_p(j,i) / exn
[1496]	1179	!
	1180	!-- Calculate fluxes
[1585]	1181	rad_net_l(j,i) = rad_net_l(j,i) + 3.0_wp * sigma_sb &
[1571]	1182	* t_surface(j,i)*4 - 4.0_wp sigma_sb &
[1551]	1183	* t_surface(j,i)*3 t_surface_p(j,i)
	1184	ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) &
	1185	- t_soil(nzb_soil,j,i))
	1186	shf_eb(j,i) = - f_shf * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
[1496]	1187
[1500]	1188	IF ( humidity ) THEN
[1571]	1189	qsws_eb(j,i) = - f_qsws * ( q_p(k+1,j,i) - q_s + dq_s_dt &
	1190	* t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
[1496]	1191
[1551]	1192	qsws_veg_eb(j,i) = - f_qsws_veg * ( q_p(k+1,j,i) - q_s &
[1571]	1193	+ dq_s_dt * t_surface(j,i) - dq_s_dt &
[1551]	1194	* t_surface_p(j,i) )
	1195
	1196	qsws_soil_eb(j,i) = - f_qsws_soil * ( q_p(k+1,j,i) - q_s &
[1571]	1197	+ dq_s_dt * t_surface(j,i) - dq_s_dt &
[1551]	1198	* t_surface_p(j,i) )
	1199
	1200	qsws_liq_eb(j,i) = - f_qsws_liq * ( q_p(k+1,j,i) - q_s &
[1571]	1201	+ dq_s_dt * t_surface(j,i) - dq_s_dt &
[1551]	1202	* t_surface_p(j,i) )
[1500]	1203	ENDIF
[1496]	1204
[1500]	1205	!
	1206	!-- Calculate the true surface resistance
[1551]	1207	IF ( qsws_eb(j,i) .EQ. 0.0_wp ) THEN
	1208	r_s(j,i) = 1.0E10_wp
[1496]	1209	ELSE
[1571]	1210	r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dt &
	1211	* t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) &
[1551]	1212	/ qsws_eb(j,i) - r_a(j,i)
[1496]	1213	ENDIF
	1214
	1215	!
	1216	!-- Calculate change in liquid water reservoir due to dew fall or
[1500]	1217	!-- evaporation of liquid water
	1218	IF ( humidity ) THEN
[1496]	1219	!
[1571]	1220	!-- If precipitation is activated, add rain water to qsws_liq_eb
	1221	!-- and qsws_soil_eb according the the vegetation coverage.
[1500]	1222	!-- precipitation_rate is given in mm.
	1223	IF ( precipitation ) THEN
[1571]	1224
	1225	!
	1226	!-- Add precipitation to liquid water reservoir, if possible.
	1227	!-- Otherwise, add the water to bare soil fraction.
	1228	IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max ) THEN
	1229	qsws_liq_eb(j,i) = qsws_liq_eb(j,i) &
	1230	+ c_veg(j,i) * precipitation_rate(j,i) &
	1231	* 0.001_wp * rho_l * l_v
	1232	ELSE
	1233	qsws_soil_eb(j,i) = qsws_soil_eb(j,i) &
	1234	+ c_veg(j,i) * precipitation_rate(j,i) &
	1235	* 0.001_wp * rho_l * l_v
	1236	ENDIF
	1237
	1238	!-- Add precipitation to bare soil according to the bare soil
	1239	!-- coverage.
	1240	qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp &
	1241	- c_veg(j,i)) * precipitation_rate(j,i) &
	1242	* 0.001_wp * rho_l * l_v
[1496]	1243	ENDIF
[1500]	1244	!
	1245	!-- If the air is saturated, check the reservoir water level
	1246	IF ( q_s .LE. q_p(k+1,j,i)) THEN
	1247	!
[1551]	1248	!-- Check if reservoir is full (avoid values > m_liq_eb_max)
	1249	!-- In that case, qsws_liq_eb goes to qsws_soil_eb. In this
	1250	!-- case qsws_veg_eb is zero anyway (because c_liq = 1),
	1251	!-- so that tend is zero and no further check is needed
	1252	IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max ) THEN
	1253	qsws_soil_eb(j,i) = qsws_soil_eb(j,i) &
	1254	+ qsws_liq_eb(j,i)
	1255	qsws_liq_eb(j,i) = 0.0_wp
[1500]	1256	ENDIF
[1496]	1257
	1258	!
[1551]	1259	!-- In case qsws_veg_eb becomes negative (unphysical behavior),
	1260	!-- let the water enter the liquid water reservoir as dew on the
[1500]	1261	!-- plant
[1551]	1262	IF ( qsws_veg_eb(j,i) .LT. 0.0_wp ) THEN
	1263	qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
	1264	qsws_veg_eb(j,i) = 0.0_wp
[1500]	1265	ENDIF
	1266	ENDIF
[1496]	1267
[1551]	1268	tend = - qsws_liq_eb(j,i) * drho_l_lv
[1496]	1269
[1551]	1270	m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend &
	1271	+ tsc(3) * tm_liq_eb_m(j,i) )
[1496]	1272
	1273	!
[1500]	1274	!-- Check if reservoir is overfull -> reduce to maximum
	1275	!-- (conservation of water is violated here)
[1551]	1276	m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
[1496]	1277
	1278	!
[1500]	1279	!-- Check if reservoir is empty (avoid values < 0.0)
	1280	!-- (conservation of water is violated here)
[1551]	1281	m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
[1496]	1282
	1283
	1284	!
[1551]	1285	!-- Calculate m_liq_eb tendencies for the next Runge-Kutta step
[1500]	1286	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1287	IF ( intermediate_timestep_count == 1 ) THEN
[1551]	1288	tm_liq_eb_m(j,i) = tend
[1500]	1289	ELSEIF ( intermediate_timestep_count < &
	1290	intermediate_timestep_count_max ) THEN
[1551]	1291	tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp &
	1292	* tm_liq_eb_m(j,i)
[1500]	1293	ENDIF
[1496]	1294	ENDIF
	1295
[1500]	1296	ENDIF
	1297
[1496]	1298	ENDDO
[1500]	1299	ENDDO
[1496]	1300
	1301	END SUBROUTINE lsm_energy_balance
	1302
	1303
	1304	!------------------------------------------------------------------------------!
	1305	! Description:
	1306	! ------------
	1307	!
[1551]	1308	! Soil model as part of the land surface model. The model predicts soil
	1309	! temperature and water content.
[1496]	1310	!------------------------------------------------------------------------------!
	1311	SUBROUTINE lsm_soil_model
	1312
	1313
	1314	IMPLICIT NONE
	1315
	1316	INTEGER(iwp) :: i !: running index
	1317	INTEGER(iwp) :: j !: running index
	1318	INTEGER(iwp) :: k !: running index
	1319
[1551]	1320	REAL(wp) :: h_vg !: Van Genuchten coef. h
[1496]	1321
[1551]	1322	REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp, & !: temp. gamma
[1496]	1323	lambda_temp, & !: temp. lambda
	1324	tend !: tendency
	1325
	1326	DO i = nxlg, nxrg
	1327	DO j = nysg, nyng
[1551]	1328	DO k = nzb_soil, nzt_soil
[1496]	1329	!
	1330	!-- Calculate volumetric heat capacity of the soil, taking into
	1331	!-- account water content
[1551]	1332	rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
	1333	+ rho_c_water * m_soil(k,j,i))
[1496]	1334
	1335	!
	1336	!-- Calculate soil heat conductivity at the center of the soil
	1337	!-- layers
[1551]	1338	ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
	1339	lambda_temp(k) = ke * (lambda_h_sat(j,i) + lambda_h_dry) + &
[1496]	1340	lambda_h_dry
	1341
	1342	ENDDO
	1343
	1344	!
	1345	!-- Calculate soil heat conductivity (lambda_h) at the _stag level
	1346	!-- using linear interpolation
[1551]	1347	DO k = nzb_soil, nzt_soil-1
[1496]	1348
	1349	lambda_h(k,j,i) = lambda_temp(k) + &
	1350	( lambda_temp(k+1) - lambda_temp(k) ) &
	1351	* 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
	1352
	1353	ENDDO
[1551]	1354	lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
[1496]	1355
	1356	!
[1551]	1357	!-- Prognostic equation for soil temperature t_soil
[1496]	1358	tend(:) = 0.0_wp
[1551]	1359	tend(0) = (1.0_wp/rho_c_total(nzb_soil,j,i)) * &
	1360	( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) &
	1361	- t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil) &
	1362	+ ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
[1496]	1363
[1551]	1364	DO k = 1, nzt_soil
	1365	tend(k) = (1.0_wp/rho_c_total(k,j,i)) &
[1496]	1366	* ( lambda_h(k,j,i) &
[1551]	1367	* ( t_soil(k+1,j,i) - t_soil(k,j,i) ) &
[1496]	1368	* ddz_soil(k) &
	1369	- lambda_h(k-1,j,i) &
[1551]	1370	* ( t_soil(k,j,i) - t_soil(k-1,j,i) ) &
[1496]	1371	* ddz_soil(k-1) &
	1372	) * ddz_soil_stag(k)
	1373	ENDDO
	1374
[1551]	1375	t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i) &
[1496]	1376	+ dt_3d * ( tsc(2) &
	1377	* tend(:) + tsc(3) &
[1551]	1378	* tt_soil_m(:,j,i) )
[1496]	1379
	1380	!
[1551]	1381	!-- Calculate t_soil tendencies for the next Runge-Kutta step
[1496]	1382	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1383	IF ( intermediate_timestep_count == 1 ) THEN
[1551]	1384	DO k = nzb_soil, nzt_soil
	1385	tt_soil_m(k,j,i) = tend(k)
[1496]	1386	ENDDO
	1387	ELSEIF ( intermediate_timestep_count < &
	1388	intermediate_timestep_count_max ) THEN
[1551]	1389	DO k = nzb_soil, nzt_soil
	1390	tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp &
	1391	* tt_soil_m(k,j,i)
[1496]	1392	ENDDO
	1393	ENDIF
	1394	ENDIF
	1395
	1396
[1551]	1397	DO k = nzb_soil, nzt_soil
	1398
[1496]	1399	!
	1400	!-- Calculate soil diffusivity at the center of the soil layers
[1551]	1401	lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat &
[1496]	1402	/ m_sat(j,i) ) * ( MAX(m_soil(k,j,i), &
[1551]	1403	m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
[1496]	1404
	1405	!
[1551]	1406	!-- Parametrization of Van Genuchten
	1407	IF ( soil_type /= 7 ) THEN
	1408	!
	1409	!-- Calculate the hydraulic conductivity after Van Genuchten
	1410	!-- (1980)
	1411	h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - &
	1412	MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i) &
	1413	/ (n_vg(j,i)-1.0_wp)) - 1.0_wp &
	1414	)**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
[1496]	1415
[1551]	1416	gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + &
	1417	(alpha_vg(j,i)h_vg)n_vg(j,i))*(1.0_wp &
	1418	-1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg &
	1419	)(n_vg(j,i)-1.0_wp))2 ) &
	1420	/ ( (1.0_wp + (alpha_vg(j,i)*h_vg &
	1421	)n_vg(j,i))((1.0_wp - 1.0_wp/n_vg(j,i)) &
	1422	*(l_vg(j,i) + 2.0_wp)) )
[1496]	1423
[1551]	1424	!
	1425	!-- Parametrization of Clapp & Hornberger
	1426	ELSE
	1427	gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i) &
	1428	/ m_sat(j,i) )*(2.0_wp b_ch + 3.0_wp)
	1429	ENDIF
	1430
[1496]	1431	ENDDO
	1432
	1433
	1434	IF ( humidity ) THEN
	1435	!
	1436	!-- Calculate soil diffusivity (lambda_w) at the _stag level
[1551]	1437	!-- using linear interpolation. To do: replace this with
	1438	!-- ECMWF-IFS Eq. 8.81
	1439	DO k = nzb_soil, nzt_soil-1
[1496]	1440
	1441	lambda_w(k,j,i) = lambda_temp(k) + &
	1442	( lambda_temp(k+1) - lambda_temp(k) ) &
[1551]	1443	* 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)
[1496]	1444	gamma_w(k,j,i) = gamma_temp(k) + &
	1445	( gamma_temp(k+1) - gamma_temp(k) ) &
[1551]	1446	* 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)
[1496]	1447
	1448	ENDDO
	1449
	1450	!
	1451	!
	1452	!-- In case of a closed bottom (= water content is conserved), set
	1453	!-- hydraulic conductivity to zero to that no water will be lost
	1454	!-- in the bottom layer.
	1455	IF ( conserve_water_content ) THEN
[1551]	1456	gamma_w(nzt_soil,j,i) = 0.0_wp
[1496]	1457	ELSE
[1551]	1458	gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil)
[1496]	1459	ENDIF
	1460
[1551]	1461	!-- The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
[1496]	1462	!-- ensures the mass conservation for water. The transpiration of
	1463	!-- plants equals the cumulative withdrawals by the roots in the
	1464	!-- soil. The scheme takes into account the availability of water
	1465	!-- in the soil layers as well as the root fraction in the
[1551]	1466	!-- respective layer. Layer with moisture below wilting point will
	1467	!-- not contribute, which reflects the preference of plants to
	1468	!-- take water from moister layers.
[1496]	1469
	1470	!
[1551]	1471	!-- Calculate the root extraction (ECMWF 7.69, modified so that the
	1472	!-- sum of root_extr = 1). The energy balance solver guarantees a
	1473	!-- positive transpiration, so that there is no need for an
	1474	!-- additional check.
	1475
[1496]	1476	m_total = 0.0_wp
[1551]	1477	DO k = nzb_soil, nzt_soil
	1478	IF ( m_soil(k,j,i) .GT. m_wilt(j,i) ) THEN
	1479	m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
	1480	ENDIF
[1496]	1481	ENDDO
	1482
[1551]	1483	DO k = nzb_soil, nzt_soil
	1484	IF ( m_soil(k,j,i) .GT. m_wilt(j,i) ) THEN
	1485	root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
	1486	ELSE
	1487	root_extr(k) = 0.0_wp
	1488	ENDIF
[1496]	1489	ENDDO
	1490
	1491	!
	1492	!-- Prognostic equation for soil water content m_soil
	1493	tend(:) = 0.0_wp
[1551]	1494	tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * ( &
	1495	m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) ) &
	1496	* ddz_soil(nzb_soil) - gamma_w(nzb_soil,j,i) - ( &
	1497	root_extr(nzb_soil) * qsws_veg_eb(j,i) &
	1498	+ qsws_soil_eb(j,i) ) * drho_l_lv ) &
	1499	* ddz_soil_stag(nzb_soil)
[1496]	1500
[1551]	1501	DO k = nzb_soil+1, nzt_soil-1
[1496]	1502	tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i) &
	1503	- m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&
	1504	- lambda_w(k-1,j,i) * (m_soil(k,j,i) - &
	1505	m_soil(k-1,j,i)) * ddz_soil(k-1) &
[1551]	1506	+ gamma_w(k-1,j,i) - (root_extr(k) &
	1507	* qsws_veg_eb(j,i) * drho_l_lv) &
[1496]	1508	) * ddz_soil_stag(k)
	1509
	1510	ENDDO
[1551]	1511	tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i) &
	1512	- lambda_w(nzt_soil-1,j,i) &
	1513	* (m_soil(nzt_soil,j,i) &
	1514	- m_soil(nzt_soil-1,j,i)) &
	1515	* ddz_soil(nzt_soil-1) &
	1516	+ gamma_w(nzt_soil-1,j,i) - ( &
	1517	root_extr(nzt_soil) &
	1518	* qsws_veg_eb(j,i) * drho_l_lv ) &
	1519	) * ddz_soil_stag(nzt_soil)
[1496]	1520
[1551]	1521	m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
[1496]	1522	+ dt_3d * ( tsc(2) * tend(:) &
	1523	+ tsc(3) * tm_soil_m(:,j,i) )
	1524
	1525	!
	1526	!-- Account for dry soils (find a better solution here!)
	1527	m_soil_p(:,j,i) = MAX(m_soil_p(:,j,i),0.0_wp)
	1528
	1529	!
	1530	!-- Calculate m_soil tendencies for the next Runge-Kutta step
	1531	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1532	IF ( intermediate_timestep_count == 1 ) THEN
[1551]	1533	DO k = nzb_soil, nzt_soil
[1496]	1534	tm_soil_m(k,j,i) = tend(k)
	1535	ENDDO
	1536	ELSEIF ( intermediate_timestep_count < &
	1537	intermediate_timestep_count_max ) THEN
[1551]	1538	DO k = nzb_soil, nzt_soil
[1496]	1539	tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp &
	1540	* tm_soil_m(k,j,i)
	1541	ENDDO
	1542	ENDIF
	1543	ENDIF
	1544
	1545	ENDIF
	1546
	1547	ENDDO
	1548	ENDDO
	1549
	1550	!
	1551	!-- Calculate surface specific humidity
	1552	IF ( humidity ) THEN
[1551]	1553	CALL calc_q_surface
[1496]	1554	ENDIF
	1555
	1556
	1557	END SUBROUTINE lsm_soil_model
	1558
	1559
	1560	!------------------------------------------------------------------------------!
	1561	! Description:
	1562	! ------------
[1551]	1563	! Calculation of specific humidity of the surface layer (surface)
[1496]	1564	!------------------------------------------------------------------------------!
[1551]	1565	SUBROUTINE calc_q_surface
[1496]	1566
	1567	IMPLICIT NONE
	1568
	1569	INTEGER :: i !: running index
	1570	INTEGER :: j !: running index
	1571	INTEGER :: k !: running index
	1572	REAL(wp) :: resistance !: aerodynamic and soil resistance term
	1573
	1574	DO i = nxlg, nxrg
	1575	DO j = nysg, nyng
	1576	k = nzb_s_inner(j,i)
	1577
	1578	!
	1579	!-- Calculate water vapour pressure at saturation
[1551]	1580	e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) &
	1581	- 273.16_wp ) / ( t_surface(j,i) &
	1582	- 35.86_wp ) )
[1496]	1583
	1584	!
	1585	!-- Calculate specific humidity at saturation
	1586	q_s = 0.622_wp * e_s / surface_pressure
	1587
	1588
	1589	resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
	1590
	1591	!
	1592	!-- Calculate specific humidity at surface
	1593	q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance) &
	1594	* q_p(k+1,j,i)
	1595
	1596	ENDDO
	1597	ENDDO
	1598
[1551]	1599	END SUBROUTINE calc_q_surface
[1496]	1600
[1551]	1601	!------------------------------------------------------------------------------!
	1602	! Description:
	1603	! ------------
	1604	! Swapping of timelevels
	1605	!------------------------------------------------------------------------------!
	1606	SUBROUTINE lsm_swap_timelevel ( mod_count )
[1496]	1607
[1551]	1608	IMPLICIT NONE
	1609
	1610	INTEGER, INTENT(IN) :: mod_count
	1611
	1612	#if defined( __nopointer )
	1613
	1614	t_surface = t_surface_p
	1615	t_soil = t_soil_p
	1616	IF ( humidity ) THEN
	1617	m_soil = m_soil_p
	1618	m_liq_eb = m_liq_eb_p
	1619	ENDIF
	1620
	1621	#else
	1622
	1623	SELECT CASE ( mod_count )
	1624
	1625	CASE ( 0 )
	1626
[1585]	1627	t_surface => t_surface_1; t_surface_p => t_surface_2
	1628	t_soil => t_soil_1; t_soil_p => t_soil_2
[1551]	1629	IF ( humidity ) THEN
[1585]	1630	m_soil => m_soil_1; m_soil_p => m_soil_2
	1631	m_liq_eb => m_liq_eb_1; m_liq_eb_p => m_liq_eb_2
[1551]	1632	ENDIF
	1633
	1634
	1635	CASE ( 1 )
	1636
[1585]	1637	t_surface => t_surface_2; t_surface_p => t_surface_1
	1638	t_soil => t_soil_2; t_soil_p => t_soil_1
[1551]	1639	IF ( humidity ) THEN
[1585]	1640	m_soil => m_soil_2; m_soil_p => m_soil_1
	1641	m_liq_eb => m_liq_eb_2; m_liq_eb_p => m_liq_eb_1
[1551]	1642	ENDIF
	1643
	1644	END SELECT
	1645	#endif
	1646
	1647	END SUBROUTINE lsm_swap_timelevel
	1648
	1649
[1496]	1650	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 |