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

Last change on this file since 2249 was 2249, checked in by sward, 7 years ago

last commit documented

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