Home

Context Navigation

land_surface_model_mod.f90 @ 4381

Last change on this file since 4381 was 4381, checked in by suehring, 5 years ago
Land-surface model: Bugfix in nested soil initialization in case no dynamic input file is present; give local error messages only onces
Property svn:keywords set to `Id`
File size: 360.1 KB

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