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

Last change on this file since 3196 was 3196, checked in by maronga, 3 years ago

set maximum value for aerodynamic resistance over horiztonal surfaces

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