Home

Context Navigation

land_surface_model_mod.f90 @ 4444

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