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

Last change on this file since 2282 was 2282, checked in by schwenkel, 4 years ago

Bugfix check saturation moisture

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