Home

Context Navigation

land_surface_model_mod.f90 @ 3014

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