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

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

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

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