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

Last change on this file since 1558 was 1556, checked in by maronga, 10 years ago

last commit documented

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