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

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

Limit aerodynamic resistance at vertial building walls; check for roughness not exceeding surface-layer height

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