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

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

Bugfixes in Makefile and land_surface_model_mod

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