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

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

Minor revision of static input file checks, bugfix in initialization of surface-fractions in LSM; minor bugfix in initialization of albedo at window-surfaces; for clearer access of albedo and emissivity introduce index for vegetation/wall, pavement/green-wall and water/window surfaces

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