Home

Context Navigation

land_surface_model_mod.f90 @ 4258

Last change on this file since 4258 was 4258, checked in by suehring, 2 years ago
Land-surface model: Revise limitation for soil moisture in case it exceeds its saturation value; Revise initialization of soil moisture and temperature in a nested run in case dynamic input information is available. This case, the soil within the child domains can be initialized separately; As part of this revision, migrate the netcdf input of soil temperature / moisture to this module, as well as the routine to inter/extrapolate soil profiles between different grids.; Plant-canopy: Check if any LAD is prescribed when plant-canopy model is applied, in order to avoid crashes in the radiation transfer model; Surface-layer fluxes: Initialization of Obukhov length also at vertical surfaces (if allocated); Urban-surface model: Add checks to ensure that relative fractions of walls, windowns and green surfaces sum-u to one; Revise message calls dealing with local checks
Property svn:keywords set to `Id`
File size: 340.8 KB

Line
1	!> @file land_surface_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! 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-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: land_surface_model_mod.f90 4258 2019-10-07 13:29:08Z suehring $
27	! - Revise limitation for soil moisture in case it exceeds its saturation
28	! value (J. Resler)
29	! - Revise initialization of soil moisture and temperature in a nested run in
30	! case dynamic input information is available. This case, the soil within
31	! the child domains can be initialized separately. (J. Resler, M. Suehring)
32	! - As part of this revision, migrate the netcdf input of soil temperature /
33	! moisture to this module, as well as the routine to inter/extrapolate soil
34	! profiles between different grids.
35	!
36	! 4251 2019-10-02 12:07:38Z maronga
37	! Bugfix: albedo_types for vegetation_type look-up table corrected.
38	!
39	! 4201 2019-08-29 15:47:27Z suehring
40	! - Limit soil moisture to its saturation moisture and give a respective
41	! warning rather than an error.
42	! - Perform checks for soil temperature only when there is no dynamic input
43	! file for the parent or possible child domains.
44	!
45	! 4194 2019-08-28 08:09:44Z suehring
46	! Apply more strict limitation of z0 over water surfaces in case it exceeds the
47	! surface-layer height, in order to avoid instabilities.
48	!
49	! 4188 2019-08-26 14:15:47Z suehring
50	! Minor adjustment in error numbers, typos corrected
51	!
52	! 4187 2019-08-26 12:43:15Z suehring
53	! Adjust message call in case of local checks
54	!
55	! 4182 2019-08-22 15:20:23Z scharf
56	! Corrected "Former revisions" section
57	!
58	! 4118 2019-07-25 16:11:45Z suehring
59	! Initialization of soil temperature and moisture via dynamic input file only
60	! for vegetation and pavement surfaces.
61	!
62	! 4110 2019-07-22 17:05:21Z suehring
63	! Relax checks for non-consistent initialization in case static or dynamic
64	! input is provided. For example, soil_temperature or deep_soil_temperature
65	! is not mandatory any more if dynamic input is available. Also, improper
66	! settings of x_type in namelist are only checked if no static file is
67	! available.
68	!
69	! 4109 2019-07-22 17:00:34Z suehring
70	! Further revision of last commit in order to avoid any side effects when
71	! albedo type is not set in namelist and default albedo type changes.
72	!
73	! 4024 2019-06-12 14:06:46Z suehring
74	! Bugfix in albedo initialization, caused crashes in rrtmg calls
75	!
76	! 3987 2019-05-22 09:52:13Z kanani
77	! Introduce alternative switch for debug output during timestepping
78	!
79	! 3964 2019-05-09 09:48:32Z suehring
80	! In a nested child domain, distinguish between soil moisture and temperature
81	! initialization from parent via dynamic input file. Further, initialize soil
82	! moisture/temperature from dynamic input file only when initialization via
83	! 'inifor' is desired.
84	!
85	! 3943 2019-05-02 09:50:41Z maronga
86	! Removed extra blank character
87	!
88	! 3941 2019-04-30 09:48:33Z suehring
89	! Check that at least one surface type is set at surface element.
90	!
91	! 3933 2019-04-25 12:33:20Z kanani
92	! Remove unused subroutine and allocation of pt_2m, this is done in surface_mod
93	! now (surfaces%pt_2m)
94	!
95	!
96	! Changes related to global restructuring of location messages and introduction
97	! of additional debug messages
98	!
99	! 3881 2019-04-10 09:31:22Z suehring
100	! Bugfix in level 3 initialization of pavement albedo type and pavement
101	! emissivity
102	!
103	! 3868 2019-04-08 11:52:36Z suehring
104	! More strict limitation of roughness length when it is in the order of the
105	! vertical grid spacing
106	!
107	! 3856 2019-04-03 11:06:59Z suehring
108	! Bugfix in lsm_init in case no surface-fractions are provided
109	!
110	! 3847 2019-04-01 14:51:44Z suehring
111	! Adjust message-call for checks that are especially carried out locally.
112	!
113	! 3832 2019-03-28 13:16:58Z raasch
114	! instrumented with openmp directives
115	!
116	! 3786 2019-03-06 16:58:03Z raasch
117	! further unused variables removed
118	!
119	! 3767 2019-02-27 08:18:02Z raasch
120	! unused variable for file index removed from rrd-subroutines parameter list
121	!
122	! 3715 2019-02-04 17:34:55Z suehring
123	! Revise check for saturation moisture
124	!
125	! 3710 2019-01-30 18:11:19Z suehring
126	! Check if soil-, water-, pavement- and vegetation types are set within a valid
127	! range.
128	!
129	! 3692 2019-01-23 14:45:49Z suehring
130	! Revise check for soil moisture higher than its saturation value
131	!
132	! 3685 2019-01-21 01:02:11Z knoop
133	! Some interface calls moved to module_interface + cleanup
134	!
135	! 3677 2019-01-17 09:07:06Z moh.hefny
136	! Removed most_method
137	!
138	! 3655 2019-01-07 16:51:22Z knoop
139	! nopointer option removed
140	!
141	! 1496 2014-12-02 17:25:50Z maronga
142	! Initial revision
143	!
144	!
145	! Description:
146	! ------------
147	!> Land surface model, consisting of a solver for the energy balance at the
148	!> surface and a multi layer soil scheme. The scheme is similar to the TESSEL
149	!> scheme implemented in the ECMWF IFS model, with modifications according to
150	!> H-TESSEL. The implementation is based on the formulation implemented in the
151	!> DALES and UCLA-LES models.
152	!>
153	!> @todo Extensive verification energy-balance solver for vertical surfaces,
154	!> e.g. parametrization of r_a
155	!> @todo Revise single land-surface processes for vertical surfaces, e.g.
156	!> treatment of humidity, etc.
157	!> @todo Consider partial absorption of the net shortwave radiation by the
158	!> skin layer.
159	!> @todo Improve surface water parameterization
160	!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
161	!> nzt_soil=3)).
162	!> @todo Implement surface runoff model (required when performing long-term LES
163	!> with considerable precipitation.
164	!> @todo Revise calculation of f2 when wilting point is non-constant in the
165	!> soil
166	!> @todo Allow for zero soil moisture (currently, it is set to wilting point)
167	!> @note No time step criterion is required as long as the soil layers do not
168	!> become too thin.
169	!> @todo Attention, pavement_subpars_1/2 are hardcoded to 8 levels, in case
170	!> more levels are used this may cause an potential bug
171	!> @todo Routine calc_q_surface required?
172	!> @todo Allow for precipitation water to enter pavements that are semi-pervious
173	!------------------------------------------------------------------------------!
174	MODULE land_surface_model_mod
175
176	USE arrays_3d, &
177	ONLY: hyp, pt, prr, q, q_p, ql, vpt, u, v, w, hyrho, exner, d_exner
178
179	USE basic_constants_and_equations_mod, &
180	ONLY: c_p, g, lv_d_cp, l_v, kappa, magnus, rho_l, r_d, r_v, rd_d_rv
181
182	USE calc_mean_profile_mod, &
183	ONLY: calc_mean_profile
184
185	USE control_parameters, &
186	ONLY: cloud_droplets, &
187	coupling_char, &
188	coupling_start_time, &
189	debug_output, debug_output_timestep, debug_string, &
190	dt_3d, &
191	end_time, humidity, intermediate_timestep_count, &
192	initializing_actions, intermediate_timestep_count_max, &
193	land_surface, max_masks, pt_surface, &
194	rho_surface, spinup, spinup_pt_mean, spinup_time, &
195	surface_pressure, timestep_scheme, tsc, &
196	time_since_reference_point
197
198	USE cpulog, &
199	ONLY: cpu_log, log_point_s
200
201	USE indices, &
202	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
203
204	USE bulk_cloud_model_mod, &
205	ONLY: bulk_cloud_model, precipitation
206
207	USE netcdf_data_input_mod, &
208	ONLY : building_type_f, &
209	char_fill, &
210	char_lod, &
211	check_existence, &
212	close_input_file, &
213	get_attribute, &
214	get_dimension_length, &
215	get_variable, &
216	init_3d, &
217	input_file_dynamic, &
218	input_pids_dynamic, &
219	input_pids_static, &
220	inquire_num_variables, &
221	inquire_variable_names, &
222	num_var_pids, &
223	open_read_file, &
224	pids_id, &
225	pavement_pars_f, &
226	pavement_subsurface_pars_f, &
227	pavement_type_f, &
228	root_area_density_lsm_f, &
229	soil_pars_f, &
230	soil_type_f, &
231	surface_fraction_f, &
232	vars_pids, &
233	vegetation_pars_f, &
234	vegetation_type_f, &
235	water_pars_f, &
236	water_type_f
237
238	USE kinds
239
240	USE pegrid
241
242	USE radiation_model_mod, &
243	ONLY: albedo, albedo_type, emissivity, force_radiation_call, &
244	radiation, radiation_scheme, unscheduled_radiation_calls
245
246	USE statistics, &
247	ONLY: hom, statistic_regions
248
249	USE surface_mod, &
250	ONLY : ind_pav_green, ind_veg_wall, ind_wat_win, &
251	surf_lsm_h, surf_lsm_v, surf_type, surface_restore_elements
252
253	IMPLICIT NONE
254
255	TYPE surf_type_lsm
256	REAL(wp), DIMENSION(:), ALLOCATABLE :: var_1d !< 1D prognostic variable
257	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_2d !< 2D prognostic variable
258	END TYPE surf_type_lsm
259
260	!
261	!-- LSM model constants
262
263	REAL(wp), PARAMETER :: &
264	b_ch = 6.04_wp, & ! Clapp & Hornberger exponent
265	lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil (W/m/K)
266	lambda_h_sm = 3.44_wp, & ! heat conductivity of the soil matrix (W/m/K)
267	lambda_h_water = 0.57_wp, & ! heat conductivity of water (W/m/K)
268	psi_sat = -0.388_wp, & ! soil matrix potential at saturation
269	rho_c_soil = 2.19E6_wp, & ! volumetric heat capacity of soil (J/m3/K)
270	rho_c_water = 4.20E6_wp, & ! volumetric heat capacity of water (J/m3/K)
271	m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m)
272
273
274	REAL(wp), DIMENSION(0:7), PARAMETER :: dz_soil_default = & ! default soil layer configuration
275	(/ 0.01_wp, 0.02_wp, 0.04_wp, &
276	0.06_wp, 0.14_wp, 0.26_wp, &
277	0.54_wp, 1.86_wp/)
278
279	REAL(wp), DIMENSION(0:3), PARAMETER :: dz_soil_ref = & ! reference four layer soil configuration used for estimating the root fractions
280	(/ 0.07_wp, 0.21_wp, 0.72_wp, &
281	1.89_wp /)
282
283	REAL(wp), DIMENSION(0:3), PARAMETER :: zs_ref = & ! reference four layer soil configuration used for estimating the root fractions
284	(/ 0.07_wp, 0.28_wp, 1.0_wp, &
285	2.89_wp /)
286
287
288	!
289	!-- LSM variables
290	CHARACTER(10) :: surface_type = 'netcdf' !< general classification. Allowed are:
291	!< 'vegetation', 'pavement', ('building'),
292	!< 'water', and 'netcdf'
293
294
295
296	INTEGER(iwp) :: nzb_soil = 0, & !< bottom of the soil model (Earth's surface)
297	nzt_soil = 7, & !< top of the soil model
298	nzt_pavement = 0, & !< top of the pavement within the soil
299	nzs = 8, & !< number of soil layers
300	pavement_depth_level = 0, & !< default NAMELIST nzt_pavement
301	pavement_type = 1, & !< default NAMELIST pavement_type
302	soil_type = 3, & !< default NAMELIST soil_type
303	vegetation_type = 2, & !< default NAMELIST vegetation_type
304	water_type = 1 !< default NAMELISt water_type
305
306
307
308	LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model
309	constant_roughness = .FALSE., & !< use fixed/dynamic roughness lengths for water surfaces
310	force_radiation_call_l = .FALSE., & !< flag to force calling of radiation routine
311	aero_resist_kray = .TRUE. !< flag to control parametrization of aerodynamic resistance at vertical surface elements
312
313	! value 9999999.9_wp -> generic available or user-defined value must be set
314	! otherwise -> no generic variable and user setting is optional
315	REAL(wp) :: alpha_vangenuchten = 9999999.9_wp, & !< NAMELIST alpha_vg
316	canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
317	c_surface = 9999999.9_wp, & !< Surface (skin) heat capacity (J/m2/K)
318	deep_soil_temperature = 9999999.9_wp, & !< Deep soil temperature (bottom boundary condition)
319	drho_l_lv, & !< (rho_l * l_v)**-1
320	field_capacity = 9999999.9_wp, & !< NAMELIST m_fc
321	f_shortwave_incoming = 9999999.9_wp, & !< NAMELIST f_sw_in
322	hydraulic_conductivity = 9999999.9_wp, & !< NAMELIST gamma_w_sat
323	ke = 0.0_wp, & !< Kersten number
324	lambda_h_sat = 0.0_wp, & !< heat conductivity for saturated soil (W/m/K)
325	lambda_surface_stable = 9999999.9_wp, & !< NAMELIST lambda_surface_s (W/m2/K)
326	lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u (W/m2/K)
327	leaf_area_index = 9999999.9_wp, & !< NAMELIST lai
328	l_vangenuchten = 9999999.9_wp, & !< NAMELIST l_vg
329	min_canopy_resistance = 9999999.9_wp, & !< NAMELIST r_canopy_min
330	min_soil_resistance = 50.0_wp, & !< NAMELIST r_soil_min
331	m_total = 0.0_wp, & !< weighted total water content of the soil (m3/m3)
332	n_vangenuchten = 9999999.9_wp, & !< NAMELIST n_vg
333	pavement_heat_capacity = 9999999.9_wp, & !< volumetric heat capacity of pavement (e.g. roads) (J/m3/K)
334	pavement_heat_conduct = 9999999.9_wp, & !< heat conductivity for pavements (e.g. roads) (W/m/K)
335	q_s = 0.0_wp, & !< saturation water vapor mixing ratio
336	residual_moisture = 9999999.9_wp, & !< NAMELIST m_res
337	rho_cp, & !< rho_surface * cp
338	rho_lv, & !< rho_ocean * l_v
339	saturation_moisture = 9999999.9_wp, & !< NAMELIST m_sat
340	skip_time_do_lsm = 0.0_wp, & !< LSM is not called before this time
341	vegetation_coverage = 9999999.9_wp, & !< NAMELIST c_veg
342	water_temperature = 9999999.9_wp, & !< water temperature
343	wilting_point = 9999999.9_wp, & !< NAMELIST m_wilt
344	z0_vegetation = 9999999.9_wp, & !< NAMELIST z0 (lsm_par)
345	z0h_vegetation = 9999999.9_wp, & !< NAMELIST z0h (lsm_par)
346	z0q_vegetation = 9999999.9_wp, & !< NAMELIST z0q (lsm_par)
347	z0_pavement = 9999999.9_wp, & !< NAMELIST z0 (lsm_par)
348	z0h_pavement = 9999999.9_wp, & !< NAMELIST z0h (lsm_par)
349	z0q_pavement = 9999999.9_wp, & !< NAMELIST z0q (lsm_par)
350	z0_water = 9999999.9_wp, & !< NAMELIST z0 (lsm_par)
351	z0h_water = 9999999.9_wp, & !< NAMELIST z0h (lsm_par)
352	z0q_water = 9999999.9_wp !< NAMELIST z0q (lsm_par)
353
354
355	REAL(wp), DIMENSION(:), ALLOCATABLE :: ddz_soil_center, & !< 1/dz_soil_center
356	ddz_soil, & !< 1/dz_soil
357	dz_soil_center, & !< soil grid spacing (center-center)
358	zs, & !< depth of the temperature/moisute levels
359	root_extr !< root extraction
360
361
362
363	REAL(wp), DIMENSION(0:20) :: root_fraction = 9999999.9_wp, & !< (NAMELIST) distribution of root surface area to the individual soil layers
364	soil_moisture = 0.0_wp, & !< NAMELIST soil moisture content (m3/m3)
365	soil_temperature = 9999999.9_wp, & !< NAMELIST soil temperature (K) +1
366	dz_soil = 9999999.9_wp, & !< (NAMELIST) soil layer depths (spacing)
367	zs_layer = 9999999.9_wp !< soil layer depths (edge)
368
369	TYPE(surf_type_lsm), POINTER :: t_soil_h, & !< Soil temperature (K), horizontal surface elements
370	t_soil_h_p, & !< Prog. soil temperature (K), horizontal surface elements
371	m_soil_h, & !< Soil moisture (m3/m3), horizontal surface elements
372	m_soil_h_p !< Prog. soil moisture (m3/m3), horizontal surface elements
373
374	TYPE(surf_type_lsm), TARGET :: t_soil_h_1, & !<
375	t_soil_h_2, & !<
376	m_soil_h_1, & !<
377	m_soil_h_2 !<
378
379	TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
380	t_soil_v, & !< Soil temperature (K), vertical surface elements
381	t_soil_v_p, & !< Prog. soil temperature (K), vertical surface elements
382	m_soil_v, & !< Soil moisture (m3/m3), vertical surface elements
383	m_soil_v_p !< Prog. soil moisture (m3/m3), vertical surface elements
384
385	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::&
386	t_soil_v_1, & !<
387	t_soil_v_2, & !<
388	m_soil_v_1, & !<
389	m_soil_v_2 !<
390
391	TYPE(surf_type_lsm), POINTER :: t_surface_h, & !< surface temperature (K), horizontal surface elements
392	t_surface_h_p, & !< progn. surface temperature (K), horizontal surface elements
393	m_liq_h, & !< liquid water reservoir (m), horizontal surface elements
394	m_liq_h_p !< progn. liquid water reservoir (m), horizontal surface elements
395
396	TYPE(surf_type_lsm), TARGET :: t_surface_h_1, & !<
397	t_surface_h_2, & !<
398	m_liq_h_1, & !<
399	m_liq_h_2 !<
400
401	TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
402	t_surface_v, & !< surface temperature (K), vertical surface elements
403	t_surface_v_p, & !< progn. surface temperature (K), vertical surface elements
404	m_liq_v, & !< liquid water reservoir (m), vertical surface elements
405	m_liq_v_p !< progn. liquid water reservoir (m), vertical surface elements
406
407	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: &
408	t_surface_v_1, & !<
409	t_surface_v_2, & !<
410	m_liq_v_1, & !<
411	m_liq_v_2 !<
412
413	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
414
415	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: t_soil_av, & !< Average of t_soil
416	m_soil_av !< Average of m_soil
417
418	TYPE(surf_type_lsm), TARGET :: tm_liq_h_m !< liquid water reservoir tendency (m), horizontal surface elements
419	TYPE(surf_type_lsm), TARGET :: tt_surface_h_m !< surface temperature tendency (K), horizontal surface elements
420	TYPE(surf_type_lsm), TARGET :: tt_soil_h_m !< t_soil storage array, horizontal surface elements
421	TYPE(surf_type_lsm), TARGET :: tm_soil_h_m !< m_soil storage array, horizontal surface elements
422
423	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: tm_liq_v_m !< liquid water reservoir tendency (m), vertical surface elements
424	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: tt_surface_v_m !< surface temperature tendency (K), vertical surface elements
425	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: tt_soil_v_m !< t_soil storage array, vertical surface elements
426	TYPE(surf_type_lsm), DIMENSION(0:3), TARGET :: tm_soil_v_m !< m_soil storage array, vertical surface elements
427
428	!
429	!-- Energy balance variables
430	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
431	c_liq_av, & !< average of c_liq
432	c_soil_av, & !< average of c_soil
433	c_veg_av, & !< average of c_veg
434	lai_av, & !< average of lai
435	qsws_liq_av, & !< average of qsws_liq
436	qsws_soil_av, & !< average of qsws_soil
437	qsws_veg_av, & !< average of qsws_veg
438	r_s_av !< average of r_s
439
440	!
441	!-- Predefined Land surface classes (vegetation_type)
442	CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
443	'user defined ', & ! 0
444	'bare soil ', & ! 1
445	'crops, mixed farming ', & ! 2
446	'short grass ', & ! 3
447	'evergreen needleleaf trees', & ! 4
448	'deciduous needleleaf trees', & ! 5
449	'evergreen broadleaf trees ', & ! 6
450	'deciduous broadleaf trees ', & ! 7
451	'tall grass ', & ! 8
452	'desert ', & ! 9
453	'tundra ', & ! 10
454	'irrigated crops ', & ! 11
455	'semidesert ', & ! 12
456	'ice caps and glaciers ', & ! 13
457	'bogs and marshes ', & ! 14
458	'evergreen shrubs ', & ! 15
459	'deciduous shrubs ', & ! 16
460	'mixed forest/woodland ', & ! 17
461	'interrupted forest ' & ! 18
462	/)
463
464	!
465	!-- Soil model classes (soil_type)
466	CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
467	'user defined', & ! 0
468	'coarse ', & ! 1
469	'medium ', & ! 2
470	'medium-fine ', & ! 3
471	'fine ', & ! 4
472	'very fine ', & ! 5
473	'organic ' & ! 6
474	/)
475
476	!
477	!-- Pavement classes
478	CHARACTER(29), DIMENSION(0:15), PARAMETER :: pavement_type_name = (/ &
479	'user defined ', & ! 0
480	'asphalt/concrete mix ', & ! 1
481	'asphalt (asphalt concrete) ', & ! 2
482	'concrete (Portland concrete) ', & ! 3
483	'sett ', & ! 4
484	'paving stones ', & ! 5
485	'cobblestone ', & ! 6
486	'metal ', & ! 7
487	'wood ', & ! 8
488	'gravel ', & ! 9
489	'fine gravel ', & ! 10
490	'pebblestone ', & ! 11
491	'woodchips ', & ! 12
492	'tartan (sports) ', & ! 13
493	'artifical turf (sports) ', & ! 14
494	'clay (sports) ' & ! 15
495	/)
496
497	!
498	!-- Water classes
499	CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
500	'user defined', & ! 0
501	'lake ', & ! 1
502	'river ', & ! 2
503	'ocean ', & ! 3
504	'pond ', & ! 4
505	'fountain ' & ! 5
506	/)
507
508	!
509	!-- Land surface parameters according to the respective classes (vegetation_type)
510	INTEGER(iwp) :: ind_v_rc_min = 0 !< index for r_canopy_min in vegetation_pars
511	INTEGER(iwp) :: ind_v_rc_lai = 1 !< index for LAI in vegetation_pars
512	INTEGER(iwp) :: ind_v_c_veg = 2 !< index for c_veg in vegetation_pars
513	INTEGER(iwp) :: ind_v_gd = 3 !< index for g_d in vegetation_pars
514	INTEGER(iwp) :: ind_v_z0 = 4 !< index for z0 in vegetation_pars
515	INTEGER(iwp) :: ind_v_z0qh = 5 !< index for z0h / z0q in vegetation_pars
516	INTEGER(iwp) :: ind_v_lambda_s = 6 !< index for lambda_s_s in vegetation_pars
517	INTEGER(iwp) :: ind_v_lambda_u = 7 !< index for lambda_s_u in vegetation_pars
518	INTEGER(iwp) :: ind_v_f_sw_in = 8 !< index for f_sw_in in vegetation_pars
519	INTEGER(iwp) :: ind_v_c_surf = 9 !< index for c_surface in vegetation_pars
520	INTEGER(iwp) :: ind_v_at = 10 !< index for albedo_type in vegetation_pars
521	INTEGER(iwp) :: ind_v_emis = 11 !< index for emissivity in vegetation_pars
522
523	INTEGER(iwp) :: ind_w_temp = 0 !< index for temperature in water_pars
524	INTEGER(iwp) :: ind_w_z0 = 1 !< index for z0 in water_pars
525	INTEGER(iwp) :: ind_w_z0h = 2 !< index for z0h in water_pars
526	INTEGER(iwp) :: ind_w_lambda_s = 3 !< index for lambda_s_s in water_pars
527	INTEGER(iwp) :: ind_w_lambda_u = 4 !< index for lambda_s_u in water_pars
528	INTEGER(iwp) :: ind_w_at = 5 !< index for albedo type in water_pars
529	INTEGER(iwp) :: ind_w_emis = 6 !< index for emissivity in water_pars
530
531	INTEGER(iwp) :: ind_p_z0 = 0 !< index for z0 in pavement_pars
532	INTEGER(iwp) :: ind_p_z0h = 1 !< index for z0h in pavement_pars
533	INTEGER(iwp) :: ind_p_at = 2 !< index for albedo type in pavement_pars
534	INTEGER(iwp) :: ind_p_emis = 3 !< index for emissivity in pavement_pars
535	INTEGER(iwp) :: ind_p_lambda_h = 0 !< index for lambda_h in pavement_subsurface_pars
536	INTEGER(iwp) :: ind_p_rho_c = 1 !< index for rho_c in pavement_pars
537	!
538	!-- Land surface parameters
539	!-- r_canopy_min, lai, c_veg, g_d z0, z0h, lambda_s_s, lambda_s_u, f_sw_in, c_surface, albedo_type, emissivity
540	REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
541	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.005_wp, 0.5E-4_wp, 0.0_wp, 0.0_wp, 0.00_wp, 0.00_wp, 17.0_wp, 0.94_wp, & ! 1
542	180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp, 0.10_wp, 0.001_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 2.0_wp, 0.95_wp, & ! 2
543	110.0_wp, 2.00_wp, 1.00_wp, 0.00_wp, 0.03_wp, 0.3E-4_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 5.0_wp, 0.95_wp, & ! 3
544	500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 6.0_wp, 0.97_wp, & ! 4
545	500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 8.0_wp, 0.97_wp, & ! 5
546	175.0_wp, 5.00_wp, 1.00_wp, 0.03_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 9.0_wp, 0.97_wp, & ! 6
547	240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 7.0_wp, 0.97_wp, & ! 7
548	100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, 0.47_wp, 0.47E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 10.0_wp, 0.97_wp, & ! 8
549	250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, 0.013_wp, 0.013E-2_wp, 15.0_wp, 15.0_wp, 0.00_wp, 0.00_wp, 11.0_wp, 0.94_wp, & ! 9
550	80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, 0.034_wp, 0.034E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 13.0_wp, 0.97_wp, & ! 10
551	180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp, 0.5_wp, 0.50E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 2.0_wp, 0.97_wp, & ! 11
552	150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, 0.17_wp, 0.17E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 12
553	0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp, 1.3E-4_wp, 58.0_wp, 58.0_wp, 0.00_wp, 0.00_wp, 14.0_wp, 0.97_wp, & ! 13
554	240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, 0.83_wp, 0.83E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 3.0_wp, 0.97_wp, & ! 14
555	225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, 0.10_wp, 0.10E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 4.0_wp, 0.97_wp, & ! 15
556	225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, 0.25_wp, 0.25E-2_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 5.0_wp, 0.97_wp, & ! 16
557	250.0_wp, 5.00_wp, 1.00_wp, 0.03_wp, 2.00_wp, 2.00_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 10.0_wp, 0.97_wp, & ! 17
558	175.0_wp, 2.50_wp, 1.00_wp, 0.03_wp, 1.10_wp, 1.10_wp, 20.0_wp, 15.0_wp, 0.03_wp, 0.00_wp, 7.0_wp, 0.97_wp & ! 18
559	/), (/ 12, 18 /) )
560
561
562	!
563	!-- Root distribution for default soil layer configuration (sum = 1)
564	!-- level 1 - level 4 according to zs_ref
565	REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ &
566	1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 1
567	0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 2
568	0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 3
569	0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & ! 4
570	0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & ! 5
571	0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & ! 6
572	0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & ! 7
573	0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & ! 8
574	1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 9
575	0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 10
576	0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 11
577	0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 12
578	0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 13
579	0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 14
580	0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 15
581	0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16
582	0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 17
583	0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 18
584	/), (/ 4, 18 /) )
585
586	!
587	!-- Soil parameters according to the following porosity classes (soil_type)
588
589	!
590	!-- Soil parameters alpha_vg, l_vg, n_vg, gamma_w_sat, m_sat, m_fc, m_wilt, m_res
591	REAL(wp), DIMENSION(0:7,1:6), PARAMETER :: soil_pars = RESHAPE( (/ &
592	3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp,& ! 1
593	3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp,& ! 2
594	0.83_wp, -0.588_wp, 1.25_wp, 0.26E-6_wp, 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp,& ! 3
595	3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp,& ! 4
596	2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp,& ! 5
597	1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp, 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6
598	/), (/ 8, 6 /) )
599
600
601	!
602	!-- TO BE FILLED
603	!-- Pavement parameters z0, z0h, albedo_type, emissivity
604	REAL(wp), DIMENSION(0:3,1:15), PARAMETER :: pavement_pars = RESHAPE( (/ &
605	5.0E-2_wp, 5.0E-4_wp, 18.0_wp, 0.97_wp, & ! 1
606	5.0E-2_wp, 5.0E-4_wp, 19.0_wp, 0.94_wp, & ! 2
607	1.0E-2_wp, 1.0E-4_wp, 20.0_wp, 0.98_wp, & ! 3
608	1.0E-2_wp, 1.0E-4_wp, 21.0_wp, 0.93_wp, & ! 4
609	1.0E-2_wp, 1.0E-4_wp, 22.0_wp, 0.97_wp, & ! 5
610	1.0E-2_wp, 1.0E-4_wp, 23.0_wp, 0.97_wp, & ! 6
611	1.0E-2_wp, 1.0E-4_wp, 24.0_wp, 0.97_wp, & ! 7
612	1.0E-2_wp, 1.0E-4_wp, 25.0_wp, 0.94_wp, & ! 8
613	1.0E-2_wp, 1.0E-4_wp, 26.0_wp, 0.98_wp, & ! 9
614	1.0E-2_wp, 1.0E-4_wp, 27.0_wp, 0.93_wp, & ! 10
615	1.0E-2_wp, 1.0E-4_wp, 28.0_wp, 0.97_wp, & ! 11
616	1.0E-2_wp, 1.0E-4_wp, 29.0_wp, 0.97_wp, & ! 12
617	1.0E-2_wp, 1.0E-4_wp, 30.0_wp, 0.97_wp, & ! 13
618	1.0E-2_wp, 1.0E-4_wp, 31.0_wp, 0.94_wp, & ! 14
619	1.0E-2_wp, 1.0E-4_wp, 32.0_wp, 0.98_wp & ! 15
620	/), (/ 4, 15 /) )
621	!
622	!-- Pavement subsurface parameters part 1: thermal conductivity (W/m/K)
623	!-- 0.0-0.01, 0.01-0.03, 0.03-0.07, 0.07-0.15, 0.15-0.30, 0.30-0.50, 0.50-1.25, 1.25-3.00
624	REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_1 = RESHAPE( (/ &
625	0.75_wp, 0.75_wp, 0.75_wp, 0.75_wp, 0.75_wp, 0.75_wp, 9999999.9_wp, 9999999.9_wp, & ! 1
626	0.75_wp, 0.75_wp, 0.75_wp, 0.75_wp, 0.75_wp, 0.75_wp, 9999999.9_wp, 9999999.9_wp, & ! 2
627	0.89_wp, 0.89_wp, 0.89_wp, 0.89_wp, 0.89_wp, 0.89_wp, 9999999.9_wp, 9999999.9_wp, & ! 3
628	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 4
629	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 5
630	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 6
631	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 7
632	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 8
633	0.70_wp, 0.70_wp, 0.70_wp, 0.70_wp, 0.70_wp, 0.70_wp, 9999999.9_wp, 9999999.9_wp, & ! 9
634	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 10
635	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 11
636	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 12
637	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 13
638	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 14
639	1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 1.00_wp, 9999999.9_wp, 9999999.9_wp & ! 15
640	/), (/ 8, 15 /) )
641
642	!
643	!-- Pavement subsurface parameters part 2: volumetric heat capacity (J/m3/K)
644	!-- 0.0-0.01, 0.01-0.03, 0.03-0.07, 0.07-0.15, 0.15-0.30, 0.30-0.50, 0.50-1.25, 1.25-3.00
645	REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_2 = RESHAPE( (/ &
646	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 1
647	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 2
648	1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 3
649	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 4
650	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 5
651	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 6
652	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 7
653	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 8
654	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 9
655	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 10
656	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 11
657	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 12
658	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 13
659	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 14
660	1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp & ! 15
661	/), (/ 8, 15 /) )
662
663	!
664	!-- TO BE FILLED
665	!-- Water parameters temperature, z0, z0h, albedo_type, emissivity,
666	REAL(wp), DIMENSION(0:6,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
667	283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 1
668	283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 2
669	283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 3
670	283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 4
671	283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp & ! 5
672	/), (/ 7, 5 /) )
673
674	SAVE
675
676
677	PRIVATE
678
679
680	!
681	!-- Public functions
682	PUBLIC lsm_boundary_condition, lsm_check_data_output, &
683	lsm_check_data_output_pr, &
684	lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,&
685	lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance, &
686	lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model, &
687	lsm_swap_timelevel, lsm_rrd_local, lsm_wrd_local
688	! !vegetat
689	!-- Public parameters, constants and initial values
690	PUBLIC aero_resist_kray, skip_time_do_lsm
691
692	!
693	!-- Public grid variables
694	PUBLIC nzb_soil, nzs, nzt_soil, zs
695
696	!
697	!-- Public prognostic variables
698	PUBLIC m_soil_h, t_soil_h
699
700	INTERFACE lsm_boundary_condition
701	MODULE PROCEDURE lsm_boundary_condition
702	END INTERFACE lsm_boundary_condition
703
704	INTERFACE lsm_check_data_output
705	MODULE PROCEDURE lsm_check_data_output
706	END INTERFACE lsm_check_data_output
707
708	INTERFACE lsm_check_data_output_pr
709	MODULE PROCEDURE lsm_check_data_output_pr
710	END INTERFACE lsm_check_data_output_pr
711
712	INTERFACE lsm_check_parameters
713	MODULE PROCEDURE lsm_check_parameters
714	END INTERFACE lsm_check_parameters
715
716	INTERFACE lsm_3d_data_averaging
717	MODULE PROCEDURE lsm_3d_data_averaging
718	END INTERFACE lsm_3d_data_averaging
719
720	INTERFACE lsm_data_output_2d
721	MODULE PROCEDURE lsm_data_output_2d
722	END INTERFACE lsm_data_output_2d
723
724	INTERFACE lsm_data_output_3d
725	MODULE PROCEDURE lsm_data_output_3d
726	END INTERFACE lsm_data_output_3d
727
728	INTERFACE lsm_define_netcdf_grid
729	MODULE PROCEDURE lsm_define_netcdf_grid
730	END INTERFACE lsm_define_netcdf_grid
731
732	INTERFACE lsm_energy_balance
733	MODULE PROCEDURE lsm_energy_balance
734	END INTERFACE lsm_energy_balance
735
736	INTERFACE lsm_header
737	MODULE PROCEDURE lsm_header
738	END INTERFACE lsm_header
739
740	INTERFACE lsm_init
741	MODULE PROCEDURE lsm_init
742	END INTERFACE lsm_init
743
744	INTERFACE lsm_init_arrays
745	MODULE PROCEDURE lsm_init_arrays
746	END INTERFACE lsm_init_arrays
747
748	INTERFACE lsm_parin
749	MODULE PROCEDURE lsm_parin
750	END INTERFACE lsm_parin
751
752	INTERFACE lsm_soil_model
753	MODULE PROCEDURE lsm_soil_model
754	END INTERFACE lsm_soil_model
755
756	INTERFACE lsm_swap_timelevel
757	MODULE PROCEDURE lsm_swap_timelevel
758	END INTERFACE lsm_swap_timelevel
759
760	INTERFACE lsm_rrd_local
761	MODULE PROCEDURE lsm_rrd_local
762	END INTERFACE lsm_rrd_local
763
764	INTERFACE lsm_wrd_local
765	MODULE PROCEDURE lsm_wrd_local
766	END INTERFACE lsm_wrd_local
767
768	CONTAINS
769
770
771	!------------------------------------------------------------------------------!
772	! Description:
773	! ------------
774	!> Set internal Neumann boundary condition at outer soil grid points
775	!> for temperature and humidity.
776	!------------------------------------------------------------------------------!
777	SUBROUTINE lsm_boundary_condition
778
779	IMPLICIT NONE
780
781	INTEGER(iwp) :: i !< grid index x-direction
782	INTEGER(iwp) :: ioff !< offset index x-direction indicating location of soil grid point
783	INTEGER(iwp) :: j !< grid index y-direction
784	INTEGER(iwp) :: joff !< offset index x-direction indicating location of soil grid point
785	INTEGER(iwp) :: k !< grid index z-direction
786	INTEGER(iwp) :: koff !< offset index x-direction indicating location of soil grid point
787	INTEGER(iwp) :: l !< running index surface-orientation
788	INTEGER(iwp) :: m !< running index surface elements
789
790	koff = surf_lsm_h%koff
791	DO m = 1, surf_lsm_h%ns
792	i = surf_lsm_h%i(m)
793	j = surf_lsm_h%j(m)
794	k = surf_lsm_h%k(m)
795	pt(k+koff,j,i) = pt(k,j,i)
796	ENDDO
797
798	DO l = 0, 3
799	ioff = surf_lsm_v(l)%ioff
800	joff = surf_lsm_v(l)%joff
801	DO m = 1, surf_lsm_v(l)%ns
802	i = surf_lsm_v(l)%i(m)
803	j = surf_lsm_v(l)%j(m)
804	k = surf_lsm_v(l)%k(m)
805	pt(k,j+joff,i+ioff) = pt(k,j,i)
806	ENDDO
807	ENDDO
808	!
809	!-- In case of humidity, set boundary conditions also for q and vpt.
810	IF ( humidity ) THEN
811	koff = surf_lsm_h%koff
812	DO m = 1, surf_lsm_h%ns
813	i = surf_lsm_h%i(m)
814	j = surf_lsm_h%j(m)
815	k = surf_lsm_h%k(m)
816	q(k+koff,j,i) = q(k,j,i)
817	vpt(k+koff,j,i) = vpt(k,j,i)
818	ENDDO
819
820	DO l = 0, 3
821	ioff = surf_lsm_v(l)%ioff
822	joff = surf_lsm_v(l)%joff
823	DO m = 1, surf_lsm_v(l)%ns
824	i = surf_lsm_v(l)%i(m)
825	j = surf_lsm_v(l)%j(m)
826	k = surf_lsm_v(l)%k(m)
827	q(k,j+joff,i+ioff) = q(k,j,i)
828	vpt(k,j+joff,i+ioff) = vpt(k,j,i)
829	ENDDO
830	ENDDO
831	ENDIF
832
833	END SUBROUTINE lsm_boundary_condition
834
835	!------------------------------------------------------------------------------!
836	! Description:
837	! ------------
838	!> Check data output for land surface model
839	!------------------------------------------------------------------------------!
840	SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
841
842
843	USE control_parameters, &
844	ONLY: data_output, message_string
845
846	IMPLICIT NONE
847
848	CHARACTER (LEN=*) :: unit !<
849	CHARACTER (LEN=*) :: var !<
850
851	INTEGER(iwp) :: i
852	INTEGER(iwp) :: ilen
853	INTEGER(iwp) :: k
854
855	SELECT CASE ( TRIM( var ) )
856
857	CASE ( 'm_soil' )
858	IF ( .NOT. land_surface ) THEN
859	message_string = 'output of "' // TRIM( var ) // '" requi' // &
860	'res land_surface = .TRUE.'
861	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
862	ENDIF
863	unit = 'm3/m3'
864
865	CASE ( 't_soil' )
866	IF ( .NOT. land_surface ) THEN
867	message_string = 'output of "' // TRIM( var ) // '" requi' // &
868	'res land_surface = .TRUE.'
869	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
870	ENDIF
871	unit = 'K'
872
873	CASE ( 'lai', 'c_liq', 'c_soil', 'c_veg', 'm_liq*', &
874	'qsws_liq', 'qsws_soil', 'qsws_veg', 'r_s' )
875	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
876	message_string = 'illegal value for data_output: "' // &
877	TRIM( var ) // '" & only 2d-horizontal ' // &
878	'cross sections are allowed for this value'
879	CALL message( 'lsm_check_data_output', 'PA0111', 1, 2, 0, 6, 0 )
880	ENDIF
881	IF ( TRIM( var ) == 'lai*' .AND. .NOT. land_surface ) THEN
882	message_string = 'output of "' // TRIM( var ) // '" requi' // &
883	'res land_surface = .TRUE.'
884	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
885	ENDIF
886	IF ( TRIM( var ) == 'c_liq*' .AND. .NOT. land_surface ) THEN
887	message_string = 'output of "' // TRIM( var ) // '" requi' // &
888	'res land_surface = .TRUE.'
889	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
890	ENDIF
891	IF ( TRIM( var ) == 'c_soil*' .AND. .NOT. land_surface ) THEN
892	message_string = 'output of "' // TRIM( var ) // '" requi' // &
893	'res land_surface = .TRUE.'
894	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
895	ENDIF
896	IF ( TRIM( var ) == 'c_veg*' .AND. .NOT. land_surface ) THEN
897	message_string = 'output of "' // TRIM( var ) // '" requi' // &
898	'res land_surface = .TRUE.'
899	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
900	ENDIF
901	IF ( TRIM( var ) == 'm_liq*' .AND. .NOT. land_surface ) THEN
902	message_string = 'output of "' // TRIM( var ) // '" requi' // &
903	'res land_surface = .TRUE.'
904	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
905	ENDIF
906	IF ( TRIM( var ) == 'qsws_liq*' .AND. .NOT. land_surface ) &
907	THEN
908	message_string = 'output of "' // TRIM( var ) // '" requi' // &
909	'res land_surface = .TRUE.'
910	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
911	ENDIF
912	IF ( TRIM( var ) == 'qsws_soil*' .AND. .NOT. land_surface ) &
913	THEN
914	message_string = 'output of "' // TRIM( var ) // '" requi' // &
915	'res land_surface = .TRUE.'
916	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
917	ENDIF
918	IF ( TRIM( var ) == 'qsws_veg*' .AND. .NOT. land_surface ) &
919	THEN
920	message_string = 'output of "' // TRIM( var ) // '" requi' // &
921	'res land_surface = .TRUE.'
922	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
923	ENDIF
924	IF ( TRIM( var ) == 'r_s*' .AND. .NOT. land_surface ) &
925	THEN
926	message_string = 'output of "' // TRIM( var ) // '" requi' // &
927	'res land_surface = .TRUE.'
928	CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
929	ENDIF
930
931	IF ( TRIM( var ) == 'lai*' ) unit = 'none'
932	IF ( TRIM( var ) == 'c_liq*' ) unit = 'none'
933	IF ( TRIM( var ) == 'c_soil*') unit = 'none'
934	IF ( TRIM( var ) == 'c_veg*' ) unit = 'none'
935	IF ( TRIM( var ) == 'm_liq*' ) unit = 'm'
936	IF ( TRIM( var ) == 'qsws_liq*' ) unit = 'W/m2'
937	IF ( TRIM( var ) == 'qsws_soil*' ) unit = 'W/m2'
938	IF ( TRIM( var ) == 'qsws_veg*' ) unit = 'W/m2'
939	IF ( TRIM( var ) == 'r_s*') unit = 's/m'
940
941	CASE DEFAULT
942	unit = 'illegal'
943
944	END SELECT
945
946
947	END SUBROUTINE lsm_check_data_output
948
949
950
951	!------------------------------------------------------------------------------!
952	! Description:
953	! ------------
954	!> Check data output of profiles for land surface model
955	!------------------------------------------------------------------------------!
956	SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
957
958	USE control_parameters, &
959	ONLY: data_output_pr, message_string
960
961	USE indices
962
963	USE profil_parameter
964
965	USE statistics
966
967	IMPLICIT NONE
968
969	CHARACTER (LEN=*) :: unit !<
970	CHARACTER (LEN=*) :: variable !<
971	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
972
973	INTEGER(iwp) :: var_count !<
974
975	SELECT CASE ( TRIM( variable ) )
976
977	CASE ( 't_soil', '#t_soil' )
978	IF ( .NOT. land_surface ) THEN
979	message_string = 'data_output_pr = ' // &
980	TRIM( data_output_pr(var_count) ) // ' is' // &
981	'not implemented for land_surface = .FALSE.'
982	CALL message( 'lsm_check_data_output_pr', 'PA0402', 1, 2, 0, 6, 0 )
983	ELSE
984	dopr_index(var_count) = 89
985	dopr_unit = 'K'
986	hom(0:nzs-1,2,89,:) = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
987	IF ( data_output_pr(var_count)(1:1) == '#' ) THEN
988	dopr_initial_index(var_count) = 90
989	hom(0:nzs-1,2,90,:) = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
990	data_output_pr(var_count) = data_output_pr(var_count)(2:)
991	ENDIF
992	unit = dopr_unit
993	ENDIF
994
995	CASE ( 'm_soil', '#m_soil' )
996	IF ( .NOT. land_surface ) THEN
997	message_string = 'data_output_pr = ' // &
998	TRIM( data_output_pr(var_count) ) // ' is' // &
999	' not implemented for land_surface = .FALSE.'
1000	CALL message( 'lsm_check_data_output_pr', 'PA0402', 1, 2, 0, 6, 0 )
1001	ELSE
1002	dopr_index(var_count) = 91
1003	dopr_unit = 'm3/m3'
1004	hom(0:nzs-1,2,91,:) = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1005	IF ( data_output_pr(var_count)(1:1) == '#' ) THEN
1006	dopr_initial_index(var_count) = 92
1007	hom(0:nzs-1,2,92,:) = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1008	data_output_pr(var_count) = data_output_pr(var_count)(2:)
1009	ENDIF
1010	unit = dopr_unit
1011	ENDIF
1012
1013
1014	CASE DEFAULT
1015	unit = 'illegal'
1016
1017	END SELECT
1018
1019
1020	END SUBROUTINE lsm_check_data_output_pr
1021
1022
1023	!------------------------------------------------------------------------------!
1024	! Description:
1025	! ------------
1026	!> Check parameters routine for land surface model
1027	!------------------------------------------------------------------------------!
1028	SUBROUTINE lsm_check_parameters
1029
1030	USE control_parameters, &
1031	ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string
1032
1033
1034	IMPLICIT NONE
1035
1036	INTEGER(iwp) :: i !< running index, x-dimension
1037	INTEGER(iwp) :: j !< running index, y-dimension
1038	INTEGER(iwp) :: k !< running index, z-dimension
1039
1040	LOGICAL :: dynamic_soil_input_parent !< flag indicating the presence of a dynamic input file for the parent
1041
1042	!
1043	!-- Check for a valid setting of surface_type. The default value is 'netcdf'.
1044	!-- In that case, the surface types are read from NetCDF file
1045	IF ( TRIM( surface_type ) /= 'vegetation' .AND. &
1046	TRIM( surface_type ) /= 'pavement' .AND. &
1047	TRIM( surface_type ) /= 'water' .AND. &
1048	TRIM( surface_type ) /= 'netcdf' ) THEN
1049	message_string = 'unknown surface type: surface_type = "' // &
1050	TRIM( surface_type ) // '"'
1051	CALL message( 'lsm_check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
1052	ENDIF
1053
1054	!
1055	!-- Dirichlet boundary conditions are required as the surface fluxes are
1056	!-- calculated from the temperature/humidity gradients in the land surface
1057	!-- model
1058	IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN
1059	message_string = 'lsm requires setting of'// &
1060	'bc_pt_b = "dirichlet" and '// &
1061	'bc_q_b = "dirichlet"'
1062	CALL message( 'lsm_check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
1063	ENDIF
1064
1065	IF ( .NOT. constant_flux_layer ) THEN
1066	message_string = 'lsm requires '// &
1067	'constant_flux_layer = .T.'
1068	CALL message( 'lsm_check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1069	ENDIF
1070
1071	IF ( .NOT. radiation ) THEN
1072	message_string = 'lsm requires '// &
1073	'the radiation model to be switched on'
1074	CALL message( 'lsm_check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1075	ENDIF
1076	!
1077	!-- Check if soil types are set within a valid range.
1078	IF ( TRIM( surface_type ) == 'vegetation' .OR. &
1079	TRIM( surface_type ) == 'pavement' ) THEN
1080	IF ( soil_type < LBOUND( soil_pars, 2 ) .AND. &
1081	soil_type > UBOUND( soil_pars, 2 ) ) THEN
1082	WRITE( message_string, * ) 'soil_type = ', soil_type, ' is out ' // &
1083	'of the valid range'
1084	CALL message( 'lsm_check_parameters', 'PA0452', 2, 2, 0, 6, 0 )
1085	ENDIF
1086	ENDIF
1087	IF ( TRIM( surface_type ) == 'netcdf' ) THEN
1088	IF ( soil_type_f%from_file ) THEN
1089	DO i = nxl, nxr
1090	DO j = nys, nyn
1091	IF ( soil_type_f%var_2d(j,i) /= soil_type_f%fill .AND. &
1092	( soil_type_f%var_2d(j,i) < LBOUND( soil_pars, 2 ) .OR. &
1093	soil_type_f%var_2d(j,i) > UBOUND( soil_pars, 2 ) ) ) THEN
1094	WRITE( message_string, * ) 'soil_type = is out of ' // &
1095	'the valid range at (j,i) = ', j, i
1096	CALL message( 'lsm_check_parameters', 'PA0452', &
1097	2, 2, myid, 6, 0 )
1098	ENDIF
1099	ENDDO
1100	ENDDO
1101	ENDIF
1102	ENDIF
1103	!
1104	!-- Check if vegetation types are set within a valid range.
1105	IF ( TRIM( surface_type ) == 'vegetation' ) THEN
1106	IF ( vegetation_type < LBOUND( vegetation_pars, 2 ) .AND. &
1107	vegetation_type > UBOUND( vegetation_pars, 2 ) ) THEN
1108	WRITE( message_string, * ) 'vegetation_type = ', vegetation_type, &
1109	' is out of the valid range'
1110	CALL message( 'lsm_check_parameters', 'PA0526', 2, 2, 0, 6, 0 )
1111	ENDIF
1112	ENDIF
1113	IF ( TRIM( surface_type ) == 'netcdf' ) THEN
1114	IF ( vegetation_type_f%from_file ) THEN
1115	DO i = nxl, nxr
1116	DO j = nys, nyn
1117	IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill .AND.&
1118	( vegetation_type_f%var(j,i) < LBOUND( vegetation_pars, 2 ) .OR.&
1119	vegetation_type_f%var(j,i) > UBOUND( vegetation_pars, 2 ) ) ) &
1120	THEN
1121	WRITE( message_string, * ) 'vegetation_type = is out of ' //&
1122	'the valid range at (j,i) = ', j, i
1123	CALL message( 'lsm_check_parameters', 'PA0526', &
1124	2, 2, myid, 6, 0 )
1125	ENDIF
1126	ENDDO
1127	ENDDO
1128	ENDIF
1129	ENDIF
1130	!
1131	!-- Check if pavement types are set within a valid range.
1132	IF ( TRIM( surface_type ) == 'pavement' ) THEN
1133	IF ( pavement_type < LBOUND( pavement_pars, 2 ) .AND. &
1134	pavement_type > UBOUND( pavement_pars, 2 ) ) THEN
1135	WRITE( message_string, * ) 'pavement_type = ', pavement_type, &
1136	' is out of the valid range'
1137	CALL message( 'lsm_check_parameters', 'PA0527', 2, 2, 0, 6, 0 )
1138	ENDIF
1139	ENDIF
1140	IF ( TRIM( surface_type ) == 'netcdf' ) THEN
1141	IF ( pavement_type_f%from_file ) THEN
1142	DO i = nxl, nxr
1143	DO j = nys, nyn
1144	IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill .AND. &
1145	( pavement_type_f%var(j,i) < LBOUND( pavement_pars, 2 ) .OR. &
1146	pavement_type_f%var(j,i) > UBOUND( pavement_pars, 2 ) ) ) THEN
1147	WRITE( message_string, * ) 'pavement_type = is out of ' // &
1148	'the valid range at (j,i) = ', j, i
1149	CALL message( 'lsm_check_parameters', 'PA0527', &
1150	2, 2, myid, 6, 0 )
1151	ENDIF
1152	ENDDO
1153	ENDDO
1154	ENDIF
1155	ENDIF
1156	!
1157	!-- Check if water types are set within a valid range.
1158	IF ( TRIM( surface_type ) == 'water' ) THEN
1159	IF ( water_type < LBOUND( water_pars, 2 ) .AND. &
1160	water_type > UBOUND( water_pars, 2 ) ) THEN
1161	WRITE( message_string, * ) 'water_type = ', water_type, &
1162	' is out of the valid range'
1163	CALL message( 'lsm_check_parameters', 'PA0528', 2, 2, 0, 6, 0 )
1164	ENDIF
1165	ENDIF
1166	IF ( TRIM( surface_type ) == 'netcdf' ) THEN
1167	IF ( water_type_f%from_file ) THEN
1168	DO i = nxl, nxr
1169	DO j = nys, nyn
1170	IF ( water_type_f%var(j,i) /= water_type_f%fill .AND. &
1171	( water_type_f%var(j,i) < LBOUND( water_pars, 2 ) .OR. &
1172	water_type_f%var(j,i) > UBOUND( water_pars, 2 ) ) ) THEN
1173	WRITE( message_string, * ) 'water_type = is out of ' // &
1174	'the valid range at (j,i) = ', j, i
1175	CALL message( 'lsm_check_parameters', 'PA0528', &
1176	2, 2, myid, 6, 0 )
1177	ENDIF
1178	ENDDO
1179	ENDDO
1180	ENDIF
1181	ENDIF
1182	!
1183	!-- Check further settings for consistency.
1184	IF ( TRIM( surface_type ) == 'vegetation' ) THEN
1185
1186	IF ( vegetation_type == 0 ) THEN
1187	IF ( min_canopy_resistance == 9999999.9_wp ) THEN
1188	message_string = 'vegetation_type = 0 (user defined)'// &
1189	'requires setting of min_canopy_resistance'// &
1190	'/= 9999999.9'
1191	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1192	ENDIF
1193
1194	IF ( leaf_area_index == 9999999.9_wp ) THEN
1195	message_string = 'vegetation_type = 0 (user_defined)'// &
1196	'requires setting of leaf_area_index'// &
1197	'/= 9999999.9'
1198	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1199	ENDIF
1200
1201	IF ( vegetation_coverage == 9999999.9_wp ) THEN
1202	message_string = 'vegetation_type = 0 (user_defined)'// &
1203	'requires setting of vegetation_coverage'// &
1204	'/= 9999999.9'
1205	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1206	ENDIF
1207
1208	IF ( canopy_resistance_coefficient == 9999999.9_wp) THEN
1209	message_string = 'vegetation_type = 0 (user_defined)'// &
1210	'requires setting of'// &
1211	'canopy_resistance_coefficient /= 9999999.9'
1212	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1213	ENDIF
1214
1215	IF ( lambda_surface_stable == 9999999.9_wp ) THEN
1216	message_string = 'vegetation_type = 0 (user_defined)'// &
1217	'requires setting of lambda_surface_stable'// &
1218	'/= 9999999.9'
1219	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1220	ENDIF
1221
1222	IF ( lambda_surface_unstable == 9999999.9_wp ) THEN
1223	message_string = 'vegetation_type = 0 (user_defined)'// &
1224	'requires setting of lambda_surface_unstable'// &
1225	'/= 9999999.9'
1226	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1227	ENDIF
1228
1229	IF ( f_shortwave_incoming == 9999999.9_wp ) THEN
1230	message_string = 'vegetation_type = 0 (user_defined)'// &
1231	'requires setting of f_shortwave_incoming'// &
1232	'/= 9999999.9'
1233	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1234	ENDIF
1235
1236	IF ( z0_vegetation == 9999999.9_wp ) THEN
1237	message_string = 'vegetation_type = 0 (user_defined)'// &
1238	'requires setting of z0_vegetation'// &
1239	'/= 9999999.9'
1240	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1241	ENDIF
1242
1243	IF ( z0h_vegetation == 9999999.9_wp ) THEN
1244	message_string = 'vegetation_type = 0 (user_defined)'// &
1245	'requires setting of z0h_vegetation'// &
1246	'/= 9999999.9'
1247	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1248	ENDIF
1249	ENDIF
1250
1251	IF ( vegetation_type == 1 ) THEN
1252	IF ( vegetation_coverage /= 9999999.9_wp .AND. vegetation_coverage &
1253	/= 0.0_wp ) THEN
1254	message_string = 'vegetation_type = 1 (bare soil)'// &
1255	' requires vegetation_coverage = 0'
1256	CALL message( 'lsm_check_parameters', 'PA0294', 1, 2, 0, 6, 0 )
1257	ENDIF
1258	ENDIF
1259
1260	ENDIF
1261
1262	IF ( TRIM( surface_type ) == 'water' ) THEN
1263
1264	IF ( water_type == 0 ) THEN
1265
1266	IF ( z0_water == 9999999.9_wp ) THEN
1267	message_string = 'water_type = 0 (user_defined)'// &
1268	'requires setting of z0_water'// &
1269	'/= 9999999.9'
1270	CALL message( 'lsm_check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
1271	ENDIF
1272
1273	IF ( z0h_water == 9999999.9_wp ) THEN
1274	message_string = 'water_type = 0 (user_defined)'// &
1275	'requires setting of z0h_water'// &
1276	'/= 9999999.9'
1277	CALL message( 'lsm_check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
1278	ENDIF
1279
1280	IF ( water_temperature == 9999999.9_wp ) THEN
1281	message_string = 'water_type = 0 (user_defined)'// &
1282	'requires setting of water_temperature'// &
1283	'/= 9999999.9'
1284	CALL message( 'lsm_check_parameters', 'PA0379', 1, 2, 0, 6, 0 )
1285	ENDIF
1286
1287	ENDIF
1288
1289	ENDIF
1290
1291	IF ( TRIM( surface_type ) == 'pavement' ) THEN
1292
1293	IF ( ANY( dz_soil /= 9999999.9_wp ) .AND. pavement_type /= 0 ) THEN
1294	message_string = 'non-default setting of dz_soil '// &
1295	'does not allow to use pavement_type /= 0)'
1296	CALL message( 'lsm_check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
1297	ENDIF
1298
1299	IF ( pavement_type == 0 ) THEN
1300
1301	IF ( z0_pavement == 9999999.9_wp ) THEN
1302	message_string = 'pavement_type = 0 (user_defined)'// &
1303	'requires setting of z0_pavement'// &
1304	'/= 9999999.9'
1305	CALL message( 'lsm_check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
1306	ENDIF
1307
1308	IF ( z0h_pavement == 9999999.9_wp ) THEN
1309	message_string = 'pavement_type = 0 (user_defined)'// &
1310	'requires setting of z0h_pavement'// &
1311	'/= 9999999.9'
1312	CALL message( 'lsm_check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
1313	ENDIF
1314
1315	IF ( pavement_heat_conduct == 9999999.9_wp ) THEN
1316	message_string = 'pavement_type = 0 (user_defined)'// &
1317	'requires setting of pavement_heat_conduct'// &
1318	'/= 9999999.9'
1319	CALL message( 'lsm_check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
1320	ENDIF
1321
1322	IF ( pavement_heat_capacity == 9999999.9_wp ) THEN
1323	message_string = 'pavement_type = 0 (user_defined)'// &
1324	'requires setting of pavement_heat_capacity'// &
1325	'/= 9999999.9'
1326	CALL message( 'lsm_check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
1327	ENDIF
1328
1329	IF ( pavement_depth_level == 0 ) THEN
1330	message_string = 'pavement_type = 0 (user_defined)'// &
1331	'requires setting of pavement_depth_level'// &
1332	'/= 0'
1333	CALL message( 'lsm_check_parameters', 'PA0474', 1, 2, 0, 6, 0 )
1334	ENDIF
1335
1336	ENDIF
1337
1338	ENDIF
1339
1340	IF ( TRIM( surface_type ) == 'netcdf' ) THEN
1341	!
1342	!-- MS: Some problme here, after calling message everythings stucks at
1343	!-- MPI_FINALIZE call.
1344	IF ( ANY( pavement_type_f%var /= pavement_type_f%fill ) .AND. &
1345	ANY( dz_soil /= 9999999.9_wp ) ) THEN
1346	message_string = 'pavement-surfaces are not allowed in ' // &
1347	'combination with a non-default setting of dz_soil'
1348	CALL message( 'lsm_check_parameters', 'PA0316', 2, 2, 0, 6, 0 )
1349	ENDIF
1350	ENDIF
1351
1352	!
1353	!-- Temporary message as long as NetCDF input is not available
1354	IF ( TRIM( surface_type ) == 'netcdf' .AND. .NOT. input_pids_static ) &
1355	THEN
1356	message_string = 'surface_type = netcdf requires static input file.'
1357	CALL message( 'lsm_check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
1358	ENDIF
1359
1360	IF ( soil_type == 0 .AND. .NOT. input_pids_static ) THEN
1361
1362	IF ( alpha_vangenuchten == 9999999.9_wp ) THEN
1363	message_string = 'soil_type = 0 (user_defined)'// &
1364	'requires setting of alpha_vangenuchten'// &
1365	'/= 9999999.9'
1366	CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1367	ENDIF
1368
1369	IF ( l_vangenuchten == 9999999.9_wp ) THEN
1370	message_string = 'soil_type = 0 (user_defined)'// &
1371	'requires setting of l_vangenuchten'// &
1372	'/= 9999999.9'
1373	CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1374	ENDIF
1375
1376	IF ( n_vangenuchten == 9999999.9_wp ) THEN
1377	message_string = 'soil_type = 0 (user_defined)'// &
1378	'requires setting of n_vangenuchten'// &
1379	'/= 9999999.9'
1380	CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1381	ENDIF
1382
1383	IF ( hydraulic_conductivity == 9999999.9_wp ) THEN
1384	message_string = 'soil_type = 0 (user_defined)'// &
1385	'requires setting of hydraulic_conductivity'// &
1386	'/= 9999999.9'
1387	CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1388	ENDIF
1389
1390	IF ( saturation_moisture == 9999999.9_wp ) THEN
1391	message_string = 'soil_type = 0 (user_defined)'// &
1392	'requires setting of saturation_moisture'// &
1393	'/= 9999999.9'
1394	CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1395	ENDIF
1396
1397	IF ( field_capacity == 9999999.9_wp ) THEN
1398	message_string = 'soil_type = 0 (user_defined)'// &
1399	'requires setting of field_capacity'// &
1400	'/= 9999999.9'
1401	CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1402	ENDIF
1403
1404	IF ( wilting_point == 9999999.9_wp ) THEN
1405	message_string = 'soil_type = 0 (user_defined)'// &
1406	'requires setting of wilting_point'// &
1407	'/= 9999999.9'
1408	CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1409	ENDIF
1410
1411	IF ( residual_moisture == 9999999.9_wp ) THEN
1412	message_string = 'soil_type = 0 (user_defined)'// &
1413	'requires setting of residual_moisture'// &
1414	'/= 9999999.9'
1415	CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1416	ENDIF
1417
1418	ENDIF
1419
1420
1421	!!! these checks are not needed for water surfaces??
1422
1423	!
1424	!-- Determine number of soil layers to be used and check whether an appropriate
1425	!-- root fraction is prescribed
1426	nzb_soil = 0
1427	nzt_soil = -1
1428	IF ( ALL( dz_soil == 9999999.9_wp ) ) THEN
1429	nzt_soil = 7
1430	dz_soil(nzb_soil:nzt_soil) = dz_soil_default
1431	ELSE
1432	DO k = 0, 19
1433	IF ( dz_soil(k) /= 9999999.9_wp ) THEN
1434	nzt_soil = nzt_soil + 1
1435	ENDIF
1436	ENDDO
1437	ENDIF
1438	nzs = nzt_soil + 1
1439
1440	!
1441	!-- Check whether valid soil temperatures are prescribed. Only check this if
1442	!-- no dynamic soil is not initialized with dynamic input.
1443	!-- In a nested case, check whether there is a dynamic input file for the
1444	!-- child (input_pids_dynamic = .T.) or one for the parent (inquire without
1445	!-- coupling_char.
1446	INQUIRE( FILE = TRIM( input_file_dynamic ), &
1447	EXIST = dynamic_soil_input_parent )
1448
1449	IF ( .NOT. input_pids_dynamic .AND. .NOT. dynamic_soil_input_parent ) THEN
1450	IF ( COUNT( soil_temperature /= 9999999.9_wp ) /= nzs ) THEN
1451	WRITE( message_string, * ) &
1452	'number of soil layers (', nzs, ') does not',&
1453	' match to the number of layers specified', &
1454	' in soil_temperature (', COUNT( &
1455	soil_temperature /= 9999999.9_wp ), ')'
1456	CALL message( 'lsm_check_parameters', 'PA0471', 1, 2, 0, 6, 0 )
1457	ENDIF
1458
1459	IF ( deep_soil_temperature == 9999999.9_wp ) THEN
1460	message_string = 'deep_soil_temperature is not set but must be'// &
1461	'/= 9999999.9'
1462	CALL message( 'lsm_check_parameters', 'PA0472', 1, 2, 0, 6, 0 )
1463	ENDIF
1464	ENDIF
1465
1466	!
1467	!-- Check whether the sum of all root fractions equals one
1468	IF ( .NOT. vegetation_type_f%from_file ) THEN
1469	IF ( vegetation_type == 0 ) THEN
1470	IF ( SUM( root_fraction(nzb_soil:nzt_soil) ) /= 1.0_wp ) THEN
1471	message_string = 'vegetation_type = 0 (user_defined)'// &
1472	'requires setting of root_fraction'// &
1473	'/= 9999999.9 and SUM(root_fraction) = 1'
1474	CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1475	ENDIF
1476	ENDIF
1477	ENDIF
1478	!
1479	!-- Calculate grid spacings. Temperature and moisture are defined at
1480	!-- the center of the soil layers, whereas gradients/fluxes are
1481	!-- defined at the edges (_layer)
1482	!
1483	!-- Allocate global 1D arrays
1484	ALLOCATE ( ddz_soil_center(nzb_soil:nzt_soil) )
1485	ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
1486	ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
1487	ALLOCATE ( zs(nzb_soil:nzt_soil+1) )
1488
1489
1490	zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
1491	zs_layer(nzb_soil) = dz_soil(nzb_soil)
1492
1493	DO k = nzb_soil+1, nzt_soil
1494	zs_layer(k) = zs_layer(k-1) + dz_soil(k)
1495	zs(k) = (zs_layer(k) + zs_layer(k-1)) * 0.5_wp
1496	ENDDO
1497
1498	dz_soil(nzt_soil+1) = zs_layer(nzt_soil) + dz_soil(nzt_soil)
1499	zs(nzt_soil+1) = zs_layer(nzt_soil) + 0.5_wp * dz_soil(nzt_soil)
1500
1501	DO k = nzb_soil, nzt_soil-1
1502	dz_soil_center(k) = zs(k+1) - zs(k)
1503	IF ( dz_soil_center(k) <= 0.0_wp ) THEN
1504	message_string = 'invalid soil layer configuration found ' // &
1505	'(dz_soil_center(k) <= 0.0)'
1506	CALL message( 'lsm_rrd_local', 'PA0140', 1, 2, 0, 6, 0 )
1507	ENDIF
1508	ENDDO
1509
1510	dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
1511
1512	ddz_soil_center = 1.0_wp / dz_soil_center
1513	ddz_soil(nzb_soil:nzt_soil) = 1.0_wp / dz_soil(nzb_soil:nzt_soil)
1514
1515
1516
1517	END SUBROUTINE lsm_check_parameters
1518
1519	!------------------------------------------------------------------------------!
1520	! Description:
1521	! ------------
1522	!> Solver for the energy balance at the surface.
1523	!------------------------------------------------------------------------------!
1524	SUBROUTINE lsm_energy_balance( horizontal, l )
1525
1526	USE pegrid
1527	USE radiation_model_mod, ONLY: rad_lw_out
1528
1529	IMPLICIT NONE
1530
1531	INTEGER(iwp) :: i !< running index
1532	INTEGER(iwp) :: i_off !< offset to determine index of surface element, seen from atmospheric grid point, for x
1533	INTEGER(iwp) :: j !< running index
1534	INTEGER(iwp) :: j_off !< offset to determine index of surface element, seen from atmospheric grid point, for y
1535	INTEGER(iwp) :: k !< running index
1536	INTEGER(iwp) :: k_off !< offset to determine index of surface element, seen from atmospheric grid point, for z
1537	INTEGER(iwp) :: ks !< running index
1538	INTEGER(iwp) :: l !< surface-facing index
1539	INTEGER(iwp) :: m !< running index concerning wall elements
1540
1541	LOGICAL :: horizontal !< Flag indicating horizontal or vertical surfaces
1542
1543	REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
1544	f1, & !< resistance correction term 1
1545	f2, & !< resistance correction term 2
1546	f3, & !< resistance correction term 3
1547	m_min, & !< minimum soil moisture
1548	e, & !< water vapour pressure
1549	e_s, & !< water vapour saturation pressure
1550	e_s_dt, & !< derivate of e_s with respect to T
1551	tend, & !< tendency
1552	dq_s_dt, & !< derivate of q_s with respect to T
1553	coef_1, & !< coef. for prognostic equation
1554	coef_2, & !< coef. for prognostic equation
1555	f_qsws, & !< factor for qsws
1556	f_qsws_veg, & !< factor for qsws_veg
1557	f_qsws_soil, & !< factor for qsws_soil
1558	f_qsws_liq, & !< factor for qsws_liq
1559	f_shf, & !< factor for shf
1560	lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
1561	lambda_surface, & !< Current value of lambda_surface (W/m2/K)
1562	m_liq_max !< maxmimum value of the liq. water reservoir
1563
1564	TYPE(surf_type_lsm), POINTER :: surf_t_surface
1565	TYPE(surf_type_lsm), POINTER :: surf_t_surface_p
1566	TYPE(surf_type_lsm), POINTER :: surf_tt_surface_m
1567	TYPE(surf_type_lsm), POINTER :: surf_m_liq
1568	TYPE(surf_type_lsm), POINTER :: surf_m_liq_p
1569	TYPE(surf_type_lsm), POINTER :: surf_tm_liq_m
1570
1571	TYPE(surf_type_lsm), POINTER :: surf_m_soil
1572	TYPE(surf_type_lsm), POINTER :: surf_t_soil
1573
1574	TYPE(surf_type), POINTER :: surf !< surface-date type variable
1575
1576
1577	IF ( debug_output_timestep ) THEN
1578	WRITE( debug_string, * ) 'lsm_energy_balance', horizontal, l
1579	CALL debug_message( debug_string, 'start' )
1580	ENDIF
1581
1582	IF ( horizontal ) THEN
1583	surf => surf_lsm_h
1584
1585	surf_t_surface => t_surface_h
1586	surf_t_surface_p => t_surface_h_p
1587	surf_tt_surface_m => tt_surface_h_m
1588	surf_m_liq => m_liq_h
1589	surf_m_liq_p => m_liq_h_p
1590	surf_tm_liq_m => tm_liq_h_m
1591	surf_m_soil => m_soil_h
1592	surf_t_soil => t_soil_h
1593	ELSE
1594	surf => surf_lsm_v(l)
1595
1596	surf_t_surface => t_surface_v(l)
1597	surf_t_surface_p => t_surface_v_p(l)
1598	surf_tt_surface_m => tt_surface_v_m(l)
1599	surf_m_liq => m_liq_v(l)
1600	surf_m_liq_p => m_liq_v_p(l)
1601	surf_tm_liq_m => tm_liq_v_m(l)
1602	surf_m_soil => m_soil_v(l)
1603	surf_t_soil => t_soil_v(l)
1604	ENDIF
1605
1606	!
1607	!-- Index offset of surface element point with respect to adjoining
1608	!-- atmospheric grid point
1609	k_off = surf%koff
1610	j_off = surf%joff
1611	i_off = surf%ioff
1612
1613	!$OMP PARALLEL PRIVATE (m, i, j, k, lambda_h_sat, ke, lambda_soil, lambda_surface, &
1614	!$OMP& c_surface_tmp, f1,m_total, f2, e_s, e, f3, m_min, m_liq_max, q_s, &
1615	!$OMP& f_qsws_veg, f_qsws_soil, f_qsws_liq, f_shf, f_qsws, e_s_dt, dq_s_dt, &
1616	!$OMP& coef_1, coef_2, tend)
1617	!$OMP DO SCHEDULE (STATIC)
1618	DO m = 1, surf%ns
1619
1620	i = surf%i(m)
1621	j = surf%j(m)
1622	k = surf%k(m)
1623
1624	!
1625	!-- Define heat conductivity between surface and soil depending on surface
1626	!-- type. For vegetation, a skin layer parameterization is used. The new
1627	!-- parameterization uses a combination of two conductivities: a constant
1628	!-- conductivity for the skin layer, and a conductivity according to the
1629	!-- uppermost soil layer. For bare soil and pavements, no skin layer is
1630	!-- applied. In these cases, the temperature is assumed to be constant
1631	!-- between the surface and the first soil layer. The heat conductivity is
1632	!-- then derived from the soil/pavement properties.
1633	!-- For water surfaces, the conductivity is already set to 1E10.
1634	!-- Moreover, the heat capacity is set. For bare soil the heat capacity is
1635	!-- the capacity of the uppermost soil layer, for pavement it is that of
1636	!-- the material involved.
1637
1638	!
1639	!-- for vegetation type surfaces, the thermal conductivity of the soil is
1640	!-- needed
1641
1642	IF ( surf%vegetation_surface(m) ) THEN
1643
1644	lambda_h_sat = lambda_h_sm*(1.0_wp - surf%m_sat(nzb_soil,m)) &
1645	lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
1646
1647	ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(nzb_soil,m) / &
1648	surf%m_sat(nzb_soil,m) ) )
1649
1650	lambda_soil = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry ) &
1651	* ddz_soil(nzb_soil) * 2.0_wp
1652
1653	!
1654	!-- When bare soil is set without a thermal conductivity (no skin layer),
1655	!-- a heat capacity is that of the soil layer, otherwise it is a
1656	!-- combination of the conductivities from the skin and the soil layer
1657	IF ( surf%lambda_surface_s(m) == 0.0_wp ) THEN
1658	surf%c_surface(m) = (rho_c_soil * (1.0_wp - surf%m_sat(nzb_soil,m))&
1659	+ rho_c_water * surf_m_soil%var_2d(nzb_soil,m) ) &
1660	* dz_soil(nzb_soil) * 0.5_wp
1661	lambda_surface = lambda_soil
1662
1663	ELSE IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))&
1664	THEN
1665	lambda_surface = surf%lambda_surface_s(m) * lambda_soil &
1666	/ ( surf%lambda_surface_s(m) + lambda_soil )
1667	ELSE
1668
1669	lambda_surface = surf%lambda_surface_u(m) * lambda_soil &
1670	/ ( surf%lambda_surface_u(m) + lambda_soil )
1671	ENDIF
1672	ELSE
1673	lambda_surface = surf%lambda_surface_s(m)
1674	ENDIF
1675
1676	!
1677	!-- Set heat capacity of the skin/surface. It is ususally zero when a skin
1678	!-- layer is used, and non-zero otherwise.
1679	c_surface_tmp = surf%c_surface(m)
1680
1681	!
1682	!-- First step: calculate aerodyamic resistance. As pt, us, ts
1683	!-- are not available for the prognostic time step, data from the last
1684	!-- time step is used here. Note that this formulation is the
1685	!-- equivalent to the ECMWF formulation using drag coefficients
1686	! IF ( bulk_cloud_model ) THEN
1687	! pt1 = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1688	! qv1 = q(k,j,i) - ql(k,j,i)
1689	! ELSEIF ( cloud_droplets ) THEN
1690	! pt1 = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1691	! qv1 = q(k,j,i)
1692	! ELSE
1693	! pt1 = pt(k,j,i)
1694	! IF ( humidity ) THEN
1695	! qv1 = q(k,j,i)
1696	! ELSE
1697	! qv1 = 0.0_wp
1698	! ENDIF
1699	! ENDIF
1700	!
1701	!-- Calculation of r_a for vertical surfaces
1702	!--
1703	!-- heat transfer coefficient for forced convection along vertical walls
1704	!-- follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
1705	!--
1706	!-- H = httc (Tsfc - Tair)
1707	!-- httc = rw * (11.8 + 4.2 * Ueff) - 4.0
1708	!--
1709	!-- rw: wall patch roughness relative to 1.0 for concrete
1710	!-- Ueff: effective wind speed
1711	!-- - 4.0 is a reduction of Rowley et al (1930) formulation based on
1712	!-- Cole and Sturrock (1977)
1713	!--
1714	!-- Ucan: Canyon wind speed
1715	!-- wstar: convective velocity
1716	!-- Qs: surface heat flux
1717	!-- zH: height of the convective layer
1718	!-- wstar = (g/TcanQszH)**(1./3.)
1719
1720	!-- Effective velocity components must always
1721	!-- be defined at scalar grid point. The wall normal component is
1722	!-- obtained by simple linear interpolation. ( An alternative would
1723	!-- be an logarithmic interpolation. )
1724	!-- A roughness lenght of 0.001 is assumed for concrete (the inverse,
1725	!-- 1000 is used in the nominator for scaling)
1726	!-- To do: detailed investigation which approach gives more reliable results!
1727	!-- Please note, in case of very small friction velocity, e.g. in little
1728	!-- holes, the resistance can become negative. For this reason, limit r_a
1729	!-- to positive values.
1730	IF ( horizontal .OR. .NOT. aero_resist_kray ) THEN
1731	surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) / &
1732	( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) )
1733	ELSE
1734	surf%r_a(m) = rho_cp / ( surf%z0(m) * 1000.0_wp &
1735	* ( 11.8_wp + 4.2_wp * &
1736	SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
1737	( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
1738	( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2, &
1739	0.01_wp ) ) &
1740	) - 4.0_wp )
1741	ENDIF
1742	!
1743	!-- Make sure that the resistance does not drop to zero for neutral
1744	!-- stratification. Also, set a maximum resistance to avoid the breakdown of
1745	!-- MOST for locations with zero wind speed
1746	IF ( surf%r_a(m) < 1.0_wp ) surf%r_a(m) = 1.0_wp
1747	IF ( surf%r_a(m) > 300.0_wp ) surf%r_a(m) = 300.0_wp
1748	!
1749	!-- Second step: calculate canopy resistance r_canopy
1750	!-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1751
1752	!-- f1: correction for incoming shortwave radiation (stomata close at
1753	!-- night)
1754	f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) / &
1755	(0.81_wp * (0.004_wp * surf%rad_sw_in(m) &
1756	+ 1.0_wp)) )
1757
1758	!
1759	!-- f2: correction for soil moisture availability to plants (the
1760	!-- integrated soil moisture must thus be considered here)
1761	!-- f2 = 0 for very dry soils
1762	m_total = 0.0_wp
1763	DO ks = nzb_soil, nzt_soil
1764	m_total = m_total + surf%root_fr(ks,m) &
1765	* MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(ks,m) )
1766	ENDDO
1767
1768	!
1769	!-- The calculation of f2 is based on only one wilting point value for all
1770	!-- soil layers. The value at k=nzb_soil is used here as a proxy but might
1771	!-- need refinement in the future.
1772	IF ( m_total > surf%m_wilt(nzb_soil,m) .AND. &
1773	m_total < surf%m_fc(nzb_soil,m) ) THEN
1774	f2 = ( m_total - surf%m_wilt(nzb_soil,m) ) / &
1775	( surf%m_fc(nzb_soil,m) - surf%m_wilt(nzb_soil,m) )
1776	ELSEIF ( m_total >= surf%m_fc(nzb_soil,m) ) THEN
1777	f2 = 1.0_wp
1778	ELSE
1779	f2 = 1.0E-20_wp
1780	ENDIF
1781
1782	!
1783	!-- Calculate water vapour pressure at saturation and convert to hPa
1784	!-- The magnus formula is limited to temperatures up to 333.15 K to
1785	!-- avoid negative values of q_s
1786	e_s = 0.01_wp * magnus( MIN(surf_t_surface%var_1d(m), 333.15_wp) )
1787
1788	!
1789	!-- f3: correction for vapour pressure deficit
1790	IF ( surf%g_d(m) /= 0.0_wp ) THEN
1791	!
1792	!-- Calculate vapour pressure
1793	e = surf%qv1(m) * surface_pressure / ( surf%qv1(m) + rd_d_rv )
1794	f3 = EXP ( - surf%g_d(m) * (e_s - e) )
1795	ELSE
1796	f3 = 1.0_wp
1797	ENDIF
1798	!
1799	!-- Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1800	!-- this calculation is obsolete, as r_canopy is not used below.
1801	!-- To do: check for very dry soil -> r_canopy goes to infinity
1802	surf%r_canopy(m) = surf%r_canopy_min(m) / &
1803	( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
1804	!
1805	!-- Third step: calculate bare soil resistance r_soil.
1806	m_min = surf%c_veg(m) * surf%m_wilt(nzb_soil,m) + &
1807	( 1.0_wp - surf%c_veg(m) ) * surf%m_res(nzb_soil,m)
1808
1809
1810	f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) / &
1811	( surf%m_fc(nzb_soil,m) - m_min )
1812	f2 = MAX( f2, 1.0E-20_wp )
1813	f2 = MIN( f2, 1.0_wp )
1814
1815	surf%r_soil(m) = surf%r_soil_min(m) / f2
1816
1817	!
1818	!-- Calculate the maximum possible liquid water amount on plants and
1819	!-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
1820	!-- assumed, while paved surfaces might hold up 1 mm of water. The
1821	!-- liquid water fraction for paved surfaces is calculated after
1822	!-- Noilhan & Planton (1989), while the ECMWF formulation is used for
1823	!-- vegetated surfaces and bare soils.
1824	IF ( surf%pavement_surface(m) ) THEN
1825	m_liq_max = m_max_depth * 5.0_wp
1826	surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
1827	ELSE
1828	m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m) &
1829	+ ( 1.0_wp - surf%c_veg(m) ) )
1830	surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq%var_1d(m) / m_liq_max )
1831	ENDIF
1832	!
1833	!-- Calculate saturation water vapor mixing ratio
1834	q_s = rd_d_rv * e_s / ( surface_pressure - e_s )
1835	!
1836	!-- In case of dewfall, set evapotranspiration to zero
1837	!-- All super-saturated water is then removed from the air
1838	IF ( humidity .AND. q_s <= surf%qv1(m) ) THEN
1839	surf%r_canopy(m) = 0.0_wp
1840	surf%r_soil(m) = 0.0_wp
1841	ENDIF
1842
1843	!
1844	!-- Calculate coefficients for the total evapotranspiration
1845	!-- In case of water surface, set vegetation and soil fluxes to zero.
1846	!-- For pavements, only evaporation of liquid water is possible.
1847	IF ( surf%water_surface(m) ) THEN
1848	f_qsws_veg = 0.0_wp
1849	f_qsws_soil = 0.0_wp
1850	f_qsws_liq = rho_lv / surf%r_a(m)
1851	ELSEIF ( surf%pavement_surface(m) ) THEN
1852	f_qsws_veg = 0.0_wp
1853	f_qsws_soil = 0.0_wp
1854	f_qsws_liq = rho_lv * surf%c_liq(m) / surf%r_a(m)
1855	ELSE
1856	f_qsws_veg = rho_lv * surf%c_veg(m) * &
1857	( 1.0_wp - surf%c_liq(m) ) / &
1858	( surf%r_a(m) + surf%r_canopy(m) )
1859	f_qsws_soil = rho_lv * (1.0_wp - surf%c_veg(m) ) / &
1860	( surf%r_a(m) + surf%r_soil(m) )
1861	f_qsws_liq = rho_lv * surf%c_veg(m) * surf%c_liq(m) / &
1862	surf%r_a(m)
1863	ENDIF
1864
1865	f_shf = rho_cp / surf%r_a(m)
1866	f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1867	!
1868	!-- Calculate derivative of q_s for Taylor series expansion
1869	e_s_dt = e_s * ( 17.62_wp / ( surf_t_surface%var_1d(m) - 29.65_wp) - &
1870	17.62_wp*( surf_t_surface%var_1d(m) - 273.15_wp) &
1871	/ ( surf_t_surface%var_1d(m) - 29.65_wp)**2 )
1872
1873	dq_s_dt = rd_d_rv * e_s_dt / ( surface_pressure - e_s_dt )
1874	!
1875	!-- Calculate net radiation radiation without longwave outgoing flux because
1876	!-- it has a dependency on surface temperature and thus enters the prognostic
1877	!-- equations directly
1878	surf%rad_net_l(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m) &
1879	+ surf%rad_lw_in(m)
1880	!
1881	!-- Calculate new skin temperature
1882	IF ( humidity ) THEN
1883	!
1884	!-- Numerator of the prognostic equation
1885	coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m) &
1886	* surf_t_surface%var_1d(m) - surf%rad_lw_out(m) &
1887	+ f_shf * surf%pt1(m) + f_qsws * ( surf%qv1(m) - q_s &
1888	+ dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface &
1889	* surf_t_soil%var_2d(nzb_soil,m)
1890
1891	!
1892	!-- Denominator of the prognostic equation
1893	coef_2 = surf%rad_lw_out_change_0(m) + f_qsws * dq_s_dt &
1894	+ lambda_surface + f_shf / exner(nzb)
1895	ELSE
1896	!
1897	!-- Numerator of the prognostic equation
1898	coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m) &
1899	* surf_t_surface%var_1d(m) - surf%rad_lw_out(m) &
1900	+ f_shf * surf%pt1(m) + lambda_surface &
1901	* surf_t_soil%var_2d(nzb_soil,m)
1902	!
1903	!-- Denominator of the prognostic equation
1904	coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exner(nzb)
1905
1906	ENDIF
1907
1908	tend = 0.0_wp
1909
1910	!
1911	!-- Implicit solution when the surface layer has no heat capacity,
1912	!-- otherwise use RK3 scheme.
1913	surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
1914	surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
1915	* dt_3d * tsc(2) )
1916
1917	!
1918	!-- Add RK3 term
1919	IF ( c_surface_tmp /= 0.0_wp ) THEN
1920
1921	surf_t_surface_p%var_1d(m) = surf_t_surface_p%var_1d(m) + dt_3d * &
1922	tsc(3) * surf_tt_surface_m%var_1d(m)
1923
1924	!
1925	!-- Calculate true tendency
1926	tend = ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) - &
1927	dt_3d * tsc(3) * surf_tt_surface_m%var_1d(m)) / (dt_3d * tsc(2))
1928	!
1929	!-- Calculate t_surface tendencies for the next Runge-Kutta step
1930	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1931	IF ( intermediate_timestep_count == 1 ) THEN
1932	surf_tt_surface_m%var_1d(m) = tend
1933	ELSEIF ( intermediate_timestep_count < &
1934	intermediate_timestep_count_max ) THEN
1935	surf_tt_surface_m%var_1d(m) = -9.5625_wp * tend + &
1936	5.3125_wp * surf_tt_surface_m%var_1d(m)
1937	ENDIF
1938	ENDIF
1939	ENDIF
1940
1941	!
1942	!-- In case of fast changes in the skin temperature, it is possible to
1943	!-- update the radiative fluxes independently from the prescribed
1944	!-- radiation call frequency. This effectively prevents oscillations,
1945	!-- especially when setting skip_time_do_radiation /= 0. The threshold
1946	!-- value of 0.2 used here is just a first guess. This method should be
1947	!-- revised in the future as tests have shown that the threshold is
1948	!-- often reached, when no oscillations would occur (causes immense
1949	!-- computing time for the radiation code).
1950	IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) ) &
1951	> 0.2_wp .AND. &
1952	unscheduled_radiation_calls ) THEN
1953	force_radiation_call_l = .TRUE.
1954	ENDIF
1955
1956	surf%pt_surface(m) = surf_t_surface_p%var_1d(m) / exner(nzb)
1957
1958	!
1959	!-- Calculate fluxes
1960	surf%rad_net_l(m) = surf%rad_net_l(m) + &
1961	surf%rad_lw_out_change_0(m) &
1962	* surf_t_surface%var_1d(m) - surf%rad_lw_out(m) &
1963	- surf%rad_lw_out_change_0(m) * surf_t_surface_p%var_1d(m)
1964
1965	surf%rad_net(m) = surf%rad_net_l(m)
1966	surf%rad_lw_out(m) = surf%rad_lw_out(m) + surf%rad_lw_out_change_0(m) * &
1967	( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
1968
1969	surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m) &
1970	- surf_t_soil%var_2d(nzb_soil,m) )
1971
1972	surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / c_p
1973
1974	!
1975	! update the 3d field of rad_lw_out array to have consistent output
1976	IF ( horizontal ) THEN
1977	IF ( radiation_scheme == 'rrtmg' ) THEN
1978	rad_lw_out(k+k_off,j+j_off,i+i_off) = surf%rad_lw_out(m)
1979	ELSE
1980	rad_lw_out(0,j+j_off,i+i_off) = surf%rad_lw_out(m)
1981	ENDIF
1982	ENDIF
1983
1984	IF ( humidity ) THEN
1985	surf%qsws(m) = - f_qsws * ( surf%qv1(m) - q_s + dq_s_dt &
1986	* surf_t_surface%var_1d(m) - dq_s_dt * &
1987	surf_t_surface_p%var_1d(m) )
1988
1989	surf%qsws_veg(m) = - f_qsws_veg * ( surf%qv1(m) - q_s &
1990	+ dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt &
1991	* surf_t_surface_p%var_1d(m) )
1992
1993	surf%qsws_soil(m) = - f_qsws_soil * ( surf%qv1(m) - q_s &
1994	+ dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt &
1995	* surf_t_surface_p%var_1d(m) )
1996
1997	surf%qsws_liq(m) = - f_qsws_liq * ( surf%qv1(m) - q_s &
1998	+ dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt &
1999	* surf_t_surface_p%var_1d(m) )
2000	ENDIF
2001
2002	!
2003	!-- Calculate the true surface resistance. ABS is used here to avoid negative
2004	!-- values that can occur for very small fluxes due to the artifical addition
2005	!-- of 1.0E-20.
2006	IF ( .NOT. humidity ) THEN
2007	surf%r_s(m) = 1.0E10_wp
2008	ELSE
2009	surf%r_s(m) = ABS(rho_lv / (f_qsws + 1.0E-20_wp) - surf%r_a(m))
2010	ENDIF
2011	!
2012	!-- Calculate change in liquid water reservoir due to dew fall or
2013	!-- evaporation of liquid water
2014	IF ( humidity ) THEN
2015	!
2016	!-- If precipitation is activated, add rain water to qsws_liq
2017	!-- and qsws_soil according the the vegetation coverage.
2018	!-- precipitation_rate is given in mm.
2019	IF ( precipitation ) THEN
2020
2021	!
2022	!-- Add precipitation to liquid water reservoir, if possible.
2023	!-- Otherwise, add the water to soil. In case of
2024	!-- pavements, the exceeding water amount is explicitly removed
2025	!-- (as fictive runoff by drainage systems)
2026	IF ( surf%pavement_surface(m) ) THEN
2027	IF ( surf_m_liq%var_1d(m) < m_liq_max ) THEN
2028	surf%qsws_liq(m) = surf%qsws_liq(m) &
2029	+ prr(k+k_off,j+j_off,i+i_off) &
2030	* hyrho(k+k_off) &
2031	* 0.001_wp * rho_l * l_v
2032	ENDIF
2033	ELSE
2034	IF ( surf_m_liq%var_1d(m) < m_liq_max ) THEN
2035	surf%qsws_liq(m) = surf%qsws_liq(m) &
2036	+ surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
2037	* hyrho(k+k_off) &
2038	* 0.001_wp * rho_l * l_v
2039	surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp - &
2040	surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
2041	* hyrho(k+k_off) &
2042	* 0.001_wp * rho_l * l_v
2043	ELSE
2044
2045	!-- Add precipitation to bare soil according to the bare soil
2046	!-- coverage.
2047	surf%qsws_soil(m) = surf%qsws_soil(m) &
2048	+ surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
2049	* hyrho(k+k_off) &
2050	* 0.001_wp * rho_l * l_v
2051
2052	ENDIF
2053	ENDIF
2054
2055	ENDIF
2056
2057	!
2058	!-- If the air is saturated, check the reservoir water level
2059	IF ( surf%qsws(m) < 0.0_wp ) THEN
2060	!
2061	!-- Check if reservoir is full (avoid values > m_liq_max)
2062	!-- In that case, qsws_liq goes to qsws_soil for pervious surfaces. In
2063	!-- this case qsws_veg is zero anyway (because c_liq = 1),
2064	!-- so that tend is zero and no further check is needed
2065	IF ( surf_m_liq%var_1d(m) == m_liq_max ) THEN
2066	IF ( .NOT. surf%pavement_surface(m)) THEN
2067	surf%qsws_soil(m) = surf%qsws_soil(m) + surf%qsws_liq(m)
2068	ENDIF
2069	surf%qsws_liq(m) = 0.0_wp
2070	ENDIF
2071
2072	!
2073	!-- In case qsws_veg becomes negative (unphysical behavior),
2074	!-- let the water enter the liquid water reservoir as dew on the
2075	!-- plant
2076	IF ( surf%qsws_veg(m) < 0.0_wp ) THEN
2077	surf%qsws_liq(m) = surf%qsws_liq(m) + surf%qsws_veg(m)
2078	surf%qsws_veg(m) = 0.0_wp
2079	ENDIF
2080	ENDIF
2081
2082	surf%qsws(m) = surf%qsws(m) / l_v
2083
2084	tend = - surf%qsws_liq(m) * drho_l_lv
2085	surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d * &
2086	( tsc(2) * tend + &
2087	tsc(3) * surf_tm_liq_m%var_1d(m) )
2088	!
2089	!-- Check if reservoir is overfull -> reduce to maximum
2090	!-- (conservation of water is violated here)
2091	surf_m_liq_p%var_1d(m) = MIN( surf_m_liq_p%var_1d(m),m_liq_max )
2092
2093	!
2094	!-- Check if reservoir is empty (avoid values < 0.0)
2095	!-- (conservation of water is violated here)
2096	surf_m_liq_p%var_1d(m) = MAX( surf_m_liq_p%var_1d(m), 0.0_wp )
2097	!
2098	!-- Calculate m_liq tendencies for the next Runge-Kutta step
2099	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2100	IF ( intermediate_timestep_count == 1 ) THEN
2101	surf_tm_liq_m%var_1d(m) = tend
2102	ELSEIF ( intermediate_timestep_count < &
2103	intermediate_timestep_count_max ) THEN
2104	surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend + &
2105	5.3125_wp * surf_tm_liq_m%var_1d(m)
2106	ENDIF
2107	ENDIF
2108
2109	ENDIF
2110
2111	ENDDO
2112	!$OMP END PARALLEL
2113
2114	!
2115	!-- Make a logical OR for all processes. Force radiation call if at
2116	!-- least one processor reached the threshold change in skin temperature
2117	IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count &
2118	== intermediate_timestep_count_max-1 ) THEN
2119	#if defined( __parallel )
2120	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2121	CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, &
2122	1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
2123	#else
2124	force_radiation_call = force_radiation_call_l
2125	#endif
2126	force_radiation_call_l = .FALSE.
2127	ENDIF
2128
2129	!
2130	!-- Calculate surface water vapor mixing ratio
2131	IF ( humidity ) THEN
2132	CALL calc_q_surface
2133	ENDIF
2134
2135	!
2136	!-- Calculate new roughness lengths (for water surfaces only)
2137	IF ( horizontal .AND. .NOT. constant_roughness ) CALL calc_z0_water_surface
2138
2139	IF ( debug_output_timestep ) THEN
2140	WRITE( debug_string, * ) 'lsm_energy_balance', horizontal, l
2141	CALL debug_message( debug_string, 'end' )
2142	ENDIF
2143
2144	CONTAINS
2145	!------------------------------------------------------------------------------!
2146	! Description:
2147	! ------------
2148	!> Calculation of mixing ratio of the skin layer (surface). It is assumend
2149	!> that the skin is always saturated.
2150	!------------------------------------------------------------------------------!
2151	SUBROUTINE calc_q_surface
2152
2153	IMPLICIT NONE
2154
2155	REAL(wp) :: e_s !< saturation water vapor pressure
2156	REAL(wp) :: q_s !< saturation mixing ratio
2157	REAL(wp) :: resistance !< aerodynamic and soil resistance term
2158
2159
2160	!$OMP PARALLEL PRIVATE (m, i, j, k, e_s, q_s, resistance)
2161	!$OMP DO SCHEDULE (STATIC)
2162	DO m = 1, surf%ns
2163
2164	i = surf%i(m)
2165	j = surf%j(m)
2166	k = surf%k(m)
2167	!
2168	!-- Calculate water vapour pressure at saturation and convert to hPa
2169	e_s = 0.01_wp * magnus( MIN(surf_t_surface_p%var_1d(m), 333.15_wp) )
2170
2171	!
2172	!-- Calculate mixing ratio at saturation
2173	q_s = rd_d_rv * e_s / ( surface_pressure - e_s )
2174
2175	resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) + 1E-5_wp )
2176
2177	!
2178	!-- Calculate mixing ratio at surface
2179	IF ( bulk_cloud_model ) THEN
2180	q(k+k_off,j+j_off,i+i_off) = resistance * q_s + &
2181	( 1.0_wp - resistance ) * &
2182	( q(k,j,i) - ql(k,j,i) )
2183	ELSE
2184	q(k+k_off,j+j_off,i+i_off) = resistance * q_s + &
2185	( 1.0_wp - resistance ) * &
2186	q(k,j,i)
2187	ENDIF
2188
2189	surf%q_surface(m) = q(k+k_off,j+j_off,i+i_off)
2190	!
2191	!-- Update virtual potential temperature
2192	surf%vpt_surface(m) = surf%pt_surface(m) * &
2193	( 1.0_wp + 0.61_wp * surf%q_surface(m) )
2194
2195
2196
2197	ENDDO
2198	!$OMP END PARALLEL
2199
2200	END SUBROUTINE calc_q_surface
2201
2202	END SUBROUTINE lsm_energy_balance
2203
2204
2205
2206	!------------------------------------------------------------------------------!
2207	! Description:
2208	! ------------
2209	!> Header output for land surface model
2210	!------------------------------------------------------------------------------!
2211	SUBROUTINE lsm_header ( io )
2212
2213
2214	IMPLICIT NONE
2215
2216	CHARACTER (LEN=86) :: t_soil_chr !< String for soil temperature profile
2217	CHARACTER (LEN=86) :: roots_chr !< String for root profile
2218	CHARACTER (LEN=86) :: vertical_index_chr !< String for the vertical index
2219	CHARACTER (LEN=86) :: m_soil_chr !< String for soil moisture
2220	CHARACTER (LEN=86) :: soil_depth_chr !< String for soil depth
2221	CHARACTER (LEN=10) :: coor_chr !< Temporary string
2222
2223	INTEGER(iwp) :: i !< Loop index over soil layers
2224
2225	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
2226
2227	t_soil_chr = ''
2228	m_soil_chr = ''
2229	soil_depth_chr = ''
2230	roots_chr = ''
2231	vertical_index_chr = ''
2232
2233	i = 1
2234	DO i = nzb_soil, nzt_soil
2235	WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
2236	t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr )
2237
2238	WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
2239	m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr )
2240
2241	WRITE (coor_chr,'(F10.2,7X)') - zs(i)
2242	soil_depth_chr = TRIM( soil_depth_chr ) // ' ' // TRIM( coor_chr )
2243
2244	WRITE (coor_chr,'(F10.2,7X)') root_fraction(i)
2245	roots_chr = TRIM( roots_chr ) // ' ' // TRIM( coor_chr )
2246
2247	WRITE (coor_chr,'(I10,7X)') i
2248	vertical_index_chr = TRIM( vertical_index_chr ) // ' ' // &
2249	TRIM( coor_chr )
2250	ENDDO
2251
2252	!
2253	!-- Write land surface model header
2254	WRITE( io, 1 )
2255	IF ( conserve_water_content ) THEN
2256	WRITE( io, 2 )
2257	ELSE
2258	WRITE( io, 3 )
2259	ENDIF
2260
2261	IF ( vegetation_type_f%from_file ) THEN
2262	WRITE( io, 5 )
2263	ELSE
2264	WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ), &
2265	TRIM (soil_type_name(soil_type) )
2266	ENDIF
2267	WRITE( io, 6 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ), &
2268	TRIM( m_soil_chr ), TRIM( roots_chr ), &
2269	TRIM( vertical_index_chr )
2270
2271	1 FORMAT (//' Land surface model information:'/ &
2272	' ------------------------------'/)
2273	2 FORMAT (' --> Soil bottom is closed (water content is conserved', &
2274	', default)')
2275	3 FORMAT (' --> Soil bottom is open (water content is not conserved)')
2276	4 FORMAT (' --> Land surface type : ',A,/ &
2277	' --> Soil porosity type : ',A)
2278	5 FORMAT (' --> Land surface type : read from file' / &
2279	' --> Soil porosity type : read from file' )
2280	6 FORMAT (/' Initial soil temperature and moisture profile:'// &
2281	' Height: ',A,' m'/ &
2282	' Temperature: ',A,' K'/ &
2283	' Moisture: ',A,' m3/m3'/ &
2284	' Root fraction: ',A,' '/ &
2285	' Grid point: ',A)
2286
2287
2288	END SUBROUTINE lsm_header
2289
2290
2291	!------------------------------------------------------------------------------!
2292	! Description:
2293	! ------------
2294	!> Initialization of the land surface model
2295	!------------------------------------------------------------------------------!
2296	SUBROUTINE lsm_init
2297
2298	USE control_parameters, &
2299	ONLY: message_string
2300
2301	USE indices, &
2302	ONLY: nx, ny, topo_min_level
2303
2304	USE pmc_handle_communicator, &
2305	ONLY: pmc_is_rootmodel
2306
2307	USE pmc_interface, &
2308	ONLY: nested_run
2309
2310	IMPLICIT NONE
2311
2312	INTEGER(iwp) :: i !< running index
2313	INTEGER(iwp) :: j !< running index
2314	INTEGER(iwp) :: k !< running index
2315	INTEGER(iwp) :: kn !< running index
2316	INTEGER(iwp) :: ko !< running index
2317	INTEGER(iwp) :: kroot !< running index
2318	INTEGER(iwp) :: kzs !< running index
2319	INTEGER(iwp) :: l !< running index surface facing
2320	INTEGER(iwp) :: m !< running index
2321	INTEGER(iwp) :: st !< soil-type index
2322	INTEGER(iwp) :: n_soil_layers_total !< temperature variable, stores the total number of soil layers + 4
2323	#if defined( __parallel )
2324	INTEGER(iwp) :: nzs_root !< number of soil layers in root domain (used in case soil data needs to be
2325	!< transferred from root model to child models)
2326
2327	REAL(wp), DIMENSION(:), ALLOCATABLE :: m_soil_root !< domain-averaged soil moisture profile in root domain
2328	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_soil_root !< domain-averaged soil temperature profile in root domain
2329	#endif
2330	REAL(wp), DIMENSION(:), ALLOCATABLE :: bound !< temporary arrays for storing index bounds
2331	REAL(wp), DIMENSION(:), ALLOCATABLE :: bound_root_fr !< temporary arrays for storing index bounds
2332	REAL(wp), DIMENSION(:), ALLOCATABLE :: pr_soil_init !< temporary array used for averaging soil profiles
2333	#if defined( __parallel )
2334	REAL(wp), DIMENSION(:), ALLOCATABLE :: z_soil_root !< vertical dimension of soil grid in root domain
2335	#endif
2336
2337	IF ( debug_output ) CALL debug_message( 'lsm_init', 'start' )
2338	!
2339	!-- If no cloud physics is used, rho_surface has not been calculated before
2340	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
2341	CALL calc_mean_profile( pt, 4 )
2342	rho_surface = hyp(nzb) / ( r_d * hom(topo_min_level+1,1,4,0) * exner(nzb) )
2343	ENDIF
2344
2345	!
2346	!-- Calculate frequently used parameters
2347	rho_cp = c_p * rho_surface
2348	rho_lv = rho_surface * l_v
2349	drho_l_lv = 1.0_wp / (rho_l * l_v)
2350
2351	!
2352	!-- Set initial values for prognostic quantities
2353	!-- Horizontal surfaces
2354	tt_surface_h_m%var_1d = 0.0_wp
2355	tt_soil_h_m%var_2d = 0.0_wp
2356	tm_soil_h_m%var_2d = 0.0_wp
2357	tm_liq_h_m%var_1d = 0.0_wp
2358	surf_lsm_h%c_liq = 0.0_wp
2359
2360	surf_lsm_h%ghf = 0.0_wp
2361
2362	surf_lsm_h%qsws_liq = 0.0_wp
2363	surf_lsm_h%qsws_soil = 0.0_wp
2364	surf_lsm_h%qsws_veg = 0.0_wp
2365
2366	surf_lsm_h%r_a = 50.0_wp
2367	surf_lsm_h%r_s = 50.0_wp
2368	surf_lsm_h%r_canopy = 0.0_wp
2369	surf_lsm_h%r_soil = 0.0_wp
2370	!
2371	!-- Do the same for vertical surfaces
2372	DO l = 0, 3
2373	tt_surface_v_m(l)%var_1d = 0.0_wp
2374	tt_soil_v_m(l)%var_2d = 0.0_wp
2375	tm_soil_v_m(l)%var_2d = 0.0_wp
2376	tm_liq_v_m(l)%var_1d = 0.0_wp
2377	surf_lsm_v(l)%c_liq = 0.0_wp
2378
2379	surf_lsm_v(l)%ghf = 0.0_wp
2380
2381	surf_lsm_v(l)%qsws_liq = 0.0_wp
2382	surf_lsm_v(l)%qsws_soil = 0.0_wp
2383	surf_lsm_v(l)%qsws_veg = 0.0_wp
2384
2385	surf_lsm_v(l)%r_a = 50.0_wp
2386	surf_lsm_v(l)%r_s = 50.0_wp
2387	surf_lsm_v(l)%r_canopy = 0.0_wp
2388	surf_lsm_v(l)%r_soil = 0.0_wp
2389	ENDDO
2390
2391	!
2392	!-- Set initial values for prognostic soil quantities
2393	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
2394	t_soil_h%var_2d = 0.0_wp
2395	m_soil_h%var_2d = 0.0_wp
2396	m_liq_h%var_1d = 0.0_wp
2397
2398	DO l = 0, 3
2399	t_soil_v(l)%var_2d = 0.0_wp
2400	m_soil_v(l)%var_2d = 0.0_wp
2401	m_liq_v(l)%var_1d = 0.0_wp
2402	ENDDO
2403	ENDIF
2404	!
2405	!-- Allocate 3D soil model arrays
2406	!-- First, for horizontal surfaces
2407	ALLOCATE ( surf_lsm_h%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2408	ALLOCATE ( surf_lsm_h%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2409	ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2410	ALLOCATE ( surf_lsm_h%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns))
2411	ALLOCATE ( surf_lsm_h%l_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2412	ALLOCATE ( surf_lsm_h%m_fc(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2413	ALLOCATE ( surf_lsm_h%m_res(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2414	ALLOCATE ( surf_lsm_h%m_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2415	ALLOCATE ( surf_lsm_h%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2416	ALLOCATE ( surf_lsm_h%n_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2417	ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2418	ALLOCATE ( surf_lsm_h%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2419	ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2420
2421	surf_lsm_h%lambda_h = 0.0_wp
2422	!
2423	!-- If required, allocate humidity-related variables for the soil model
2424	IF ( humidity ) THEN
2425	ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2426	ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2427
2428	surf_lsm_h%lambda_w = 0.0_wp
2429	ENDIF
2430	!
2431	!-- For vertical surfaces
2432	DO l = 0, 3
2433	ALLOCATE ( surf_lsm_v(l)%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2434	ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2435	ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2436	ALLOCATE ( surf_lsm_v(l)%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns))
2437	ALLOCATE ( surf_lsm_v(l)%l_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2438	ALLOCATE ( surf_lsm_v(l)%m_fc(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2439	ALLOCATE ( surf_lsm_v(l)%m_res(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2440	ALLOCATE ( surf_lsm_v(l)%m_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2441	ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2442	ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2443	ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2444	ALLOCATE ( surf_lsm_v(l)%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2445	ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2446
2447	surf_lsm_v(l)%lambda_h = 0.0_wp
2448
2449	!
2450	!-- If required, allocate humidity-related variables for the soil model
2451	IF ( humidity ) THEN
2452	ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2453	ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2454
2455	surf_lsm_v(l)%lambda_w = 0.0_wp
2456	ENDIF
2457	ENDDO
2458	!
2459	!-- Allocate albedo type and emissivity for vegetation, water and pavement
2460	!-- fraction.
2461	!-- Set default values at each surface element.
2462	ALLOCATE ( surf_lsm_h%albedo_type(0:2,1:surf_lsm_h%ns) )
2463	ALLOCATE ( surf_lsm_h%emissivity(0:2,1:surf_lsm_h%ns) )
2464	!
2465	!-- Initialize albedo type according to its default type, in order to set values
2466	!-- independent on default albedo_type in radiation model.
2467	surf_lsm_h%albedo_type(ind_veg_wall,:) = &
2468	INT( vegetation_pars(ind_v_at,vegetation_type) )
2469	surf_lsm_h%albedo_type(ind_wat_win,:) = &
2470	INT( water_pars(ind_w_at,water_type) )
2471	surf_lsm_h%albedo_type(ind_pav_green,:) = &
2472	INT( pavement_pars(ind_p_at,pavement_type) )
2473	surf_lsm_h%emissivity = emissivity
2474	DO l = 0, 3
2475	ALLOCATE ( surf_lsm_v(l)%albedo_type(0:2,1:surf_lsm_v(l)%ns) )
2476	ALLOCATE ( surf_lsm_v(l)%emissivity(0:2,1:surf_lsm_v(l)%ns) )
2477	!
2478	!-- Initialize albedo type according to its default type, in order to
2479	!-- set values independent on default albedo_type in radiation model.
2480	surf_lsm_v(l)%albedo_type(ind_veg_wall,:) = &
2481	INT( vegetation_pars(ind_v_at,vegetation_type) )
2482	surf_lsm_v(l)%albedo_type(ind_wat_win,:) = &
2483	INT( water_pars(ind_w_at,water_type) )
2484	surf_lsm_v(l)%albedo_type(ind_pav_green,:) = &
2485	INT( pavement_pars(ind_p_at,pavement_type) )
2486	surf_lsm_v(l)%emissivity = emissivity
2487	ENDDO
2488	!
2489	!-- Allocate arrays for relative surface fraction.
2490	!-- 0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction
2491	ALLOCATE( surf_lsm_h%frac(0:2,1:surf_lsm_h%ns) )
2492	surf_lsm_h%frac = 0.0_wp
2493	DO l = 0, 3
2494	ALLOCATE( surf_lsm_v(l)%frac(0:2,1:surf_lsm_v(l)%ns) )
2495	surf_lsm_v(l)%frac = 0.0_wp
2496	ENDDO
2497	!
2498	!-- For vertical walls only - allocate special flag indicating if any building is on
2499	!-- top of any natural surfaces. Used for initialization only.
2500	DO l = 0, 3
2501	ALLOCATE( surf_lsm_v(l)%building_covered(1:surf_lsm_v(l)%ns) )
2502	ENDDO
2503	!
2504	!-- Allocate arrays for the respective types and their names on the surface
2505	!-- elements. This will be required to treat deposition of chemical species.
2506	ALLOCATE( surf_lsm_h%pavement_type(1:surf_lsm_h%ns) )
2507	ALLOCATE( surf_lsm_h%vegetation_type(1:surf_lsm_h%ns) )
2508	ALLOCATE( surf_lsm_h%water_type(1:surf_lsm_h%ns) )
2509
2510	surf_lsm_h%pavement_type = 0
2511	surf_lsm_h%vegetation_type = 0
2512	surf_lsm_h%water_type = 0
2513
2514	ALLOCATE( surf_lsm_h%pavement_type_name(1:surf_lsm_h%ns) )
2515	ALLOCATE( surf_lsm_h%vegetation_type_name(1:surf_lsm_h%ns) )
2516	ALLOCATE( surf_lsm_h%water_type_name(1:surf_lsm_h%ns) )
2517
2518	surf_lsm_h%pavement_type_name = 'none'
2519	surf_lsm_h%vegetation_type_name = 'none'
2520	surf_lsm_h%water_type_name = 'none'
2521
2522	DO l = 0, 3
2523	ALLOCATE( surf_lsm_v(l)%pavement_type(1:surf_lsm_v(l)%ns) )
2524	ALLOCATE( surf_lsm_v(l)%vegetation_type(1:surf_lsm_v(l)%ns) )
2525	ALLOCATE( surf_lsm_v(l)%water_type(1:surf_lsm_v(l)%ns) )
2526
2527	surf_lsm_v(l)%pavement_type = 0
2528	surf_lsm_v(l)%vegetation_type = 0
2529	surf_lsm_v(l)%water_type = 0
2530
2531	ALLOCATE( surf_lsm_v(l)%pavement_type_name(1:surf_lsm_v(l)%ns) )
2532	ALLOCATE( surf_lsm_v(l)%vegetation_type_name(1:surf_lsm_v(l)%ns) )
2533	ALLOCATE( surf_lsm_v(l)%water_type_name(1:surf_lsm_v(l)%ns) )
2534
2535	surf_lsm_v(l)%pavement_type_name = 'none'
2536	surf_lsm_v(l)%vegetation_type_name = 'none'
2537	surf_lsm_v(l)%water_type_name = 'none'
2538	ENDDO
2539
2540	!
2541	!-- Set flag parameter for the prescribed surface type depending on user
2542	!-- input. Set surface fraction to 1 for the respective type.
2543	SELECT CASE ( TRIM( surface_type ) )
2544
2545	CASE ( 'vegetation' )
2546
2547	surf_lsm_h%vegetation_surface = .TRUE.
2548	surf_lsm_h%frac(ind_veg_wall,:) = 1.0_wp
2549	DO l = 0, 3
2550	surf_lsm_v(l)%vegetation_surface = .TRUE.
2551	surf_lsm_v(l)%frac(ind_veg_wall,:) = 1.0_wp
2552	ENDDO
2553
2554	CASE ( 'water' )
2555
2556	surf_lsm_h%water_surface = .TRUE.
2557	surf_lsm_h%frac(ind_wat_win,:) = 1.0_wp
2558	!
2559	!-- Note, vertical water surface does not really make sense.
2560	DO l = 0, 3
2561	surf_lsm_v(l)%water_surface = .TRUE.
2562	surf_lsm_v(l)%frac(ind_wat_win,:) = 1.0_wp
2563	ENDDO
2564
2565	CASE ( 'pavement' )
2566
2567	surf_lsm_h%pavement_surface = .TRUE.
2568	surf_lsm_h%frac(ind_pav_green,:) = 1.0_wp
2569	DO l = 0, 3
2570	surf_lsm_v(l)%pavement_surface = .TRUE.
2571	surf_lsm_v(l)%frac(ind_pav_green,:) = 1.0_wp
2572	ENDDO
2573
2574	CASE ( 'netcdf' )
2575
2576	DO m = 1, surf_lsm_h%ns
2577	i = surf_lsm_h%i(m)
2578	j = surf_lsm_h%j(m)
2579	IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
2580	surf_lsm_h%vegetation_surface(m) = .TRUE.
2581	IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill ) &
2582	surf_lsm_h%pavement_surface(m) = .TRUE.
2583	IF ( water_type_f%var(j,i) /= water_type_f%fill ) &
2584	surf_lsm_h%water_surface(m) = .TRUE.
2585	!
2586	!-- Check if at least one type is set.
2587	IF ( .NOT. surf_lsm_h%vegetation_surface(m) .AND. &
2588	.NOT. surf_lsm_h%pavement_surface(m) .AND. &
2589	.NOT. surf_lsm_h%water_surface(m) ) THEN
2590	WRITE( message_string, * ) 'Horizontal surface element ' // &
2591	' at i, j = ', i, j, &
2592	' is neither a vegetation, ' // &
2593	'pavement, nor a water surface.'
2594	CALL message( 'land_surface_model_mod', 'PA0619', &
2595	2, 2, myid, 6, 0 )
2596	ENDIF
2597
2598	ENDDO
2599	!
2600	!-- For vertical surfaces some special checks and treatment are
2601	!-- required for correct initialization.
2602	DO l = 0, 3
2603	DO m = 1, surf_lsm_v(l)%ns
2604	!
2605	!-- Only for vertical surfaces. Check if at the grid point where
2606	!-- the wall is defined (i+ioff, j+joff) is any building.
2607	!-- This case, no natural surfaces properties will be defined at
2608	!-- at this grid point, leading to problems in the initialization.
2609	!-- To overcome this, define a special flag which
2610	!-- indicates that a building is defined at the wall grid point
2611	!-- and take the surface properties from the adjoining grid
2612	!-- point, i.e. without offset values.
2613	!-- Further, there can occur a special case where elevation
2614	!-- changes are larger than building heights. This case, (j,i)
2615	!-- and (j+joff,i+ioff) grid points may be both covered by
2616	!-- buildings, but vertical, but vertically natural walls may
2617	!-- be located between the buildings. This case, it is not
2618	!-- guaranteed that information about natural surface types is
2619	!-- given, neither at (j,i) nor at (j+joff,i+ioff), again leading
2620	!-- to non-initialized surface properties.
2621	surf_lsm_v(l)%building_covered(m) = .FALSE.
2622	!
2623	!-- Wall grid point is building-covered. This case, set
2624	!-- flag indicating that surface properties are initialized
2625	!-- from neighboring reference grid point, which is not
2626	!-- building_covered.
2627	IF ( building_type_f%from_file ) THEN
2628	i = surf_lsm_v(l)%i(m)
2629	j = surf_lsm_v(l)%j(m)
2630	IF ( building_type_f%var(j+surf_lsm_v(l)%joff, &
2631	i+surf_lsm_v(l)%ioff) /= &
2632	building_type_f%fill ) &
2633	surf_lsm_v(l)%building_covered(m) = .TRUE.
2634	!
2635	!-- Wall grid point as well as neighboring reference grid
2636	!-- point are both building-covered. This case, surface
2637	!-- properties are not necessarily defined (not covered by
2638	!-- checks for static input file) at this surface. Hence,
2639	!-- initialize surface properties by simply setting
2640	!-- vegetation_type_f to bare-soil bulk parametrization.
2641	!-- soil_type_f as well as surface_fractions_f will be set
2642	!-- also.
2643	IF ( building_type_f%var(j+surf_lsm_v(l)%joff, &
2644	i+surf_lsm_v(l)%ioff) /= &
2645	building_type_f%fill .AND. &
2646	building_type_f%var(j,i) /= building_type_f%fill ) &
2647	THEN
2648	vegetation_type_f%var(j,i) = 1 ! bare soil
2649	soil_type_f%var_2d(j,i) = 1
2650	!
2651	!-- Set surface_fraction if provided in static input,
2652	!-- else, in case no tiles are used, this will be done
2653	!-- on basis of the prescribed types (vegetation/pavement/
2654	!-- water_type).
2655	IF ( surface_fraction_f%from_file ) THEN
2656	surface_fraction_f%frac(ind_veg_wall,j,i) = 1.0_wp
2657	surface_fraction_f%frac(ind_pav_green,j,i) = 0.0_wp
2658	surface_fraction_f%frac(ind_wat_win,j,i) = 0.0_wp
2659	ENDIF
2660	ENDIF
2661
2662	ENDIF
2663	!
2664	!-- Normally proceed with setting surface types.
2665	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
2666	surf_lsm_v(l)%building_covered(m) )
2667	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
2668	surf_lsm_v(l)%building_covered(m) )
2669	IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
2670	surf_lsm_v(l)%vegetation_surface(m) = .TRUE.
2671	IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill ) &
2672	surf_lsm_v(l)%pavement_surface(m) = .TRUE.
2673	IF ( water_type_f%var(j,i) /= water_type_f%fill ) &
2674	surf_lsm_v(l)%water_surface(m) = .TRUE.
2675	!
2676	!-- Check if at least one type is set.
2677	IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m) .AND. &
2678	.NOT. surf_lsm_v(l)%pavement_surface(m) .AND. &
2679	.NOT. surf_lsm_v(l)%water_surface(m) ) THEN
2680	WRITE( message_string, * ) 'Vertical surface element ' //&
2681	' at i, j = ', i, j, &
2682	' is neither a vegetation, ' // &
2683	'pavement, nor a water surface.'
2684	CALL message( 'land_surface_model_mod', 'PA0619', &
2685	2, 2, myid, 6, 0 )
2686	ENDIF
2687	ENDDO
2688	ENDDO
2689
2690	END SELECT
2691	!
2692	!-- In case of netcdf input file, further initialize surface fractions.
2693	!-- At the moment only 1 surface is given at a location, so that the fraction
2694	!-- is either 0 or 1. This will be revised later. If surface fraction
2695	!-- is not given in static input file, relative fractions will be derived
2696	!-- from given surface type. In this case, only 1 type is given at a certain
2697	!-- location (already checked).
2698	IF ( input_pids_static .AND. surface_fraction_f%from_file ) THEN
2699	DO m = 1, surf_lsm_h%ns
2700	i = surf_lsm_h%i(m)
2701	j = surf_lsm_h%j(m)
2702	!
2703	!-- 0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction
2704	surf_lsm_h%frac(ind_veg_wall,m) = &
2705	surface_fraction_f%frac(ind_veg_wall,j,i)
2706	surf_lsm_h%frac(ind_pav_green,m) = &
2707	surface_fraction_f%frac(ind_pav_green,j,i)
2708	surf_lsm_h%frac(ind_wat_win,m) = &
2709	surface_fraction_f%frac(ind_wat_win,j,i)
2710
2711	ENDDO
2712	DO l = 0, 3
2713	DO m = 1, surf_lsm_v(l)%ns
2714	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
2715	surf_lsm_v(l)%building_covered(m) )
2716	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
2717	surf_lsm_v(l)%building_covered(m) )
2718	!
2719	!-- 0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction
2720	surf_lsm_v(l)%frac(ind_veg_wall,m) = &
2721	surface_fraction_f%frac(ind_veg_wall,j,i)
2722	surf_lsm_v(l)%frac(ind_pav_green,m) = &
2723	surface_fraction_f%frac(ind_pav_green,j,i)
2724	surf_lsm_v(l)%frac(ind_wat_win,m) = &
2725	surface_fraction_f%frac(ind_wat_win,j,i)
2726
2727	ENDDO
2728	ENDDO
2729	ELSEIF ( input_pids_static .AND. .NOT. surface_fraction_f%from_file ) &
2730	THEN
2731	DO m = 1, surf_lsm_h%ns
2732	i = surf_lsm_h%i(m)
2733	j = surf_lsm_h%j(m)
2734
2735	IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
2736	surf_lsm_h%frac(ind_veg_wall,m) = 1.0_wp
2737	IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill ) &
2738	surf_lsm_h%frac(ind_pav_green,m) = 1.0_wp
2739	IF ( water_type_f%var(j,i) /= water_type_f%fill ) &
2740	surf_lsm_h%frac(ind_wat_win,m) = 1.0_wp
2741	ENDDO
2742	DO l = 0, 3
2743	DO m = 1, surf_lsm_v(l)%ns
2744	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
2745	surf_lsm_v(l)%building_covered(m) )
2746	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
2747	surf_lsm_v(l)%building_covered(m) )
2748
2749	IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
2750	surf_lsm_v(l)%frac(ind_veg_wall,m) = 1.0_wp
2751	IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill ) &
2752	surf_lsm_v(l)%frac(ind_pav_green,m) = 1.0_wp
2753	IF ( water_type_f%var(j,i) /= water_type_f%fill ) &
2754	surf_lsm_v(l)%frac(ind_wat_win,m) = 1.0_wp
2755	ENDDO
2756	ENDDO
2757	ENDIF
2758	!
2759	!-- Level 1, initialization of soil parameters.
2760	!-- It is possible to overwrite each parameter by setting the respecticy
2761	!-- NAMELIST variable to a value /= 9999999.9.
2762	IF ( soil_type /= 0 ) THEN
2763
2764	IF ( alpha_vangenuchten == 9999999.9_wp ) THEN
2765	alpha_vangenuchten = soil_pars(0,soil_type)
2766	ENDIF
2767
2768	IF ( l_vangenuchten == 9999999.9_wp ) THEN
2769	l_vangenuchten = soil_pars(1,soil_type)
2770	ENDIF
2771
2772	IF ( n_vangenuchten == 9999999.9_wp ) THEN
2773	n_vangenuchten = soil_pars(2,soil_type)
2774	ENDIF
2775
2776	IF ( hydraulic_conductivity == 9999999.9_wp ) THEN
2777	hydraulic_conductivity = soil_pars(3,soil_type)
2778	ENDIF
2779
2780	IF ( saturation_moisture == 9999999.9_wp ) THEN
2781	saturation_moisture = soil_pars(4,soil_type)
2782	ENDIF
2783
2784	IF ( field_capacity == 9999999.9_wp ) THEN
2785	field_capacity = soil_pars(5,soil_type)
2786	ENDIF
2787
2788	IF ( wilting_point == 9999999.9_wp ) THEN
2789	wilting_point = soil_pars(6,soil_type)
2790	ENDIF
2791
2792	IF ( residual_moisture == 9999999.9_wp ) THEN
2793	residual_moisture = soil_pars(7,soil_type)
2794	ENDIF
2795
2796	ENDIF
2797	!
2798	!-- Map values to the respective 2D/3D arrays
2799	!-- Horizontal surfaces
2800	surf_lsm_h%alpha_vg = alpha_vangenuchten
2801	surf_lsm_h%l_vg = l_vangenuchten
2802	surf_lsm_h%n_vg = n_vangenuchten
2803	surf_lsm_h%gamma_w_sat = hydraulic_conductivity
2804	surf_lsm_h%m_sat = saturation_moisture
2805	surf_lsm_h%m_fc = field_capacity
2806	surf_lsm_h%m_wilt = wilting_point
2807	surf_lsm_h%m_res = residual_moisture
2808	surf_lsm_h%r_soil_min = min_soil_resistance
2809	!
2810	!-- Vertical surfaces
2811	DO l = 0, 3
2812	surf_lsm_v(l)%alpha_vg = alpha_vangenuchten
2813	surf_lsm_v(l)%l_vg = l_vangenuchten
2814	surf_lsm_v(l)%n_vg = n_vangenuchten
2815	surf_lsm_v(l)%gamma_w_sat = hydraulic_conductivity
2816	surf_lsm_v(l)%m_sat = saturation_moisture
2817	surf_lsm_v(l)%m_fc = field_capacity
2818	surf_lsm_v(l)%m_wilt = wilting_point
2819	surf_lsm_v(l)%m_res = residual_moisture
2820	surf_lsm_v(l)%r_soil_min = min_soil_resistance
2821	ENDDO
2822	!
2823	!-- Level 2, initialization of soil parameters via soil_type read from file.
2824	!-- Soil parameters are initialized for each (y,x)-grid point
2825	!-- individually using default paramter settings according to the given
2826	!-- soil type.
2827	IF ( soil_type_f%from_file ) THEN
2828	!
2829	!-- Level of detail = 1, i.e. a homogeneous soil distribution along the
2830	!-- vertical dimension is assumed.
2831	IF ( soil_type_f%lod == 1 ) THEN
2832	!
2833	!-- Horizontal surfaces
2834	DO m = 1, surf_lsm_h%ns
2835	i = surf_lsm_h%i(m)
2836	j = surf_lsm_h%j(m)
2837
2838	st = soil_type_f%var_2d(j,i)
2839	IF ( st /= soil_type_f%fill ) THEN
2840	surf_lsm_h%alpha_vg(:,m) = soil_pars(0,st)
2841	surf_lsm_h%l_vg(:,m) = soil_pars(1,st)
2842	surf_lsm_h%n_vg(:,m) = soil_pars(2,st)
2843	surf_lsm_h%gamma_w_sat(:,m) = soil_pars(3,st)
2844	surf_lsm_h%m_sat(:,m) = soil_pars(4,st)
2845	surf_lsm_h%m_fc(:,m) = soil_pars(5,st)
2846	surf_lsm_h%m_wilt(:,m) = soil_pars(6,st)
2847	surf_lsm_h%m_res(:,m) = soil_pars(7,st)
2848	ENDIF
2849	ENDDO
2850	!
2851	!-- Vertical surfaces ( assumes the soil type given at respective (x,y)
2852	DO l = 0, 3
2853	DO m = 1, surf_lsm_v(l)%ns
2854	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
2855	surf_lsm_v(l)%building_covered(m) )
2856	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
2857	surf_lsm_v(l)%building_covered(m) )
2858
2859	st = soil_type_f%var_2d(j,i)
2860	IF ( st /= soil_type_f%fill ) THEN
2861	surf_lsm_v(l)%alpha_vg(:,m) = soil_pars(0,st)
2862	surf_lsm_v(l)%l_vg(:,m) = soil_pars(1,st)
2863	surf_lsm_v(l)%n_vg(:,m) = soil_pars(2,st)
2864	surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars(3,st)
2865	surf_lsm_v(l)%m_sat(:,m) = soil_pars(4,st)
2866	surf_lsm_v(l)%m_fc(:,m) = soil_pars(5,st)
2867	surf_lsm_v(l)%m_wilt(:,m) = soil_pars(6,st)
2868	surf_lsm_v(l)%m_res(:,m) = soil_pars(7,st)
2869	ENDIF
2870	ENDDO
2871	ENDDO
2872	!
2873	!-- Level of detail = 2, i.e. soil type and thus the soil parameters
2874	!-- can be heterogeneous along the vertical dimension.
2875	ELSE
2876	!
2877	!-- Horizontal surfaces
2878	DO m = 1, surf_lsm_h%ns
2879	i = surf_lsm_h%i(m)
2880	j = surf_lsm_h%j(m)
2881
2882	DO k = nzb_soil, nzt_soil
2883	st = soil_type_f%var_3d(k,j,i)
2884	IF ( st /= soil_type_f%fill ) THEN
2885	surf_lsm_h%alpha_vg(k,m) = soil_pars(0,st)
2886	surf_lsm_h%l_vg(k,m) = soil_pars(1,st)
2887	surf_lsm_h%n_vg(k,m) = soil_pars(2,st)
2888	surf_lsm_h%gamma_w_sat(k,m) = soil_pars(3,st)
2889	surf_lsm_h%m_sat(k,m) = soil_pars(4,st)
2890	surf_lsm_h%m_fc(k,m) = soil_pars(5,st)
2891	surf_lsm_h%m_wilt(k,m) = soil_pars(6,st)
2892	surf_lsm_h%m_res(k,m) = soil_pars(7,st)
2893	ENDIF
2894	ENDDO
2895	ENDDO
2896	!
2897	!-- Vertical surfaces ( assumes the soil type given at respective (x,y)
2898	DO l = 0, 3
2899	DO m = 1, surf_lsm_v(l)%ns
2900	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
2901	surf_lsm_v(l)%building_covered(m) )
2902	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
2903	surf_lsm_v(l)%building_covered(m) )
2904
2905	DO k = nzb_soil, nzt_soil
2906	st = soil_type_f%var_3d(k,j,i)
2907	IF ( st /= soil_type_f%fill ) THEN
2908	surf_lsm_v(l)%alpha_vg(k,m) = soil_pars(0,st)
2909	surf_lsm_v(l)%l_vg(k,m) = soil_pars(1,st)
2910	surf_lsm_v(l)%n_vg(k,m) = soil_pars(2,st)
2911	surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars(3,st)
2912	surf_lsm_v(l)%m_sat(k,m) = soil_pars(4,st)
2913	surf_lsm_v(l)%m_fc(k,m) = soil_pars(5,st)
2914	surf_lsm_v(l)%m_wilt(k,m) = soil_pars(6,st)
2915	surf_lsm_v(l)%m_res(k,m) = soil_pars(7,st)
2916	ENDIF
2917	ENDDO
2918	ENDDO
2919	ENDDO
2920	ENDIF
2921	ENDIF
2922	!
2923	!-- Level 3, initialization of single soil parameters at single z,x,y
2924	!-- position via soil_pars read from file.
2925	IF ( soil_pars_f%from_file ) THEN
2926	!
2927	!-- Level of detail = 1, i.e. a homogeneous vertical distribution of soil
2928	!-- parameters is assumed.
2929	!-- Horizontal surfaces
2930	IF ( soil_pars_f%lod == 1 ) THEN
2931	!
2932	!-- Horizontal surfaces
2933	DO m = 1, surf_lsm_h%ns
2934	i = surf_lsm_h%i(m)
2935	j = surf_lsm_h%j(m)
2936
2937	IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill ) &
2938	surf_lsm_h%alpha_vg(:,m) = soil_pars_f%pars_xy(0,j,i)
2939	IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill ) &
2940	surf_lsm_h%l_vg(:,m) = soil_pars_f%pars_xy(1,j,i)
2941	IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill ) &
2942	surf_lsm_h%n_vg(:,m) = soil_pars_f%pars_xy(2,j,i)
2943	IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill ) &
2944	surf_lsm_h%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
2945	IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill ) &
2946	surf_lsm_h%m_sat(:,m) = soil_pars_f%pars_xy(4,j,i)
2947	IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill ) &
2948	surf_lsm_h%m_fc(:,m) = soil_pars_f%pars_xy(5,j,i)
2949	IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill ) &
2950	surf_lsm_h%m_wilt(:,m) = soil_pars_f%pars_xy(6,j,i)
2951	IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill ) &
2952	surf_lsm_h%m_res(:,m) = soil_pars_f%pars_xy(7,j,i)
2953
2954	ENDDO
2955	!
2956	!-- Vertical surfaces
2957	DO l = 0, 3
2958	DO m = 1, surf_lsm_v(l)%ns
2959	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
2960	surf_lsm_v(l)%building_covered(m) )
2961	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
2962	surf_lsm_v(l)%building_covered(m) )
2963
2964	IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill ) &
2965	surf_lsm_v(l)%alpha_vg(:,m) = soil_pars_f%pars_xy(0,j,i)
2966	IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill ) &
2967	surf_lsm_v(l)%l_vg(:,m) = soil_pars_f%pars_xy(1,j,i)
2968	IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill ) &
2969	surf_lsm_v(l)%n_vg(:,m) = soil_pars_f%pars_xy(2,j,i)
2970	IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill ) &
2971	surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
2972	IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill ) &
2973	surf_lsm_v(l)%m_sat(:,m) = soil_pars_f%pars_xy(4,j,i)
2974	IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill ) &
2975	surf_lsm_v(l)%m_fc(:,m) = soil_pars_f%pars_xy(5,j,i)
2976	IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill ) &
2977	surf_lsm_v(l)%m_wilt(:,m) = soil_pars_f%pars_xy(6,j,i)
2978	IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill ) &
2979	surf_lsm_v(l)%m_res(:,m) = soil_pars_f%pars_xy(7,j,i)
2980
2981	ENDDO
2982	ENDDO
2983	!
2984	!-- Level of detail = 2, i.e. soil parameters can be set at each soil
2985	!-- layer individually.
2986	ELSE
2987	!
2988	!-- Horizontal surfaces
2989	DO m = 1, surf_lsm_h%ns
2990	i = surf_lsm_h%i(m)
2991	j = surf_lsm_h%j(m)
2992
2993	DO k = nzb_soil, nzt_soil
2994	IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill ) &
2995	surf_lsm_h%alpha_vg(k,m) = soil_pars_f%pars_xyz(0,k,j,i)
2996	IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill ) &
2997	surf_lsm_h%l_vg(k,m) = soil_pars_f%pars_xyz(1,k,j,i)
2998	IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill ) &
2999	surf_lsm_h%n_vg(k,m) = soil_pars_f%pars_xyz(2,k,j,i)
3000	IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill ) &
3001	surf_lsm_h%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
3002	IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill ) &
3003	surf_lsm_h%m_sat(k,m) = soil_pars_f%pars_xyz(4,k,j,i)
3004	IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill ) &
3005	surf_lsm_h%m_fc(k,m) = soil_pars_f%pars_xyz(5,k,j,i)
3006	IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill ) &
3007	surf_lsm_h%m_wilt(k,m) = soil_pars_f%pars_xyz(6,k,j,i)
3008	IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill ) &
3009	surf_lsm_h%m_res(k,m) = soil_pars_f%pars_xyz(7,k,j,i)
3010	ENDDO
3011
3012	ENDDO
3013	!
3014	!-- Vertical surfaces
3015	DO l = 0, 3
3016	DO m = 1, surf_lsm_v(l)%ns
3017	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
3018	surf_lsm_v(l)%building_covered(m) )
3019	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
3020	surf_lsm_v(l)%building_covered(m) )
3021
3022	DO k = nzb_soil, nzt_soil
3023	IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill ) &
3024	surf_lsm_v(l)%alpha_vg(k,m) = soil_pars_f%pars_xyz(0,k,j,i)
3025	IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill ) &
3026	surf_lsm_v(l)%l_vg(k,m) = soil_pars_f%pars_xyz(1,k,j,i)
3027	IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill ) &
3028	surf_lsm_v(l)%n_vg(k,m) = soil_pars_f%pars_xyz(2,k,j,i)
3029	IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill ) &
3030	surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
3031	IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill ) &
3032	surf_lsm_v(l)%m_sat(k,m) = soil_pars_f%pars_xyz(4,k,j,i)
3033	IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill ) &
3034	surf_lsm_v(l)%m_fc(k,m) = soil_pars_f%pars_xyz(5,k,j,i)
3035	IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill ) &
3036	surf_lsm_v(l)%m_wilt(k,m) = soil_pars_f%pars_xyz(6,k,j,i)
3037	IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill ) &
3038	surf_lsm_v(l)%m_res(k,m) = soil_pars_f%pars_xyz(7,k,j,i)
3039	ENDDO
3040
3041	ENDDO
3042	ENDDO
3043
3044	ENDIF
3045	ENDIF
3046
3047	!
3048	!-- Level 1, initialization of vegetation parameters. A horizontally
3049	!-- homogeneous distribution is assumed here.
3050	IF ( vegetation_type /= 0 ) THEN
3051
3052	IF ( min_canopy_resistance == 9999999.9_wp ) THEN
3053	min_canopy_resistance = vegetation_pars(ind_v_rc_min,vegetation_type)
3054	ENDIF
3055
3056	IF ( leaf_area_index == 9999999.9_wp ) THEN
3057	leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type)
3058	ENDIF
3059
3060	IF ( vegetation_coverage == 9999999.9_wp ) THEN
3061	vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type)
3062	ENDIF
3063
3064	IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN
3065	canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type)
3066	ENDIF
3067
3068	IF ( z0_vegetation == 9999999.9_wp ) THEN
3069	z0_vegetation = vegetation_pars(ind_v_z0,vegetation_type)
3070	ENDIF
3071
3072	IF ( z0h_vegetation == 9999999.9_wp ) THEN
3073	z0h_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
3074	ENDIF
3075
3076	IF ( z0q_vegetation == 9999999.9_wp ) THEN
3077	z0q_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
3078	ENDIF
3079
3080	IF ( lambda_surface_stable == 9999999.9_wp ) THEN
3081	lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type)
3082	ENDIF
3083
3084	IF ( lambda_surface_unstable == 9999999.9_wp ) THEN
3085	lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type)
3086	ENDIF
3087
3088	IF ( f_shortwave_incoming == 9999999.9_wp ) THEN
3089	f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type)
3090	ENDIF
3091
3092	IF ( c_surface == 9999999.9_wp ) THEN
3093	c_surface = vegetation_pars(ind_v_c_surf,vegetation_type)
3094	ENDIF
3095
3096	IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN
3097	albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type))
3098	ENDIF
3099
3100	IF ( emissivity == 9999999.9_wp ) THEN
3101	emissivity = vegetation_pars(ind_v_emis,vegetation_type)
3102	ENDIF
3103
3104	ENDIF
3105	!
3106	!-- Map values onto horizontal elemements
3107	DO m = 1, surf_lsm_h%ns
3108	IF ( surf_lsm_h%vegetation_surface(m) ) THEN
3109	surf_lsm_h%r_canopy_min(m) = min_canopy_resistance
3110	surf_lsm_h%lai(m) = leaf_area_index
3111	surf_lsm_h%c_veg(m) = vegetation_coverage
3112	surf_lsm_h%g_d(m) = canopy_resistance_coefficient
3113	surf_lsm_h%z0(m) = z0_vegetation
3114	surf_lsm_h%z0h(m) = z0h_vegetation
3115	surf_lsm_h%z0q(m) = z0q_vegetation
3116	surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
3117	surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
3118	surf_lsm_h%f_sw_in(m) = f_shortwave_incoming
3119	surf_lsm_h%c_surface(m) = c_surface
3120	surf_lsm_h%albedo_type(ind_veg_wall,m) = albedo_type
3121	surf_lsm_h%emissivity(ind_veg_wall,m) = emissivity
3122
3123	surf_lsm_h%vegetation_type(m) = vegetation_type
3124	surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
3125	ELSE
3126	surf_lsm_h%lai(m) = 0.0_wp
3127	surf_lsm_h%c_veg(m) = 0.0_wp
3128	surf_lsm_h%g_d(m) = 0.0_wp
3129	ENDIF
3130
3131	ENDDO
3132	!
3133	!-- Map values onto vertical elements, even though this does not make
3134	!-- much sense.
3135	DO l = 0, 3
3136	DO m = 1, surf_lsm_v(l)%ns
3137	IF ( surf_lsm_v(l)%vegetation_surface(m) ) THEN
3138	surf_lsm_v(l)%r_canopy_min(m) = min_canopy_resistance
3139	surf_lsm_v(l)%lai(m) = leaf_area_index
3140	surf_lsm_v(l)%c_veg(m) = vegetation_coverage
3141	surf_lsm_v(l)%g_d(m) = canopy_resistance_coefficient
3142	surf_lsm_v(l)%z0(m) = z0_vegetation
3143	surf_lsm_v(l)%z0h(m) = z0h_vegetation
3144	surf_lsm_v(l)%z0q(m) = z0q_vegetation
3145	surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
3146	surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
3147	surf_lsm_v(l)%f_sw_in(m) = f_shortwave_incoming
3148	surf_lsm_v(l)%c_surface(m) = c_surface
3149	surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = albedo_type
3150	surf_lsm_v(l)%emissivity(ind_veg_wall,m) = emissivity
3151
3152	surf_lsm_v(l)%vegetation_type(m) = vegetation_type
3153	surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
3154	ELSE
3155	surf_lsm_v(l)%lai(m) = 0.0_wp
3156	surf_lsm_v(l)%c_veg(m) = 0.0_wp
3157	surf_lsm_v(l)%g_d(m) = 0.0_wp
3158	ENDIF
3159	ENDDO
3160	ENDDO
3161
3162	!
3163	!-- Level 2, initialization of vegation parameters via vegetation_type read
3164	!-- from file. Vegetation parameters are initialized for each (y,x)-grid point
3165	!-- individually using default paramter settings according to the given
3166	!-- vegetation type.
3167	IF ( vegetation_type_f%from_file ) THEN
3168	!
3169	!-- Horizontal surfaces
3170	DO m = 1, surf_lsm_h%ns
3171	i = surf_lsm_h%i(m)
3172	j = surf_lsm_h%j(m)
3173
3174	st = vegetation_type_f%var(j,i)
3175	IF ( st /= vegetation_type_f%fill .AND. st /= 0 ) THEN
3176	surf_lsm_h%r_canopy_min(m) = vegetation_pars(ind_v_rc_min,st)
3177	surf_lsm_h%lai(m) = vegetation_pars(ind_v_rc_lai,st)
3178	surf_lsm_h%c_veg(m) = vegetation_pars(ind_v_c_veg,st)
3179	surf_lsm_h%g_d(m) = vegetation_pars(ind_v_gd,st)
3180	surf_lsm_h%z0(m) = vegetation_pars(ind_v_z0,st)
3181	surf_lsm_h%z0h(m) = vegetation_pars(ind_v_z0qh,st)
3182	surf_lsm_h%z0q(m) = vegetation_pars(ind_v_z0qh,st)
3183	surf_lsm_h%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
3184	surf_lsm_h%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
3185	surf_lsm_h%f_sw_in(m) = vegetation_pars(ind_v_f_sw_in,st)
3186	surf_lsm_h%c_surface(m) = vegetation_pars(ind_v_c_surf,st)
3187	surf_lsm_h%albedo_type(ind_veg_wall,m) = INT( vegetation_pars(ind_v_at,st) )
3188	surf_lsm_h%emissivity(ind_veg_wall,m) = vegetation_pars(ind_v_emis,st)
3189
3190	surf_lsm_h%vegetation_type(m) = st
3191	surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(st)
3192	ENDIF
3193	ENDDO
3194	!
3195	!-- Vertical surfaces
3196	DO l = 0, 3
3197	DO m = 1, surf_lsm_v(l)%ns
3198	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
3199	surf_lsm_v(l)%building_covered(m) )
3200	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
3201	surf_lsm_v(l)%building_covered(m) )
3202
3203	st = vegetation_type_f%var(j,i)
3204	IF ( st /= vegetation_type_f%fill .AND. st /= 0 ) THEN
3205	surf_lsm_v(l)%r_canopy_min(m) = vegetation_pars(ind_v_rc_min,st)
3206	surf_lsm_v(l)%lai(m) = vegetation_pars(ind_v_rc_lai,st)
3207	surf_lsm_v(l)%c_veg(m) = vegetation_pars(ind_v_c_veg,st)
3208	surf_lsm_v(l)%g_d(m) = vegetation_pars(ind_v_gd,st)
3209	surf_lsm_v(l)%z0(m) = vegetation_pars(ind_v_z0,st)
3210	surf_lsm_v(l)%z0h(m) = vegetation_pars(ind_v_z0qh,st)
3211	surf_lsm_v(l)%z0q(m) = vegetation_pars(ind_v_z0qh,st)
3212	surf_lsm_v(l)%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
3213	surf_lsm_v(l)%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
3214	surf_lsm_v(l)%f_sw_in(m) = vegetation_pars(ind_v_f_sw_in,st)
3215	surf_lsm_v(l)%c_surface(m) = vegetation_pars(ind_v_c_surf,st)
3216	surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = INT( vegetation_pars(ind_v_at,st) )
3217	surf_lsm_v(l)%emissivity(ind_veg_wall,m) = vegetation_pars(ind_v_emis,st)
3218
3219	surf_lsm_v(l)%vegetation_type(m) = st
3220	surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(st)
3221	ENDIF
3222	ENDDO
3223	ENDDO
3224	ENDIF
3225	!
3226	!-- Level 3, initialization of vegation parameters at single (x,y)
3227	!-- position via vegetation_pars read from file.
3228	IF ( vegetation_pars_f%from_file ) THEN
3229	!
3230	!-- Horizontal surfaces
3231	DO m = 1, surf_lsm_h%ns
3232
3233	i = surf_lsm_h%i(m)
3234	j = surf_lsm_h%j(m)
3235	!
3236	!-- If surface element is not a vegetation surface and any value in
3237	!-- vegetation_pars is given, neglect this information and give an
3238	!-- informative message that this value will not be used.
3239	IF ( .NOT. surf_lsm_h%vegetation_surface(m) .AND. &
3240	ANY( vegetation_pars_f%pars_xy(:,j,i) /= &
3241	vegetation_pars_f%fill ) ) THEN
3242	WRITE( message_string, * ) &
3243	'surface element at grid point (j,i) = (', &
3244	j, i, ') is not a vegetation surface, ', &
3245	'so that information given in ', &
3246	'vegetation_pars at this point is neglected.'
3247	CALL message( 'land_surface_model_mod', 'PA0436', 0, 0, myid, 6, 0 )
3248	ELSE
3249
3250	IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /= &
3251	vegetation_pars_f%fill ) &
3252	surf_lsm_h%r_canopy_min(m) = &
3253	vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
3254	IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /= &
3255	vegetation_pars_f%fill ) &
3256	surf_lsm_h%lai(m) = &
3257	vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
3258	IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /= &
3259	vegetation_pars_f%fill ) &
3260	surf_lsm_h%c_veg(m) = &
3261	vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
3262	IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /= &
3263	vegetation_pars_f%fill ) &
3264	surf_lsm_h%g_d(m) = &
3265	vegetation_pars_f%pars_xy(ind_v_gd,j,i)
3266	IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /= &
3267	vegetation_pars_f%fill ) &
3268	surf_lsm_h%z0(m) = &
3269	vegetation_pars_f%pars_xy(ind_v_z0,j,i)
3270	IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /= &
3271	vegetation_pars_f%fill ) THEN
3272	surf_lsm_h%z0h(m) = &
3273	vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3274	surf_lsm_h%z0q(m) = &
3275	vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3276	ENDIF
3277	IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /= &
3278	vegetation_pars_f%fill ) &
3279	surf_lsm_h%lambda_surface_s(m) = &
3280	vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
3281	IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /= &
3282	vegetation_pars_f%fill ) &
3283	surf_lsm_h%lambda_surface_u(m) = &
3284	vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
3285	IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /= &
3286	vegetation_pars_f%fill ) &
3287	surf_lsm_h%f_sw_in(m) = &
3288	vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
3289	IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /= &
3290	vegetation_pars_f%fill ) &
3291	surf_lsm_h%c_surface(m) = &
3292	vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
3293	IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /= &
3294	vegetation_pars_f%fill ) &
3295	surf_lsm_h%albedo_type(ind_veg_wall,m) = &
3296	INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
3297	IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /= &
3298	vegetation_pars_f%fill ) &
3299	surf_lsm_h%emissivity(ind_veg_wall,m) = &
3300	vegetation_pars_f%pars_xy(ind_v_emis,j,i)
3301	ENDIF
3302	ENDDO
3303	!
3304	!-- Vertical surfaces
3305	DO l = 0, 3
3306	DO m = 1, surf_lsm_v(l)%ns
3307	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
3308	surf_lsm_v(l)%building_covered(m) )
3309	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
3310	surf_lsm_v(l)%building_covered(m) )
3311	!
3312	!-- If surface element is not a vegetation surface and any value in
3313	!-- vegetation_pars is given, neglect this information and give an
3314	!-- informative message that this value will not be used.
3315	IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m) .AND. &
3316	ANY( vegetation_pars_f%pars_xy(:,j,i) /= &
3317	vegetation_pars_f%fill ) ) THEN
3318	WRITE( message_string, * ) &
3319	'surface element at grid point (j,i) = (', &
3320	j, i, ') is not a vegetation surface, ', &
3321	'so that information given in ', &
3322	'vegetation_pars at this point is neglected.'
3323	CALL message( 'land_surface_model_mod', 'PA0436', 0, 0, myid, 6, 0 )
3324	ELSE
3325
3326	IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /= &
3327	vegetation_pars_f%fill ) &
3328	surf_lsm_v(l)%r_canopy_min(m) = &
3329	vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
3330	IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /= &
3331	vegetation_pars_f%fill ) &
3332	surf_lsm_v(l)%lai(m) = &
3333	vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
3334	IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /= &
3335	vegetation_pars_f%fill ) &
3336	surf_lsm_v(l)%c_veg(m) = &
3337	vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
3338	IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /= &
3339	vegetation_pars_f%fill ) &
3340	surf_lsm_v(l)%g_d(m) = &
3341	vegetation_pars_f%pars_xy(ind_v_gd,j,i)
3342	IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /= &
3343	vegetation_pars_f%fill ) &
3344	surf_lsm_v(l)%z0(m) = &
3345	vegetation_pars_f%pars_xy(ind_v_z0,j,i)
3346	IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /= &
3347	vegetation_pars_f%fill ) THEN
3348	surf_lsm_v(l)%z0h(m) = &
3349	vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3350	surf_lsm_v(l)%z0q(m) = &
3351	vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3352	ENDIF
3353	IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /= &
3354	vegetation_pars_f%fill ) &
3355	surf_lsm_v(l)%lambda_surface_s(m) = &
3356	vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
3357	IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /= &
3358	vegetation_pars_f%fill ) &
3359	surf_lsm_v(l)%lambda_surface_u(m) = &
3360	vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
3361	IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /= &
3362	vegetation_pars_f%fill ) &
3363	surf_lsm_v(l)%f_sw_in(m) = &
3364	vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
3365	IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /= &
3366	vegetation_pars_f%fill ) &
3367	surf_lsm_v(l)%c_surface(m) = &
3368	vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
3369	IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /= &
3370	vegetation_pars_f%fill ) &
3371	surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = &
3372	INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
3373	IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /= &
3374	vegetation_pars_f%fill ) &
3375	surf_lsm_v(l)%emissivity(ind_veg_wall,m) = &
3376	vegetation_pars_f%pars_xy(ind_v_emis,j,i)
3377	ENDIF
3378
3379	ENDDO
3380	ENDDO
3381	ENDIF
3382
3383	!
3384	!-- Level 1, initialization of water parameters. A horizontally
3385	!-- homogeneous distribution is assumed here.
3386	IF ( water_type /= 0 ) THEN
3387
3388	IF ( water_temperature == 9999999.9_wp ) THEN
3389	water_temperature = water_pars(ind_w_temp,water_type)
3390	ENDIF
3391
3392	IF ( z0_water == 9999999.9_wp ) THEN
3393	z0_water = water_pars(ind_w_z0,water_type)
3394	ENDIF
3395
3396	IF ( z0h_water == 9999999.9_wp ) THEN
3397	z0h_water = water_pars(ind_w_z0h,water_type)
3398	ENDIF
3399
3400	IF ( z0q_water == 9999999.9_wp ) THEN
3401	z0q_water = water_pars(ind_w_z0h,water_type)
3402	ENDIF
3403
3404	IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN
3405	albedo_type = INT(water_pars(ind_w_at,water_type))
3406	ENDIF
3407
3408	IF ( emissivity == 9999999.9_wp ) THEN
3409	emissivity = water_pars(ind_w_emis,water_type)
3410	ENDIF
3411
3412	ENDIF
3413	!
3414	!-- Map values onto horizontal elemements
3415	DO m = 1, surf_lsm_h%ns
3416	IF ( surf_lsm_h%water_surface(m) ) THEN
3417	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) &
3418	t_soil_h%var_2d(:,m) = water_temperature
3419	surf_lsm_h%z0(m) = z0_water
3420	surf_lsm_h%z0h(m) = z0h_water
3421	surf_lsm_h%z0q(m) = z0q_water
3422	surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
3423	surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp
3424	surf_lsm_h%c_surface(m) = 0.0_wp
3425	surf_lsm_h%albedo_type(ind_wat_win,m) = albedo_type
3426	surf_lsm_h%emissivity(ind_wat_win,m) = emissivity
3427
3428	surf_lsm_h%water_type(m) = water_type
3429	surf_lsm_h%water_type_name(m) = water_type_name(water_type)
3430	ENDIF
3431	ENDDO
3432	!
3433	!-- Map values onto vertical elements, even though this does not make
3434	!-- much sense.
3435	DO l = 0, 3
3436	DO m = 1, surf_lsm_v(l)%ns
3437	IF ( surf_lsm_v(l)%water_surface(m) ) THEN
3438	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) &
3439	t_soil_v(l)%var_2d(:,m) = water_temperature
3440	surf_lsm_v(l)%z0(m) = z0_water
3441	surf_lsm_v(l)%z0h(m) = z0h_water
3442	surf_lsm_v(l)%z0q(m) = z0q_water
3443	surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
3444	surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp
3445	surf_lsm_v(l)%c_surface(m) = 0.0_wp
3446	surf_lsm_v(l)%albedo_type(ind_wat_win,m) = albedo_type
3447	surf_lsm_v(l)%emissivity(ind_wat_win,m) = emissivity
3448
3449	surf_lsm_v(l)%water_type(m) = water_type
3450	surf_lsm_v(l)%water_type_name(m) = water_type_name(water_type)
3451	ENDIF
3452	ENDDO
3453	ENDDO
3454	!
3455	!
3456	!-- Level 2, initialization of water parameters via water_type read
3457	!-- from file. Water surfaces are initialized for each (y,x)-grid point
3458	!-- individually using default paramter settings according to the given
3459	!-- water type.
3460	!-- Note, parameter 3/4 of water_pars are albedo and emissivity,
3461	!-- whereas paramter 3/4 of water_pars_f are heat conductivities!
3462	IF ( water_type_f%from_file ) THEN
3463	!
3464	!-- Horizontal surfaces
3465	DO m = 1, surf_lsm_h%ns
3466	i = surf_lsm_h%i(m)
3467	j = surf_lsm_h%j(m)
3468
3469	st = water_type_f%var(j,i)
3470	IF ( st /= water_type_f%fill .AND. st /= 0 ) THEN
3471	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) &
3472	t_soil_h%var_2d(:,m) = water_pars(ind_w_temp,st)
3473	surf_lsm_h%z0(m) = water_pars(ind_w_z0,st)
3474	surf_lsm_h%z0h(m) = water_pars(ind_w_z0h,st)
3475	surf_lsm_h%z0q(m) = water_pars(ind_w_z0h,st)
3476	surf_lsm_h%lambda_surface_s(m) = water_pars(ind_w_lambda_s,st)
3477	surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st)
3478	surf_lsm_h%c_surface(m) = 0.0_wp
3479	surf_lsm_h%albedo_type(ind_wat_win,m) = INT( water_pars(ind_w_at,st) )
3480	surf_lsm_h%emissivity(ind_wat_win,m) = water_pars(ind_w_emis,st)
3481
3482	surf_lsm_h%water_type(m) = st
3483	surf_lsm_h%water_type_name(m) = water_type_name(st)
3484	ENDIF
3485	ENDDO
3486	!
3487	!-- Vertical surfaces
3488	DO l = 0, 3
3489	DO m = 1, surf_lsm_v(l)%ns
3490	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
3491	surf_lsm_v(l)%building_covered(m) )
3492	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
3493	surf_lsm_v(l)%building_covered(m) )
3494
3495	st = water_type_f%var(j,i)
3496	IF ( st /= water_type_f%fill .AND. st /= 0 ) THEN
3497	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) &
3498	t_soil_v(l)%var_2d(:,m) = water_pars(ind_w_temp,st)
3499	surf_lsm_v(l)%z0(m) = water_pars(ind_w_z0,st)
3500	surf_lsm_v(l)%z0h(m) = water_pars(ind_w_z0h,st)
3501	surf_lsm_v(l)%z0q(m) = water_pars(ind_w_z0h,st)
3502	surf_lsm_v(l)%lambda_surface_s(m) = &
3503	water_pars(ind_w_lambda_s,st)
3504	surf_lsm_v(l)%lambda_surface_u(m) = &
3505	water_pars(ind_w_lambda_u,st)
3506	surf_lsm_v(l)%c_surface(m) = 0.0_wp
3507	surf_lsm_v(l)%albedo_type(ind_wat_win,m) = &
3508	INT( water_pars(ind_w_at,st) )
3509	surf_lsm_v(l)%emissivity(ind_wat_win,m) = &
3510	water_pars(ind_w_emis,st)
3511
3512	surf_lsm_v(l)%water_type(m) = st
3513	surf_lsm_v(l)%water_type_name(m) = water_type_name(st)
3514	ENDIF
3515	ENDDO
3516	ENDDO
3517	ENDIF
3518
3519	!
3520	!-- Level 3, initialization of water parameters at single (x,y)
3521	!-- position via water_pars read from file.
3522	IF ( water_pars_f%from_file ) THEN
3523	!
3524	!-- Horizontal surfaces
3525	DO m = 1, surf_lsm_h%ns
3526	i = surf_lsm_h%i(m)
3527	j = surf_lsm_h%j(m)
3528	!
3529	!-- If surface element is not a water surface and any value in
3530	!-- water_pars is given, neglect this information and give an
3531	!-- informative message that this value will not be used.
3532	IF ( .NOT. surf_lsm_h%water_surface(m) .AND. &
3533	ANY( water_pars_f%pars_xy(:,j,i) /= water_pars_f%fill ) ) THEN
3534	WRITE( message_string, * ) &
3535	'surface element at grid point (j,i) = (', &
3536	j, i, ') is not a water surface, ', &
3537	'so that information given in ', &
3538	'water_pars at this point is neglected.'
3539	CALL message( 'land_surface_model_mod', 'PA0645', 0, 0, myid, 6, 0 )
3540	ELSE
3541	IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /= &
3542	water_pars_f%fill .AND. &
3543	TRIM( initializing_actions ) /= 'read_restart_data' ) &
3544	t_soil_h%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
3545
3546	IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= water_pars_f%fill ) &
3547	surf_lsm_h%z0(m) = water_pars_f%pars_xy(ind_w_z0,j,i)
3548
3549	IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= water_pars_f%fill )&
3550	THEN
3551	surf_lsm_h%z0h(m) = water_pars_f%pars_xy(ind_w_z0h,j,i)
3552	surf_lsm_h%z0q(m) = water_pars_f%pars_xy(ind_w_z0h,j,i)
3553	ENDIF
3554	IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /= &
3555	water_pars_f%fill ) &
3556	surf_lsm_h%lambda_surface_s(m) = &
3557	water_pars_f%pars_xy(ind_w_lambda_s,j,i)
3558
3559	IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /= &
3560	water_pars_f%fill ) &
3561	surf_lsm_h%lambda_surface_u(m) = &
3562	water_pars_f%pars_xy(ind_w_lambda_u,j,i)
3563
3564	IF ( water_pars_f%pars_xy(ind_w_at,j,i) /= &
3565	water_pars_f%fill ) &
3566	surf_lsm_h%albedo_type(ind_wat_win,m) = &
3567	INT( water_pars_f%pars_xy(ind_w_at,j,i) )
3568
3569	IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /= &
3570	water_pars_f%fill ) &
3571	surf_lsm_h%emissivity(ind_wat_win,m) = &
3572	water_pars_f%pars_xy(ind_w_emis,j,i)
3573	ENDIF
3574	ENDDO
3575	!
3576	!-- Vertical surfaces
3577	DO l = 0, 3
3578	DO m = 1, surf_lsm_v(l)%ns
3579	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
3580	surf_lsm_v(l)%building_covered(m) )
3581	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
3582	surf_lsm_v(l)%building_covered(m) )
3583	!
3584	!-- If surface element is not a water surface and any value in
3585	!-- water_pars is given, neglect this information and give an
3586	!-- informative message that this value will not be used.
3587	IF ( .NOT. surf_lsm_v(l)%water_surface(m) .AND. &
3588	ANY( water_pars_f%pars_xy(:,j,i) /= &
3589	water_pars_f%fill ) ) THEN
3590	WRITE( message_string, * ) &
3591	'surface element at grid point (j,i) = (', &
3592	j, i, ') is not a water surface, ', &
3593	'so that information given in ', &
3594	'water_pars at this point is neglected.'
3595	CALL message( 'land_surface_model_mod', 'PA0645', &
3596	0, 0, myid, 6, 0 )
3597	ELSE
3598
3599	IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /= &
3600	water_pars_f%fill .AND. &
3601	TRIM( initializing_actions ) /= 'read_restart_data' ) &
3602	t_soil_v(l)%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
3603
3604	IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= &
3605	water_pars_f%fill ) &
3606	surf_lsm_v(l)%z0(m) = water_pars_f%pars_xy(ind_w_z0,j,i)
3607
3608	IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= &
3609	water_pars_f%fill ) THEN
3610	surf_lsm_v(l)%z0h(m) = water_pars_f%pars_xy(ind_w_z0h,j,i)
3611	surf_lsm_v(l)%z0q(m) = water_pars_f%pars_xy(ind_w_z0h,j,i)
3612	ENDIF
3613
3614	IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /= &
3615	water_pars_f%fill ) &
3616	surf_lsm_v(l)%lambda_surface_s(m) = &
3617	water_pars_f%pars_xy(ind_w_lambda_s,j,i)
3618
3619	IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /= &
3620	water_pars_f%fill ) &
3621	surf_lsm_v(l)%lambda_surface_u(m) = &
3622	water_pars_f%pars_xy(ind_w_lambda_u,j,i)
3623
3624	IF ( water_pars_f%pars_xy(ind_w_at,j,i) /= &
3625	water_pars_f%fill ) &
3626	surf_lsm_v(l)%albedo_type(ind_wat_win,m) = &
3627	INT( water_pars_f%pars_xy(ind_w_at,j,i) )
3628
3629	IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /= &
3630	water_pars_f%fill ) &
3631	surf_lsm_v(l)%emissivity(ind_wat_win,m) = &
3632	water_pars_f%pars_xy(ind_w_emis,j,i)
3633	ENDIF
3634	ENDDO
3635	ENDDO
3636
3637	ENDIF
3638	!
3639	!-- Initialize pavement-type surfaces, level 1
3640	IF ( pavement_type /= 0 ) THEN
3641
3642	!
3643	!-- When a pavement_type is used, overwrite a possible setting of
3644	!-- the pavement depth as it is already defined by the pavement type
3645	pavement_depth_level = 0
3646
3647	IF ( z0_pavement == 9999999.9_wp ) THEN
3648	z0_pavement = pavement_pars(ind_p_z0,pavement_type)
3649	ENDIF
3650
3651	IF ( z0h_pavement == 9999999.9_wp ) THEN
3652	z0h_pavement = pavement_pars(ind_p_z0h,pavement_type)
3653	ENDIF
3654
3655	IF ( z0q_pavement == 9999999.9_wp ) THEN
3656	z0q_pavement = pavement_pars(ind_p_z0h,pavement_type)
3657	ENDIF
3658
3659	IF ( pavement_heat_conduct == 9999999.9_wp ) THEN
3660	pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
3661	ENDIF
3662
3663	IF ( pavement_heat_capacity == 9999999.9_wp ) THEN
3664	pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
3665	ENDIF
3666
3667	IF ( albedo_type == 9999999 .AND. albedo == 9999999.9_wp ) THEN
3668	albedo_type = INT(pavement_pars(ind_p_at,pavement_type))
3669	ENDIF
3670
3671	IF ( emissivity == 9999999.9_wp ) THEN
3672	emissivity = pavement_pars(ind_p_emis,pavement_type)
3673	ENDIF
3674
3675	!
3676	!-- If the depth level of the pavement is not set, determine it from
3677	!-- lookup table.
3678	IF ( pavement_depth_level == 0 ) THEN
3679	DO k = nzb_soil, nzt_soil
3680	IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &
3681	.OR. pavement_subsurface_pars_2(k,pavement_type) == 9999999.9_wp)&
3682	THEN
3683	nzt_pavement = k-1
3684	EXIT
3685	ENDIF
3686	ENDDO
3687	ELSE
3688	nzt_pavement = pavement_depth_level
3689	ENDIF
3690
3691	ENDIF
3692	!
3693	!-- Level 1 initialization of pavement type surfaces. Horizontally
3694	!-- homogeneous characteristics are assumed
3695	surf_lsm_h%nzt_pavement = pavement_depth_level
3696	DO m = 1, surf_lsm_h%ns
3697	IF ( surf_lsm_h%pavement_surface(m) ) THEN
3698	surf_lsm_h%nzt_pavement(m) = nzt_pavement
3699	surf_lsm_h%z0(m) = z0_pavement
3700	surf_lsm_h%z0h(m) = z0h_pavement
3701	surf_lsm_h%z0q(m) = z0q_pavement
3702	surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct &
3703	* ddz_soil(nzb_soil) &
3704	* 2.0_wp
3705	surf_lsm_h%lambda_surface_u(m) = pavement_heat_conduct &
3706	* ddz_soil(nzb_soil) &
3707	* 2.0_wp
3708	surf_lsm_h%c_surface(m) = pavement_heat_capacity &
3709	* dz_soil(nzb_soil) &
3710	* 0.25_wp
3711
3712	surf_lsm_h%albedo_type(ind_pav_green,m) = albedo_type
3713	surf_lsm_h%emissivity(ind_pav_green,m) = emissivity
3714
3715	surf_lsm_h%pavement_type(m) = pavement_type
3716	surf_lsm_h%pavement_type_name(m) = pavement_type_name(pavement_type)
3717
3718	IF ( pavement_type /= 0 ) THEN
3719	DO k = nzb_soil, surf_lsm_h%nzt_pavement(m)
3720	surf_lsm_h%lambda_h_def(k,m) = &
3721	pavement_subsurface_pars_1(k,pavement_type)
3722	surf_lsm_h%rho_c_total_def(k,m) = &
3723	pavement_subsurface_pars_2(k,pavement_type)
3724	ENDDO
3725	ELSE
3726	surf_lsm_h%lambda_h_def(:,m) = pavement_heat_conduct
3727	surf_lsm_h%rho_c_total_def(:,m) = pavement_heat_capacity
3728	ENDIF
3729	ENDIF
3730	ENDDO
3731
3732	DO l = 0, 3
3733	surf_lsm_v(l)%nzt_pavement = pavement_depth_level
3734	DO m = 1, surf_lsm_v(l)%ns
3735	IF ( surf_lsm_v(l)%pavement_surface(m) ) THEN
3736	surf_lsm_v(l)%nzt_pavement(m) = nzt_pavement
3737	surf_lsm_v(l)%z0(m) = z0_pavement
3738	surf_lsm_v(l)%z0h(m) = z0h_pavement
3739	surf_lsm_v(l)%z0q(m) = z0q_pavement
3740	surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct &
3741	* ddz_soil(nzb_soil) &
3742	* 2.0_wp
3743	surf_lsm_v(l)%lambda_surface_u(m) = pavement_heat_conduct &
3744	* ddz_soil(nzb_soil) &
3745	* 2.0_wp
3746	surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity &
3747	* dz_soil(nzb_soil) &
3748	* 0.25_wp
3749
3750	surf_lsm_v(l)%albedo_type(ind_pav_green,m) = albedo_type
3751	surf_lsm_v(l)%emissivity(ind_pav_green,m) = emissivity
3752
3753	surf_lsm_v(l)%pavement_type(m) = pavement_type
3754	surf_lsm_v(l)%pavement_type_name(m) = pavement_type_name(pavement_type)
3755
3756	IF ( pavement_type /= 0 ) THEN
3757	DO k = nzb_soil, surf_lsm_v(l)%nzt_pavement(m)
3758	surf_lsm_v(l)%lambda_h_def(k,m) = &
3759	pavement_subsurface_pars_1(k,pavement_type)
3760	surf_lsm_v(l)%rho_c_total_def(k,m) = &
3761	pavement_subsurface_pars_2(k,pavement_type)
3762	ENDDO
3763	ELSE
3764	surf_lsm_v(l)%lambda_h_def(:,m) = pavement_heat_conduct
3765	surf_lsm_v(l)%rho_c_total_def(:,m) = pavement_heat_capacity
3766	ENDIF
3767	ENDIF
3768	ENDDO
3769	ENDDO
3770	!
3771	!-- Level 2 initialization of pavement type surfaces via pavement_type read
3772	!-- from file. Pavement surfaces are initialized for each (y,x)-grid point
3773	!-- individually.
3774	IF ( pavement_type_f%from_file ) THEN
3775	!
3776	!-- Horizontal surfaces
3777	DO m = 1, surf_lsm_h%ns
3778	i = surf_lsm_h%i(m)
3779	j = surf_lsm_h%j(m)
3780
3781	st = pavement_type_f%var(j,i)
3782	IF ( st /= pavement_type_f%fill .AND. st /= 0 ) THEN
3783	!
3784	!-- Determine deepmost index of pavement layer
3785	DO k = nzb_soil, nzt_soil
3786	IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp &
3787	.OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp) &
3788	THEN
3789	surf_lsm_h%nzt_pavement(m) = k-1
3790	EXIT
3791	ENDIF
3792	ENDDO
3793
3794	surf_lsm_h%z0(m) = pavement_pars(ind_p_z0,st)
3795	surf_lsm_h%z0h(m) = pavement_pars(ind_p_z0h,st)
3796	surf_lsm_h%z0q(m) = pavement_pars(ind_p_z0h,st)
3797
3798	surf_lsm_h%lambda_surface_s(m) = &
3799	pavement_subsurface_pars_1(0,st) &
3800	* ddz_soil(nzb_soil) &
3801	* 2.0_wp
3802	surf_lsm_h%lambda_surface_u(m) = &
3803	pavement_subsurface_pars_1(0,st) &
3804	* ddz_soil(nzb_soil) &
3805	* 2.0_wp
3806	surf_lsm_h%c_surface(m) = &
3807	pavement_subsurface_pars_2(0,st)&
3808	* dz_soil(nzb_soil) &
3809	* 0.25_wp
3810	surf_lsm_h%albedo_type(ind_pav_green,m) = INT( pavement_pars(ind_p_at,st) )
3811	surf_lsm_h%emissivity(ind_pav_green,m) = pavement_pars(ind_p_emis,st)
3812
3813	surf_lsm_h%pavement_type(m) = st
3814	surf_lsm_h%pavement_type_name(m) = pavement_type_name(st)
3815
3816	DO k = nzb_soil, surf_lsm_h%nzt_pavement(m)
3817	surf_lsm_h%lambda_h_def(k,m) = &
3818	pavement_subsurface_pars_1(k,pavement_type)
3819	surf_lsm_h%rho_c_total_def(k,m) = &
3820	pavement_subsurface_pars_2(k,pavement_type)
3821	ENDDO
3822	ENDIF
3823	ENDDO
3824	!
3825	!-- Vertical surfaces
3826	DO l = 0, 3
3827	DO m = 1, surf_lsm_v(l)%ns
3828	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
3829	surf_lsm_v(l)%building_covered(m) )
3830	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
3831	surf_lsm_v(l)%building_covered(m) )
3832
3833	st = pavement_type_f%var(j,i)
3834	IF ( st /= pavement_type_f%fill .AND. st /= 0 ) THEN
3835	!
3836	!-- Determine deepmost index of pavement layer
3837	DO k = nzb_soil, nzt_soil
3838	IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp &
3839	.OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp) &
3840	THEN
3841	surf_lsm_v(l)%nzt_pavement(m) = k-1
3842	EXIT
3843	ENDIF
3844	ENDDO
3845
3846	surf_lsm_v(l)%z0(m) = pavement_pars(ind_p_z0,st)
3847	surf_lsm_v(l)%z0h(m) = pavement_pars(ind_p_z0h,st)
3848	surf_lsm_v(l)%z0q(m) = pavement_pars(ind_p_z0h,st)
3849
3850	surf_lsm_v(l)%lambda_surface_s(m) = &
3851	pavement_subsurface_pars_1(0,st) &
3852	* ddz_soil(nzb_soil) &
3853	* 2.0_wp
3854	surf_lsm_v(l)%lambda_surface_u(m) = &
3855	pavement_subsurface_pars_1(0,st) &
3856	* ddz_soil(nzb_soil) &
3857	* 2.0_wp
3858
3859	surf_lsm_v(l)%c_surface(m) = &
3860	pavement_subsurface_pars_2(0,st) &
3861	* dz_soil(nzb_soil) &
3862	* 0.25_wp
3863	surf_lsm_v(l)%albedo_type(ind_pav_green,m) = &
3864	INT( pavement_pars(ind_p_at,st) )
3865	surf_lsm_v(l)%emissivity(ind_pav_green,m) = &
3866	pavement_pars(ind_p_emis,st)
3867
3868	surf_lsm_v(l)%pavement_type(m) = st
3869	surf_lsm_v(l)%pavement_type_name(m) = pavement_type_name(st)
3870
3871	DO k = nzb_soil, surf_lsm_v(l)%nzt_pavement(m)
3872	surf_lsm_v(l)%lambda_h_def(k,m) = &
3873	pavement_subsurface_pars_1(k,pavement_type)
3874	surf_lsm_v(l)%rho_c_total_def(k,m) = &
3875	pavement_subsurface_pars_2(k,pavement_type)
3876	ENDDO
3877	ENDIF
3878	ENDDO
3879	ENDDO
3880	ENDIF
3881	!
3882	!-- Level 3, initialization of pavement parameters at single (x,y)
3883	!-- position via pavement_pars read from file.
3884	IF ( pavement_pars_f%from_file ) THEN
3885	!
3886	!-- Horizontal surfaces
3887	DO m = 1, surf_lsm_h%ns
3888	i = surf_lsm_h%i(m)
3889	j = surf_lsm_h%j(m)
3890	!
3891	!-- If surface element is not a pavement surface and any value in
3892	!-- pavement_pars is given, neglect this information and give an
3893	!-- informative message that this value will not be used.
3894	IF ( .NOT. surf_lsm_h%pavement_surface(m) .AND. &
3895	ANY( pavement_pars_f%pars_xy(:,j,i) /= &
3896	pavement_pars_f%fill ) ) THEN
3897	WRITE( message_string, * ) &
3898	'surface element at grid point (j,i) = (', &
3899	j, i, ') is not a pavement surface, ', &
3900	'so that information given in ', &
3901	'pavement_pars at this point is neglected.'
3902	CALL message( 'land_surface_model_mod', 'PA0647', 0, 0, myid, 6, 0 )
3903	ELSE
3904	IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /= &
3905	pavement_pars_f%fill ) &
3906	surf_lsm_h%z0(m) = pavement_pars_f%pars_xy(ind_p_z0,j,i)
3907	IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /= &
3908	pavement_pars_f%fill ) THEN
3909	surf_lsm_h%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
3910	surf_lsm_h%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
3911	ENDIF
3912	IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i) &
3913	/= pavement_subsurface_pars_f%fill ) THEN
3914	surf_lsm_h%lambda_surface_s(m) = &
3915	pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3916	* ddz_soil(nzb_soil) &
3917	* 2.0_wp
3918	surf_lsm_h%lambda_surface_u(m) = &
3919	pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3920	* ddz_soil(nzb_soil) &
3921	* 2.0_wp
3922	ENDIF
3923	IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) /= &
3924	pavement_subsurface_pars_f%fill ) THEN
3925	surf_lsm_h%c_surface(m) = &
3926	pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) &
3927	* dz_soil(nzb_soil) &
3928	* 0.25_wp
3929	ENDIF
3930	IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /= &
3931	pavement_pars_f%fill ) &
3932	surf_lsm_h%albedo_type(ind_pav_green,m) = &
3933	INT( pavement_pars_f%pars_xy(ind_p_at,j,i) )
3934	IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /= &
3935	pavement_pars_f%fill ) &
3936	surf_lsm_h%emissivity(ind_pav_green,m) = &
3937	pavement_pars_f%pars_xy(ind_p_emis,j,i)
3938	ENDIF
3939
3940	ENDDO
3941	!
3942	!-- Vertical surfaces
3943	DO l = 0, 3
3944	DO m = 1, surf_lsm_v(l)%ns
3945	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
3946	surf_lsm_v(l)%building_covered(m) )
3947	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
3948	surf_lsm_v(l)%building_covered(m) )
3949	!
3950	!-- If surface element is not a pavement surface and any value in
3951	!-- pavement_pars is given, neglect this information and give an
3952	!-- informative message that this value will not be used.
3953	IF ( .NOT. surf_lsm_v(l)%pavement_surface(m) .AND. &
3954	ANY( pavement_pars_f%pars_xy(:,j,i) /= &
3955	pavement_pars_f%fill ) ) THEN
3956	WRITE( message_string, * ) &
3957	'surface element at grid point (j,i) = (', &
3958	j, i, ') is not a pavement surface, ', &
3959	'so that information given in ', &
3960	'pavement_pars at this point is neglected.'
3961	CALL message( 'land_surface_model_mod', 'PA0647', 0, 0, myid, 6, 0 )
3962	ELSE
3963
3964	IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /= &
3965	pavement_pars_f%fill ) &
3966	surf_lsm_v(l)%z0(m) = pavement_pars_f%pars_xy(ind_p_z0,j,i)
3967	IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /= &
3968	pavement_pars_f%fill ) THEN
3969	surf_lsm_v(l)%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
3970	surf_lsm_v(l)%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
3971	ENDIF
3972	IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3973	/= pavement_subsurface_pars_f%fill ) THEN
3974	surf_lsm_v(l)%lambda_surface_s(m) = &
3975	pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3976	* ddz_soil(nzb_soil) &
3977	* 2.0_wp
3978	surf_lsm_v(l)%lambda_surface_u(m) = &
3979	pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3980	* ddz_soil(nzb_soil) &
3981	* 2.0_wp
3982	ENDIF
3983	IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) &
3984	/= pavement_subsurface_pars_f%fill ) THEN
3985	surf_lsm_v(l)%c_surface(m) = &
3986	pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)&
3987	* dz_soil(nzb_soil) &
3988	* 0.25_wp
3989	ENDIF
3990	IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /= &
3991	pavement_pars_f%fill ) &
3992	surf_lsm_v(l)%albedo_type(ind_pav_green,m) = &
3993	INT( pavement_pars_f%pars_xy(ind_p_at,j,i) )
3994
3995	IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /= &
3996	pavement_pars_f%fill ) &
3997	surf_lsm_v(l)%emissivity(ind_pav_green,m) = &
3998	pavement_pars_f%pars_xy(ind_p_emis,j,i)
3999	ENDIF
4000	ENDDO
4001	ENDDO
4002	ENDIF
4003	!
4004	!-- Moreover, for grid points which are flagged with pavement-type 0 or whre
4005	!-- pavement_subsurface_pars_f is provided, soil heat conductivity and
4006	!-- capacity are initialized with parameters given in
4007	!-- pavement_subsurface_pars read from file.
4008	IF ( pavement_subsurface_pars_f%from_file ) THEN
4009	!
4010	!-- Set pavement depth to nzt_soil. Please note, this is just a
4011	!-- workaround at the moment.
4012	DO m = 1, surf_lsm_h%ns
4013	IF ( surf_lsm_h%pavement_surface(m) ) THEN
4014
4015	i = surf_lsm_h%i(m)
4016	j = surf_lsm_h%j(m)
4017
4018	surf_lsm_h%nzt_pavement(m) = nzt_soil
4019
4020	DO k = nzb_soil, nzt_soil
4021	surf_lsm_h%lambda_h_def(k,m) = &
4022	pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i)
4023	surf_lsm_h%rho_c_total_def(k,m) = &
4024	pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i)
4025	ENDDO
4026
4027	ENDIF
4028	ENDDO
4029	DO l = 0, 3
4030	DO m = 1, surf_lsm_v(l)%ns
4031	IF ( surf_lsm_v(l)%pavement_surface(m) ) THEN
4032
4033	i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff, &
4034	surf_lsm_v(l)%building_covered(m) )
4035	j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff, &
4036	surf_lsm_v(l)%building_covered(m) )
4037
4038	surf_lsm_v(l)%nzt_pavement(m) = nzt_soil
4039
4040	DO k = nzb_soil, nzt_soil
4041	surf_lsm_v(l)%lambda_h_def(k,m) = &
4042	pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i)
4043	surf_lsm_v(l)%rho_c_total_def(k,m) = &
4044	pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i)
4045	ENDDO
4046
4047	ENDIF
4048	ENDDO
4049	ENDDO
4050	ENDIF
4051
4052	!
4053	!-- Initial run actions
4054	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
4055	!
4056	!-- First, initialize soil temperature and moisture.
4057	!-- According to the initialization for surface and soil parameters,
4058	!-- initialize soil moisture and temperature via a level approach. This
4059	!-- is to assure that all surface elements are initialized, even if
4060	!-- data provided from input file contains fill values at some locations.
4061	!-- Level 1, initialization via profiles given in parameter file
4062	DO m = 1, surf_lsm_h%ns
4063	IF ( surf_lsm_h%vegetation_surface(m) .OR. &
4064	surf_lsm_h%pavement_surface(m) ) THEN
4065	DO k = nzb_soil, nzt_soil
4066	t_soil_h%var_2d(k,m) = soil_temperature(k)
4067	m_soil_h%var_2d(k,m) = soil_moisture(k)
4068	ENDDO
4069	t_soil_h%var_2d(nzt_soil+1,m) = deep_soil_temperature
4070	ENDIF
4071	ENDDO
4072	DO l = 0, 3
4073	DO m = 1, surf_lsm_v(l)%ns
4074	IF ( surf_lsm_v(l)%vegetation_surface(m) .OR. &
4075	surf_lsm_v(l)%pavement_surface(m) ) THEN
4076	DO k = nzb_soil, nzt_soil
4077	t_soil_v(l)%var_2d(k,m) = soil_temperature(k)
4078	m_soil_v(l)%var_2d(k,m) = soil_moisture(k)
4079	ENDDO
4080	t_soil_v(l)%var_2d(nzt_soil+1,m) = deep_soil_temperature
4081	ENDIF
4082	ENDDO
4083	ENDDO
4084	!
4085	!-- Level 2 initialization of the soil, read soil properties from
4086	!-- dynamic input file.
4087	IF ( input_pids_dynamic ) THEN
4088	!
4089	!-- CPU measurement
4090	CALL cpu_log( log_point_s(85), 'NetCDF input init', 'start' )
4091	#if defined ( __netcdf )
4092	!
4093	!-- Open file in read-only mode
4094	CALL open_read_file( TRIM( input_file_dynamic ) // &
4095	TRIM( coupling_char ), pids_id )
4096	!
4097	!-- Inquire all variable names
4098	CALL inquire_num_variables( pids_id, num_var_pids )
4099	!
4100	!-- Allocate memory to store variable names.
4101	ALLOCATE( vars_pids(1:num_var_pids) )
4102	CALL inquire_variable_names( pids_id, vars_pids )
4103	!
4104	!-- Read vertical dimension for soil depth.
4105	IF ( check_existence( vars_pids, 'zsoil' ) ) &
4106	CALL get_dimension_length( pids_id, init_3d%nzs, 'zsoil' )
4107	!
4108	!-- Read also the horizontal dimensions required for soil initialization.
4109	!-- Please note, in case of non-nested runs or in case of root domain,
4110	!-- these data is already available, but will be read again for the sake
4111	!-- of clearness.
4112	CALL get_dimension_length( pids_id, init_3d%nx, 'x' )
4113	CALL get_dimension_length( pids_id, init_3d%ny, 'y' )
4114	!
4115	!-- Check for correct horizontal and vertical dimension. Please note,
4116	!-- in case of non-nested runs or in case of root domain, these checks
4117	!-- are already performed
4118	IF ( init_3d%nx-1 /= nx .OR. init_3d%ny-1 /= ny ) THEN
4119	message_string = 'Number of horizontal grid points in '// &
4120	'dynamic input file does not match ' // &
4121	'the number of numeric grid points.'
4122	CALL message( 'lsm_init', 'PA0543', 1, 2, 0, 6, 0 )
4123	ENDIF
4124	!
4125	!-- Read vertical dimensions. Later, these are required for eventual
4126	!-- inter- and extrapolations of the initialization data.
4127	IF ( check_existence( vars_pids, 'zsoil' ) ) THEN
4128	ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
4129	CALL get_variable( pids_id, 'zsoil', init_3d%z_soil )
4130	ENDIF
4131	!
4132	!-- Read initial data for soil moisture
4133	IF ( check_existence( vars_pids, 'init_soil_m' ) ) THEN
4134	!
4135	!-- Read attributes for the fill value and level-of-detail
4136	CALL get_attribute( pids_id, char_fill, &
4137	init_3d%fill_msoil, &
4138	.FALSE., 'init_soil_m' )
4139	CALL get_attribute( pids_id, char_lod, &
4140	init_3d%lod_msoil, &
4141	.FALSE., 'init_soil_m' )
4142	!
4143	!-- level-of-detail 1 - read initialization profile
4144	IF ( init_3d%lod_msoil == 1 ) THEN
4145	ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
4146
4147	CALL get_variable( pids_id, 'init_soil_m', &
4148	init_3d%msoil_1d(0:init_3d%nzs-1) )
4149	!
4150	!-- level-of-detail 2 - read 3D initialization data
4151	ELSEIF ( init_3d%lod_msoil == 2 ) THEN
4152	ALLOCATE ( init_3d%msoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
4153
4154	CALL get_variable( pids_id, 'init_soil_m', &
4155	init_3d%msoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
4156	nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
4157
4158	ENDIF
4159	init_3d%from_file_msoil = .TRUE.
4160	ENDIF
4161	!
4162	!-- Read soil temperature
4163	IF ( check_existence( vars_pids, 'init_soil_t' ) ) THEN
4164	!
4165	!-- Read attributes for the fill value and level-of-detail
4166	CALL get_attribute( pids_id, char_fill, &
4167	init_3d%fill_tsoil, &
4168	.FALSE., 'init_soil_t' )
4169	CALL get_attribute( pids_id, char_lod, &
4170	init_3d%lod_tsoil, &
4171	.FALSE., 'init_soil_t' )
4172	!
4173	!-- level-of-detail 1 - read initialization profile
4174	IF ( init_3d%lod_tsoil == 1 ) THEN
4175	ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
4176
4177	CALL get_variable( pids_id, 'init_soil_t', &
4178	init_3d%tsoil_1d(0:init_3d%nzs-1) )
4179
4180	!
4181	!-- level-of-detail 2 - read 3D initialization data
4182	ELSEIF ( init_3d%lod_tsoil == 2 ) THEN
4183	ALLOCATE ( init_3d%tsoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
4184
4185	CALL get_variable( pids_id, 'init_soil_t', &
4186	init_3d%tsoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
4187	nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
4188	ENDIF
4189	init_3d%from_file_tsoil = .TRUE.
4190	ENDIF
4191	!
4192	!-- Close input file
4193	CALL close_input_file( pids_id )
4194	#endif
4195	!
4196	!-- End of CPU measurement
4197	CALL cpu_log( log_point_s(85), 'NetCDF input init', 'stop' )
4198	ENDIF
4199	!
4200	!-- In case no dynamic input is available for a child domain but the
4201	!-- parent is initialized with dynamic input file, the different soil
4202	!-- states can lead to significant discrepancies in the atmospheric
4203	!-- surface forcing. For this reason, the child domain is initialized with
4204	!-- domain-averaged soil profiles from the root domain, even if
4205	!-- no initialization with inifor is set. Note, as long as a dynamic
4206	!-- input file with soil information is available for the child domain,
4207	!-- the input file information will be used.
4208	IF ( nested_run ) THEN
4209	#if defined( __parallel )
4210	!
4211	!-- In case of a nested run, first average the soil profiles in the
4212	!-- root domain.
4213	IF ( pmc_is_rootmodel() ) THEN
4214	!
4215	!-- Child domains will be only initialized with horizontally
4216	!-- averaged soil profiles in parent domain (for sake of simplicity).
4217	!-- If required, average soil data on root parent domain before the
4218	!-- soil profiles are distributed onto the child domains.
4219	!-- Start with soil moisture.
4220	IF ( init_3d%from_file_msoil .AND. init_3d%lod_msoil == 2 ) &
4221	THEN
4222	ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
4223	DO k = 0, init_3d%nzs-1
4224	pr_soil_init(k) = SUM( init_3d%msoil_3d(k,nys:nyn,nxl:nxr) )
4225	ENDDO
4226	!
4227	!-- Allocate 1D array for soil-moisture profile (will not be
4228	!-- allocated in lod==2 case).
4229	ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
4230	init_3d%msoil_1d = 0.0_wp
4231	CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%msoil_1d(0), &
4232	SIZE(pr_soil_init), &
4233	MPI_REAL, MPI_SUM, comm2d, ierr )
4234
4235	init_3d%msoil_1d = init_3d%msoil_1d / &
4236	REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
4237	DEALLOCATE( pr_soil_init )
4238	ENDIF
4239	!
4240	!-- Proceed with soil temperature.
4241	IF ( init_3d%from_file_tsoil .AND. init_3d%lod_tsoil == 2 ) &
4242	THEN
4243	ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
4244
4245	DO k = 0, init_3d%nzs-1
4246	pr_soil_init(k) = SUM( init_3d%tsoil_3d(k,nys:nyn,nxl:nxr) )
4247	ENDDO
4248	!
4249	!-- Allocate 1D array for soil-temperature profile (will not be
4250	!-- allocated in lod==2 case).
4251	ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
4252	init_3d%tsoil_1d = 0.0_wp
4253	CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%tsoil_1d(0), &
4254	SIZE(pr_soil_init), &
4255	MPI_REAL, MPI_SUM, comm2d, ierr )
4256	init_3d%tsoil_1d = init_3d%tsoil_1d / &
4257	REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
4258	DEALLOCATE( pr_soil_init )
4259
4260	ENDIF
4261	ENDIF
4262	!
4263	!-- Broadcast number of soil layers in root model to all childs.
4264	!-- Note, only process 0 in COMM_WORLD is sending.
4265	IF ( pmc_is_rootmodel() ) nzs_root = init_3d%nzs
4266
4267	CALL MPI_BCAST( nzs_root, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
4268	!
4269	!-- Allocate dummy arrays for soil moisture and temperature profiles
4270	!-- on all domains.
4271	ALLOCATE( z_soil_root(1:nzs_root) )
4272	ALLOCATE( m_soil_root(0:nzs_root-1) )
4273	ALLOCATE( t_soil_root(0:nzs_root-1) )
4274	!
4275	!-- Distribute the mean soil profiles to all child domains.
4276	IF ( pmc_is_rootmodel() ) THEN
4277	z_soil_root = init_3d%z_soil
4278	m_soil_root = init_3d%msoil_1d
4279	t_soil_root = init_3d%tsoil_1d
4280	ENDIF
4281
4282	CALL MPI_BCAST( z_soil_root, SIZE( z_soil_root ), &
4283	MPI_REAL, 0, MPI_COMM_WORLD, ierr )
4284	CALL MPI_BCAST( m_soil_root, SIZE( m_soil_root ), &
4285	MPI_REAL, 0, MPI_COMM_WORLD, ierr )
4286	CALL MPI_BCAST( t_soil_root, SIZE( t_soil_root ), &
4287	MPI_REAL, 0, MPI_COMM_WORLD, ierr )
4288	!
4289	!-- In the following, the child domains decide whether they take
4290	!-- the information from the root domain or not.
4291	IF ( .NOT. pmc_is_rootmodel() ) THEN
4292	!
4293	!-- If soil moisture or temperature isn't in dynamic input file for
4294	!-- the child, take the information provided from the root model.
4295	!-- Start with z-dimension
4296	IF ( .NOT. init_3d%from_file_msoil .OR. &
4297	.NOT. init_3d%from_file_msoil ) THEN
4298	init_3d%nzs = nzs_root
4299	ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
4300	init_3d%z_soil(1:init_3d%nzs) = z_soil_root
4301	ENDIF
4302	!
4303	!-- Take soil moisture. Note, control flags from_file... and LoD
4304	!-- need to be set.
4305	IF ( .NOT. init_3d%from_file_msoil ) THEN
4306	ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
4307	init_3d%lod_msoil = 1
4308	init_3d%from_file_msoil = .TRUE.
4309
4310	init_3d%msoil_1d = m_soil_root
4311	ENDIF
4312	!
4313	!-- Take soil temperature. Note, control flags from_file... and LoD
4314	!-- need to be set.
4315	IF ( .NOT. init_3d%from_file_tsoil ) THEN
4316	ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
4317	init_3d%lod_tsoil = 1
4318	init_3d%from_file_tsoil = .TRUE.
4319
4320	init_3d%tsoil_1d = t_soil_root
4321	ENDIF
4322	ENDIF
4323
4324	DEALLOCATE( z_soil_root )
4325	DEALLOCATE( m_soil_root )
4326	DEALLOCATE( t_soil_root )
4327	ENDIF
4328	#endif
4329	!
4330	!-- Proceed with Level 2 initialization.
4331	IF ( init_3d%from_file_msoil ) THEN
4332
4333	IF ( init_3d%lod_msoil == 1 ) THEN
4334	DO m = 1, surf_lsm_h%ns
4335	IF ( surf_lsm_h%vegetation_surface(m) .OR. &
4336	surf_lsm_h%pavement_surface(m) ) THEN
4337
4338	CALL interpolate_soil_profile( &
4339	m_soil_h%var_2d(nzb_soil:nzt_soil,m), &
4340	init_3d%msoil_1d(:), &
4341	zs(nzb_soil:nzt_soil), init_3d%z_soil, &
4342	nzb_soil, nzt_soil, &
4343	nzb_soil, init_3d%nzs-1 )
4344	ENDIF
4345	ENDDO
4346	DO l = 0, 3
4347	DO m = 1, surf_lsm_v(l)%ns
4348	IF ( surf_lsm_v(l)%vegetation_surface(m) .OR. &
4349	surf_lsm_v(l)%pavement_surface(m) ) THEN
4350	CALL interpolate_soil_profile( &
4351	m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
4352	init_3d%msoil_1d(:), &
4353	zs(nzb_soil:nzt_soil), init_3d%z_soil, &
4354	nzb_soil, nzt_soil, &
4355	nzb_soil, init_3d%nzs-1 )
4356	ENDIF
4357	ENDDO
4358	ENDDO
4359	ELSE
4360
4361	DO m = 1, surf_lsm_h%ns
4362	IF ( surf_lsm_h%vegetation_surface(m) .OR. &
4363	surf_lsm_h%pavement_surface(m) ) THEN
4364	i = surf_lsm_h%i(m)
4365	j = surf_lsm_h%j(m)
4366
4367	IF ( init_3d%msoil_3d(0,j,i) /= init_3d%fill_msoil ) &
4368	CALL interpolate_soil_profile( &
4369	m_soil_h%var_2d(nzb_soil:nzt_soil,m), &
4370	init_3d%msoil_3d(:,j,i), &
4371	zs(nzb_soil:nzt_soil), init_3d%z_soil, &
4372	nzb_soil, nzt_soil, &
4373	nzb_soil, init_3d%nzs-1 )
4374	ENDIF
4375	ENDDO
4376	DO l = 0, 3
4377	DO m = 1, surf_lsm_v(l)%ns
4378	IF ( surf_lsm_v(l)%vegetation_surface(m) .OR. &
4379	surf_lsm_v(l)%pavement_surface(m) ) THEN
4380	!
4381	!-- Note, in contrast to the static input data the dynamic
4382	!-- input do not need to be checked whether a grid point
4383	!-- is building covered. This is because soil data in the
4384	!-- dynamic input is provided for the whole domain.
4385	i = surf_lsm_v(l)%i(m)
4386	j = surf_lsm_v(l)%j(m)
4387
4388	IF ( init_3d%msoil_3d(0,j,i) /= init_3d%fill_msoil ) &
4389	CALL interpolate_soil_profile( &
4390	m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
4391	init_3d%msoil_3d(:,j,i), &
4392	zs(nzb_soil:nzt_soil), init_3d%z_soil, &
4393	nzb_soil, nzt_soil, &
4394	nzb_soil, init_3d%nzs-1 )
4395	ENDIF
4396	ENDDO
4397	ENDDO
4398	ENDIF
4399	ENDIF
4400	!
4401	!-- Soil temperature
4402	IF ( init_3d%from_file_tsoil ) THEN
4403
4404	IF ( init_3d%lod_tsoil == 1 ) THEN ! change to 1 if provided correctly by INIFOR
4405	DO m = 1, surf_lsm_h%ns
4406	IF ( surf_lsm_h%vegetation_surface(m) .OR. &
4407	surf_lsm_h%pavement_surface(m) ) THEN
4408	CALL interpolate_soil_profile( &
4409	t_soil_h%var_2d(nzb_soil:nzt_soil,m), &
4410	init_3d%tsoil_1d(:), &
4411	zs(nzb_soil:nzt_soil), init_3d%z_soil, &
4412	nzb_soil, nzt_soil, &
4413	nzb_soil, init_3d%nzs-1 )
4414	!
4415	!-- Set boundary condition, i.e. deep soil temperature
4416	t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m)
4417	ENDIF
4418	ENDDO
4419	DO l = 0, 3
4420	DO m = 1, surf_lsm_v(l)%ns
4421	IF ( surf_lsm_v(l)%vegetation_surface(m) .OR. &
4422	surf_lsm_v(l)%pavement_surface(m) ) THEN
4423	CALL interpolate_soil_profile( &
4424	t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
4425	init_3d%tsoil_1d(:), &
4426	zs(nzb_soil:nzt_soil), init_3d%z_soil, &
4427	nzb_soil, nzt_soil, &
4428	nzb_soil, init_3d%nzs-1 )
4429	!
4430	!-- Set boundary condition, i.e. deep soil temperature
4431	t_soil_v(l)%var_2d(nzt_soil+1,m) = &
4432	t_soil_v(l)%var_2d(nzt_soil,m)
4433	ENDIF
4434	ENDDO
4435	ENDDO
4436	ELSE
4437
4438	DO m = 1, surf_lsm_h%ns
4439	IF ( surf_lsm_h%vegetation_surface(m) .OR. &
4440	surf_lsm_h%pavement_surface(m) ) THEN
4441	i = surf_lsm_h%i(m)
4442	j = surf_lsm_h%j(m)
4443
4444	IF ( init_3d%tsoil_3d(0,j,i) /= init_3d%fill_tsoil ) &
4445	CALL interpolate_soil_profile( &
4446	t_soil_h%var_2d(nzb_soil:nzt_soil,m), &
4447	init_3d%tsoil_3d(:,j,i), &
4448	zs(nzb_soil:nzt_soil), init_3d%z_soil, &
4449	nzb_soil, nzt_soil, &
4450	nzb_soil, init_3d%nzs-1 )
4451	!
4452	!-- Set boundary condition, i.e. deep soil temperature
4453	t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m)
4454	ENDIF
4455	ENDDO
4456	DO l = 0, 3
4457	DO m = 1, surf_lsm_v(l)%ns
4458	IF ( surf_lsm_v(l)%vegetation_surface(m) .OR. &
4459	surf_lsm_v(l)%pavement_surface(m) ) THEN
4460	!
4461	!-- Note, in contrast to the static input data the dynamic
4462	!-- input do not need to be checked whether a grid point
4463	!-- is building covered. This is because soil data in the
4464	!-- dynamic input is provided for the whole domain.
4465	i = surf_lsm_v(l)%i(m)
4466	j = surf_lsm_v(l)%j(m)
4467
4468	IF ( init_3d%tsoil_3d(0,j,i) /= init_3d%fill_tsoil ) &
4469	CALL interpolate_soil_profile( &
4470	t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
4471	init_3d%tsoil_3d(:,j,i), &
4472	zs(nzb_soil:nzt_soil), init_3d%z_soil, &
4473	nzb_soil, nzt_soil, &
4474	nzb_soil, init_3d%nzs-1 )
4475	!
4476	!-- Set boundary condition, i.e. deep soil temperature
4477	t_soil_v(l)%var_2d(nzt_soil+1,m) = &
4478	t_soil_v(l)%var_2d(nzt_soil,m)
4479	ENDIF
4480	ENDDO
4481	ENDDO
4482	ENDIF
4483	ENDIF
4484	!
4485	!-- After soil moisture and temperature are finally initialized, check
4486	!-- if soil moisture is higher than its saturation value. If this would
4487	!-- be the case, the soil model parametrization will produce floating
4488	!-- point errors. Hence, limit the soil moisture to its saturation value
4489	!-- and give a warning.
4490	DO m = 1, surf_lsm_h%ns
4491	IF ( surf_lsm_h%vegetation_surface(m) .OR. &
4492	surf_lsm_h%pavement_surface(m) ) THEN
4493	DO k = nzb_soil, nzt_soil
4494	IF ( m_soil_h%var_2d(k,m) > surf_lsm_h%m_sat(k,m) ) THEN
4495	m_soil_h%var_2d(k,m) = surf_lsm_h%m_sat(k,m)
4496	WRITE( message_string, * ) 'soil moisture is higher '// &
4497	'than its saturation value at (k,j,i) ', k, &
4498	surf_lsm_h%i(m), surf_lsm_h%j(m), ' and is ' // &
4499	'thus limited to this value to maintain stability.'
4500	CALL message( 'lsm_init', 'PA0458', 0, 1, myid, 6, 0 )
4501	ENDIF
4502	ENDDO
4503	ENDIF
4504	ENDDO
4505	DO l = 0, 3
4506	DO m = 1, surf_lsm_v(l)%ns
4507	IF ( surf_lsm_v(l)%vegetation_surface(m) .OR. &
4508	surf_lsm_v(l)%pavement_surface(m) ) THEN
4509	DO k = nzb_soil, nzt_soil
4510	IF ( m_soil_v(l)%var_2d(k,m) > surf_lsm_v(l)%m_sat(k,m) )&
4511	THEN
4512	m_soil_v(l)%var_2d(k,m) = surf_lsm_v(l)%m_sat(k,m)
4513	WRITE( message_string, * ) &
4514	'soil moisture is higher '// &
4515	'than its saturation value at (k,j,i) ', k, &
4516	surf_lsm_v(l)%i(m), surf_lsm_v(l)%j(m), &
4517	' and is ' // &
4518	'thus limited to this value to maintain stability.'
4519	CALL message( 'lsm_init', 'PA0458', 0, 1, myid, 6, 0 )
4520	ENDIF
4521	ENDDO
4522	ENDIF
4523	ENDDO
4524	ENDDO
4525
4526	!
4527	!-- Further initialization
4528	DO m = 1, surf_lsm_h%ns
4529
4530	i = surf_lsm_h%i(m)
4531	j = surf_lsm_h%j(m)
4532	k = surf_lsm_h%k(m)
4533	!
4534	!-- Initialize surface temperature with soil temperature in the uppermost
4535	!-- uppermost layer
4536	t_surface_h%var_1d(m) = t_soil_h%var_2d(nzb_soil,m)
4537	surf_lsm_h%pt_surface(m) = t_soil_h%var_2d(nzb_soil,m) / exner(nzb)
4538
4539	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
4540	surf_lsm_h%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
4541	ELSE
4542	surf_lsm_h%pt1(m) = pt(k,j,i)
4543	ENDIF
4544	!
4545	!-- Assure that r_a cannot be zero at model start
4546	IF ( surf_lsm_h%pt1(m) == surf_lsm_h%pt_surface(m) ) &
4547	surf_lsm_h%pt1(m) = surf_lsm_h%pt1(m) + 1.0E-20_wp
4548
4549	surf_lsm_h%us(m) = 0.1_wp
4550	surf_lsm_h%ts(m) = ( surf_lsm_h%pt1(m) - surf_lsm_h%pt_surface(m) )&
4551	/ surf_lsm_h%r_a(m)
4552	surf_lsm_h%shf(m) = - surf_lsm_h%us(m) * surf_lsm_h%ts(m) &
4553	* rho_surface
4554	ENDDO
4555	!
4556	!-- Vertical surfaces
4557	DO l = 0, 3
4558	DO m = 1, surf_lsm_v(l)%ns
4559	i = surf_lsm_v(l)%i(m)
4560	j = surf_lsm_v(l)%j(m)
4561	k = surf_lsm_v(l)%k(m)
4562	!
4563	!-- Initialize surface temperature with soil temperature in the uppermost
4564	!-- uppermost layer
4565	t_surface_v(l)%var_1d(m) = t_soil_v(l)%var_2d(nzb_soil,m)
4566	surf_lsm_v(l)%pt_surface(m) = t_soil_v(l)%var_2d(nzb_soil,m) / exner(nzb)
4567
4568	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
4569	surf_lsm_v(l)%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
4570	ELSE
4571	surf_lsm_v(l)%pt1(m) = pt(k,j,i)
4572	ENDIF
4573
4574	!
4575	!-- Assure that r_a cannot be zero at model start
4576	IF ( surf_lsm_v(l)%pt1(m) == surf_lsm_v(l)%pt_surface(m) ) &
4577	surf_lsm_v(l)%pt1(m) = surf_lsm_v(l)%pt1(m) + 1.0E-20_wp
4578	!
4579	!-- Set artifical values for ts and us so that r_a has its initial value
4580	!-- for the first time step. Only for interior core domain, not for ghost points
4581	surf_lsm_v(l)%us(m) = 0.1_wp
4582	surf_lsm_v(l)%ts(m) = ( surf_lsm_v(l)%pt1(m) - surf_lsm_v(l)%pt_surface(m) ) /&
4583	surf_lsm_v(l)%r_a(m)