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

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

some improvements in land surface model

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