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

Last change on this file since 1956 was 1950, checked in by maronga, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 109.9 KB
Line 
1!> @file land_surface_model_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: land_surface_model_mod.f90 1950 2016-06-17 07:21:31Z hellstea $
26!
27! 1949 2016-06-17 07:19:16Z maronga
28! Bugfix: calculation of qsws_soil_eb with precipitation = .TRUE. gave
29! qsws_soil_eb = 0 due to a typo
30!
31! 1856 2016-04-13 12:56:17Z maronga
32! Bugfix: for water surfaces, the initial water surface temperature is set equal
33! to the intital skin temperature. Moreover, the minimum value of r_a is now
34! 1.0 to avoid too large fluxes at the first model time step
35!
36! 1849 2016-04-08 11:33:18Z hoffmann
37! prr moved to arrays_3d
38!
39! 1826 2016-04-07 12:01:39Z maronga
40! Cleanup after modularization
41!
42! 1817 2016-04-06 15:44:20Z maronga
43! Added interface for lsm_init_arrays. Added subroutines for check_parameters,
44! header, and parin. Renamed some subroutines.
45!
46! 1788 2016-03-10 11:01:04Z maronga
47! Bugfix: calculate lambda_surface based on temperature gradient between skin
48! layer and soil layer instead of Obukhov length
49! Changed: moved calculation of surface specific humidity to energy balance solver
50! New: water surfaces are available by using a fixed sea surface temperature.
51! The roughness lengths are calculated dynamically using the Charnock
52! parameterization. This involves the new roughness length for moisture z0q.
53! New: modified solution of the energy balance solver and soil model for
54! paved surfaces (i.e. asphalt concrete).
55! Syntax layout improved.
56! Changed: parameter dewfall removed.
57!
58! 1783 2016-03-06 18:36:17Z raasch
59! netcdf variables moved to netcdf module
60!
61! 1757 2016-02-22 15:49:32Z maronga
62! Bugfix: set tm_soil_m to zero after allocation. Added parameter
63! unscheduled_radiation_calls to control calls of the radiation model based on
64! the skin temperature change during one time step (preliminary version). Set
65! qsws_soil_eb to zero at model start (previously set to qsws_eb). Removed MAX
66! function as it cannot be vectorized.
67!
68! 1709 2015-11-04 14:47:01Z maronga
69! Renamed pt_1 and qv_1 to pt1 and qv1.
70! Bugfix: set initial values for t_surface_p in case of restart runs
71! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
72! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
73! Added todo action
74!
75! 1697 2015-10-28 17:14:10Z raasch
76! bugfix: misplaced cpp-directive
77!
78! 1695 2015-10-27 10:03:11Z maronga
79! Bugfix: REAL constants provided with KIND-attribute in call of
80! Replaced rif with ol
81!
82! 1691 2015-10-26 16:17:44Z maronga
83! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes:
84! Soil temperatures are now defined at the edges of the layers, calculation of
85! shb_eb corrected, prognostic equation for skin temperature corrected. Surface
86! fluxes are now directly transfered to atmosphere
87!
88! 1682 2015-10-07 23:56:08Z knoop
89! Code annotations made doxygen readable
90!
91! 1590 2015-05-08 13:56:27Z maronga
92! Bugfix: definition of character strings requires same length for all elements
93!
94! 1585 2015-04-30 07:05:52Z maronga
95! Modifications for RRTMG. Changed tables to PARAMETER type.
96!
97! 1571 2015-03-12 16:12:49Z maronga
98! Removed upper-case variable names. Corrected distribution of precipitation to
99! the liquid water reservoir and the bare soil fractions.
100!
101! 1555 2015-03-04 17:44:27Z maronga
102! Added output of r_a and r_s
103!
104! 1553 2015-03-03 17:33:54Z maronga
105! Improved better treatment of roughness lengths. Added default soil temperature
106! profile
107!
108! 1551 2015-03-03 14:18:16Z maronga
109! Flux calculation is now done in prandtl_fluxes. Added support for data output.
110! Vertical indices have been replaced. Restart runs are now possible. Some
111! variables have beem renamed. Bugfix in the prognostic equation for the surface
112! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
113! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
114! the hydraulic conductivity. Bugfix for root fraction and extraction
115! calculation
116!
117! intrinsic function MAX and MIN
118!
119! 1500 2014-12-03 17:42:41Z maronga
120! Corrected calculation of aerodynamic resistance (r_a).
121! Precipitation is now added to liquid water reservoir using LE_liq.
122! Added support for dry runs.
123!
124! 1496 2014-12-02 17:25:50Z maronga
125! Initial revision
126!
127!
128! Description:
129! ------------
130!> Land surface model, consisting of a solver for the energy balance at the
131!> surface and a four layer soil scheme. The scheme is similar to the TESSEL
132!> scheme implemented in the ECMWF IFS model, with modifications according to
133!> H-TESSEL. The implementation is based on the formulation implemented in the
134!> DALES and UCLA-LES models.
135!>
136!> @todo Consider partial absorption of the net shortwave radiation by the
137!>       skin layer.
138!> @todo Improve surface water parameterization
139!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
140!>       nzt_soil=3)).
141!> @todo Implement surface runoff model (required when performing long-term LES
142!>       with considerable precipitation.
143!> @todo Fix crashes with radiation_scheme == 'constant'
144!>
145!> @note No time step criterion is required as long as the soil layers do not
146!>       become too thin.
147!------------------------------------------------------------------------------!
148 MODULE land_surface_model_mod
149 
150    USE arrays_3d,                                                             &
151        ONLY:  hyp, ol, pt, pt_p, prr, q, q_p, ql, qsws, shf, ts, us, vpt, z0, &
152               z0h, z0q
153
154    USE cloud_parameters,                                                      &
155        ONLY:  cp, hyrho, l_d_cp, l_d_r, l_v, pt_d_t, rho_l, r_d, r_v
156
157    USE control_parameters,                                                    &
158        ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
159               initializing_actions, intermediate_timestep_count_max,          &
160               max_masks, precipitation, pt_surface,                           &
161               rho_surface, roughness_length, surface_pressure,                &
162               timestep_scheme, tsc, z0h_factor, time_since_reference_point
163
164    USE indices,                                                               &
165        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 
166
167    USE kinds
168
169    USE pegrid
170
171    USE radiation_model_mod,                                                   &
172        ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
173               rad_lw_out, rad_lw_out_change_0, sigma_sb,                      &
174               unscheduled_radiation_calls
175       
176#if defined ( __rrtmg )
177    USE radiation_model_mod,                                                   &
178        ONLY:  rrtm_idrv
179#endif
180
181    USE statistics,                                                            &
182        ONLY:  hom, statistic_regions
183
184    IMPLICIT NONE
185
186!
187!-- LSM model constants
188    INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !< bottom of the soil model (to be switched)
189                               nzt_soil = 3, & !< top of the soil model (to be switched)
190                               nzs = 4         !< number of soil layers (fixed for now)
191
192    REAL(wp), PARAMETER ::                     &
193              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
194              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
195              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
196              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
197              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
198              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
199              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
200              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
201
202
203!
204!-- LSM variables
205    INTEGER(iwp) :: veg_type  = 2, & !< NAMELIST veg_type_2d
206                    soil_type = 3    !< NAMELIST soil_type_2d
207
208    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  soil_type_2d, &  !< soil type, 0: user-defined, 1-7: generic (see list)
209                                                  veg_type_2d      !< vegetation type, 0: user-defined, 1-19: generic (see list)
210
211    LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface,     & !< flag parameter for water surfaces (classes 14+15)
212                                            pave_surface,      & !< flag parameter for pavements (asphalt etc.) (class 20)
213                                            building_surface     !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
214
215    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
216               force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls
217               land_surface = .FALSE.              !< flag parameter indicating wheather the lsm is used
218
219!   value 9999999.9_wp -> generic available or user-defined value must be set
220!   otherwise -> no generic variable and user setting is optional
221    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
222                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
223                c_surface   = 20000.0_wp,               & !< Surface (skin) heat capacity
224                drho_l_lv,                              & !< (rho_l * l_v)**-1
225                exn,                                    & !< value of the Exner function
226                e_s = 0.0_wp,                           & !< saturation water vapour pressure
227                field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
228                f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
229                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
230                ke = 0.0_wp,                            & !< Kersten number
231                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil
232                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s
233                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u
234                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
235                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
236                min_canopy_resistance = 9999999.9_wp,   & !< NAMELIST r_canopy_min
237                min_soil_resistance = 50.0_wp,          & !< NAMELIST r_soil_min
238                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
239                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
240                pave_depth = 9999999.9_wp,              & !< depth of the pavement
241                pave_heat_capacity = 1.94E6_wp,         & !< volumetric heat capacity of pavement (e.g. roads)
242                pave_heat_conductivity = 1.00_wp,       & !< heat conductivity for pavements (e.g. roads)
243                q_s = 0.0_wp,                           & !< saturation specific humidity
244                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
245                rho_cp,                                 & !< rho_surface * cp
246                rho_lv,                                 & !< rho * l_v
247                rd_d_rv,                                & !< r_d / r_v
248                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
249                skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
250                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
251                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
252                z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
253                z0h_eb = 9999999.9_wp,                  & !< NAMELIST z0h (lsm_par)
254                z0q_eb = 9999999.9_wp                     !< NAMELIST z0q (lsm_par)
255
256    REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
257              ddz_soil_stag,                  & !< 1/dz_soil_stag
258              dz_soil_stag,                   & !< soil grid spacing (center-center)
259              root_extr = 0.0_wp,             & !< root extraction
260              root_fraction = (/9999999.9_wp, 9999999.9_wp,    &
261                                9999999.9_wp, 9999999.9_wp /), & !< distribution of root surface area to the individual soil layers
262              zs = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),    & !< soil layer depths (m)
263              soil_moisture = 0.0_wp          !< soil moisture content (m3/m3)
264
265    REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) ::   &
266              soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp,  283.0_wp,    & !< soil temperature (K)
267                                   283.0_wp /),                                &                                   
268              ddz_soil,                                                        & !< 1/dz_soil
269              dz_soil                                                            !< soil grid spacing (edge-edge)
270
271#if defined( __nopointer )
272    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !< surface temperature (K)
273                                                     t_surface_p, & !< progn. surface temperature (K)
274                                                     m_liq_eb,    & !< liquid water reservoir (m)
275                                                     m_liq_eb_av, & !< liquid water reservoir (m)
276                                                     m_liq_eb_p     !< progn. liquid water reservoir (m)
277#else
278    REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
279                                         t_surface_p,    & 
280                                         m_liq_eb,       & 
281                                         m_liq_eb_p
282
283    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
284                                                     m_liq_eb_av,              &
285                                                     m_liq_eb_1, m_liq_eb_2
286#endif
287
288!
289!-- Temporal tendencies for time stepping           
290    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !< surface temperature tendency (K)
291                                             tm_liq_eb_m      !< liquid water reservoir tendency (m)
292
293!
294!-- Energy balance variables               
295    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
296              alpha_vg,         & !< coef. of Van Genuchten
297              c_liq,            & !< liquid water coverage (of vegetated area)
298              c_liq_av,         & !< average of c_liq
299              c_soil_av,        & !< average of c_soil
300              c_veg,            & !< vegetation coverage
301              c_veg_av,         & !< average of c_veg
302              f_sw_in,          & !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
303              ghf_eb,           & !< ground heat flux
304              ghf_eb_av,        & !< average of ghf_eb
305              gamma_w_sat,      & !< hydraulic conductivity at saturation
306              g_d,              & !< coefficient for dependence of r_canopy on water vapour pressure deficit
307              lai,              & !< leaf area index
308              lai_av,           & !< average of lai
309              lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type)
310              lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type)
311              l_vg,             & !< coef. of Van Genuchten
312              m_fc,             & !< soil moisture at field capacity (m3/m3)
313              m_res,            & !< residual soil moisture
314              m_sat,            & !< saturation soil moisture (m3/m3)
315              m_wilt,           & !< soil moisture at permanent wilting point (m3/m3)
316              n_vg,             & !< coef. Van Genuchten 
317              qsws_eb,          & !< surface flux of latent heat (total)
318              qsws_eb_av,       & !< average of qsws_eb
319              qsws_liq_eb,      & !< surface flux of latent heat (liquid water portion)
320              qsws_liq_eb_av,   & !< average of qsws_liq_eb
321              qsws_soil_eb,     & !< surface flux of latent heat (soil portion)
322              qsws_soil_eb_av,  & !< average of qsws_soil_eb
323              qsws_veg_eb,      & !< surface flux of latent heat (vegetation portion)
324              qsws_veg_eb_av,   & !< average of qsws_veg_eb
325              rad_net_l,        & !< local copy of rad_net (net radiation at surface)
326              r_a,              & !< aerodynamic resistance
327              r_a_av,           & !< average of r_a
328              r_canopy,         & !< canopy resistance
329              r_soil,           & !< soil resistance
330              r_soil_min,       & !< minimum soil resistance
331              r_s,              & !< total surface resistance (combination of r_soil and r_canopy)
332              r_s_av,           & !< average of r_s
333              r_canopy_min,     & !< minimum canopy (stomatal) resistance
334              shf_eb,           & !< surface flux of sensible heat
335              shf_eb_av           !< average of shf_eb
336
337
338    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
339              lambda_h, &   !< heat conductivity of soil (W/m/K)                           
340              lambda_w, &   !< hydraulic diffusivity of soil (?)
341              gamma_w,  &   !< hydraulic conductivity of soil (W/m/K)
342              rho_c_total   !< volumetric heat capacity of the actual soil matrix (?)
343
344#if defined( __nopointer )
345    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
346              t_soil,    & !< Soil temperature (K)
347              t_soil_av, & !< Average of t_soil
348              t_soil_p,  & !< Prog. soil temperature (K)
349              m_soil,    & !< Soil moisture (m3/m3)
350              m_soil_av, & !< Average of m_soil
351              m_soil_p     !< Prog. soil moisture (m3/m3)
352#else
353    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
354              t_soil, t_soil_p, &
355              m_soil, m_soil_p   
356
357    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
358              t_soil_av, t_soil_1, t_soil_2,                                   &
359              m_soil_av, m_soil_1, m_soil_2
360#endif
361
362
363    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
364              tt_soil_m, & !< t_soil storage array
365              tm_soil_m, & !< m_soil storage array
366              root_fr      !< root fraction (sum=1)
367
368
369!
370!-- Predefined Land surface classes (veg_type)
371    CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ &
372                                   'user defined              ',    & ! 0
373                                   'crops, mixed farming      ',    & !  1
374                                   'short grass               ',    & !  2
375                                   'evergreen needleleaf trees',    & !  3
376                                   'deciduous needleleaf trees',    & !  4
377                                   'evergreen broadleaf trees ',    & !  5
378                                   'deciduous broadleaf trees ',    & !  6
379                                   'tall grass                ',    & !  7
380                                   'desert                    ',    & !  8
381                                   'tundra                    ',    & !  9
382                                   'irrigated crops           ',    & ! 10
383                                   'semidesert                ',    & ! 11
384                                   'ice caps and glaciers     ',    & ! 12
385                                   'bogs and marshes          ',    & ! 13
386                                   'inland water              ',    & ! 14
387                                   'ocean                     ',    & ! 15
388                                   'evergreen shrubs          ',    & ! 16
389                                   'deciduous shrubs          ',    & ! 17
390                                   'mixed forest/woodland     ',    & ! 18
391                                   'interrupted forest        ',    & ! 19
392                                   'pavements/roads           '     & ! 20
393                                                                 /)
394
395!
396!-- Soil model classes (soil_type)
397    CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
398                                   'user defined',                  & ! 0
399                                   'coarse      ',                  & ! 1
400                                   'medium      ',                  & ! 2
401                                   'medium-fine ',                  & ! 3
402                                   'fine        ',                  & ! 4
403                                   'very fine   ',                  & ! 5
404                                   'organic     ',                  & ! 6
405                                   'loamy (CH)  '                   & ! 7
406                                                                 /)
407!
408!-- Land surface parameters according to the respective classes (veg_type)
409
410!
411!-- Land surface parameters I
412!--                          r_canopy_min,     lai,   c_veg,     g_d
413    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ &
414                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
415                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
416                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
417                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
418                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
419                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
420                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
421                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
422                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
423                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
424                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
425                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
426                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
427                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
428                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
429                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
430                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
431                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
432                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,  & ! 19
433                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp   & ! 20
434                                 /), (/ 4, 20 /) )
435
436!
437!-- Land surface parameters II          z0,         z0h,         z0q
438    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ & 
439                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & !  1
440                                   0.20_wp,  0.20E-2_wp,  0.20E-2_wp,       & !  2
441                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  3
442                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  4
443                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  5
444                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  6
445                                   0.47_wp,  0.47E-2_wp,  0.47E-2_wp,       & !  7
446                                  0.013_wp, 0.013E-2_wp, 0.013E-2_wp,       & !  8
447                                  0.034_wp, 0.034E-2_wp, 0.034E-2_wp,       & !  9
448                                    0.5_wp,  0.50E-2_wp,  0.50E-2_wp,       & ! 10
449                                   0.17_wp,  0.17E-2_wp,  0.17E-2_wp,       & ! 11
450                                 1.3E-3_wp,   1.3E-4_wp,   1.3E-4_wp,       & ! 12
451                                   0.83_wp,  0.83E-2_wp,  0.83E-2_wp,       & ! 13
452                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 14
453                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 15
454                                   0.10_wp,  0.10E-2_wp,  0.10E-2_wp,       & ! 16
455                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & ! 17
456                                   2.00_wp,  2.00E-2_wp,  2.00E-2_wp,       & ! 18
457                                   1.10_wp,  1.10E-2_wp,  1.10E-2_wp,       & ! 19
458                                 1.0E-4_wp,   1.0E-5_wp,   1.0E-5_wp        & ! 20
459                                 /), (/ 3, 20 /) )
460
461!
462!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
463    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ &
464                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
465                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
466                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  3
467                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  4
468                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  5
469                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  6
470                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  7
471                                      15.0_wp,       15.0_wp, 0.00_wp,     & !  8
472                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  9
473                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 10
474                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
475                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
476                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13
477                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 14
478                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 15
479                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
480                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
481                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
482                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 19
483                                       0.0_wp,        0.0_wp, 0.00_wp      & ! 20
484                                      /), (/ 3, 20 /) )
485
486!
487!-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
488    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ &
489                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
490                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
491                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  3
492                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  4
493                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  5
494                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  6
495                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  7
496                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  8
497                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & !  9
498                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 10
499                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 11
500                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 12
501                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 13
502                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 14
503                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 15
504                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
505                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17
506                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
507                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 19
508                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp             & ! 20
509                                 /), (/ 4, 20 /) )
510
511!
512!-- Soil parameters according to the following porosity classes (soil_type)
513
514!
515!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
516    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
517                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
518                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
519                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
520                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
521                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
522                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
523                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
524                                 /), (/ 4, 7 /) )
525
526!
527!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
528    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
529                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
530                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
531                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
532                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
533                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
534                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
535                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
536                                 /), (/ 4, 7 /) )
537
538
539    SAVE
540
541
542    PRIVATE
543
544   
545!
546!-- Public functions
547    PUBLIC lsm_check_data_output, lsm_check_data_output_pr,                    &
548           lsm_check_parameters, lsm_energy_balance, lsm_header, lsm_init,     &
549           lsm_init_arrays, lsm_parin, lsm_soil_model, lsm_swap_timelevel 
550!
551!-- Public parameters, constants and initial values
552    PUBLIC land_surface, skip_time_do_lsm
553
554!
555!-- Public grid variables
556    PUBLIC nzb_soil, nzs, nzt_soil, zs
557
558!
559!-- Public 2D output variables
560    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
561           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
562           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
563           r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
564
565!
566!-- Public prognostic variables
567    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
568
569
570    INTERFACE lsm_check_data_output
571       MODULE PROCEDURE lsm_check_data_output
572    END INTERFACE lsm_check_data_output
573   
574    INTERFACE lsm_check_data_output_pr
575       MODULE PROCEDURE lsm_check_data_output_pr
576    END INTERFACE lsm_check_data_output_pr
577   
578    INTERFACE lsm_check_parameters
579       MODULE PROCEDURE lsm_check_parameters
580    END INTERFACE lsm_check_parameters
581   
582    INTERFACE lsm_energy_balance
583       MODULE PROCEDURE lsm_energy_balance
584    END INTERFACE lsm_energy_balance
585
586    INTERFACE lsm_header
587       MODULE PROCEDURE lsm_header
588    END INTERFACE lsm_header
589   
590    INTERFACE lsm_init
591       MODULE PROCEDURE lsm_init
592    END INTERFACE lsm_init
593
594    INTERFACE lsm_init_arrays
595       MODULE PROCEDURE lsm_init_arrays
596    END INTERFACE lsm_init_arrays
597   
598    INTERFACE lsm_parin
599       MODULE PROCEDURE lsm_parin
600    END INTERFACE lsm_parin
601   
602    INTERFACE lsm_soil_model
603       MODULE PROCEDURE lsm_soil_model
604    END INTERFACE lsm_soil_model
605
606    INTERFACE lsm_swap_timelevel
607       MODULE PROCEDURE lsm_swap_timelevel
608    END INTERFACE lsm_swap_timelevel
609
610 CONTAINS
611
612!------------------------------------------------------------------------------!
613! Description:
614! ------------
615!> Check data output for land surface model
616!------------------------------------------------------------------------------!
617    SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
618 
619 
620       USE control_parameters,                                                 &
621           ONLY: data_output, message_string
622
623       IMPLICIT NONE
624
625       CHARACTER (LEN=*) ::  unit     !<
626       CHARACTER (LEN=*) ::  var      !<
627
628       INTEGER(iwp) :: i
629       INTEGER(iwp) :: ilen   
630       INTEGER(iwp) :: k
631
632       SELECT CASE ( TRIM( var ) )
633
634          CASE ( 'm_soil' )
635             IF (  .NOT.  land_surface )  THEN
636                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
637                         'res land_surface = .TRUE.'
638                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
639             ENDIF
640             unit = 'm3/m3'
641             
642          CASE ( 't_soil' )
643             IF (  .NOT.  land_surface )  THEN
644                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
645                         'res land_surface = .TRUE.'
646                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
647             ENDIF
648             unit = 'K'   
649             
650          CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',&
651                 'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',  &
652                 'r_a*', 'r_s*', 'shf_eb*' )
653             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
654                message_string = 'illegal value for data_output: "' //         &
655                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
656                                 'cross sections are allowed for this value'
657                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
658             ENDIF
659             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
660                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
661                                 'res land_surface = .TRUE.'
662                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
663             ENDIF
664             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
665                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
666                                 'res land_surface = .TRUE.'
667                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
668             ENDIF
669             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
670                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
671                                 'res land_surface = .TRUE.'
672                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
673             ENDIF
674             IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
675                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
676                                 'res land_surface = .TRUE.'
677                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
678             ENDIF
679             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
680                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
681                                 'res land_surface = .TRUE.'
682                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
683             ENDIF
684             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
685                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
686                                 'res land_surface = .TRUE.'
687                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
688             ENDIF
689             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
690                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
691                                 'res land_surface = .TRUE.'
692                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
693             ENDIF
694             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
695             THEN
696                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
697                                 'res land_surface = .TRUE.'
698                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
699             ENDIF
700             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
701             THEN
702                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
703                                 'res land_surface = .TRUE.'
704                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
705             ENDIF
706             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
707             THEN
708                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
709                                 'res land_surface = .TRUE.'
710                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
711             ENDIF
712             IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
713             THEN
714                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
715                                 'res land_surface = .TRUE.'
716                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
717             ENDIF
718             IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
719             THEN
720                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
721                                 'res land_surface = .TRUE.'
722                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
723             ENDIF
724
725             IF ( TRIM( var ) == 'lai*'   )  unit = 'none' 
726             IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
727             IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
728             IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
729             IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
730             IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
731             IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
732             IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
733             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
734             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
735             IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
736             IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
737             IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
738             
739          CASE DEFAULT
740             unit = 'illegal'
741
742       END SELECT
743
744
745    END SUBROUTINE lsm_check_data_output
746
747 SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
748 
749       USE control_parameters,                                                 &
750           ONLY: data_output_pr, message_string
751
752       USE indices
753
754       USE profil_parameter
755
756       USE statistics
757
758       IMPLICIT NONE
759   
760       CHARACTER (LEN=*) ::  unit      !<
761       CHARACTER (LEN=*) ::  variable  !<
762       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
763 
764       INTEGER(iwp) ::  user_pr_index !<
765       INTEGER(iwp) ::  var_count     !<
766
767       SELECT CASE ( TRIM( variable ) )
768       
769          CASE ( 't_soil', '#t_soil' )
770             IF (  .NOT.  land_surface )  THEN
771                message_string = 'data_output_pr = ' //                        &
772                                 TRIM( data_output_pr(var_count) ) // ' is' // &
773                                 'not implemented for land_surface = .FALSE.'
774                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
775             ELSE
776                dopr_index(var_count) = 89
777                dopr_unit     = 'K'
778                hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
779                IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
780                   dopr_initial_index(var_count) = 90
781                   hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
782                   data_output_pr(var_count)     = data_output_pr(var_count)(2:)
783                ENDIF
784                unit = dopr_unit
785             ENDIF
786
787          CASE ( 'm_soil', '#m_soil' )
788             IF (  .NOT.  land_surface )  THEN
789                message_string = 'data_output_pr = ' //                        &
790                                 TRIM( data_output_pr(var_count) ) // ' is' // &
791                                 ' not implemented for land_surface = .FALSE.'
792                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
793             ELSE
794                dopr_index(var_count) = 91
795                dopr_unit     = 'm3/m3'
796                hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
797                IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
798                   dopr_initial_index(var_count) = 92
799                   hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
800                   data_output_pr(var_count)     = data_output_pr(var_count)(2:)
801                ENDIF
802                 unit = dopr_unit
803             ENDIF
804
805
806          CASE DEFAULT
807             unit = 'illegal'
808
809       END SELECT
810
811
812    END SUBROUTINE lsm_check_data_output_pr
813 
814 
815!------------------------------------------------------------------------------!
816! Description:
817! ------------
818!> Check parameters routine for land surface model
819!------------------------------------------------------------------------------!
820    SUBROUTINE lsm_check_parameters
821
822       USE control_parameters,                                                 &
823           ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string,         &
824                 most_method, topography
825                 
826       USE radiation_model_mod,                                                &
827           ONLY: radiation
828   
829   
830       IMPLICIT NONE
831
832 
833!
834!--    Dirichlet boundary conditions are required as the surface fluxes are
835!--    calculated from the temperature/humidity gradients in the land surface
836!--    model
837       IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
838          message_string = 'lsm requires setting of'//                         &
839                           'bc_pt_b = "dirichlet" and '//                      &
840                           'bc_q_b  = "dirichlet"'
841          CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
842       ENDIF
843
844       IF (  .NOT.  constant_flux_layer )  THEN
845          message_string = 'lsm requires '//                                   &
846                           'constant_flux_layer = .T.'
847          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
848       ENDIF
849
850       IF ( topography /= 'flat' )  THEN
851          message_string = 'lsm cannot be used ' //                            & 
852                           'in combination with  topography /= "flat"'
853          CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
854       ENDIF
855
856       IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
857              most_method == 'lookup' )  THEN
858           WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
859                                      'allowed in combination with ',          &
860                                      'most_method = ', most_method
861          CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
862       ENDIF
863
864       IF ( veg_type == 0 )  THEN
865          IF ( SUM( root_fraction ) /= 1.0_wp )  THEN
866             message_string = 'veg_type = 0 (user_defined)'//                  &
867                              'requires setting of root_fraction(0:3)'//       &
868                              '/= 9999999.9 and SUM(root_fraction) = 1'
869             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
870          ENDIF
871 
872          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
873             message_string = 'veg_type = 0 (user defined)'//                  &
874                              'requires setting of min_canopy_resistance'//    &
875                              '/= 9999999.9'
876             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
877          ENDIF
878
879          IF ( leaf_area_index == 9999999.9_wp )  THEN
880             message_string = 'veg_type = 0 (user_defined)'//                  &
881                              'requires setting of leaf_area_index'//          &
882                              '/= 9999999.9'
883             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
884          ENDIF
885
886          IF ( vegetation_coverage == 9999999.9_wp )  THEN
887             message_string = 'veg_type = 0 (user_defined)'//                  &
888                              'requires setting of vegetation_coverage'//      &
889                              '/= 9999999.9'
890             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
891          ENDIF
892
893          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
894             message_string = 'veg_type = 0 (user_defined)'//                  &
895                              'requires setting of'//                          &
896                              'canopy_resistance_coefficient /= 9999999.9'
897             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
898          ENDIF
899
900          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
901             message_string = 'veg_type = 0 (user_defined)'//                  &
902                              'requires setting of lambda_surface_stable'//    &
903                              '/= 9999999.9'
904             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
905          ENDIF
906
907          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
908             message_string = 'veg_type = 0 (user_defined)'//                  &
909                              'requires setting of lambda_surface_unstable'//  &
910                              '/= 9999999.9'
911             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
912          ENDIF
913
914          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
915             message_string = 'veg_type = 0 (user_defined)'//                  &
916                              'requires setting of f_shortwave_incoming'//     &
917                              '/= 9999999.9'
918             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
919          ENDIF
920
921          IF ( z0_eb == 9999999.9_wp )  THEN
922             message_string = 'veg_type = 0 (user_defined)'//                  &
923                              'requires setting of z0_eb'//                    &
924                              '/= 9999999.9'
925             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
926          ENDIF
927
928          IF ( z0h_eb == 9999999.9_wp )  THEN
929             message_string = 'veg_type = 0 (user_defined)'//                  &
930                              'requires setting of z0h_eb'//                   &
931                              '/= 9999999.9'
932             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
933          ENDIF
934
935
936       ENDIF
937
938       IF ( soil_type == 0 )  THEN
939
940          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
941             message_string = 'soil_type = 0 (user_defined)'//                 &
942                              'requires setting of alpha_vangenuchten'//       &
943                              '/= 9999999.9'
944             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
945          ENDIF
946
947          IF ( l_vangenuchten == 9999999.9_wp )  THEN
948             message_string = 'soil_type = 0 (user_defined)'//                 &
949                              'requires setting of l_vangenuchten'//           &
950                              '/= 9999999.9'
951             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
952          ENDIF
953
954          IF ( n_vangenuchten == 9999999.9_wp )  THEN
955             message_string = 'soil_type = 0 (user_defined)'//                 &
956                              'requires setting of n_vangenuchten'//           &
957                              '/= 9999999.9'
958             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
959          ENDIF
960
961          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
962             message_string = 'soil_type = 0 (user_defined)'//                 &
963                              'requires setting of hydraulic_conductivity'//   &
964                              '/= 9999999.9'
965             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
966          ENDIF
967
968          IF ( saturation_moisture == 9999999.9_wp )  THEN
969             message_string = 'soil_type = 0 (user_defined)'//                 &
970                              'requires setting of saturation_moisture'//      &
971                              '/= 9999999.9'
972             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
973          ENDIF
974
975          IF ( field_capacity == 9999999.9_wp )  THEN
976             message_string = 'soil_type = 0 (user_defined)'//                 &
977                              'requires setting of field_capacity'//           &
978                              '/= 9999999.9'
979             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
980          ENDIF
981
982          IF ( wilting_point == 9999999.9_wp )  THEN
983             message_string = 'soil_type = 0 (user_defined)'//                 &
984                              'requires setting of wilting_point'//            &
985                              '/= 9999999.9'
986             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
987          ENDIF
988
989          IF ( residual_moisture == 9999999.9_wp )  THEN
990             message_string = 'soil_type = 0 (user_defined)'//                 &
991                              'requires setting of residual_moisture'//        &
992                              '/= 9999999.9'
993             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
994          ENDIF
995
996       ENDIF
997
998       IF (  .NOT.  radiation )  THEN
999          message_string = 'lsm requires '//                                   &
1000                           'radiation = .T.'
1001          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1002       ENDIF
1003       
1004       
1005    END SUBROUTINE lsm_check_parameters
1006 
1007!------------------------------------------------------------------------------!
1008! Description:
1009! ------------
1010!> Solver for the energy balance at the surface.
1011!------------------------------------------------------------------------------!
1012    SUBROUTINE lsm_energy_balance
1013
1014
1015       IMPLICIT NONE
1016
1017       INTEGER(iwp) ::  i         !< running index
1018       INTEGER(iwp) ::  j         !< running index
1019       INTEGER(iwp) ::  k, ks     !< running index
1020
1021       REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
1022                   f1,          & !< resistance correction term 1
1023                   f2,          & !< resistance correction term 2
1024                   f3,          & !< resistance correction term 3
1025                   m_min,       & !< minimum soil moisture
1026                   e,           & !< water vapour pressure
1027                   e_s,         & !< water vapour saturation pressure
1028                   e_s_dt,      & !< derivate of e_s with respect to T
1029                   tend,        & !< tendency
1030                   dq_s_dt,     & !< derivate of q_s with respect to T
1031                   coef_1,      & !< coef. for prognostic equation
1032                   coef_2,      & !< coef. for prognostic equation
1033                   f_qsws,      & !< factor for qsws_eb
1034                   f_qsws_veg,  & !< factor for qsws_veg_eb
1035                   f_qsws_soil, & !< factor for qsws_soil_eb
1036                   f_qsws_liq,  & !< factor for qsws_liq_eb
1037                   f_shf,       & !< factor for shf_eb
1038                   lambda_surface, & !< Current value of lambda_surface
1039                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
1040                   pt1,         & !< potential temperature at first grid level
1041                   qv1            !< specific humidity at first grid level
1042
1043!
1044!--    Calculate the exner function for the current time step
1045       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
1046
1047       DO  i = nxlg, nxrg
1048          DO  j = nysg, nyng
1049             k = nzb_s_inner(j,i)
1050
1051!
1052!--          Set lambda_surface according to stratification between skin layer and soil
1053             IF (  .NOT.  pave_surface(j,i) )  THEN
1054
1055                c_surface_tmp = c_surface
1056
1057                IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
1058                   lambda_surface = lambda_surface_s(j,i)
1059                ELSE
1060                   lambda_surface = lambda_surface_u(j,i)
1061                ENDIF
1062             ELSE
1063
1064                c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
1065                lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
1066
1067             ENDIF
1068
1069!
1070!--          First step: calculate aerodyamic resistance. As pt, us, ts
1071!--          are not available for the prognostic time step, data from the last
1072!--          time step is used here. Note that this formulation is the
1073!--          equivalent to the ECMWF formulation using drag coefficients
1074             IF ( cloud_physics )  THEN
1075                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
1076                qv1 = q(k+1,j,i) - ql(k+1,j,i)
1077             ELSE
1078                pt1 = pt(k+1,j,i)
1079                qv1 = q(k+1,j,i)
1080             ENDIF
1081
1082             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
1083
1084!
1085!--          Make sure that the resistance does not drop to zero for neutral
1086!--          stratification
1087             IF ( ABS(r_a(j,i)) < 1.0_wp )  r_a(j,i) = 1.0_wp
1088
1089!
1090!--          Second step: calculate canopy resistance r_canopy
1091!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1092 
1093!--          f1: correction for incoming shortwave radiation (stomata close at
1094!--          night)
1095             IF ( radiation_scheme /= 'constant' )  THEN
1096                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
1097                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
1098                               + 1.0_wp)) )
1099             ELSE
1100                f1 = 1.0_wp
1101             ENDIF
1102
1103
1104!
1105!--          f2: correction for soil moisture availability to plants (the
1106!--          integrated soil moisture must thus be considered here)
1107!--          f2 = 0 for very dry soils
1108             m_total = 0.0_wp
1109             DO  ks = nzb_soil, nzt_soil
1110                 m_total = m_total + root_fr(ks,j,i)                           &
1111                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
1112             ENDDO 
1113
1114             IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
1115                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
1116             ELSEIF ( m_total >= m_fc(j,i) )  THEN
1117                f2 = 1.0_wp
1118             ELSE
1119                f2 = 1.0E-20_wp
1120             ENDIF
1121
1122!
1123!--          Calculate water vapour pressure at saturation
1124             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1125                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
1126
1127!
1128!--          f3: correction for vapour pressure deficit
1129             IF ( g_d(j,i) /= 0.0_wp )  THEN
1130!
1131!--             Calculate vapour pressure
1132                e  = qv1 * surface_pressure / 0.622_wp
1133                f3 = EXP ( -g_d(j,i) * (e_s - e) )
1134             ELSE
1135                f3 = 1.0_wp
1136             ENDIF
1137
1138!
1139!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1140!--          this calculation is obsolete, as r_canopy is not used below.
1141!--          To do: check for very dry soil -> r_canopy goes to infinity
1142             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
1143                                             + 1.0E-20_wp)
1144
1145!
1146!--          Third step: calculate bare soil resistance r_soil. The Clapp &
1147!--          Hornberger parametrization does not consider c_veg.
1148             IF ( soil_type_2d(j,i) /= 7 )  THEN
1149                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
1150                        m_res(j,i)
1151             ELSE
1152                m_min = m_wilt(j,i)
1153             ENDIF
1154
1155             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
1156             f2 = MAX(f2,1.0E-20_wp)
1157             f2 = MIN(f2,1.0_wp)
1158
1159             r_soil(j,i) = r_soil_min(j,i) / f2
1160
1161!
1162!--          Calculate the maximum possible liquid water amount on plants and
1163!--          bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
1164!--          assumed, while paved surfaces might hold up 1 mm of water. The
1165!--          liquid water fraction for paved surfaces is calculated after
1166!--          Noilhan & Planton (1989), while the ECMWF formulation is used for
1167!--          vegetated surfaces and bare soils.
1168             IF ( pave_surface(j,i) )  THEN
1169                m_liq_eb_max = m_max_depth * 5.0_wp
1170                c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
1171             ELSE
1172                m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)           &
1173                               + (1.0_wp - c_veg(j,i)) )
1174                c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
1175             ENDIF
1176
1177!
1178!--          Calculate saturation specific humidity
1179             q_s = 0.622_wp * e_s / surface_pressure
1180
1181!
1182!--          In case of dewfall, set evapotranspiration to zero
1183!--          All super-saturated water is then removed from the air
1184             IF ( humidity  .AND.  q_s <= qv1 )  THEN
1185                r_canopy(j,i) = 0.0_wp
1186                r_soil(j,i)   = 0.0_wp
1187             ENDIF
1188
1189!
1190!--          Calculate coefficients for the total evapotranspiration
1191!--          In case of water surface, set vegetation and soil fluxes to zero.
1192!--          For pavements, only evaporation of liquid water is possible.
1193             IF ( water_surface(j,i) )  THEN
1194                f_qsws_veg  = 0.0_wp
1195                f_qsws_soil = 0.0_wp
1196                f_qsws_liq  = rho_lv / r_a(j,i)
1197             ELSEIF ( pave_surface (j,i) )  THEN
1198                f_qsws_veg  = 0.0_wp
1199                f_qsws_soil = 0.0_wp
1200                f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
1201             ELSE
1202                f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/     &
1203                              (r_a(j,i) + r_canopy(j,i))
1204                f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +     &
1205                                                          r_soil(j,i))
1206                f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
1207             ENDIF
1208!
1209!--          If soil moisture is below wilting point, plants do no longer
1210!--          transpirate.
1211!              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
1212!                 f_qsws_veg = 0.0_wp
1213!              ENDIF
1214
1215             f_shf  = rho_cp / r_a(j,i)
1216             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1217
1218!
1219!--          Calculate derivative of q_s for Taylor series expansion
1220             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
1221                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
1222                              / (t_surface(j,i) - 35.86_wp)**2 )
1223
1224             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
1225
1226!
1227!--          Add LW up so that it can be removed in prognostic equation
1228             rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
1229
1230!
1231!--          Calculate new skin temperature
1232             IF ( humidity )  THEN
1233#if defined ( __rrtmg )
1234!
1235!--             Numerator of the prognostic equation
1236                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
1237                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
1238                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
1239                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
1240                         * t_soil(nzb_soil,j,i)
1241
1242!
1243!--             Denominator of the prognostic equation
1244                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
1245                         + lambda_surface + f_shf / exn
1246#else
1247
1248!
1249!--             Numerator of the prognostic equation
1250                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
1251                         * t_surface(j,i) ** 4                                 &
1252                         + f_shf * pt1 + f_qsws * ( qv1                        &
1253                         - q_s + dq_s_dt * t_surface(j,i) )                    &
1254                         + lambda_surface * t_soil(nzb_soil,j,i)
1255
1256!
1257!--             Denominator of the prognostic equation
1258                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
1259                         * dq_s_dt + lambda_surface + f_shf / exn
1260 
1261#endif
1262             ELSE
1263
1264#if defined ( __rrtmg )
1265!
1266!--             Numerator of the prognostic equation
1267                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
1268                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
1269                         + f_shf * pt1  + lambda_surface                       &
1270                         * t_soil(nzb_soil,j,i)
1271
1272!
1273!--             Denominator of the prognostic equation
1274                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
1275#else
1276
1277!
1278!--             Numerator of the prognostic equation
1279                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
1280                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
1281                          + lambda_surface * t_soil(nzb_soil,j,i)
1282
1283!
1284!--             Denominator of the prognostic equation
1285                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
1286                         + lambda_surface + f_shf / exn
1287#endif
1288             ENDIF
1289
1290             tend = 0.0_wp
1291
1292!
1293!--          Implicit solution when the surface layer has no heat capacity,
1294!--          otherwise use RK3 scheme.
1295             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *    &
1296                                t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
1297                                * dt_3d * tsc(2) ) 
1298
1299!
1300!--          Add RK3 term
1301             IF ( c_surface_tmp /= 0.0_wp )  THEN
1302
1303                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
1304                                   * tt_surface_m(j,i)
1305
1306!
1307!--             Calculate true tendency
1308                tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
1309                       * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
1310!
1311!--             Calculate t_surface tendencies for the next Runge-Kutta step
1312                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1313                   IF ( intermediate_timestep_count == 1 )  THEN
1314                      tt_surface_m(j,i) = tend
1315                   ELSEIF ( intermediate_timestep_count <                      &
1316                            intermediate_timestep_count_max )  THEN
1317                      tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
1318                                          * tt_surface_m(j,i)
1319                   ENDIF
1320                ENDIF
1321             ENDIF
1322
1323!
1324!--          In case of fast changes in the skin temperature, it is possible to
1325!--          update the radiative fluxes independently from the prescribed
1326!--          radiation call frequency. This effectively prevents oscillations,
1327!--          especially when setting skip_time_do_radiation /= 0. The threshold
1328!--          value of 0.2 used here is just a first guess. This method should be
1329!--          revised in the future as tests have shown that the threshold is
1330!--          often reached, when no oscillations would occur (causes immense
1331!--          computing time for the radiation code).
1332             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.     &
1333                  unscheduled_radiation_calls )  THEN
1334                force_radiation_call_l = .TRUE.
1335             ENDIF
1336
1337             pt(k,j,i) = t_surface_p(j,i) / exn
1338
1339!
1340!--          Calculate fluxes
1341#if defined ( __rrtmg )
1342             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
1343                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
1344                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
1345
1346             IF ( rrtm_idrv == 1 )  THEN
1347                rad_net(j,i) = rad_net_l(j,i)
1348                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
1349                                      + rad_lw_out_change_0(j,i)               &
1350                                      * ( t_surface_p(j,i) - t_surface(j,i) )
1351             ENDIF
1352#else
1353             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb             &
1354                                * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
1355                                * t_surface(j,i)**3 * t_surface_p(j,i)
1356#endif
1357
1358             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
1359                              - t_soil(nzb_soil,j,i))
1360
1361             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
1362
1363             shf(j,i) = shf_eb(j,i) / rho_cp
1364
1365             IF ( humidity )  THEN
1366                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
1367                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
1368
1369                qsws(j,i) = qsws_eb(j,i) / rho_lv
1370
1371                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
1372                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1373                                    * t_surface_p(j,i) )
1374
1375                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
1376                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1377                                    * t_surface_p(j,i) )
1378
1379                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
1380                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1381                                    * t_surface_p(j,i) )
1382             ENDIF
1383
1384!
1385!--          Calculate the true surface resistance
1386             IF ( qsws_eb(j,i) == 0.0_wp )  THEN
1387                r_s(j,i) = 1.0E10_wp
1388             ELSE
1389                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                    &
1390                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
1391                           / qsws_eb(j,i) - r_a(j,i)
1392             ENDIF
1393
1394!
1395!--          Calculate change in liquid water reservoir due to dew fall or
1396!--          evaporation of liquid water
1397             IF ( humidity )  THEN
1398!
1399!--             If precipitation is activated, add rain water to qsws_liq_eb
1400!--             and qsws_soil_eb according the the vegetation coverage.
1401!--             precipitation_rate is given in mm.
1402                IF ( precipitation )  THEN
1403
1404!
1405!--                Add precipitation to liquid water reservoir, if possible.
1406!--                Otherwise, add the water to soil. In case of
1407!--                pavements, the exceeding water amount is implicitely removed
1408!--                as runoff as qsws_soil_eb is then not used in the soil model
1409                   IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
1410                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
1411                                       + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
1412                                       * 0.001_wp * rho_l * l_v
1413                   ELSE
1414                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1415                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
1416                                        * 0.001_wp * rho_l * l_v
1417                   ENDIF
1418
1419!--                Add precipitation to bare soil according to the bare soil
1420!--                coverage.
1421                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp             &
1422                                       - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
1423                                       * 0.001_wp * rho_l * l_v
1424                ENDIF
1425
1426!
1427!--             If the air is saturated, check the reservoir water level
1428                IF ( qsws_eb(j,i) < 0.0_wp )  THEN
1429
1430!
1431!--                Check if reservoir is full (avoid values > m_liq_eb_max)
1432!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
1433!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
1434!--                so that tend is zero and no further check is needed
1435                   IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
1436                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1437                                           + qsws_liq_eb(j,i)
1438
1439                      qsws_liq_eb(j,i)  = 0.0_wp
1440                   ENDIF
1441
1442!
1443!--                In case qsws_veg_eb becomes negative (unphysical behavior),
1444!--                let the water enter the liquid water reservoir as dew on the
1445!--                plant
1446                   IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
1447                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
1448                      qsws_veg_eb(j,i) = 0.0_wp
1449                   ENDIF
1450                ENDIF                   
1451 
1452                tend = - qsws_liq_eb(j,i) * drho_l_lv
1453
1454                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
1455                                                   + tsc(3) * tm_liq_eb_m(j,i) )
1456
1457!
1458!--             Check if reservoir is overfull -> reduce to maximum
1459!--             (conservation of water is violated here)
1460                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
1461
1462!
1463!--             Check if reservoir is empty (avoid values < 0.0)
1464!--             (conservation of water is violated here)
1465                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
1466
1467
1468!
1469!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
1470                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1471                   IF ( intermediate_timestep_count == 1 )  THEN
1472                      tm_liq_eb_m(j,i) = tend
1473                   ELSEIF ( intermediate_timestep_count <                      &
1474                            intermediate_timestep_count_max )  THEN
1475                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
1476                                      * tm_liq_eb_m(j,i)
1477                   ENDIF
1478                ENDIF
1479
1480             ENDIF
1481
1482          ENDDO
1483       ENDDO
1484
1485!
1486!--    Make a logical OR for all processes. Force radiation call if at
1487!--    least one processor reached the threshold change in skin temperature
1488       IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count    &
1489            == intermediate_timestep_count_max-1 )  THEN
1490#if defined( __parallel )
1491          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1492          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
1493                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
1494#else
1495          force_radiation_call = force_radiation_call_l
1496#endif
1497          force_radiation_call_l = .FALSE.
1498       ENDIF
1499
1500!
1501!--    Calculate surface specific humidity
1502       IF ( humidity )  THEN
1503          CALL calc_q_surface
1504       ENDIF
1505
1506!
1507!--    Calculate new roughness lengths (for water surfaces only)
1508       CALL calc_z0_water_surface
1509
1510
1511    END SUBROUTINE lsm_energy_balance
1512
1513!------------------------------------------------------------------------------!
1514! Description:
1515! ------------
1516!> Header output for land surface model
1517!------------------------------------------------------------------------------!
1518    SUBROUTINE lsm_header ( io )
1519
1520
1521       IMPLICIT NONE
1522
1523       CHARACTER (LEN=86) ::  t_soil_chr          !< String for soil temperature profile
1524       CHARACTER (LEN=86) ::  roots_chr           !< String for root profile
1525       CHARACTER (LEN=86) ::  vertical_index_chr  !< String for the vertical index
1526       CHARACTER (LEN=86) ::  m_soil_chr          !< String for soil moisture
1527       CHARACTER (LEN=86) ::  soil_depth_chr      !< String for soil depth
1528       CHARACTER (LEN=10) ::  coor_chr            !< Temporary string
1529   
1530       INTEGER(iwp) ::  i                         !< Loop index over soil layers
1531 
1532       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1533 
1534       t_soil_chr = ''
1535       m_soil_chr    = ''
1536       soil_depth_chr  = '' 
1537       roots_chr        = '' 
1538       vertical_index_chr   = ''
1539
1540       i = 1
1541       DO i = nzb_soil, nzt_soil
1542          WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
1543          t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr )
1544
1545          WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
1546          m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr )
1547
1548          WRITE (coor_chr,'(F10.2,7X)')  - zs(i)
1549          soil_depth_chr = TRIM( soil_depth_chr ) // ' '  // TRIM( coor_chr )
1550
1551          WRITE (coor_chr,'(F10.2,7X)')  root_fraction(i)
1552          roots_chr = TRIM( roots_chr ) // ' '  // TRIM( coor_chr )
1553
1554          WRITE (coor_chr,'(I10,7X)')  i
1555          vertical_index_chr = TRIM( vertical_index_chr ) // ' '  //           &
1556                               TRIM( coor_chr )
1557       ENDDO
1558
1559!
1560!--    Write land surface model header
1561       WRITE( io,  1 )
1562       IF ( conserve_water_content )  THEN
1563          WRITE( io, 2 )
1564       ELSE
1565          WRITE( io, 3 )
1566       ENDIF
1567
1568       WRITE( io, 4 ) TRIM( veg_type_name(veg_type) ),                         &
1569                        TRIM (soil_type_name(soil_type) )
1570       WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
1571                        TRIM( m_soil_chr ), TRIM( roots_chr ),                 &
1572                        TRIM( vertical_index_chr )
1573
15741   FORMAT (//' Land surface model information:'/                              &
1575              ' ------------------------------'/)
15762   FORMAT ('    --> Soil bottom is closed (water content is conserved',       &
1577            ', default)')
15783   FORMAT ('    --> Soil bottom is open (water content is not conserved)')         
15794   FORMAT ('    --> Land surface type  : ',A,/                                &
1580            '    --> Soil porosity type : ',A)
15815   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
1582            '       Height:        ',A,'  m'/                                  &
1583            '       Temperature:   ',A,'  K'/                                  &
1584            '       Moisture:      ',A,'  m**3/m**3'/                          &
1585            '       Root fraction: ',A,'  '/                                   &
1586            '       Grid point:    ',A)
1587
1588    END SUBROUTINE lsm_header
1589
1590
1591!------------------------------------------------------------------------------!
1592! Description:
1593! ------------
1594!> Initialization of the land surface model
1595!------------------------------------------------------------------------------!
1596    SUBROUTINE lsm_init
1597   
1598
1599       IMPLICIT NONE
1600
1601       INTEGER(iwp) ::  i !< running index
1602       INTEGER(iwp) ::  j !< running index
1603       INTEGER(iwp) ::  k !< running index
1604
1605       REAL(wp) :: pt1   !< potential temperature at first grid level
1606
1607
1608!
1609!--    Calculate Exner function
1610       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
1611
1612
1613!
1614!--    If no cloud physics is used, rho_surface has not been calculated before
1615       IF (  .NOT.  cloud_physics )  THEN
1616          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
1617       ENDIF
1618
1619!
1620!--    Calculate frequently used parameters
1621       rho_cp    = cp * rho_surface
1622       rd_d_rv   = r_d / r_v
1623       rho_lv    = rho_surface * l_v
1624       drho_l_lv = 1.0_wp / (rho_l * l_v)
1625
1626!
1627!--    Set inital values for prognostic quantities
1628       tt_surface_m = 0.0_wp
1629       tt_soil_m    = 0.0_wp
1630       tm_soil_m    = 0.0_wp
1631       tm_liq_eb_m  = 0.0_wp
1632       c_liq        = 0.0_wp
1633
1634       ghf_eb = 0.0_wp
1635       shf_eb = rho_cp * shf
1636
1637       IF ( humidity )  THEN
1638          qsws_eb = rho_lv * qsws
1639       ELSE
1640          qsws_eb = 0.0_wp
1641       ENDIF
1642
1643       qsws_liq_eb  = 0.0_wp
1644       qsws_soil_eb = 0.0_wp
1645       qsws_veg_eb  = 0.0_wp
1646
1647       r_a        = 50.0_wp
1648       r_s        = 50.0_wp
1649       r_canopy   = 0.0_wp
1650       r_soil     = 0.0_wp
1651
1652!
1653!--    Allocate 3D soil model arrays
1654       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1655       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1656       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1657
1658       lambda_h = 0.0_wp
1659!
1660!--    If required, allocate humidity-related variables for the soil model
1661       IF ( humidity )  THEN
1662          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1663          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
1664
1665          lambda_w = 0.0_wp 
1666       ENDIF
1667
1668!
1669!--    Calculate grid spacings. Temperature and moisture are defined at
1670!--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
1671!--    at the centers
1672       dz_soil(nzb_soil) = zs(nzb_soil)
1673
1674       DO  k = nzb_soil+1, nzt_soil
1675          dz_soil(k) = zs(k) - zs(k-1)
1676       ENDDO
1677       dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
1678
1679       DO  k = nzb_soil, nzt_soil-1
1680          dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
1681       ENDDO
1682       dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
1683
1684       ddz_soil      = 1.0_wp / dz_soil
1685       ddz_soil_stag = 1.0_wp / dz_soil_stag
1686
1687!
1688!--    Initialize standard soil types. It is possible to overwrite each
1689!--    parameter by setting the respecticy NAMELIST variable to a
1690!--    value /= 9999999.9.
1691       IF ( soil_type /= 0 )  THEN 
1692 
1693          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
1694             alpha_vangenuchten = soil_pars(0,soil_type)
1695          ENDIF
1696
1697          IF ( l_vangenuchten == 9999999.9_wp )  THEN
1698             l_vangenuchten = soil_pars(1,soil_type)
1699          ENDIF
1700
1701          IF ( n_vangenuchten == 9999999.9_wp )  THEN
1702             n_vangenuchten = soil_pars(2,soil_type)           
1703          ENDIF
1704
1705          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
1706             hydraulic_conductivity = soil_pars(3,soil_type)           
1707          ENDIF
1708
1709          IF ( saturation_moisture == 9999999.9_wp )  THEN
1710             saturation_moisture = m_soil_pars(0,soil_type)           
1711          ENDIF
1712
1713          IF ( field_capacity == 9999999.9_wp )  THEN
1714             field_capacity = m_soil_pars(1,soil_type)           
1715          ENDIF
1716
1717          IF ( wilting_point == 9999999.9_wp )  THEN
1718             wilting_point = m_soil_pars(2,soil_type)           
1719          ENDIF
1720
1721          IF ( residual_moisture == 9999999.9_wp )  THEN
1722             residual_moisture = m_soil_pars(3,soil_type)       
1723          ENDIF
1724
1725       ENDIF   
1726
1727!
1728!--    Map values to the respective 2D arrays
1729       alpha_vg      = alpha_vangenuchten
1730       l_vg          = l_vangenuchten
1731       n_vg          = n_vangenuchten 
1732       gamma_w_sat   = hydraulic_conductivity
1733       m_sat         = saturation_moisture
1734       m_fc          = field_capacity
1735       m_wilt        = wilting_point
1736       m_res         = residual_moisture
1737       r_soil_min    = min_soil_resistance
1738
1739!
1740!--    Initial run actions
1741       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1742
1743          t_soil    = 0.0_wp
1744          m_liq_eb  = 0.0_wp
1745          m_soil    = 0.0_wp
1746
1747!
1748!--       Map user settings of T and q for each soil layer
1749!--       (make sure that the soil moisture does not drop below the permanent
1750!--       wilting point) -> problems with devision by zero)
1751          DO  k = nzb_soil, nzt_soil
1752             t_soil(k,:,:)    = soil_temperature(k)
1753             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
1754             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
1755          ENDDO
1756          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
1757
1758!
1759!--       Calculate surface temperature
1760          t_surface   = pt_surface * exn
1761
1762!
1763!--       Set artifical values for ts and us so that r_a has its initial value
1764!--       for the first time step
1765          DO  i = nxlg, nxrg
1766             DO  j = nysg, nyng
1767                k = nzb_s_inner(j,i)
1768
1769                IF ( cloud_physics )  THEN
1770                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
1771                ELSE
1772                   pt1 = pt(k+1,j,i)
1773                ENDIF
1774
1775!
1776!--             Assure that r_a cannot be zero at model start
1777                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
1778
1779                us(j,i)  = 0.1_wp
1780                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
1781                shf(j,i) = - us(j,i) * ts(j,i)
1782             ENDDO
1783          ENDDO
1784
1785!
1786!--    Actions for restart runs
1787       ELSE
1788
1789          DO  i = nxlg, nxrg
1790             DO  j = nysg, nyng
1791                k = nzb_s_inner(j,i)               
1792                t_surface(j,i) = pt(k,j,i) * exn
1793             ENDDO
1794          ENDDO
1795
1796       ENDIF
1797
1798       DO  k = nzb_soil, nzt_soil
1799          root_fr(k,:,:) = root_fraction(k)
1800       ENDDO
1801
1802       IF ( veg_type /= 0 )  THEN
1803          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
1804             min_canopy_resistance = veg_pars(0,veg_type)
1805          ENDIF
1806          IF ( leaf_area_index == 9999999.9_wp )  THEN
1807             leaf_area_index = veg_pars(1,veg_type)         
1808          ENDIF
1809          IF ( vegetation_coverage == 9999999.9_wp )  THEN
1810             vegetation_coverage = veg_pars(2,veg_type)     
1811          ENDIF
1812          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
1813              canopy_resistance_coefficient= veg_pars(3,veg_type)     
1814          ENDIF
1815          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
1816             lambda_surface_stable = surface_pars(0,veg_type)         
1817          ENDIF
1818          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
1819             lambda_surface_unstable = surface_pars(1,veg_type)       
1820          ENDIF
1821          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
1822             f_shortwave_incoming = surface_pars(2,veg_type)       
1823          ENDIF
1824          IF ( z0_eb == 9999999.9_wp )  THEN
1825             roughness_length = roughness_par(0,veg_type) 
1826             z0_eb            = roughness_par(0,veg_type) 
1827          ENDIF
1828          IF ( z0h_eb == 9999999.9_wp )  THEN
1829             z0h_eb = roughness_par(1,veg_type)
1830          ENDIF
1831          IF ( z0q_eb == 9999999.9_wp )  THEN
1832             z0q_eb = roughness_par(2,veg_type)
1833          ENDIF
1834          z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp )
1835
1836          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
1837             DO  k = nzb_soil, nzt_soil
1838                root_fr(k,:,:) = root_distribution(k,veg_type)
1839                root_fraction(k) = root_distribution(k,veg_type)
1840             ENDDO
1841          ENDIF
1842
1843       ELSE
1844
1845          IF ( z0_eb == 9999999.9_wp )  THEN
1846             z0_eb = roughness_length
1847          ENDIF
1848          IF ( z0h_eb == 9999999.9_wp )  THEN
1849             z0h_eb = z0_eb * z0h_factor
1850          ENDIF
1851          IF ( z0q_eb == 9999999.9_wp )  THEN
1852             z0q_eb = z0_eb * z0h_factor
1853          ENDIF
1854
1855       ENDIF
1856
1857!
1858!--    For surfaces covered with pavement, set depth of the pavement (with dry
1859!--    soil below). The depth must be greater than the first soil layer depth
1860       IF ( veg_type == 20 )  THEN
1861          IF ( pave_depth == 9999999.9_wp )  THEN
1862             pave_depth = zs(nzb_soil) 
1863          ELSE
1864             pave_depth = MAX( zs(nzb_soil), pave_depth )
1865          ENDIF
1866       ENDIF
1867
1868!
1869!--    Map vegetation and soil types to 2D array to allow for heterogeneous
1870!--    surfaces via user interface see below
1871       veg_type_2d = veg_type
1872       soil_type_2d = soil_type
1873
1874!
1875!--    Map vegetation parameters to the respective 2D arrays
1876       r_canopy_min         = min_canopy_resistance
1877       lai                  = leaf_area_index
1878       c_veg                = vegetation_coverage
1879       g_d                  = canopy_resistance_coefficient
1880       lambda_surface_s     = lambda_surface_stable
1881       lambda_surface_u     = lambda_surface_unstable
1882       f_sw_in              = f_shortwave_incoming
1883       z0                   = z0_eb
1884       z0h                  = z0h_eb
1885       z0q                  = z0q_eb
1886
1887!
1888!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
1889       CALL user_init_land_surface
1890
1891!
1892!--    Set flag parameter if vegetation type was set to a water surface. Also
1893!--    set temperature to a constant value in all "soil" layers.
1894       DO  i = nxlg, nxrg
1895          DO  j = nysg, nyng
1896             IF ( veg_type_2d(j,i) == 14  .OR.  veg_type_2d(j,i) == 15 )  THEN
1897                water_surface(j,i) = .TRUE.
1898                t_soil(:,j,i) = t_surface(j,i)
1899             ELSEIF ( veg_type_2d(j,i) == 20 )  THEN
1900                pave_surface(j,i) = .TRUE.
1901                m_soil(:,j,i) = 0.0_wp
1902             ENDIF
1903
1904          ENDDO
1905       ENDDO
1906
1907!
1908!--    Calculate new roughness lengths (for water surfaces only)
1909       CALL calc_z0_water_surface
1910
1911       t_soil_p    = t_soil
1912       m_soil_p    = m_soil
1913       m_liq_eb_p  = m_liq_eb
1914       t_surface_p = t_surface
1915
1916
1917
1918!--    Store initial profiles of t_soil and m_soil (assuming they are
1919!--    horizontally homogeneous on this PE)
1920       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
1921                                                nysg,nxlg), 2,                 &
1922                                                statistic_regions+1 )
1923       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
1924                                                nysg,nxlg), 2,                 &
1925                                                statistic_regions+1 )
1926
1927    END SUBROUTINE lsm_init
1928
1929
1930!------------------------------------------------------------------------------!
1931! Description:
1932! ------------
1933!> Allocate land surface model arrays and define pointers
1934!------------------------------------------------------------------------------!
1935    SUBROUTINE lsm_init_arrays
1936   
1937
1938       IMPLICIT NONE
1939
1940!
1941!--    Allocate surface and soil temperature / humidity
1942#if defined( __nopointer )
1943       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
1944       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
1945       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1946       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1947       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
1948       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
1949       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
1950       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
1951#else
1952       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
1953       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
1954       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1955       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1956       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
1957       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
1958       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
1959       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
1960#endif
1961
1962!
1963!--    Allocate intermediate timestep arrays
1964       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
1965       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1966       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
1967       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1968
1969!
1970!--    Allocate 2D vegetation model arrays
1971       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
1972       ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) )
1973       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
1974       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
1975       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
1976       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
1977       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
1978       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
1979       ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
1980       ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
1981       ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
1982       ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
1983       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
1984       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
1985       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
1986       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
1987       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
1988       ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) )
1989       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
1990       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
1991       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
1992       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
1993       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
1994       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
1995       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
1996       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
1997       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
1998       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
1999       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
2000       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
2001       ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) )
2002       ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) )
2003       ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) )
2004
2005#if ! defined( __nopointer )
2006!
2007!--    Initial assignment of the pointers
2008       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
2009       t_surface => t_surface_1; t_surface_p => t_surface_2
2010       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
2011       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
2012#endif
2013
2014
2015    END SUBROUTINE lsm_init_arrays
2016
2017
2018!------------------------------------------------------------------------------!
2019! Description:
2020! ------------
2021!> Parin for &lsmpar for land surface model
2022!------------------------------------------------------------------------------!
2023    SUBROUTINE lsm_parin
2024
2025
2026       IMPLICIT NONE
2027
2028       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
2029       
2030       NAMELIST /lsm_par/         alpha_vangenuchten, c_surface,               &
2031                                  canopy_resistance_coefficient,               &
2032                                  conserve_water_content,                      &
2033                                  f_shortwave_incoming, field_capacity,        & 
2034                                  hydraulic_conductivity,                      &
2035                                  lambda_surface_stable,                       &
2036                                  lambda_surface_unstable, leaf_area_index,    &
2037                                  l_vangenuchten, min_canopy_resistance,       &
2038                                  min_soil_resistance, n_vangenuchten,         &
2039                                  pave_depth, pave_heat_capacity,              &
2040                                  pave_heat_conductivity,                      &
2041                                  residual_moisture, root_fraction,            &
2042                                  saturation_moisture, skip_time_do_lsm,       &
2043                                  soil_moisture, soil_temperature, soil_type,  &
2044                                  vegetation_coverage, veg_type, wilting_point,& 
2045                                  zs, z0_eb, z0h_eb, z0q_eb
2046       
2047       line = ' '
2048       
2049!
2050!--    Try to find land surface model package
2051       REWIND ( 11 )
2052       line = ' '
2053       DO   WHILE ( INDEX( line, '&lsm_par' ) == 0 )
2054          READ ( 11, '(A)', END=10 )  line
2055       ENDDO
2056       BACKSPACE ( 11 )
2057
2058!
2059!--    Read user-defined namelist
2060       READ ( 11, lsm_par )
2061
2062!
2063!--    Set flag that indicates that the land surface model is switched on
2064       land_surface = .TRUE.
2065
2066 10    CONTINUE
2067       
2068
2069    END SUBROUTINE lsm_parin
2070
2071
2072!------------------------------------------------------------------------------!
2073! Description:
2074! ------------
2075!> Soil model as part of the land surface model. The model predicts soil
2076!> temperature and water content.
2077!------------------------------------------------------------------------------!
2078    SUBROUTINE lsm_soil_model
2079
2080
2081       IMPLICIT NONE
2082
2083       INTEGER(iwp) ::  i   !< running index
2084       INTEGER(iwp) ::  j   !< running index
2085       INTEGER(iwp) ::  k   !< running index
2086
2087       REAL(wp)     :: h_vg !< Van Genuchten coef. h
2088
2089       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
2090                                                 lambda_temp, & !< temp. lambda
2091                                                 tend           !< tendency
2092
2093       DO  i = nxlg, nxrg
2094          DO  j = nysg, nyng
2095
2096             IF ( pave_surface(j,i) )  THEN
2097                rho_c_total(nzb_soil,j,i) = pave_heat_capacity
2098                lambda_temp(nzb_soil)     = pave_heat_conductivity
2099             ENDIF
2100
2101             IF (  .NOT.  water_surface(j,i) )  THEN
2102                DO  k = nzb_soil, nzt_soil
2103
2104
2105                   IF ( pave_surface(j,i)  .AND.  zs(k) <= pave_depth )  THEN
2106                   
2107                      rho_c_total(k,j,i) = pave_heat_capacity
2108                      lambda_temp(k)     = pave_heat_conductivity   
2109
2110                   ELSE           
2111!
2112!--                   Calculate volumetric heat capacity of the soil, taking
2113!--                   into account water content
2114                      rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
2115                                           + rho_c_water * m_soil(k,j,i))
2116
2117!
2118!--                   Calculate soil heat conductivity at the center of the soil
2119!--                   layers
2120                      lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *    &
2121                                     lambda_h_water ** m_soil(k,j,i)
2122
2123                      ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i)             &
2124                                                     / m_sat(j,i)))
2125
2126                      lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
2127                                       lambda_h_dry
2128                   ENDIF
2129
2130                ENDDO
2131
2132!
2133!--             Calculate soil heat conductivity (lambda_h) at the _stag level
2134!--             using linear interpolation. For pavement surface, the
2135!--             true pavement depth is considered
2136                DO  k = nzb_soil, nzt_soil-1
2137                   IF ( pave_surface(j,i)  .AND.  zs(k)   < pave_depth         &
2138                                           .AND.  zs(k+1) > pave_depth )  THEN
2139                      lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1)  &
2140                                        * lambda_temp(k)                       &
2141                                        + ( 1.0_wp - ( pave_depth - zs(k) )    &
2142                                        / dz_soil(k+1) ) * lambda_temp(k+1)
2143                   ELSE
2144                      lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
2145                                        * 0.5_wp
2146                   ENDIF
2147                ENDDO
2148                lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
2149
2150
2151
2152
2153!
2154!--             Prognostic equation for soil temperature t_soil
2155                tend(:) = 0.0_wp
2156
2157                tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *          &
2158                          ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)  &
2159                            - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)    &
2160                            + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
2161
2162                DO  k = nzb_soil+1, nzt_soil
2163                   tend(k) = (1.0_wp/rho_c_total(k,j,i))                       &
2164                             * (   lambda_h(k,j,i)                             &
2165                                 * ( t_soil(k+1,j,i) - t_soil(k,j,i) )         &
2166                                 * ddz_soil(k+1)                               &
2167                                 - lambda_h(k-1,j,i)                           &
2168                                 * ( t_soil(k,j,i) - t_soil(k-1,j,i) )         &
2169                                 * ddz_soil(k)                                 &
2170                               ) * ddz_soil_stag(k)
2171
2172                ENDDO
2173
2174                t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)&
2175                                                  + dt_3d * ( tsc(2)           &
2176                                                  * tend(nzb_soil:nzt_soil)    & 
2177                                                  + tsc(3)                     &
2178                                                  * tt_soil_m(:,j,i) )   
2179
2180!
2181!--             Calculate t_soil tendencies for the next Runge-Kutta step
2182                IF ( timestep_scheme(1:5) == 'runge' )  THEN
2183                   IF ( intermediate_timestep_count == 1 )  THEN
2184                      DO  k = nzb_soil, nzt_soil
2185                         tt_soil_m(k,j,i) = tend(k)
2186                      ENDDO
2187                   ELSEIF ( intermediate_timestep_count <                      &
2188                            intermediate_timestep_count_max )  THEN
2189                      DO  k = nzb_soil, nzt_soil
2190                         tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
2191                                         * tt_soil_m(k,j,i)
2192                      ENDDO
2193                   ENDIF
2194                ENDIF
2195
2196
2197                DO  k = nzb_soil, nzt_soil
2198
2199!
2200!--                Calculate soil diffusivity at the center of the soil layers
2201                   lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat       &
2202                                     / m_sat(j,i) ) * ( MAX( m_soil(k,j,i),    &
2203                                     m_wilt(j,i) ) / m_sat(j,i) )**(           &
2204                                     b_ch + 2.0_wp )
2205
2206!
2207!--                Parametrization of Van Genuchten
2208                   IF ( soil_type /= 7 )  THEN
2209!
2210!--                   Calculate the hydraulic conductivity after Van Genuchten
2211!--                   (1980)
2212                      h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -    &
2213                                 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**(      &
2214                                 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp  &
2215                             )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i)
2216
2217
2218                      gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +         &
2219                                      ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**(  &
2220                                      1.0_wp - 1.0_wp / n_vg(j,i) ) - (        &
2221                                      alpha_vg(j,i) * h_vg )**( n_vg(j,i)      &
2222                                      - 1.0_wp) )**2 )                         &
2223                                      / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg    &
2224                                      )**n_vg(j,i) )**( ( 1.0_wp  - 1.0_wp     &
2225                                      / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) )
2226
2227!
2228!--                Parametrization of Clapp & Hornberger
2229                   ELSE
2230                      gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i)       &
2231                                      / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
2232                   ENDIF
2233
2234                ENDDO
2235
2236!
2237!--             Prognostic equation for soil moisture content. Only performed,
2238!--             when humidity is enabled in the atmosphere and the surface type
2239!--             is not pavement (implies dry soil below).
2240                IF ( humidity  .AND.  .NOT.  pave_surface(j,i) )  THEN
2241!
2242!--                Calculate soil diffusivity (lambda_w) at the _stag level
2243!--                using linear interpolation. To do: replace this with
2244!--                ECMWF-IFS Eq. 8.81
2245                   DO  k = nzb_soil, nzt_soil-1
2246                     
2247                      lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
2248                                        * 0.5_wp
2249                      gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )    &
2250                                        * 0.5_wp
2251
2252                   ENDDO
2253
2254!
2255!
2256!--                In case of a closed bottom (= water content is conserved),
2257!--                set hydraulic conductivity to zero to that no water will be
2258!--                lost in the bottom layer.
2259                   IF ( conserve_water_content )  THEN
2260                      gamma_w(nzt_soil,j,i) = 0.0_wp
2261                   ELSE
2262                      gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
2263                   ENDIF     
2264
2265!--                The root extraction (= root_extr * qsws_veg_eb / (rho_l     
2266!--                * l_v)) ensures the mass conservation for water. The         
2267!--                transpiration of plants equals the cumulative withdrawals by
2268!--                the roots in the soil. The scheme takes into account the
2269!--                availability of water in the soil layers as well as the root
2270!--                fraction in the respective layer. Layer with moisture below
2271!--                wilting point will not contribute, which reflects the
2272!--                preference of plants to take water from moister layers.
2273
2274!
2275!--                Calculate the root extraction (ECMWF 7.69, the sum of
2276!--                root_extr = 1). The energy balance solver guarantees a
2277!--                positive transpiration, so that there is no need for an
2278!--                additional check.
2279                   DO  k = nzb_soil, nzt_soil
2280                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
2281                          m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
2282                       ENDIF
2283                   ENDDO 
2284
2285                   IF ( m_total > 0.0_wp )  THEN
2286                      DO  k = nzb_soil, nzt_soil
2287                         IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
2288                            root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)      &
2289                                                            / m_total
2290                         ELSE
2291                            root_extr(k) = 0.0_wp
2292                         ENDIF
2293                      ENDDO
2294                   ENDIF
2295
2296!
2297!--                Prognostic equation for soil water content m_soil.
2298                   tend(:) = 0.0_wp
2299
2300                   tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (               &
2301                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
2302                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
2303                               root_extr(nzb_soil) * qsws_veg_eb(j,i)          &
2304                               + qsws_soil_eb(j,i) ) * drho_l_lv )             &
2305                               * ddz_soil_stag(nzb_soil)
2306
2307                   DO  k = nzb_soil+1, nzt_soil-1
2308                      tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)          &
2309                                - m_soil(k,j,i) ) * ddz_soil(k+1)              &
2310                                - gamma_w(k,j,i)                               &
2311                                - lambda_w(k-1,j,i) * (m_soil(k,j,i) -         &
2312                                m_soil(k-1,j,i)) * ddz_soil(k)                 &
2313                                + gamma_w(k-1,j,i) - (root_extr(k)             &
2314                                * qsws_veg_eb(j,i) * drho_l_lv)                &
2315                                ) * ddz_soil_stag(k)
2316
2317                   ENDDO
2318                   tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                  &
2319                                           - lambda_w(nzt_soil-1,j,i)          &
2320                                           * (m_soil(nzt_soil,j,i)             &
2321                                           - m_soil(nzt_soil-1,j,i))           &
2322                                           * ddz_soil(nzt_soil)                &
2323                                           + gamma_w(nzt_soil-1,j,i) - (       &
2324                                             root_extr(nzt_soil)               &
2325                                           * qsws_veg_eb(j,i) * drho_l_lv  )   &
2326                                     ) * ddz_soil_stag(nzt_soil)             
2327
2328                   m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
2329                                                   + dt_3d * ( tsc(2) * tend(:)   &
2330                                                   + tsc(3) * tm_soil_m(:,j,i) )   
2331   
2332!
2333!--                Account for dry soils (find a better solution here!)
2334                   DO  k = nzb_soil, nzt_soil
2335                      IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
2336                   ENDDO
2337
2338!
2339!--                Calculate m_soil tendencies for the next Runge-Kutta step
2340                   IF ( timestep_scheme(1:5) == 'runge' )  THEN
2341                      IF ( intermediate_timestep_count == 1 )  THEN
2342                         DO  k = nzb_soil, nzt_soil
2343                            tm_soil_m(k,j,i) = tend(k)
2344                         ENDDO
2345                      ELSEIF ( intermediate_timestep_count <                   &
2346                               intermediate_timestep_count_max )  THEN
2347                         DO  k = nzb_soil, nzt_soil
2348                            tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp&
2349                                     * tm_soil_m(k,j,i)
2350                         ENDDO
2351                      ENDIF
2352                   ENDIF
2353
2354                ENDIF
2355
2356             ENDIF
2357
2358          ENDDO
2359       ENDDO
2360
2361    END SUBROUTINE lsm_soil_model
2362
2363 
2364!------------------------------------------------------------------------------!
2365! Description:
2366! ------------
2367!> Swapping of timelevels
2368!------------------------------------------------------------------------------!
2369    SUBROUTINE lsm_swap_timelevel ( mod_count )
2370
2371       IMPLICIT NONE
2372
2373       INTEGER, INTENT(IN) :: mod_count
2374
2375#if defined( __nopointer )
2376
2377       t_surface    = t_surface_p
2378       t_soil       = t_soil_p
2379       IF ( humidity )  THEN
2380          m_soil    = m_soil_p
2381          m_liq_eb  = m_liq_eb_p
2382       ENDIF
2383
2384#else
2385   
2386       SELECT CASE ( mod_count )
2387
2388          CASE ( 0 )
2389
2390             t_surface  => t_surface_1; t_surface_p  => t_surface_2
2391             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
2392             IF ( humidity )  THEN
2393                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
2394                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
2395             ENDIF
2396
2397
2398          CASE ( 1 )
2399
2400             t_surface  => t_surface_2; t_surface_p  => t_surface_1
2401             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
2402             IF ( humidity )  THEN
2403                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
2404                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
2405             ENDIF
2406
2407       END SELECT
2408#endif
2409
2410    END SUBROUTINE lsm_swap_timelevel
2411
2412
2413!------------------------------------------------------------------------------!
2414! Description:
2415! ------------
2416!> Calculation of roughness length for open water (lakes, ocean). The
2417!> parameterization follows Charnock (1955). Two different implementations
2418!> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al.
2419!> 2012)
2420!------------------------------------------------------------------------------!
2421    SUBROUTINE calc_z0_water_surface
2422
2423       USE control_parameters,                                                 &
2424           ONLY: g, kappa, molecular_viscosity
2425
2426       IMPLICIT NONE
2427
2428       INTEGER :: i  !< running index
2429       INTEGER :: j  !< running index
2430
2431       REAL(wp), PARAMETER :: alpha_ch  = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF
2432!       REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number)
2433!       REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization
2434!       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
2435
2436
2437       DO  i = nxlg, nxrg   
2438          DO  j = nysg, nyng
2439             IF ( water_surface(j,i) )  THEN
2440
2441!
2442!--             Disabled: FLake parameterization. Ideally, the Charnock
2443!--             coefficient should depend on the water depth and the fetch
2444!--             length
2445!                re_0 = z0(j,i) * us(j,i) / molecular_viscosity
2446!       
2447!                z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
2448!                              alpha_ch * us(j,i) / g )
2449!
2450!                z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
2451!                z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
2452
2453!
2454!--              Set minimum roughness length for u* > 0.2
2455!                IF ( us(j,i) > 0.2_wp )  THEN
2456!                   z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
2457!                   z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
2458!                ENDIF
2459
2460!
2461!--             ECMWF IFS model parameterization after Beljaars (1994). At low
2462!--             wind speed, the sea surface becomes aerodynamically smooth and
2463!--             the roughness scales with the viscosity. At high wind speed, the
2464!--             Charnock relation is used.
2465                z0(j,i) =   ( 0.11_wp * molecular_viscosity / us(j,i) )        &
2466                          + ( alpha_ch * us(j,i)**2 / g )
2467
2468                z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
2469                z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
2470
2471             ENDIF
2472          ENDDO
2473       ENDDO
2474
2475    END SUBROUTINE calc_z0_water_surface
2476
2477
2478!------------------------------------------------------------------------------!
2479! Description:
2480! ------------
2481!> Calculation of specific humidity of the skin layer (surface). It is assumend
2482!> that the skin is always saturated.
2483!------------------------------------------------------------------------------!
2484    SUBROUTINE calc_q_surface
2485
2486       IMPLICIT NONE
2487
2488       INTEGER :: i              !< running index
2489       INTEGER :: j              !< running index
2490       INTEGER :: k              !< running index
2491
2492       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
2493
2494       DO  i = nxlg, nxrg   
2495          DO  j = nysg, nyng
2496             k = nzb_s_inner(j,i)
2497
2498!
2499!--          Calculate water vapour pressure at saturation
2500             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i)   &
2501                   - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) )
2502
2503!
2504!--          Calculate specific humidity at saturation
2505             q_s = 0.622_wp * e_s / surface_pressure
2506
2507             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
2508
2509!
2510!--          Calculate specific humidity at surface
2511             IF ( cloud_physics )  THEN
2512                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
2513                             * ( q(k+1,j,i) - ql(k+1,j,i) )
2514             ELSE
2515                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
2516                             * q(k+1,j,i)
2517             ENDIF
2518
2519!
2520!--          Update virtual potential temperature
2521             vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
2522
2523          ENDDO
2524       ENDDO
2525
2526    END SUBROUTINE calc_q_surface
2527
2528
2529 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.