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

Last change on this file since 3768 was 3767, checked in by raasch, 2 years ago

unused variables removed from rrd-subroutines parameter list

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