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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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