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

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

bugfix for last commit

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