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

Last change on this file since 3049 was 3049, checked in by Giersch, 3 years ago

Revision history corrected

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