Home

Context Navigation

land_surface_model_mod.f90 @ 2938

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