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

Last change on this file since 4429 was 4429, checked in by raasch, 19 months ago

serial (non-MPI) test case added, several bugfixes for the serial mode

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