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

Last change on this file since 4339 was 4339, checked in by suehring, 16 months ago

Bugfix, character strenght too short, caused crash on NEC

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