Home

Context Navigation

land_surface_model_mod.f90 @ 2798

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