Home

Context Navigation

land_surface_model_mod.f90 @ 4313

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