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

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

bugfix for previous revision

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