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

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

last commit documented

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