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

Last change on this file since 3026 was 3026, checked in by schwenkel, 3 years ago

Changed the name specific humidity to mixing ratio

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