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

Last change on this file since 3004 was 3004, checked in by Giersch, 3 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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