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

Last change on this file since 2307 was 2307, checked in by suehring, 4 years ago

Bugfix,

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