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

Last change on this file since 2512 was 2512, checked in by raasch, 4 years ago

upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny; no output if redundant ghost layer data to NetCDF files

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