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

Last change on this file since 1585 was 1585, checked in by maronga, 6 years ago

Added support for RRTMG radiation code

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