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

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

last commit documented

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