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

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

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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