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

Last change on this file since 3045 was 3045, checked in by Giersch, 3 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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