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

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

Bugfix in referencing buildings on orography top; minor bugfix for last LSM commit

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