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

Last change on this file since 3246 was 3246, checked in by sward, 3 years ago

Added error handling for wrong input parameters

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