Home

Context Navigation

land_surface_model_mod.f90 @ 3142

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