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

Last change on this file since 2881 was 2881, checked in by maronga, 7 years ago

bugfix in land surface model and new option in spinup

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