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

Last change on this file since 2805 was 2805, checked in by suehring, 4 years ago

Bugfix in re-mapping surface-element data in case of restarts; bugfix in initialization of water-surface roughness; bugfix - uninitialized resistance at urban-type surfaces

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