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

Last change on this file since 3241 was 3241, checked in by raasch, 3 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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