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

Last change on this file since 2573 was 2573, checked in by scharf, 4 years ago

commit documented

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