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

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

added new spinup mechanism for surface/radiation models

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