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

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

Adjustments according new topography and surface-modelling concept implemented

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