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

Last change on this file since 3127 was 3127, checked in by maronga, 3 years ago

bugfix for high vegetation surface temperatures

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