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

Last change on this file since 1572 was 1572, checked in by maronga, 7 years ago

last commit documented

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