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

Last change on this file since 2247 was 2247, checked in by sward, 4 years ago

last commit documented

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