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

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

improvements for spinup mechanism

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