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

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

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

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