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

Last change on this file since 3233 was 3233, checked in by schwenkel, 3 years ago

Adapted for the use of cloud_droplets

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