Home

Context Navigation

land_surface_model_mod.f90 @ 4495

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