source: palm/trunk/SOURCE/land_surface_model.f90 @ 1683

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 71.5 KB
Line 
1!> @file land_surface_model.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-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: land_surface_model.f90 1683 2015-10-07 23:57:51Z knoop $
26! Bugfix: REAL constants provided with KIND-attribute in call of
27! 1682 2015-10-07 23:56:08Z knoop
28! Code annotations made doxygen readable
29!
30! 1590 2015-05-08 13:56:27Z maronga
31! Bugfix: definition of character strings requires same length for all elements
32!
33! 1585 2015-04-30 07:05:52Z maronga
34! Modifications for RRTMG. Changed tables to PARAMETER type.
35!
36! 1571 2015-03-12 16:12:49Z maronga
37! Removed upper-case variable names. Corrected distribution of precipitation to
38! the liquid water reservoir and the bare soil fractions.
39!
40! 1555 2015-03-04 17:44:27Z maronga
41! Added output of r_a and r_s
42!
43! 1553 2015-03-03 17:33:54Z maronga
44! Improved better treatment of roughness lengths. Added default soil temperature
45! profile
46!
47! 1551 2015-03-03 14:18:16Z maronga
48! Flux calculation is now done in prandtl_fluxes. Added support for data output.
49! Vertical indices have been replaced. Restart runs are now possible. Some
50! variables have beem renamed. Bugfix in the prognostic equation for the surface
51! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
52! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
53! the hydraulic conductivity. Bugfix for root fraction and extraction
54! calculation
55!
56! intrinsic function MAX and MIN
57!
58! 1500 2014-12-03 17:42:41Z maronga
59! Corrected calculation of aerodynamic resistance (r_a).
60! Precipitation is now added to liquid water reservoir using LE_liq.
61! Added support for dry runs.
62!
63! 1496 2014-12-02 17:25:50Z maronga
64! Initial revision
65!
66!
67! Description:
68! ------------
69!> Land surface model, consisting of a solver for the energy balance at the
70!> surface and a four layer soil scheme. The scheme is similar to the TESSEL
71!> scheme implemented in the ECMWF IFS model, with modifications according to
72!> H-TESSEL. The implementation is based on the formulation implemented in the
73!> DALES and UCLA-LES models.
74!>
75!> @todo Check dewfall parametrization for fog simulations.
76!> @todo Consider partial absorption of the net shortwave radiation by the surface layer.
77!> @todo Allow for water surfaces, check performance for bare soils.
78!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3)).
79!> @todo Implement surface runoff model (required when performing long-term LES with
80!>       considerable precipitation.
81!>
82!> @note  No time step criterion is required as long as the soil layers do not become
83!>        too thin.
84!------------------------------------------------------------------------------!
85 MODULE land_surface_model_mod
86 
87     USE arrays_3d,                                                            &
88         ONLY:  pt, pt_p, q_p, qsws, rif, shf, ts, us, z0, z0h
89
90     USE cloud_parameters,                                                     &
91         ONLY:  cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
92
93     USE control_parameters,                                                   &
94         ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,   &
95                initializing_actions, intermediate_timestep_count_max,         &
96                max_masks, precipitation, pt_surface, rho_surface,             &
97                roughness_length, surface_pressure, timestep_scheme, tsc,      &
98                z0h_factor, time_since_reference_point
99
100     USE indices,                                                              &
101         ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner
102
103     USE kinds
104
105    USE netcdf_control,                                                        &
106        ONLY:  dots_label, dots_num, dots_unit
107
108     USE radiation_model_mod,                                                  &
109         ONLY:  radiation_scheme, rad_net, rad_sw_in, sigma_sb
110
111     USE statistics,                                                           &
112         ONLY:  hom, statistic_regions
113
114    IMPLICIT NONE
115
116!
117!-- LSM model constants
118    INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !< bottom of the soil model (to be switched)
119                               nzt_soil = 3, & !< top of the soil model (to be switched)
120                               nzs = 4         !< number of soil layers (fixed for now)
121
122    INTEGER(iwp) :: dots_soil = 0  !< starting index for timeseries output
123
124    INTEGER(iwp), DIMENSION(0:1) :: id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,  &
125                                    id_dim_zs_3d, id_var_zs_xy,                & 
126                                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d
127                                   
128    INTEGER(iwp), DIMENSION(1:max_masks,0:1) :: id_dim_zs_mask, id_var_zs_mask
129
130    REAL(wp), PARAMETER ::                     &
131              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
132              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
133              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
134              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
135              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
136              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
137              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
138              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
139
140
141!
142!-- LSM variables
143    INTEGER(iwp) :: veg_type  = 2, & !< vegetation type, 0: user-defined, 1-19: generic (see list)
144                    soil_type = 3    !< soil type, 0: user-defined, 1-6: generic (see list)
145
146    LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model
147               dewfall = .TRUE.,                & !< allow/inhibit dewfall
148               land_surface = .FALSE.             !< flag parameter indicating wheather the lsm is used
149
150!   value 9999999.9_wp -> generic available or user-defined value must be set
151!   otherwise -> no generic variable and user setting is optional
152    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
153                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
154                c_surface   = 20000.0_wp,               & !< Surface (skin) heat capacity
155                drho_l_lv,                              & !< (rho_l * l_v)**-1
156                exn,                                    & !< value of the Exner function
157                e_s = 0.0_wp,                           & !< saturation water vapour pressure
158                field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
159                f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
160                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
161                ke = 0.0_wp,                            & !< Kersten number
162                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s
163                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u
164                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
165                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
166                min_canopy_resistance = 9999999.9_wp,   & !< NAMELIST r_canopy_min
167                min_soil_resistance = 50.0_wp,          & !< NAMELIST r_soil_min
168                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
169                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
170                q_s = 0.0_wp,                           & !< saturation specific humidity
171                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
172                rho_cp,                                 & !< rho_surface * cp
173                rho_lv,                                 & !< rho * l_v
174                rd_d_rv,                                & !< r_d / r_v
175                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
176                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
177                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
178                z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
179                z0h_eb = 9999999.9_wp                    !< NAMELIST z0h (lsm_par)
180
181    REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
182              ddz_soil,                     & !< 1/dz_soil
183              ddz_soil_stag,                & !< 1/dz_soil_stag
184              dz_soil,                      & !< soil grid spacing (center-center)
185              dz_soil_stag,                 & !< soil grid spacing (edge-edge)
186              root_extr = 0.0_wp,           & !< root extraction
187              root_fraction = (/9999999.9_wp, 9999999.9_wp,    &
188                                9999999.9_wp, 9999999.9_wp /), & !< distribution of root surface area to the individual soil layers
189              zs = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),    & !< soil layer depths (m)
190              soil_moisture = 0.0_wp          !< soil moisture content (m3/m3)
191
192    REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) ::   &
193              soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp,  283.0_wp,    &
194                                   283.0_wp /) !< soil temperature (K)
195
196#if defined( __nopointer )
197    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !< surface temperature (K)
198                                                     t_surface_p, & !< progn. surface temperature (K)
199                                                     m_liq_eb,    & !< liquid water reservoir (m)
200                                                     m_liq_eb_av, & !< liquid water reservoir (m)
201                                                     m_liq_eb_p     !< progn. liquid water reservoir (m)
202#else
203    REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
204                                         t_surface_p,    & 
205                                         m_liq_eb,       & 
206                                         m_liq_eb_p
207
208    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
209                                                     m_liq_eb_av,              &
210                                                     m_liq_eb_1, m_liq_eb_2
211#endif
212
213!
214!-- Temporal tendencies for time stepping           
215    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !< surface temperature tendency (K)
216                                             tm_liq_eb_m      !< liquid water reservoir tendency (m)
217
218!
219!-- Energy balance variables               
220    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
221              alpha_vg,         & !< coef. of Van Genuchten
222              c_liq,            & !< liquid water coverage (of vegetated area)
223              c_liq_av,         & !< average of c_liq
224              c_soil_av,        & !< average of c_soil
225              c_veg,            & !< vegetation coverage
226              c_veg_av,         & !< average of c_veg
227              f_sw_in,          & !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
228              ghf_eb,           & !< ground heat flux
229              ghf_eb_av,        & !< average of ghf_eb
230              gamma_w_sat,      & !< hydraulic conductivity at saturation
231              g_d,              & !< coefficient for dependence of r_canopy on water vapour pressure deficit
232              lai,              & !< leaf area index
233              lai_av,           & !< average of lai
234              lambda_h_sat,     & !< heat conductivity for dry soil
235              lambda_surface_s,    & !< coupling between surface and soil (depends on vegetation type)
236              lambda_surface_u,    & !< coupling between surface and soil (depends on vegetation type)
237              l_vg,             & !< coef. of Van Genuchten
238              m_fc,             & !< soil moisture at field capacity (m3/m3)
239              m_res,            & !< residual soil moisture
240              m_sat,            & !< saturation soil moisture (m3/m3)
241              m_wilt,           & !< soil moisture at permanent wilting point (m3/m3)
242              n_vg,             & !< coef. Van Genuchten 
243              qsws_eb,          & !< surface flux of latent heat (total)
244              qsws_eb_av,       & !< average of qsws_eb
245              qsws_liq_eb,      & !< surface flux of latent heat (liquid water portion)
246              qsws_liq_eb_av,   & !< average of qsws_liq_eb
247              qsws_soil_eb,     & !< surface flux of latent heat (soil portion)
248              qsws_soil_eb_av,  & !< average of qsws_soil_eb
249              qsws_veg_eb,      & !< surface flux of latent heat (vegetation portion)
250              qsws_veg_eb_av,   & !< average of qsws_veg_eb
251              rad_net_l,        & !< local copy of rad_net (net radiation at surface)
252              r_a,              & !< aerodynamic resistance
253              r_a_av,           & !< avergae of r_a
254              r_canopy,         & !< canopy resistance
255              r_soil,           & !< soil resitance
256              r_soil_min,       & !< minimum soil resistance
257              r_s,              & !< total surface resistance (combination of r_soil and r_canopy)
258              r_s_av,           & !< avergae of r_s
259              r_canopy_min,     & !< minimum canopy (stomatal) resistance
260              shf_eb,           & !< surface flux of sensible heat
261              shf_eb_av           !< average of shf_eb
262
263    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
264              lambda_h, &   !< heat conductivity of soil (?)                           
265              lambda_w, &   !< hydraulic diffusivity of soil (?)
266              gamma_w,  &   !< hydraulic conductivity of soil (?)
267              rho_c_total   !< volumetric heat capacity of the actual soil matrix (?)
268
269#if defined( __nopointer )
270    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
271              t_soil,    & !< Soil temperature (K)
272              t_soil_av, & !< Average of t_soil
273              t_soil_p,  & !< Prog. soil temperature (K)
274              m_soil,    & !< Soil moisture (m3/m3)
275              m_soil_av, & !< Average of m_soil
276              m_soil_p     !< Prog. soil moisture (m3/m3)
277#else
278    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
279              t_soil, t_soil_p, &
280              m_soil, m_soil_p   
281
282    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
283              t_soil_av, t_soil_1, t_soil_2,                                   &
284              m_soil_av, m_soil_1, m_soil_2
285#endif
286
287
288    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
289              tt_soil_m, & !< t_soil storage array
290              tm_soil_m, & !< m_soil storage array
291              root_fr      !< root fraction (sum=1)
292
293
294!
295!-- Predefined Land surface classes (veg_type)
296    CHARACTER(26), DIMENSION(0:19), PARAMETER :: veg_type_name = (/ &
297                                   'user defined              ',    & ! 0
298                                   'crops, mixed farming      ',    & !  1
299                                   'short grass               ',    & !  2
300                                   'evergreen needleleaf trees',    & !  3
301                                   'deciduous needleleaf trees',    & !  4
302                                   'evergreen broadleaf trees ',    & !  5
303                                   'deciduous broadleaf trees ',    & !  6
304                                   'tall grass                ',    & !  7
305                                   'desert                    ',    & !  8
306                                   'tundra                    ',    & !  9
307                                   'irrigated crops           ',    & ! 10
308                                   'semidesert                ',    & ! 11
309                                   'ice caps and glaciers     ',    & ! 12
310                                   'bogs and marshes          ',    & ! 13
311                                   'inland water              ',    & ! 14
312                                   'ocean                     ',    & ! 15
313                                   'evergreen shrubs          ',    & ! 16
314                                   'deciduous shrubs          ',    & ! 17
315                                   'mixed forest/woodland     ',    & ! 18
316                                   'interrupted forest        '     & ! 19
317                                                                 /)
318
319!
320!-- Soil model classes (soil_type)
321    CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
322                                   'user defined',                  & ! 0
323                                   'coarse      ',                  & ! 1
324                                   'medium      ',                  & ! 2
325                                   'medium-fine ',                  & ! 3
326                                   'fine        ',                  & ! 4
327                                   'very fine   ',                  & ! 5
328                                   'organic     ',                  & ! 6
329                                   'loamy (CH)  '                   & ! 7
330                                                                 /)
331!
332!-- Land surface parameters according to the respective classes (veg_type)
333
334!
335!-- Land surface parameters I
336!--                          r_canopy_min,     lai,   c_veg,     g_d
337    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: veg_pars = RESHAPE( (/ &
338                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
339                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
340                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
341                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
342                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
343                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
344                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
345                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
346                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
347                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
348                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
349                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
350                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
351                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
352                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
353                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
354                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
355                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
356                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp   & ! 19
357                                 /), (/ 4, 19 /) )
358
359!
360!-- Land surface parameters II          z0,         z0h
361    REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ & 
362                                   0.25_wp,  0.25E-2_wp,                    & !  1
363                                   0.20_wp,  0.20E-2_wp,                    & !  2
364                                   2.00_wp,     2.00_wp,                    & !  3
365                                   2.00_wp,     2.00_wp,                    & !  4
366                                   2.00_wp,     2.00_wp,                    & !  5
367                                   2.00_wp,     2.00_wp,                    & !  6
368                                   0.47_wp,  0.47E-2_wp,                    & !  7
369                                  0.013_wp, 0.013E-2_wp,                    & !  8
370                                  0.034_wp, 0.034E-2_wp,                    & !  9
371                                    0.5_wp,  0.50E-2_wp,                    & ! 10
372                                   0.17_wp,  0.17E-2_wp,                    & ! 11
373                                 1.3E-3_wp,   1.3E-4_wp,                    & ! 12
374                                   0.83_wp,  0.83E-2_wp,                    & ! 13
375                                   0.00_wp,  0.00E-2_wp,                    & ! 14
376                                   0.00_wp,  0.00E-2_wp,                    & ! 15
377                                   0.10_wp,  0.10E-2_wp,                    & ! 16
378                                   0.25_wp,  0.25E-2_wp,                    & ! 17
379                                   2.00_wp,  2.00E-2_wp,                    & ! 18
380                                   1.10_wp,  1.10E-2_wp                     & ! 19
381                                 /), (/ 2, 19 /) )
382
383!
384!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
385    REAL(wp), DIMENSION(0:2,1:19), PARAMETER :: surface_pars = RESHAPE( (/ &
386                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
387                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
388                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  3
389                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  4
390                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  5
391                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  6
392                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  7
393                                      15.0_wp,       15.0_wp, 0.00_wp,     & !  8
394                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  9
395                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 10
396                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
397                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
398                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13
399                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 14
400                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 15
401                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
402                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
403                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
404                                      20.0_wp,       15.0_wp, 0.03_wp      & ! 19
405                                      /), (/ 3, 19 /) )
406
407!
408!-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
409    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: root_distribution = RESHAPE( (/ &
410                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
411                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
412                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  3
413                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  4
414                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  5
415                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  6
416                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  7
417                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  8
418                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & !  9
419                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 10
420                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 11
421                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 12
422                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 13
423                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 14
424                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 15
425                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
426                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17
427                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
428                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 19
429                                 /), (/ 4, 19 /) )
430
431!
432!-- Soil parameters according to the following porosity classes (soil_type)
433
434!
435!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
436    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
437                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
438                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
439                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
440                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
441                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
442                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
443                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
444                                 /), (/ 4, 7 /) )
445
446!
447!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
448    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
449                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
450                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
451                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
452                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
453                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
454                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
455                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
456                                 /), (/ 4, 7 /) )
457
458
459    SAVE
460
461
462    PRIVATE
463
464
465!
466!-- Public parameters, constants and initial values
467    PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
468           conserve_water_content, dewfall, field_capacity,                    &
469           f_shortwave_incoming, hydraulic_conductivity, init_lsm,             &
470           init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable,    &
471           land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model,  &
472           lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance,          &
473           min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp,     &
474           rho_lv, root_fraction, saturation_moisture, soil_moisture,          &
475           soil_temperature, soil_type, soil_type_name, vegetation_coverage,   &
476           veg_type, veg_type_name, wilting_point, z0_eb, z0h_eb
477
478!
479!-- Public grid and NetCDF variables
480    PUBLIC dots_soil, id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,                &
481           id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy, id_var_zs_xz,           &
482           id_var_zs_yz, id_var_zs_3d, id_var_zs_mask, nzb_soil, nzs, nzt_soil,&
483           zs
484
485
486!
487!-- Public 2D output variables
488    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
489           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
490           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
491           r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
492
493
494!
495!-- Public prognostic variables
496    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
497
498    INTERFACE init_lsm
499       MODULE PROCEDURE init_lsm
500    END INTERFACE init_lsm
501
502    INTERFACE lsm_energy_balance
503       MODULE PROCEDURE lsm_energy_balance
504    END INTERFACE lsm_energy_balance
505
506    INTERFACE lsm_soil_model
507       MODULE PROCEDURE lsm_soil_model
508    END INTERFACE lsm_soil_model
509
510    INTERFACE lsm_swap_timelevel
511       MODULE PROCEDURE lsm_swap_timelevel
512    END INTERFACE lsm_swap_timelevel
513
514 CONTAINS
515
516
517!------------------------------------------------------------------------------!
518! Description:
519! ------------
520!> Allocate land surface model arrays and define pointers
521!------------------------------------------------------------------------------!
522    SUBROUTINE init_lsm_arrays
523   
524
525       IMPLICIT NONE
526
527!
528!--    Allocate surface and soil temperature / humidity
529#if defined( __nopointer )
530       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
531       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
532       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
533       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
534       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
535       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
536       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
537       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
538#else
539       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
540       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
541       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
542       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
543       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
544       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
545       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
546       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
547#endif
548
549!
550!--    Allocate intermediate timestep arrays
551       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
552       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
553       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
554       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
555
556!
557!--    Allocate 2D vegetation model arrays
558       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
559       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
560       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
561       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
562       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
563       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
564       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
565       ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
566       ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
567       ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
568       ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
569       ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
570       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
571       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
572       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
573       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
574       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
575       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
576       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
577       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
578       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
579       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
580       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
581       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
582       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
583       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
584       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
585       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
586       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
587
588#if ! defined( __nopointer )
589!
590!--    Initial assignment of the pointers
591       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
592       t_surface => t_surface_1; t_surface_p => t_surface_2
593       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
594       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
595#endif
596
597
598    END SUBROUTINE init_lsm_arrays
599
600!------------------------------------------------------------------------------!
601! Description:
602! ------------
603!> Initialization of the land surface model
604!------------------------------------------------------------------------------!
605    SUBROUTINE init_lsm
606   
607
608       IMPLICIT NONE
609
610       INTEGER(iwp) ::  i !< running index
611       INTEGER(iwp) ::  j !< running index
612       INTEGER(iwp) ::  k !< running index
613
614
615!
616!--    Calculate Exner function
617       exn       = ( surface_pressure / 1000.0_wp )**0.286_wp
618
619
620!
621!--    If no cloud physics is used, rho_surface has not been calculated before
622       IF ( .NOT. cloud_physics )  THEN
623          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
624       ENDIF
625
626!
627!--    Calculate frequently used parameters
628       rho_cp    = cp * rho_surface
629       rd_d_rv   = r_d / r_v
630       rho_lv    = rho_surface * l_v
631       drho_l_lv = 1.0_wp / (rho_l * l_v)
632
633!
634!--    Set inital values for prognostic quantities
635       tt_surface_m = 0.0_wp
636       tt_soil_m    = 0.0_wp
637       tm_liq_eb_m  = 0.0_wp
638       c_liq        = 0.0_wp
639
640       ghf_eb = 0.0_wp
641       shf_eb = rho_cp * shf
642
643       IF ( humidity )  THEN
644          qsws_eb = rho_l * l_v * qsws
645       ELSE
646          qsws_eb = 0.0_wp
647       ENDIF
648
649       qsws_liq_eb  = 0.0_wp
650       qsws_soil_eb = qsws_eb
651       qsws_veg_eb  = 0.0_wp
652
653       r_a        = 50.0_wp
654       r_s        = 50.0_wp
655       r_canopy   = 0.0_wp
656       r_soil     = 0.0_wp
657
658!
659!--    Allocate 3D soil model arrays
660       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
661       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
662       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
663
664       lambda_h = 0.0_wp
665!
666!--    If required, allocate humidity-related variables for the soil model
667       IF ( humidity )  THEN
668          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
669          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
670
671          lambda_w = 0.0_wp 
672       ENDIF
673
674!
675!--    Calculate grid spacings. Temperature and moisture are defined at
676!--    the center of the soil layers, whereas gradients/fluxes are defined
677!--    at the edges (_stag)
678       dz_soil_stag(nzb_soil) = zs(nzb_soil)
679
680       DO k = nzb_soil+1, nzt_soil
681          dz_soil_stag(k) = zs(k) - zs(k-1)
682       ENDDO
683
684       DO k = nzb_soil, nzt_soil-1
685          dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k))
686       ENDDO
687       dz_soil(nzt_soil) = dz_soil_stag(nzt_soil)
688
689       ddz_soil      = 1.0_wp / dz_soil
690       ddz_soil_stag = 1.0_wp / dz_soil_stag
691
692!
693!--    Initialize standard soil types. It is possible to overwrite each
694!--    parameter by setting the respecticy NAMELIST variable to a
695!--    value /= 9999999.9.
696       IF ( soil_type /= 0 )  THEN 
697 
698          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
699             alpha_vangenuchten = soil_pars(0,soil_type)
700          ENDIF
701
702          IF ( l_vangenuchten == 9999999.9_wp )  THEN
703             l_vangenuchten = soil_pars(1,soil_type)
704          ENDIF
705
706          IF ( n_vangenuchten == 9999999.9_wp )  THEN
707             n_vangenuchten = soil_pars(2,soil_type)           
708          ENDIF
709
710          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
711             hydraulic_conductivity = soil_pars(3,soil_type)           
712          ENDIF
713
714          IF ( saturation_moisture == 9999999.9_wp )  THEN
715             saturation_moisture = m_soil_pars(0,soil_type)           
716          ENDIF
717
718          IF ( field_capacity == 9999999.9_wp )  THEN
719             field_capacity = m_soil_pars(1,soil_type)           
720          ENDIF
721
722          IF ( wilting_point == 9999999.9_wp )  THEN
723             wilting_point = m_soil_pars(2,soil_type)           
724          ENDIF
725
726          IF ( residual_moisture == 9999999.9_wp )  THEN
727             residual_moisture = m_soil_pars(3,soil_type)       
728          ENDIF
729
730       ENDIF   
731
732       alpha_vg      = alpha_vangenuchten
733       l_vg          = l_vangenuchten
734       n_vg          = n_vangenuchten
735       gamma_w_sat   = hydraulic_conductivity
736       m_sat         = saturation_moisture
737       m_fc          = field_capacity
738       m_wilt        = wilting_point
739       m_res         = residual_moisture
740       r_soil_min    = min_soil_resistance
741
742!
743!--    Initial run actions
744       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
745
746          t_soil    = 0.0_wp
747          m_liq_eb  = 0.0_wp
748          m_soil    = 0.0_wp
749
750!
751!--       Map user settings of T and q for each soil layer
752!--       (make sure that the soil moisture does not drop below the permanent
753!--       wilting point) -> problems with devision by zero)
754          DO k = nzb_soil, nzt_soil
755             t_soil(k,:,:)    = soil_temperature(k)
756             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
757             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
758          ENDDO
759          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
760
761!
762!--       Calculate surface temperature
763          t_surface = pt_surface * exn
764
765!
766!--       Set artifical values for ts and us so that r_a has its initial value for
767!--       the first time step
768          DO  i = nxlg, nxrg
769             DO  j = nysg, nyng
770                k = nzb_s_inner(j,i)
771
772!
773!--             Assure that r_a cannot be zero at model start
774                IF ( pt(k+1,j,i) == pt(k,j,i) )  pt(k+1,j,i) = pt(k+1,j,i)     &
775                                                 + 1.0E-10_wp
776
777                us(j,i) = 0.1_wp
778                ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
779                shf(j,i) = - us(j,i) * ts(j,i)
780             ENDDO
781          ENDDO
782
783!
784!--    Actions for restart runs
785       ELSE
786
787          DO  i = nxlg, nxrg
788             DO  j = nysg, nyng
789                k = nzb_s_inner(j,i)               
790                t_surface(j,i) = pt(k,j,i) * exn
791             ENDDO
792          ENDDO
793
794       ENDIF
795
796!
797!--    Calculate saturation soil heat conductivity
798       lambda_h_sat(:,:) = lambda_h_sm ** (1.0_wp - m_sat(:,:)) *              &
799                           lambda_h_water ** m_sat(:,:)
800
801
802
803
804       DO k = nzb_soil, nzt_soil
805          root_fr(k,:,:) = root_fraction(k)
806       ENDDO
807
808       IF ( veg_type /= 0 )  THEN
809          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
810             min_canopy_resistance = veg_pars(0,veg_type)
811          ENDIF
812          IF ( leaf_area_index == 9999999.9_wp )  THEN
813             leaf_area_index = veg_pars(1,veg_type)         
814          ENDIF
815          IF ( vegetation_coverage == 9999999.9_wp )  THEN
816             vegetation_coverage = veg_pars(2,veg_type)     
817          ENDIF
818          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
819              canopy_resistance_coefficient= veg_pars(3,veg_type)     
820          ENDIF
821          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
822             lambda_surface_stable = surface_pars(0,veg_type)         
823          ENDIF
824          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
825             lambda_surface_unstable = surface_pars(1,veg_type)       
826          ENDIF
827          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
828             f_shortwave_incoming = surface_pars(2,veg_type)       
829          ENDIF
830          IF ( z0_eb == 9999999.9_wp )  THEN
831             roughness_length = roughness_par(0,veg_type) 
832             z0_eb            = roughness_par(0,veg_type) 
833          ENDIF
834          IF ( z0h_eb == 9999999.9_wp )  THEN
835             z0h_eb = roughness_par(1,veg_type)
836          ENDIF
837          z0h_factor = z0h_eb / z0_eb
838
839          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
840             DO k = nzb_soil, nzt_soil
841                root_fr(k,:,:) = root_distribution(k,veg_type)
842                root_fraction(k) = root_distribution(k,veg_type)
843             ENDDO
844          ENDIF
845
846       ELSE
847
848          IF ( z0_eb == 9999999.9_wp )  THEN
849             z0_eb = roughness_length
850          ENDIF
851          IF ( z0h_eb == 9999999.9_wp )  THEN
852             z0h_eb = z0_eb * z0h_factor
853          ENDIF
854
855       ENDIF
856
857!
858!--    Initialize vegetation
859       r_canopy_min         = min_canopy_resistance
860       lai                  = leaf_area_index
861       c_veg                = vegetation_coverage
862       g_d                  = canopy_resistance_coefficient
863       lambda_surface_s     = lambda_surface_stable
864       lambda_surface_u     = lambda_surface_unstable
865       f_sw_in              = f_shortwave_incoming
866       z0                   = z0_eb
867       z0h                  = z0h_eb
868
869!
870!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
871       CALL user_init_land_surface
872
873
874       t_soil_p    = t_soil
875       m_soil_p    = m_soil
876       m_liq_eb_p  = m_liq_eb
877
878!--    Store initial profiles of t_soil and m_soil (assuming they are
879!--    horizontally homogeneous on this PE)
880       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
881                                                nysg,nxlg), 2,                 &
882                                                statistic_regions+1 )
883       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
884                                                nysg,nxlg), 2,                 &
885                                                statistic_regions+1 )
886
887!
888!--    Calculate humidity at the surface
889       IF ( humidity )  THEN
890          CALL calc_q_surface
891       ENDIF
892
893
894
895!
896!--    Add timeseries for land surface model
897
898       dots_label(dots_num+1) = "ghf_eb"
899       dots_label(dots_num+2) = "shf_eb"
900       dots_label(dots_num+3) = "qsws_eb"
901       dots_label(dots_num+4) = "qsws_liq_eb"
902       dots_label(dots_num+5) = "qsws_soil_eb"
903       dots_label(dots_num+6) = "qsws_veg_eb"
904       dots_unit(dots_num+1:dots_num+6) = "W/m2"
905
906       dots_label(dots_num+7) = "r_a"
907       dots_label(dots_num+8) = "r_s"
908       dots_unit(dots_num+7:dots_num+8) = "s/m"
909
910       dots_soil = dots_num + 1
911       dots_num  = dots_num + 8
912
913       RETURN
914
915    END SUBROUTINE init_lsm
916
917
918
919!------------------------------------------------------------------------------!
920! Description:
921! ------------
922!> Solver for the energy balance at the surface.
923!> @note Surface fluxes are calculated in the land surface model, but these are
924!>       not used in the atmospheric code. The fluxes are calculated afterwards in
925!>       prandtl_fluxes using the surface values of temperature and humidity as
926!>       provided by the land surface model. In this way, the fluxes in the land
927!>       surface model are not equal to the ones calculated in prandtl_fluxes
928!------------------------------------------------------------------------------!
929    SUBROUTINE lsm_energy_balance
930
931
932       IMPLICIT NONE
933
934       INTEGER(iwp) ::  i         !< running index
935       INTEGER(iwp) ::  j         !< running index
936       INTEGER(iwp) ::  k, ks     !< running index
937
938       REAL(wp) :: f1,          & !< resistance correction term 1
939                   f2,          & !< resistance correction term 2
940                   f3,          & !< resistance correction term 3
941                   m_min,       & !< minimum soil moisture
942                   T_1,         & !< actual temperature at first grid point
943                   e,           & !< water vapour pressure
944                   e_s,         & !< water vapour saturation pressure
945                   e_s_dt,      & !< derivate of e_s with respect to T
946                   tend,        & !< tendency
947                   dq_s_dt,     & !< derivate of q_s with respect to T
948                   coef_1,      & !< coef. for prognostic equation
949                   coef_2,      & !< coef. for prognostic equation
950                   f_qsws,      & !< factor for qsws_eb
951                   f_qsws_veg,  & !< factor for qsws_veg_eb
952                   f_qsws_soil, & !< factor for qsws_soil_eb
953                   f_qsws_liq,  & !< factor for qsws_liq_eb
954                   f_shf,       & !< factor for shf_eb
955                   lambda_surface, & !< Current value of lambda_surface
956                   m_liq_eb_max   !< maxmimum value of the liq. water reservoir
957
958
959!
960!--    Calculate the exner function for the current time step
961       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
962
963       
964
965       DO i = nxlg, nxrg
966          DO j = nysg, nyng
967             k = nzb_s_inner(j,i)
968
969!
970!--          Set lambda_surface according to stratification
971             IF ( rif(j,i) >= 0.0_wp )  THEN
972                lambda_surface = lambda_surface_s(j,i)
973             ELSE
974                lambda_surface = lambda_surface_u(j,i)
975             ENDIF
976
977!
978!--          First step: calculate aerodyamic resistance. As pt, us, ts
979!--          are not available for the prognostic time step, data from the last
980!--          time step is used here. Note that this formulation is the
981!--          equivalent to the ECMWF formulation using drag coefficients
982             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) +       &
983                         1.0E-20_wp)
984
985!
986!--          Second step: calculate canopy resistance r_canopy
987!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
988 
989!--          f1: correction for incoming shortwave radiation (stomata close at
990!--          night)
991             IF ( radiation_scheme /= 'constant' )  THEN
992                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /&
993                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)        &
994                               + 1.0_wp)) )
995             ELSE
996                f1 = 1.0_wp
997             ENDIF
998
999!
1000!--          f2: correction for soil moisture availability to plants (the
1001!--          integrated soil moisture must thus be considered here)
1002!--          f2 = 0 for very dry soils
1003             m_total = 0.0_wp
1004             DO ks = nzb_soil, nzt_soil
1005                 m_total = m_total + root_fr(ks,j,i)                           &
1006                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
1007             ENDDO
1008
1009             IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT. m_fc(j,i) )  THEN
1010                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
1011             ELSEIF ( m_total .GE. m_fc(j,i) )  THEN
1012                f2 = 1.0_wp
1013             ELSE
1014                f2 = 1.0E-20_wp
1015             ENDIF
1016
1017!
1018!--          Calculate water vapour pressure at saturation
1019             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1020                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
1021
1022!
1023!--          f3: correction for vapour pressure deficit
1024             IF ( g_d(j,i) /= 0.0_wp )  THEN
1025!
1026!--             Calculate vapour pressure
1027                e  = q_p(k+1,j,i) * surface_pressure / 0.622_wp
1028                f3 = EXP ( -g_d(j,i) * (e_s - e) )
1029             ELSE
1030                f3 = 1.0_wp
1031             ENDIF
1032
1033!
1034!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1035!--          this calculation is obsolete, as r_canopy is not used below.
1036!--          To do: check for very dry soil -> r_canopy goes to infinity
1037             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
1038                                             + 1.0E-20_wp)
1039
1040!
1041!--          Third step: calculate bare soil resistance r_soil. The Clapp &
1042!--          Hornberger parametrization does not consider c_veg.
1043             IF ( soil_type /= 7 )  THEN
1044                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
1045                        m_res(j,i)
1046             ELSE
1047                m_min = m_wilt(j,i)
1048             ENDIF
1049
1050             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
1051             f2 = MAX(f2,1.0E-20_wp)
1052             f2 = MIN(f2,1.0_wp)
1053
1054             r_soil(j,i) = r_soil_min(j,i) / f2
1055
1056!
1057!--          Calculate fraction of liquid water reservoir
1058             m_liq_eb_max = m_max_depth * lai(j,i)
1059             c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
1060             
1061
1062!
1063!--          Calculate saturation specific humidity
1064             q_s = 0.622_wp * e_s / surface_pressure
1065
1066!
1067!--          In case of dewfall, set evapotranspiration to zero
1068!--          All super-saturated water is then removed from the air
1069             IF ( humidity .AND. dewfall .AND. q_s .LE. q_p(k+1,j,i) )  THEN
1070                r_canopy(j,i) = 0.0_wp
1071                r_soil(j,i)   = 0.0_wp
1072             ENDIF
1073
1074!
1075!--          Calculate coefficients for the total evapotranspiration
1076             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
1077                           (r_a(j,i) + r_canopy(j,i))
1078
1079             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
1080                                                          r_soil(j,i))
1081             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
1082
1083
1084!
1085!--          If soil moisture is below wilting point, plants do no longer
1086!--          transpirate.
1087!              IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
1088!                 f_qsws_veg = 0.0_wp
1089!              ENDIF
1090
1091!
1092!--          If dewfall is deactivated, vegetation, soil and liquid water
1093!--          reservoir are not allowed to take up water from the super-saturated
1094!--          air.
1095             IF ( humidity )  THEN
1096                IF ( q_s .LE. q_p(k+1,j,i) )  THEN
1097                   IF ( .NOT. dewfall )  THEN
1098                      f_qsws_veg  = 0.0_wp
1099                      f_qsws_soil = 0.0_wp
1100                      f_qsws_liq  = 0.0_wp
1101                   ENDIF
1102                ENDIF
1103             ENDIF
1104
1105             f_shf  = rho_cp / r_a(j,i)
1106             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1107
1108!
1109!--          Calculate derivative of q_s for Taylor series expansion
1110             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
1111                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
1112                              / (t_surface(j,i) - 35.86_wp)**2 )
1113
1114             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
1115
1116             T_1 = pt_p(k+1,j,i) * exn
1117
1118!
1119!--          Add LW up so that it can be removed in prognostic equation
1120             rad_net_l(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
1121
1122             IF ( humidity )  THEN
1123
1124!
1125!--             Numerator of the prognostic equation
1126                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
1127                         + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s   &
1128                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
1129                         * t_soil(nzb_soil,j,i)
1130
1131!
1132!--             Denominator of the prognostic equation
1133                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
1134                         * dq_s_dt + lambda_surface + f_shf / exn
1135
1136             ELSE
1137
1138!
1139!--             Numerator of the prognostic equation
1140                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
1141                         + f_shf / exn * T_1 + lambda_surface                  &
1142                         * t_soil(nzb_soil,j,i)
1143
1144!
1145!--             Denominator of the prognostic equation
1146                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
1147                         + lambda_surface + f_shf / exn
1148
1149             ENDIF
1150
1151             tend = 0.0_wp
1152
1153!
1154!--          Implicit solution when the surface layer has no heat capacity,
1155!--          otherwise use RK3 scheme.
1156             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface *        &
1157                                t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
1158                                * tsc(2) ) 
1159!
1160!--          Add RK3 term
1161             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
1162                                * tt_surface_m(j,i)
1163
1164!
1165!--          Calculate true tendency
1166             tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
1167                    * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
1168!
1169!--          Calculate t_surface tendencies for the next Runge-Kutta step
1170             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1171                IF ( intermediate_timestep_count == 1 )  THEN
1172                   tt_surface_m(j,i) = tend
1173                ELSEIF ( intermediate_timestep_count <                         &
1174                         intermediate_timestep_count_max )  THEN
1175                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
1176                                       * tt_surface_m(j,i)
1177                ENDIF
1178             ENDIF
1179
1180             pt_p(k,j,i) = t_surface_p(j,i) / exn
1181!
1182!--          Calculate fluxes
1183             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb                 &
1184                              * t_surface(j,i)**4 - 4.0_wp * sigma_sb          &
1185                              * t_surface(j,i)**3 * t_surface_p(j,i)
1186             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
1187                              - t_soil(nzb_soil,j,i))
1188             shf_eb(j,i)    = - f_shf  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
1189
1190             IF ( humidity )  THEN
1191                qsws_eb(j,i)  = - f_qsws    * ( q_p(k+1,j,i) - q_s + dq_s_dt   &
1192                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
1193
1194                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( q_p(k+1,j,i) - q_s       &
1195                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1196                                    * t_surface_p(j,i) )
1197
1198                qsws_soil_eb(j,i) = - f_qsws_soil * ( q_p(k+1,j,i) - q_s       &
1199                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1200                                    * t_surface_p(j,i) )
1201
1202                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( q_p(k+1,j,i) - q_s       &
1203                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1204                                    * t_surface_p(j,i) )
1205             ENDIF
1206
1207!
1208!--          Calculate the true surface resistance
1209             IF ( qsws_eb(j,i) .EQ. 0.0_wp )  THEN
1210                r_s(j,i) = 1.0E10_wp
1211             ELSE
1212                r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dt           &
1213                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
1214                           / qsws_eb(j,i) - r_a(j,i)
1215             ENDIF
1216
1217!
1218!--          Calculate change in liquid water reservoir due to dew fall or
1219!--          evaporation of liquid water
1220             IF ( humidity )  THEN
1221!
1222!--             If precipitation is activated, add rain water to qsws_liq_eb
1223!--             and qsws_soil_eb according the the vegetation coverage.
1224!--             precipitation_rate is given in mm.
1225                IF ( precipitation )  THEN
1226
1227!
1228!--                Add precipitation to liquid water reservoir, if possible.
1229!--                Otherwise, add the water to bare soil fraction.
1230                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
1231                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
1232                                       + c_veg(j,i) * precipitation_rate(j,i)  &
1233                                       * 0.001_wp * rho_l * l_v
1234                   ELSE
1235                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1236                                        + c_veg(j,i) * precipitation_rate(j,i) &
1237                                        * 0.001_wp * rho_l * l_v
1238                   ENDIF
1239
1240!--                Add precipitation to bare soil according to the bare soil
1241!--                coverage.
1242                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp             &
1243                                       - c_veg(j,i)) * precipitation_rate(j,i) &
1244                                       * 0.001_wp * rho_l * l_v
1245                ENDIF
1246!
1247!--             If the air is saturated, check the reservoir water level
1248                IF ( q_s .LE. q_p(k+1,j,i))  THEN
1249!
1250!--                Check if reservoir is full (avoid values > m_liq_eb_max)
1251!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
1252!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
1253!--                so that tend is zero and no further check is needed
1254                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
1255                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1256                                           + qsws_liq_eb(j,i)
1257                      qsws_liq_eb(j,i)  = 0.0_wp
1258                   ENDIF
1259
1260!
1261!--                In case qsws_veg_eb becomes negative (unphysical behavior),
1262!--                let the water enter the liquid water reservoir as dew on the
1263!--                plant
1264                   IF ( qsws_veg_eb(j,i) .LT. 0.0_wp )  THEN
1265                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
1266                      qsws_veg_eb(j,i) = 0.0_wp
1267                   ENDIF
1268                ENDIF                   
1269 
1270                tend = - qsws_liq_eb(j,i) * drho_l_lv
1271
1272                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
1273                                                   + tsc(3) * tm_liq_eb_m(j,i) )
1274
1275!
1276!--             Check if reservoir is overfull -> reduce to maximum
1277!--             (conservation of water is violated here)
1278                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
1279
1280!
1281!--             Check if reservoir is empty (avoid values < 0.0)
1282!--             (conservation of water is violated here)
1283                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
1284
1285
1286!
1287!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
1288                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1289                   IF ( intermediate_timestep_count == 1 )  THEN
1290                      tm_liq_eb_m(j,i) = tend
1291                   ELSEIF ( intermediate_timestep_count <                      &
1292                            intermediate_timestep_count_max )  THEN
1293                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
1294                                      * tm_liq_eb_m(j,i)
1295                   ENDIF
1296                ENDIF
1297
1298             ENDIF
1299
1300          ENDDO
1301       ENDDO
1302
1303    END SUBROUTINE lsm_energy_balance
1304
1305
1306!------------------------------------------------------------------------------!
1307! Description:
1308! ------------
1309!> Soil model as part of the land surface model. The model predicts soil
1310!> temperature and water content.
1311!------------------------------------------------------------------------------!
1312    SUBROUTINE lsm_soil_model
1313
1314
1315       IMPLICIT NONE
1316
1317       INTEGER(iwp) ::  i   !< running index
1318       INTEGER(iwp) ::  j   !< running index
1319       INTEGER(iwp) ::  k   !< running index
1320
1321       REAL(wp)     :: h_vg !< Van Genuchten coef. h
1322
1323       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
1324                                               lambda_temp, & !< temp. lambda
1325                                               tend           !< tendency
1326
1327       DO i = nxlg, nxrg   
1328          DO j = nysg, nyng
1329             DO k = nzb_soil, nzt_soil
1330!
1331!--             Calculate volumetric heat capacity of the soil, taking into
1332!--             account water content
1333                rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i))       &
1334                                     + rho_c_water * m_soil(k,j,i))
1335
1336!
1337!--             Calculate soil heat conductivity at the center of the soil
1338!--             layers
1339                ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
1340                lambda_temp(k) = ke * (lambda_h_sat(j,i) + lambda_h_dry) +     &
1341                                 lambda_h_dry
1342
1343             ENDDO
1344
1345!
1346!--          Calculate soil heat conductivity (lambda_h) at the _stag level
1347!--          using linear interpolation
1348             DO k = nzb_soil, nzt_soil-1
1349                 
1350                lambda_h(k,j,i) = lambda_temp(k) +                             &
1351                                  ( lambda_temp(k+1) - lambda_temp(k) )        &
1352                                  * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
1353
1354             ENDDO
1355             lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
1356
1357!
1358!--          Prognostic equation for soil temperature t_soil
1359             tend(:) = 0.0_wp
1360             tend(0) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *                    &
1361                       ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)     &
1362                         - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil)         &
1363                         + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
1364
1365             DO  k = 1, nzt_soil
1366                tend(k) = (1.0_wp/rho_c_total(k,j,i))                          &
1367                          * (   lambda_h(k,j,i)                                &
1368                              * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
1369                              * ddz_soil(k)                                    &
1370                              - lambda_h(k-1,j,i)                              &
1371                              * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
1372                              * ddz_soil(k-1)                                  &
1373                            ) * ddz_soil_stag(k)
1374             ENDDO
1375
1376             t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)   &
1377                                             + dt_3d * ( tsc(2)                &
1378                                             * tend(:) + tsc(3)                &
1379                                             * tt_soil_m(:,j,i) )   
1380
1381!
1382!--          Calculate t_soil tendencies for the next Runge-Kutta step
1383             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1384                IF ( intermediate_timestep_count == 1 )  THEN
1385                   DO  k = nzb_soil, nzt_soil
1386                      tt_soil_m(k,j,i) = tend(k)
1387                   ENDDO
1388                ELSEIF ( intermediate_timestep_count <                         &
1389                         intermediate_timestep_count_max )  THEN
1390                   DO  k = nzb_soil, nzt_soil
1391                      tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
1392                                         * tt_soil_m(k,j,i)
1393                   ENDDO
1394                ENDIF
1395             ENDIF
1396
1397
1398             DO k = nzb_soil, nzt_soil
1399
1400!
1401!--             Calculate soil diffusivity at the center of the soil layers
1402                lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat          &
1403                                  / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
1404                                  m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
1405
1406!
1407!--             Parametrization of Van Genuchten
1408                IF ( soil_type /= 7 )  THEN
1409!
1410!--                Calculate the hydraulic conductivity after Van Genuchten
1411!--                (1980)
1412                   h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -       &
1413                              MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i)   &
1414                              / (n_vg(j,i)-1.0_wp)) - 1.0_wp                   &
1415                          )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
1416
1417                   gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +            &
1418                                   (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp   &
1419                                   -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg    &
1420                                   )**(n_vg(j,i)-1.0_wp))**2 )                 &
1421                                   / ( (1.0_wp + (alpha_vg(j,i)*h_vg           &
1422                                   )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) &
1423                                   *(l_vg(j,i) + 2.0_wp)) )
1424
1425!
1426!--             Parametrization of Clapp & Hornberger
1427                ELSE
1428                   gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i)           &
1429                                   / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
1430                ENDIF
1431
1432             ENDDO
1433
1434
1435             IF ( humidity )  THEN
1436!
1437!--             Calculate soil diffusivity (lambda_w) at the _stag level
1438!--             using linear interpolation. To do: replace this with
1439!--             ECMWF-IFS Eq. 8.81
1440                DO k = nzb_soil, nzt_soil-1
1441                     
1442                   lambda_w(k,j,i) = lambda_temp(k) +                          &
1443                                     ( lambda_temp(k+1) - lambda_temp(k) )     &
1444                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)
1445                   gamma_w(k,j,i)  = gamma_temp(k) +                           &
1446                                     ( gamma_temp(k+1) - gamma_temp(k) )       &
1447                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)                 
1448
1449                ENDDO
1450
1451!
1452!
1453!--             In case of a closed bottom (= water content is conserved), set
1454!--             hydraulic conductivity to zero to that no water will be lost
1455!--             in the bottom layer.
1456                IF ( conserve_water_content )  THEN
1457                   gamma_w(nzt_soil,j,i) = 0.0_wp
1458                ELSE
1459                   gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil)
1460                ENDIF     
1461
1462!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
1463!--             ensures the mass conservation for water. The transpiration of
1464!--             plants equals the cumulative withdrawals by the roots in the
1465!--             soil. The scheme takes into account the availability of water
1466!--             in the soil layers as well as the root fraction in the
1467!--             respective layer. Layer with moisture below wilting point will
1468!--             not contribute, which reflects the preference of plants to
1469!--             take water from moister layers.
1470
1471!
1472!--             Calculate the root extraction (ECMWF 7.69, modified so that the
1473!--             sum of root_extr = 1). The energy balance solver guarantees a
1474!--             positive transpiration, so that there is no need for an
1475!--             additional check.
1476
1477                m_total = 0.0_wp
1478                DO k = nzb_soil, nzt_soil
1479                    IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
1480                       m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
1481                    ENDIF
1482                ENDDO 
1483
1484                DO k = nzb_soil, nzt_soil
1485                   IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
1486                      root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
1487                   ELSE
1488                      root_extr(k) = 0.0_wp
1489                   ENDIF
1490                ENDDO
1491
1492!
1493!--             Prognostic equation for soil water content m_soil
1494                tend(:) = 0.0_wp
1495                tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (                  &
1496                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
1497                            * ddz_soil(nzb_soil) - gamma_w(nzb_soil,j,i) - (   &
1498                            root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
1499                            + qsws_soil_eb(j,i) ) * drho_l_lv )                &
1500                            * ddz_soil_stag(nzb_soil)
1501
1502                DO  k = nzb_soil+1, nzt_soil-1
1503                   tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
1504                               - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&
1505                               - lambda_w(k-1,j,i) * (m_soil(k,j,i) -          &
1506                               m_soil(k-1,j,i)) * ddz_soil(k-1)                &
1507                               + gamma_w(k-1,j,i) - (root_extr(k)              &
1508                               * qsws_veg_eb(j,i) * drho_l_lv)                 &
1509                             ) * ddz_soil_stag(k)
1510
1511                ENDDO
1512                tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                     &
1513                                        - lambda_w(nzt_soil-1,j,i)             &
1514                                        * (m_soil(nzt_soil,j,i)                &
1515                                        - m_soil(nzt_soil-1,j,i))              &
1516                                        * ddz_soil(nzt_soil-1)                 &
1517                                        + gamma_w(nzt_soil-1,j,i) - (          &
1518                                          root_extr(nzt_soil)                  &
1519                                        * qsws_veg_eb(j,i) * drho_l_lv  )      &
1520                                      ) * ddz_soil_stag(nzt_soil)             
1521
1522                m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
1523                                                + dt_3d * ( tsc(2) * tend(:)   &
1524                                                + tsc(3) * tm_soil_m(:,j,i) )   
1525
1526!
1527!--             Account for dry soils (find a better solution here!)
1528                m_soil_p(:,j,i) = MAX(m_soil_p(:,j,i),0.0_wp)
1529
1530!
1531!--             Calculate m_soil tendencies for the next Runge-Kutta step
1532                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1533                   IF ( intermediate_timestep_count == 1 )  THEN
1534                      DO  k = nzb_soil, nzt_soil
1535                         tm_soil_m(k,j,i) = tend(k)
1536                      ENDDO
1537                   ELSEIF ( intermediate_timestep_count <                      &
1538                            intermediate_timestep_count_max )  THEN
1539                      DO  k = nzb_soil, nzt_soil
1540                         tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
1541                                     * tm_soil_m(k,j,i)
1542                      ENDDO
1543                   ENDIF
1544                ENDIF
1545
1546             ENDIF
1547
1548          ENDDO
1549       ENDDO
1550
1551!
1552!--    Calculate surface specific humidity
1553       IF ( humidity )  THEN
1554          CALL calc_q_surface
1555       ENDIF
1556
1557
1558    END SUBROUTINE lsm_soil_model
1559
1560
1561!------------------------------------------------------------------------------!
1562! Description:
1563! ------------
1564!> Calculation of specific humidity of the surface layer (surface)
1565!------------------------------------------------------------------------------!
1566    SUBROUTINE calc_q_surface
1567
1568       IMPLICIT NONE
1569
1570       INTEGER :: i              !< running index
1571       INTEGER :: j              !< running index
1572       INTEGER :: k              !< running index
1573       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
1574
1575       DO i = nxlg, nxrg   
1576          DO j = nysg, nyng
1577             k = nzb_s_inner(j,i)
1578
1579!
1580!--          Calculate water vapour pressure at saturation
1581             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1582                                              - 273.16_wp ) / ( t_surface(j,i) &
1583                                              - 35.86_wp ) )
1584
1585!
1586!--          Calculate specific humidity at saturation
1587             q_s = 0.622_wp * e_s / surface_pressure
1588
1589
1590             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
1591
1592!
1593!--          Calculate specific humidity at surface
1594             q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance)             &
1595                          * q_p(k+1,j,i)
1596
1597          ENDDO
1598       ENDDO
1599
1600    END SUBROUTINE calc_q_surface
1601
1602!------------------------------------------------------------------------------!
1603! Description:
1604! ------------
1605!> Swapping of timelevels
1606!------------------------------------------------------------------------------!
1607    SUBROUTINE lsm_swap_timelevel ( mod_count )
1608
1609       IMPLICIT NONE
1610
1611       INTEGER, INTENT(IN) :: mod_count
1612
1613#if defined( __nopointer )
1614
1615       t_surface    = t_surface_p
1616       t_soil       = t_soil_p
1617       IF ( humidity )  THEN
1618          m_soil    = m_soil_p
1619          m_liq_eb  = m_liq_eb_p
1620       ENDIF
1621
1622#else
1623   
1624       SELECT CASE ( mod_count )
1625
1626          CASE ( 0 )
1627
1628             t_surface  => t_surface_1; t_surface_p  => t_surface_2
1629             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
1630             IF ( humidity )  THEN
1631                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
1632                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
1633             ENDIF
1634
1635
1636          CASE ( 1 )
1637
1638             t_surface  => t_surface_2; t_surface_p  => t_surface_1
1639             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
1640             IF ( humidity )  THEN
1641                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
1642                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
1643             ENDIF
1644
1645       END SELECT
1646#endif
1647
1648    END SUBROUTINE lsm_swap_timelevel
1649
1650
1651 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.