Home

Context Navigation

land_surface_model_mod.f90 @ 4312

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