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

Last change on this file since 2504 was 2504, checked in by maronga, 4 years ago

lsm now allows for moist soil and roots below paved surfaces

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