Home

Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

land_surface_model_mod.f90

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