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

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

Speed-up NetCDF input; Revise NetCDF-input routines and remove input via io-blocks; Temporarily revoke renaming of input variables in dynamic driver; More detailed error messages created; Bugfix in mapping 3D buildings; Bugfix in land-surface model at pavement surfaces; Bugfix in initialization with inifor

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