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

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

revision of root fraction calculation in land surface model

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