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

Last change on this file since 4360 was 4360, checked in by suehring, 19 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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