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

Last change on this file since 4296 was 4296, checked in by maronga, 17 months ago

improved treatment of interception reservoir in land surface model

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