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

Last change on this file since 4441 was 4441, checked in by suehring, 3 years ago

Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better vectorization in the radiation interactions; Set back turbulent length scale to 8 x grid spacing in the parametrized mode for the synthetic turbulence generator (was accidentally changed in last commit)

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