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

Last change on this file since 4381 was 4381, checked in by suehring, 17 months ago

Land-surface model: Bugfix in nested soil initialization in case no dynamic input file is present; give local error messages only onces

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