source: palm/trunk/SOURCE/land_surface_model_mod.f90 @ 3215

Last change on this file since 3215 was 3215, checked in by suehring, 3 years ago

changes for commit 3209 documented

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