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

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

bugfix for last commit

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