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

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

modified calculation of surface resistance in land surface scheme

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