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

Last change on this file since 4442 was 4442, checked in by suehring, 4 years ago

last commit documented

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