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

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

land surface model released

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