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

Last change on this file since 1533 was 1514, checked in by heinze, 10 years ago

last commmit documented

  • Property svn:keywords set to Id
File size: 53.7 KB
RevLine 
[1496]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!
[1514]23!
[1496]24! Former revisions:
25! -----------------
26! $Id: land_surface_model.f90 1514 2014-12-19 09:14:55Z hoffmann $
[1514]27! Bugfix: REAL constants provided with KIND-attribute in call of
28! intrinsic function MAX and MIN
[1496]29!
[1501]30! 1500 2014-12-03 17:42:41Z maronga
31! Corrected calculation of aerodynamic resistance (r_a).
32! Precipitation is now added to liquid water reservoir using LE_liq.
33! Added support for dry runs.
34!
[1497]35! 1496 2014-12-02 17:25:50Z maronga
36! Initial revision
37!
[1496]38!
39! Description:
40! ------------
41! Land surface model, consisting of a solver for the energy balance at the
42! surface and a four layer soil scheme. The scheme is similar to the TESSEL
43! scheme implemented in the ECMWF IFS model, with modifications according to
44! H-TESSEL. The implementation is based on the formulation implemented in the
45! DALES model.
[1500]46!
47! To do list:
48! -----------
49! - Add support for binary I/O support
50! - Add support for lsm data output
51! - Check for time step criterion
52! - Check use with RK-2 and Euler time-stepping
53! - Adaption for use with cloud physics (liquid water potential temperature)
54! - Check reaction of plants at wilting point and at atmospheric saturation
55! - Consider partial absorption of the net shortwave radiation by the skin layer
56! - Allow for water surfaces, check performance for bare soils
[1496]57!------------------------------------------------------------------------------!
58     USE arrays_3d,                                                            &
59         ONLY:  pt, pt_p, q, q_p, qsws, rif, shf, ts, us, z0, z0h
60
61     USE cloud_parameters,                                                     &
[1500]62         ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
[1496]63
64     USE control_parameters,                                                   &
65         ONLY: dt_3d, humidity, intermediate_timestep_count,                   &
[1500]66               intermediate_timestep_count_max, precipitation, pt_surface,     &
67               rho_surface, surface_pressure, timestep_scheme, tsc
[1496]68
69     USE indices,                                                              &
70         ONLY: nxlg, nxrg, nyng, nysg, nzb_s_inner 
71
72     USE kinds
73
74     USE radiation_model_mod,                                                  &
75         ONLY: Rn, SW_in, sigma_SB
76
77
78    IMPLICIT NONE
79
80!
81!-- LSM model constants
82    INTEGER(iwp), PARAMETER :: soil_layers = 4 !: number of soil layers (fixed for now)
83
84    REAL(wp), PARAMETER ::                     &
85              b_CH               = 6.04_wp,    & ! Clapp & Hornberger exponent
86              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil
87              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
88              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
89              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
90              rhoC_soil          = 2.19E6_wp,  & ! volumetric heat capacity of soil
91              rhoC_water         = 4.20E6_wp,  & ! volumetric heat capacity of water
92              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
93
94
95!
96!-- LSM variables
97    INTEGER(iwp) :: veg_type  = 2, & !: vegetation type, 0: user-defined, 1-19: generic (see list)
98                    soil_type = 3    !: soil type, 0: user-defined, 1-6: generic (see list)
99
100    LOGICAL :: conserve_water_content = .TRUE., & !: open or closed bottom surface for the soil model
101               land_surface = .FALSE.             !: flag parameter indicating wheather the lsm is used
102
103!   value 9999999.9_wp -> generic available or user-defined value must be set
104!   otherwise -> no generic variable and user setting is optional
105    REAL(wp) :: alpha_VanGenuchten = 0.0_wp,            & !: NAMELIST alpha_VG
106                canopy_resistance_coefficient = 0.0_wp, & !: NAMELIST gD
107                C_skin   = 20000.0_wp,                  & !: Skin heat capacity
108                drho_l_lv,                              & !: (rho_l * l_v)**-1
109                exn,                                    & !: value of the Exner function
110                e_s = 0.0_wp,                           & !: saturation water vapour pressure
111                field_capacity = 0.0_wp,                & !: NAMELIST m_fc
112                f_shortwave_incoming = 9999999.9_wp,    & !: NAMELIST f_SW_in
113                hydraulic_conductivity = 0.0_wp,        & !: NAMELIST gamma_w_sat
114                Ke = 0.0_wp,                            & !: Kersten number
115                lambda_skin_stable = 9999999.9_wp,      & !: NAMELIST lambda_skin_s
116                lambda_skin_unstable = 9999999.9_wp,    & !: NAMELIST lambda_skin_u
117                leaf_area_index = 9999999.9_wp,         & !: NAMELIST LAI
118                l_VanGenuchten = 0.0_WP,                & !: NAMELIST l_VG
119                min_canopy_resistance = 110.0_wp,       & !: NAMELIST r_s_min
120                m_total = 0.0_wp,                       & !: weighed total water content of the soil (m3/m3)
121                n_VanGenuchten = 0.0_WP,                & !: NAMELIST n_VG
122                q_s = 0.0_wp,                           & !: saturation specific humidity
123                residual_moisture = 0.0_wp,             & !: NAMELIST m_res
124                rho_cp,                                 & !: rho_surface * cp
125                rho_lv,                                 & !: rho * l_v
126                rd_d_rv,                                & !: r_d / r_v
127                saturation_moisture = 0.0_wp,           & !: NAMELIST m_sat
128                vegetation_coverage = 9999999.9_wp,     & !: NAMELIST c_veg
129                wilting_point = 0.0_wp                    !: NAMELIST m_wilt
130
131    REAL(wp), DIMENSION(0:soil_layers-1) :: &
132              ddz_soil,                     & !: 1/dz_soil
133              ddz_soil_stag,                & !: 1/dz_soil_stag
134              dz_soil,                      & !: soil grid spacing (center-center)
135              dz_soil_stag,                 & !: soil grid spacing (edge-edge)
136              root_extr = 0.0_wp,           & !: root extraction
137              root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers
138              soil_level = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),   & !: soil layer depths (m)
[1500]139              soil_moisture = 0.0_wp          !: soil moisture content (m3/m3)
[1496]140
141    REAL(wp), DIMENSION(0:soil_layers) ::   &
142              soil_temperature = 9999999.9_wp !: soil temperature (K)
143
144#if defined( __nopointer )
145    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0,    & !: skin temperature (K)
146                                                     T_0_p,  & !: progn. skin temperature (K)
147                                                     m_liq,  & !: liquid water reservoir (m)
148                                                     m_liq_p   !: progn. liquid water reservoir (m)
149#else
150    REAL(wp), DIMENSION(:,:), POINTER :: T_0,   &
151                                         T_0_p, & 
152                                         m_liq, & 
153                                         m_liq_p
154
155    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0_1, T_0_2,    &
156                                                     m_liq_1, m_liq_2
157#endif
158
159!
160!-- Temporal tendencies for time stepping           
161    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tT_0_m,  & !: skin temperature tendency (K)
162                                             tm_liq_m   !: liquid water reservoir tendency (m)
163
164!
165!-- Energy balance variables               
166    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
167              alpha_VG,      & !: coef. of Van Genuchten
168              c_liq,         & !: liquid water coverage (of vegetated area)
169              c_veg,         & !: vegetation coverage   
170              f_SW_in,       & !: ?
171              G,             & !: surface soil heat flux
172              H,             & !: surface flux of sensible heat
173              gamma_w_sat,   & !: hydraulic conductivity at saturation
174              gD,            & !: coefficient for dependence of r_canopy on water vapour pressure deficit
175              LAI,           & !: leaf area index
176              LE,            & !: surface flux of latent heat (total)
177              LE_veg,        & !: surface flux of latent heat (vegetation portion)
178              LE_soil,       & !: surface flux of latent heat (soil portion)
179              LE_liq,        & !: surface flux of latent heat (liquid water portion)
180              lambda_h_sat,  & !: heat conductivity for dry soil
181              lambda_skin_s, & !: coupling between skin and soil (depends on vegetation type)
182              lambda_skin_u, & !: coupling between skin and soil (depends on vegetation type)
183              l_VG,          & !: coef. of Van Genuchten
184              m_fc,          & !: soil moisture at field capacity (m3/m3)
185              m_res,         & !: residual soil moisture
186              m_sat,         & !: saturation soil moisture (m3/m3)
187              m_wilt,        & !: soil moisture at permanent wilting point (m3/m3)
188              n_VG,          & !: coef. Van Genuchten 
189              r_a,           & !: aerodynamic resistance
190              r_canopy,      & !: canopy resistance
191              r_soil,        & !: soil resitance
192              r_soil_min,    & !: minimum soil resistance
193              r_s,           & !: total surface resistance (combination of r_soil and r_canopy)         
194              r_s_min          !: minimum canopy resistance
195
196    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
197              lambda_h, &   !: heat conductivity of soil (?)                           
198              lambda_w, &   !: hydraulic diffusivity of soil (?)
199              gamma_w,  &   !: hydraulic conductivity of soil (?)
200              rhoC_total    !: volumetric heat capacity of the actual soil matrix (?)
201
202#if defined( __nopointer )
203    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
204              T_soil,    & !: Soil temperature (K)
205              T_soil_p,  & !: Prog. soil temperature (K)
206              m_soil,    & !: Soil moisture (m3/m3)
207              m_soil_p     !: Prog. soil moisture (m3/m3)
208#else
209    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
210              T_soil, T_soil_p, &
211              m_soil, m_soil_p   
212
213    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
214              T_soil_1, T_soil_2, &
215              m_soil_1, m_soil_2
216
217
218#endif
219
220
221    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
222              tT_soil_m, & !: T_soil storage array
223              tm_soil_m, & !: m_soil storage array
224              root_fr      !: root fraction (sum=1)
225
226!
227!--  Land surface parameters according to the following classes (veg_type)
228!--  (0 user defined)
229!--  1 crops, mixed farming
230!--  2 short grass
231!--  3 evergreen needleleaf trees
232!--  4 deciduous needleleaf trees
233!--  5 evergreen broadleaf trees
234!--  6 deciduous broadleaf trees
235!--  7 tall grass
236!--  8 desert
237!--  9 tundra
238!-- 10 irrigated crops
239!-- 11 semidesert
240!-- 12 ice caps and glaciers
241!-- 13 bogs and marshes
242!-- 14 inland water
243!-- 15 ocean
244!-- 16 evergreen shrubs
245!-- 17 deciduous shrubs
246!-- 18 mixed forest/woodland
247!-- 19 interrupted forest
248
249!
250!-- Land surface parameters I     r_s_min,     LAI,   c_veg,      gD
251    REAL(wp), DIMENSION(0:3,1:19) :: veg_pars = RESHAPE( (/           &
252                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & !  1
253                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & !  2
254                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  3
255                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  4
256                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  5
257                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & !  6
258                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & !  7
259                                 250.0_wp, 0.50_wp, 0.00_wp, 0.00_wp, & !  8
260                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & !  9
261                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10
262                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, & ! 11
263                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
264                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, & ! 13
265                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
266                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
267                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, & ! 16
268                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17
269                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18
270                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp  & ! 19
271                                 /), (/ 4, 19 /) )
272
273!
274!-- Land surface parameters II          z0,         z0h
275    REAL(wp), DIMENSION(0:1,1:19) :: roughness_par = RESHAPE( (/ & 
276                                   0.25_wp,  0.25E-2_wp,         & !  1
277                                   0.20_wp,  0.20E-2_wp,         & !  2
278                                   2.00_wp,     2.00_wp,         & !  3
279                                   2.00_wp,     2.00_wp,         & !  4
280                                   2.00_wp,     2.00_wp,         & !  5
281                                   2.00_wp,     2.00_wp,         & !  6
282                                   0.47_wp,  0.47E-2_wp,         & !  7
283                                  0.013_wp, 0.013E-2_wp,         & !  8
284                                  0.034_wp, 0.034E-2_wp,         & !  9
285                                    0.5_wp,  0.50E-2_wp,         & ! 10
286                                   0.17_wp,  0.17E-2_wp,         & ! 11
287                                 1.3E-3_wp,   1.3E-4_wp,         & ! 12
288                                   0.83_wp,  0.83E-2_wp,         & ! 13
289                                   0.00_wp,  0.00E-2_wp,         & ! 14
290                                   0.00_wp,  0.00E-2_wp,         & ! 15
291                                   0.10_wp,  0.10E-2_wp,         & ! 16
292                                   0.25_wp,  0.25E-2_wp,         & ! 17
293                                   2.00_wp,  2.00E-2_wp,         & ! 18
294                                   1.10_wp,  1.10E-2_wp          & ! 19
295                                 /), (/ 2, 19 /) )
296
297!
298!-- Land surface parameters III lambda_skin_s, lambda_skin_u, f_SW_in
299    REAL(wp), DIMENSION(0:2,1:19) :: skin_pars = RESHAPE( (/           &
300                                      10.0_wp,       10.0_wp, 0.05_wp, & !  1
301                                      10.0_wp,       10.0_wp, 0.05_wp, & !  2
302                                      20.0_wp,       15.0_wp, 0.03_wp, & !  3
303                                      20.0_wp,       15.0_wp, 0.03_wp, & !  4
304                                      20.0_wp,       15.0_wp, 0.03_wp, & !  5
305                                      20.0_wp,       15.0_wp, 0.03_wp, & !  6
306                                      10.0_wp,       10.0_wp, 0.05_wp, & !  7
307                                      15.0_wp,       15.0_wp, 0.00_wp, & !  8
308                                      10.0_wp,       10.0_wp, 0.05_wp, & !  9
309                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 10
310                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 11
311                                      58.0_wp,       58.0_wp, 0.00_wp, & ! 12
312                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 13
313                                    1.0E20_wp,     1.0E20_wp, 0.00_wp, & ! 14
314                                    1.0E20_wp,     1.0E20_wp, 0.00_wp, & ! 15
315                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 16
316                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 17
317                                      20.0_wp,       15.0_wp, 0.03_wp, & ! 18
318                                      20.0_wp,       15.0_wp, 0.03_wp  & ! 19
319                                      /), (/ 3, 19 /) )
320
321!
322!-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
323    REAL(wp), DIMENSION(0:3,1:19) :: root_distribution = RESHAPE( (/ &
324                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & !  1
325                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & !  2
326                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & !  3
327                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & !  4
328                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & !  5
329                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & !  6
330                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & !  7
331                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & !  8
332                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & !  9
333                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 10
334                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 11
335                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
336                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 13
337                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
338                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
339                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16
340                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17
341                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18
342                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp  & ! 19
343                                 /), (/ 4, 19 /) )
344
345!
346!-- Soil parameters according to the following porosity classes (soil_type)
347!-- (0 user defined)
348!-- 1 coarse
349!-- 2 medium
350!-- 3 medium-fine
351!-- 4 fine
352!-- 5 very fine
353!-- 6 organic
354!
355!-- Soil parameters I           alpha_VG,      l_VG,    n_VG, gamma_w_sat
356    REAL(wp), DIMENSION(0:3,1:6) :: soil_pars = RESHAPE( (/                &
357                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
358                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
359                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
360                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
361                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
362                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp  & ! 6
363                                 /), (/ 4, 6 /) )
364
365!
366!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
367    REAL(wp), DIMENSION(0:3,1:6) :: m_soil_pars = RESHAPE( (/            &
368                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
369                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
370                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
371                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
372                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
373                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp  & ! 6
374                                 /), (/ 4, 6 /) )
375
376
377    SAVE
378
379
380    PRIVATE
381
382
383    PUBLIC alpha_VanGenuchten, C_skin, canopy_resistance_coefficient,          &
384           conserve_water_content,      field_capacity, f_shortwave_incoming,  &
385           hydraulic_conductivity, init_lsm, lambda_skin_stable,               &
386           lambda_skin_unstable, land_surface, leaf_area_index,                &
387           lsm_energy_balance, lsm_soil_model, l_VanGenuchten,                 &
388           min_canopy_resistance, n_VanGenuchten, residual_moisture,           &
389           root_fraction, saturation_moisture, soil_level, soil_moisture,      &
390           soil_temperature, soil_type, vegetation_coverage, veg_type,         &
391           wilting_point
392
393#if defined( __nopointer )
394    PUBLIC m_liq, m_liq_p, m_soil, m_soil_p, T_0, T_0_p, T_soil, T_soil_p
395#else
396    PUBLIC m_liq, m_liq_1, m_liq_2, m_liq_p, m_soil, m_soil_1, m_soil_2,       &
397           m_soil_p, T_0, T_0_1, T_0_2, T_0_p, T_soil, T_soil_1, T_soil_2,     &
398           T_soil_p
399#endif
400
401
402    INTERFACE init_lsm
403       MODULE PROCEDURE init_lsm
404    END INTERFACE init_lsm
405
406    INTERFACE lsm_energy_balance
407       MODULE PROCEDURE lsm_energy_balance
408    END INTERFACE lsm_energy_balance
409
410    INTERFACE lsm_soil_model
411       MODULE PROCEDURE lsm_soil_model
412    END INTERFACE lsm_soil_model
413
414
415 CONTAINS
416
417
418!------------------------------------------------------------------------------!
419! Description:
420! ------------
421!-- Initialization of the land surface model
422!------------------------------------------------------------------------------!
423    SUBROUTINE init_lsm
424   
425
426       IMPLICIT NONE
427
428       INTEGER(iwp) ::  i !: running index
429       INTEGER(iwp) ::  j !: running index
430       INTEGER(iwp) ::  k !: running index
431
432
433!
434!--    Calculate frequently used parameters
435       rho_cp    = cp * rho_surface
436       rd_d_rv   = r_d / r_v
437       rho_lv    = rho_surface * l_v
438       drho_l_lv = 1.0 / (rho_l * l_v)
439
440!
441!--    Allocate skin and soil temperature / humidity
442#if defined( __nopointer )
443       ALLOCATE ( T_0(nysg:nyng,nxlg:nxrg) )
444       ALLOCATE ( T_0_p(nysg:nyng,nxlg:nxrg) )
445#else
446       ALLOCATE ( T_0_1(nysg:nyng,nxlg:nxrg) )
447       ALLOCATE ( T_0_2(nysg:nyng,nxlg:nxrg) )
448#endif
449
450       ALLOCATE ( tT_0_m(nysg:nyng,nxlg:nxrg) )
451
452#if defined( __nopointer )
453       ALLOCATE ( T_soil(0:soil_layers,nysg:nyng,nxlg:nxrg) )
454       ALLOCATE ( T_soil_p(0:soil_layers,nysg:nyng,nxlg:nxrg) )
455#else
456       ALLOCATE ( T_soil_1(0:soil_layers,nysg:nyng,nxlg:nxrg) )
457       ALLOCATE ( T_soil_2(0:soil_layers,nysg:nyng,nxlg:nxrg) )
458#endif
459
460       ALLOCATE ( tT_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
461
462#if defined( __nopointer )
463       ALLOCATE ( m_liq(nysg:nyng,nxlg:nxrg) )
464       ALLOCATE ( m_liq_p(nysg:nyng,nxlg:nxrg) )
465#else
466       ALLOCATE ( m_liq_1(nysg:nyng,nxlg:nxrg) )
467       ALLOCATE ( m_liq_2(nysg:nyng,nxlg:nxrg) )
468#endif
469
470       ALLOCATE ( tm_liq_m(nysg:nyng,nxlg:nxrg) )
471
472#if defined( __nopointer )
473       ALLOCATE ( m_soil(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
474       ALLOCATE ( m_soil_p(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
475#else
476       ALLOCATE ( m_soil_1(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
477       ALLOCATE ( m_soil_2(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
478#endif
479
480       ALLOCATE ( tm_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
481
482
483#if ! defined( __nopointer )
484!
485!--    Initial assignment of the pointers
486       T_soil => T_soil_1; T_soil_p => T_soil_2
487       T_0 => T_0_1; T_0_p => T_0_2
488       m_soil => m_soil_1; m_soil_p => m_soil_2
489       m_liq => m_liq_1; m_liq_p => m_liq_2
490#endif
491
492       T_0    = 0.0_wp
493       T_0_p  = 0.0_wp
494       tT_0_m = 0.0_wp
495
496       T_soil    = 0.0_wp
497       T_soil_p  = 0.0_wp
498       tT_soil_m = 0.0_wp
499
500       m_liq    = 0.0_wp
501       m_liq_p  = 0.0_wp
502       tm_liq_m = 0.0_wp
503
504       m_soil    = 0.0_wp
505       m_soil_p  = 0.0_wp
506       tm_soil_m = 0.0_wp
507
508!
509!--    Allocate 2D vegetation model arrays
510       ALLOCATE ( alpha_VG(nysg:nyng,nxlg:nxrg) )
511       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
512       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
513       ALLOCATE ( f_SW_in(nysg:nyng,nxlg:nxrg) )
514       ALLOCATE ( G(nysg:nyng,nxlg:nxrg) )
515       ALLOCATE ( H(nysg:nyng,nxlg:nxrg) )
516       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
517       ALLOCATE ( gD(nysg:nyng,nxlg:nxrg) )
518       ALLOCATE ( LAI(nysg:nyng,nxlg:nxrg) )
519       ALLOCATE ( LE(nysg:nyng,nxlg:nxrg) )
520       ALLOCATE ( LE_veg(nysg:nyng,nxlg:nxrg) )
521       ALLOCATE ( LE_soil(nysg:nyng,nxlg:nxrg) )
522       ALLOCATE ( LE_liq(nysg:nyng,nxlg:nxrg) )
523       ALLOCATE ( l_VG(nysg:nyng,nxlg:nxrg) )
524       ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
525       ALLOCATE ( lambda_skin_u(nysg:nyng,nxlg:nxrg) )
526       ALLOCATE ( lambda_skin_s(nysg:nyng,nxlg:nxrg) )
527       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
528       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
529       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
530       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
531       ALLOCATE ( n_VG(nysg:nyng,nxlg:nxrg) )
532       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
533       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
534       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
535       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
536       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
537       ALLOCATE ( r_s_min(nysg:nyng,nxlg:nxrg) )
538
539!
540!--    Set initial and default values
541       c_liq   = 0.0_wp
542       c_veg   = 0.0_wp
543       f_SW_in = 0.05_wp
544       gD      = 0.0_wp
545       LAI     = 0.0_wp
546       lambda_skin_u = 10.0_wp
547       lambda_skin_s = 10.0_wp
548
549
550       G       = 0.0_wp
551       H       = rho_cp * shf
[1500]552
553       IF ( humidity )  THEN
554          LE = rho_l * l_v * qsws
555       ELSE
556          LE = 0.0_wp
557       ENDIF
558
[1496]559       LE_veg  = 0.0_wp
560       LE_soil = LE
561       LE_liq  = 0.0_wp
562
563       r_a        = 50.0_wp
564       r_canopy   = 0.0_wp
565       r_soil     = 0.0_wp
566       r_soil_min = 50.0_wp
567       r_s        = 110.0_wp
568       r_s_min    = min_canopy_resistance
569
570!
571!--    Allocate 3D soil model arrays
572       ALLOCATE ( root_fr(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
573       ALLOCATE ( lambda_h(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
574       ALLOCATE ( rhoC_total(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
575
576       lambda_h = 0.0_wp
577!
578!--    If required, allocate humidity-related variables for the soil model
579       IF ( humidity )  THEN
580          ALLOCATE ( lambda_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
581          ALLOCATE ( gamma_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )   
582
583          lambda_w = 0.0_wp 
584       ENDIF
585
586!
587!--    Calculate grid spacings. Temperature and moisture are defined at
588!--    the center of the soil layers, whereas gradients/fluxes are defined
589!--    at the edges (_stag)
590       dz_soil_stag(0) = soil_level(0)
591
592       DO k = 1, soil_layers-1
593          dz_soil_stag(k) = soil_level(k) - soil_level(k-1)
594       ENDDO
595
596       DO k = 0, soil_layers-2
597          dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k))
598       ENDDO
599       dz_soil(soil_layers-1) = dz_soil_stag(soil_layers-1)
600
601       ddz_soil      = 1.0 / dz_soil
602       ddz_soil_stag = 1.0 / dz_soil_stag
603!
604!--    Initialize soil
605       IF ( soil_type .NE. 0 )  THEN   
606          alpha_VG    = soil_pars(0,soil_type)
607          l_VG        = soil_pars(1,soil_type)
608          n_VG        = soil_pars(2,soil_type)   
609          gamma_w_sat = soil_pars(3,soil_type) 
610          m_sat       = m_soil_pars(0,soil_type)
611          m_fc        = m_soil_pars(1,soil_type)   
612          m_wilt      = m_soil_pars(2,soil_type) 
613          m_res       = m_soil_pars(3,soil_type)
614       ELSE
615          alpha_VG    = alpha_VanGenuchten
616          l_VG        = l_VanGenuchten
617          n_VG        = n_VanGenuchten 
618          gamma_w_sat = hydraulic_conductivity
619          m_sat       = saturation_moisture
620          m_fc        = field_capacity
621          m_wilt      = wilting_point
622          m_res       = residual_moisture
623       ENDIF   
624
625!
626!--    Map user settings of T and q for each soil layer
627!--    (make sure that the soil moisture does not drop below the permanent
628!--    wilting point) -> problems with devision by zero)
629       DO k = 0, soil_layers-1
630          T_soil(k,:,:)  = soil_temperature(k)
631          m_soil(k,:,:)  = MAX(soil_moisture(k),m_wilt(:,:))
632       ENDDO
633       T_soil(soil_layers,:,:) = soil_temperature(soil_layers)
634
635
636       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
637       T_0  = pt_surface * exn
638
639       T_soil_p = T_soil
640       m_soil_p = m_soil
641
642!
643!--    Calculate saturation soil heat conductivity
644       lambda_h_sat(:,:) = lambda_h_sm ** (1.0_wp - m_sat(:,:)) *              &
645                           lambda_h_water ** m_sat(:,:)
646
647!
648!--    Initialize vegetation
649       IF ( veg_type .NE. 0 )  THEN
650
651          r_s_min              = veg_pars(0,veg_type)
652          LAI                  = veg_pars(1,veg_type)
653          c_veg                = veg_pars(2,veg_type)
654          gD                   = veg_pars(3,veg_type)
655          lambda_skin_s        = skin_pars(0,veg_type)
656          lambda_skin_u        = skin_pars(1,veg_type)
657          f_SW_in              = skin_pars(2,veg_type)
658          z0                   = roughness_par(0,veg_type)
659          z0h                  = roughness_par(1,veg_type)
660
661
662          DO k = 0, soil_layers-1
663             root_fr(k,:,:) = root_distribution(k,veg_type)
664          ENDDO
665
666       ELSE
667
668          DO k = 0, soil_layers-1
669             root_fr(k,:,:) = root_fraction(k)
670          ENDDO
671
672       ENDIF
673
674!
675!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
676       CALL user_init_land_surface
677
678!
679!--    Set artifical values for ts and us so that r_a has its initial value for
680!--    the first time step
681       DO  i = nxlg, nxrg
682          DO  j = nysg, nyng
683             k = nzb_s_inner(j,i)
[1500]684!
685!--          Assure that r_a cannot be zero at model start
686             IF ( pt(k+1,j,i) == pt(k,j,i) )  pt(k+1,j,i) = pt(k+1,j,i) + 1.0E-10_wp
687
[1496]688             us(j,i) = 0.1_wp
689             ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
690             shf(j,i) = - us(j,i) * ts(j,i)
691          ENDDO
692       ENDDO
693
694!
695!--    Calculate humidity at the surface
696       IF ( humidity )  THEN
697          CALL calc_q0
698       ENDIF
699
700       RETURN
701
702    END SUBROUTINE init_lsm
703
704
705
706!------------------------------------------------------------------------------!
707! Description:
708! ------------
709!
710!------------------------------------------------------------------------------!
711    SUBROUTINE lsm_energy_balance
712
713
714       IMPLICIT NONE
715
716       INTEGER(iwp) ::  i         !: running index
717       INTEGER(iwp) ::  j         !: running index
718       INTEGER(iwp) ::  k, ks     !: running index
719
720       REAL(wp) :: f1,          & !: resistance correction term 1
721                   f2,          & !: resistance correction term 2
722                   f3,          & !: resistance correction term 3
723                   m_min,       & !: minimum soil moisture
724                   T_1,         & !: actual temperature at first grid point
725                   e,           & !: water vapour pressure
726                   e_s,         & !: water vapour saturation pressure
727                   e_s_dT,      & !: derivate of e_s with respect to T
728                   tend,        & !: tendency
729                   dq_s_dT,     & !: derivate of q_s with respect to T
730                   coef_1,      & !: coef. for prognostic equation
731                   coef_2,      & !: coef. for prognostic equation
732                   f_LE,        & !: factor for LE
733                   f_LE_veg,    & !: factor for LE_veg
734                   f_LE_soil,   & !: factor for LE_soil
735                   f_LE_liq,    & !: factor for LE_liq
736                   f_H,         & !: factor for H
737                   lambda_skin, & !: Current value of lambda_skin
738                   m_liq_max      !: maxmimum value of the liquid water reservoir
739
740!
741!--    Calculate the exner function for the current time step
742       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
743
744
745       DO i = nxlg, nxrg
746          DO j = nysg, nyng
[1500]747             k = nzb_s_inner(j,i)
[1496]748
749!
750!--          Set lambda_skin according to stratification
751             IF ( rif(j,i) >= 0.0_wp )  THEN
752                lambda_skin = lambda_skin_s(j,i)
753             ELSE
754                lambda_skin = lambda_skin_u(j,i)
755             ENDIF
[1500]756
[1496]757!
[1500]758!--          First step: calculate aerodyamic resistance. As pt, us, ts
759!--          are not available for the prognostic time step, data from the last
760!--          time step is used here. Note that this formulation is the
761!--          equivalent to the ECMWF formulation using drag coefficients
762             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20)
[1496]763
764!
765!--          Second step: calculate canopy resistance r_canopy
766!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
767 
768!--          f1: correction for incoming shortwave radiation
769             f1 = MIN(1.0_wp, ( 0.004_wp * SW_in(j,i) + 0.05_wp ) /     &
770                              (0.81_wp * (0.004_wp * SW_in(j,i) + 1.0_wp) ) )
771
772!
773!--          f2: correction for soil moisture f2=0 for very dry soil
774             m_total = 0.0_wp
775             DO ks = 0, soil_layers-1
776                 m_total = m_total + root_fr(ks,j,i) * m_soil(ks,j,i)
777             ENDDO 
778
779             IF (  m_total .GT. m_wilt(j,i) .AND. m_total .LE. m_fc(j,i) )  THEN
780                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
781             ELSE
782                f2 = 1.0E-20_wp
783             ENDIF
784
785!
786!--          Calculate water vapour pressure at saturation
787!--          (T_0 should be replaced by liquid water temp?!)
788             e_s = 0.01 * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) - 273.16_wp )&
789                                           / ( T_0(j,i) - 35.86_wp ) )
790
791!
792!--          f3: correction for vapour pressure deficit
793             IF ( gD(j,i) .NE. 0.0_wp )  THEN
794!
795!--             Calculate vapour pressure
796                e  = q_p(k+1,j,i) * surface_pressure / 0.622
797                f3 = EXP ( -gD(j,i) * (e_s - e) )
798             ELSE
799                f3 = 1.0_wp
800             ENDIF
801
802!
803!--          To do: check for very dry soil -> r_canopy goes to infinity
804             r_canopy(j,i)  = r_s_min(j,i) / (LAI(j,i) * f1 * f2 * f3 + 1.0E-20)
805
806!
807!--          Third step: calculate bare soil resistance r_soil
808             m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *        &
809                     m_res(j,i)
810
811             f2 = ( m_soil(0,j,i) - m_min ) / ( m_fc(j,i) - m_min )
[1513]812             f2 = MAX(f2,1.0E-20_wp)
[1496]813
814             r_soil(j,i) = r_soil_min(j,i) / f2
815
816!
817!--          Calculate fraction of liquid water reservoir
818             m_liq_max = m_max_depth * LAI(j,i)
[1513]819             c_liq(j,i) = MIN(1.0_wp, m_liq(j,i)/m_liq_max)
[1496]820
821             q_s = 0.622_wp * e_s / surface_pressure
[1500]822
823!
824!--          In case of dew fall, set resistances to zero.
825!--          To do: what does that physically reasoning behind this?
826             IF ( humidity )  THEN
827                IF ( q_s .LE. q_p(k+1,j,i) )  THEN
828                   r_canopy(j,i) = 0.0_wp
829                   r_soil(j,i) = 0.0_wp
830                ENDIF
[1496]831             ENDIF 
832
833
834!
835!--          Calculate coefficients for the total evapotranspiration
836             f_LE_veg  = rho_lv * c_veg(j,i) * (1.0 - c_liq(j,i)) / (r_a(j,i)  &
837                                                + r_canopy(j,i))
838             f_LE_soil = rho_lv * (1.0 - c_veg(j,i)) / (r_a(j,i) + r_soil(j,i))
839             f_LE_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
840
841
[1500]842!
843!--          If soil moisture is below wilting point, plants do no longer
844!--          transpirate.
845             IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
846                f_LE_veg = 0.0
847             ENDIF
[1496]848
849             f_H  = rho_cp / r_a(j,i)
850             f_LE = f_LE_veg + f_LE_soil + f_LE_liq
851       
852!
853!--          Calculate derivative of q_s for Taylor series expansion
854             e_s_dT = e_s * ( 17.269_wp / (T_0(j,i) - 35.86_wp) -              &
855                              17.269_wp*(T_0(j,i) - 273.16_wp) / (T_0(j,i)     &
856                              - 35.86_wp)**2 )
857
858             dq_s_dT = 0.622_wp * e_s_dT / surface_pressure
859
860             T_1 = pt_p(k+1,j,i) * exn
861
862!
863!--          Add LW up so that it can be removed in prognostic equation
864             Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4
865
[1500]866             IF ( humidity )  THEN
867
[1496]868!
[1500]869!--             Numerator of the prognostic equation
870                coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H     &
871                         / exn * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT   &
872                         * T_0(j,i) ) + lambda_skin * T_soil(0,j,i)
[1496]873
874!
[1500]875!--             Denominator of the prognostic equation
876                coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT    &
877                         + lambda_skin + f_H / exn
[1496]878
[1500]879             ELSE
880
881!
882!--             Numerator of the prognostic equation
883                coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H /   &
884                         exn * T_1 + lambda_skin * T_soil(0,j,i)
885
886!
887!--             Denominator of the prognostic equation
888                coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3                     &
889                         + lambda_skin + f_H / exn
890
891             ENDIF
892
[1496]893             tend = 0.0_wp
894
895!
896!--          Implicit solution when the skin layer has no heat capacity,
897!--          otherwise use RK3 scheme.
898             T_0_p(j,i) = ( coef_1 * dt_3d * tsc(2) + C_skin * T_0(j,i) ) /    &
899                          ( C_skin + coef_2 * dt_3d * tsc(2) ) 
900
901!
902!--          Add RK3 term
903             T_0_p(j,i) = T_0_p(j,i) + dt_3d * tsc(3) * tT_soil_m(0,j,i)
904
905!
906!--          Calculate true tendency
907             tend = (T_0_p(j,i) - T_0(j,i) - tsc(3) * tT_0_m(j,i)) / (dt_3d    &
908                      * tsc(2))
909
910!
911!--          Calculate T_0 tendencies for the next Runge-Kutta step
912             IF ( timestep_scheme(1:5) == 'runge' )  THEN
913                IF ( intermediate_timestep_count == 1 )  THEN
914                   tT_0_m(j,i) = tend
915                ELSEIF ( intermediate_timestep_count <                         &
916                         intermediate_timestep_count_max )  THEN
917                   tT_0_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tT_0_m(j,i)
918                ENDIF
919             ENDIF
920
921             pt_p(k,j,i) = T_0_p(j,i) / exn
922!
923!--          Calculate fluxes
924             Rn(j,i)        = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i)**4        &
925                              - 4.0_wp * sigma_SB * T_0(j,i)**3 * T_0_p(j,i)
926             G(j,i)         = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i))
927             H(j,i)         = - f_H  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
928
[1500]929             IF ( humidity )  THEN
930                LE(j,i)        = - f_LE      * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
931                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
[1496]932
[1500]933                LE_veg(j,i)    = - f_LE_veg  * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
934                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
935                LE_soil(j,i)   = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
936                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
937                LE_liq(j,i)    = - f_LE_liq  * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
938                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
939             ENDIF
[1496]940
941!              IF ( i == 1 .AND. j == 1 )  THEN
942!                 PRINT*, "Rn", Rn(j,i)
943!                  PRINT*, "H", H(j,i)
944!                 PRINT*, "LE", LE(j,i)
945!                 PRINT*, "LE_liq", LE_liq(j,i)
946!                 PRINT*, "LE_veg", LE_veg(j,i)
947!                 PRINT*, "LE_soil", LE_soil(j,i)
948!                 PRINT*, "G", G(j,i)
949!              ENDIF
950
[1500]951!
952!--          Calculate the true surface resistance
[1496]953             IF ( LE(j,i) .EQ. 0.0 )  THEN
954                r_s(j,i) = 1.0E10
955             ELSE
956                r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dT * T_0(j,i)&
957                           - dq_s_dT * T_0_p(j,i) ) / LE(j,i) - r_a(j,i)
958             ENDIF
959
960!
[1500]961!--          Calculate fluxes in the atmosphere
962             shf(j,i) = H(j,i) / rho_cp
963
964!
[1496]965!--          Calculate change in liquid water reservoir due to dew fall or
[1500]966!--          evaporation of liquid water
967             IF ( humidity )  THEN
[1496]968!
[1500]969!--             If precipitation is activated, add rain water to LE_liq.
970!--             precipitation_rate is given in mm.
971                IF ( precipitation )  THEN
972                   LE_liq(j,i) = LE_liq(j,i) + precipitation_rate(j,i)         &
973                                               * 0.001_wp * rho_l * l_v
[1496]974                ENDIF
[1500]975!
976!--             If the air is saturated, check the reservoir water level
977                IF ( q_s .LE. q_p(k+1,j,i))  THEN
978!
979!--                Check if reservoir is full (avoid values > m_liq_max)
980!--                In that case, LE_liq goes to LE_soil. In this case
981!--                LE_veg is zero anyway (because c_liq = 1), so that tend is
982!--                zero and no further check is needed
983                   IF ( m_liq(j,i) .EQ. m_liq_max )  THEN
984                      LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i)
985                      LE_liq(j,i) = 0.0_wp
986                   ENDIF
[1496]987
988!
[1500]989!--                In case LE_veg becomes negative (unphysical behavior), let
990!--                the water enter the liquid water reservoir as dew on the
991!--                plant
992                   IF ( LE_veg(j,i) .LT. 0.0_wp )  THEN
993                      LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)
994                      LE_veg(j,i) = 0.0_wp
995                   ENDIF
996                ENDIF                   
[1496]997 
[1500]998                tend = - LE_liq(j,i) * drho_l_lv
[1496]999
[1500]1000                m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend            &
[1496]1001                                                   + tsc(3) * tm_liq_m(j,i) )
1002
1003!
[1500]1004!--             Check if reservoir is overfull -> reduce to maximum
1005!--             (conservation of water is violated here)
1006                m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max)
[1496]1007
1008!
[1500]1009!--             Check if reservoir is empty (avoid values < 0.0)
1010!--             (conservation of water is violated here)
1011                m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp)
[1496]1012
1013
1014!
[1500]1015!--             Calculate m_liq tendencies for the next Runge-Kutta step
1016                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1017                   IF ( intermediate_timestep_count == 1 )  THEN
1018                      tm_liq_m(j,i) = tend
1019                   ELSEIF ( intermediate_timestep_count <                      &
1020                            intermediate_timestep_count_max )  THEN
1021                      tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp            &
1022                                      * tm_liq_m(j,i)
1023                   ENDIF
[1496]1024                ENDIF
1025
1026!
[1500]1027!--             Calculate moisture flux in the atmosphere
1028                qsws(j,i) = LE(j,i) / rho_lv
[1496]1029
[1500]1030             ENDIF
1031
[1496]1032          ENDDO
[1500]1033       ENDDO
[1496]1034
1035
1036
1037    END SUBROUTINE lsm_energy_balance
1038
1039
1040!------------------------------------------------------------------------------!
1041! Description:
1042! ------------
1043!
1044!------------------------------------------------------------------------------!
1045    SUBROUTINE lsm_soil_model
1046
1047
1048       IMPLICIT NONE
1049
1050       INTEGER(iwp) ::  i   !: running index
1051       INTEGER(iwp) ::  j   !: running index
1052       INTEGER(iwp) ::  k   !: running index
1053
1054       REAL(wp)     :: h_VG !: Van Genuchten coef. h
1055
1056       REAL(wp), DIMENSION(0:soil_layers-1) :: gamma_temp,  & !: temp. gamma
1057                                               lambda_temp, & !: temp. lambda
1058                                               tend           !: tendency
1059
1060       DO i = nxlg, nxrg   
1061          DO j = nysg, nyng
1062             DO k = 0, soil_layers-1
1063!
1064!--             Calculate volumetric heat capacity of the soil, taking into
1065!--             account water content
1066                rhoC_total(k,j,i) = (rhoC_soil * (1.0 - m_sat(j,i))            &
1067                                     + rhoC_water * m_soil(k,j,i))
1068
1069!
1070!--             Calculate soil heat conductivity at the center of the soil
1071!--             layers
[1513]1072                Ke = 1.0 + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
[1496]1073                lambda_temp(k) = Ke * (lambda_h_sat(j,i) + lambda_h_dry) +     &
1074                                 lambda_h_dry
1075
1076             ENDDO
1077
1078!
1079!--          Calculate soil heat conductivity (lambda_h) at the _stag level
1080!--          using linear interpolation
1081             DO k = 0, soil_layers-2
1082                 
1083                lambda_h(k,j,i) = lambda_temp(k) +                             &
1084                                  ( lambda_temp(k+1) - lambda_temp(k) )        &
1085                                  * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
1086
1087             ENDDO
1088             lambda_h(soil_layers-1,j,i) = lambda_temp(soil_layers-1)
1089
1090!
1091!--          Prognostic equation for soil temperature T_soil
1092             tend(:) = 0.0_wp
1093             tend(0) = (1.0/rhoC_total(0,j,i)) *                               &
1094                       ( lambda_h(0,j,i) * ( T_soil(1,j,i) - T_soil(0,j,i) )   &
1095                         * ddz_soil(0) + G(j,i) ) * ddz_soil_stag(0)
1096
1097             DO  k = 1, soil_layers-1
1098                tend(k) = (1.0/rhoC_total(k,j,i))                              &
1099                          * (   lambda_h(k,j,i)                                &
1100                              * ( T_soil(k+1,j,i) - T_soil(k,j,i) )            &
1101                              * ddz_soil(k)                                    &
1102                              - lambda_h(k-1,j,i)                              &
1103                              * ( T_soil(k,j,i) - T_soil(k-1,j,i) )            &
1104                              * ddz_soil(k-1)                                  &
1105                            ) * ddz_soil_stag(k)
1106             ENDDO
1107
1108             T_soil_p(0:soil_layers-1,j,i) = T_soil(0:soil_layers-1,j,i)       &
1109                                             + dt_3d * ( tsc(2)                &
1110                                             * tend(:) + tsc(3)                &
1111                                             * tT_soil_m(:,j,i) )   
1112
1113!
1114!--          Calculate T_soil tendencies for the next Runge-Kutta step
1115             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1116                IF ( intermediate_timestep_count == 1 )  THEN
1117                   DO  k = 0, soil_layers-1
1118                      tT_soil_m(k,j,i) = tend(k)
1119                   ENDDO
1120                ELSEIF ( intermediate_timestep_count <                         &
1121                         intermediate_timestep_count_max )  THEN
1122                   DO  k = 0, soil_layers-1
1123                      tT_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
1124                                         * tT_soil_m(k,j,i)
1125                   ENDDO
1126                ENDIF
1127             ENDIF
1128
1129
1130             DO k = 0, soil_layers-1
1131!
1132!--             Calculate soil diffusivity at the center of the soil layers
1133                lambda_temp(k) = (- b_CH * gamma_w_sat(j,i) * psi_sat          &
1134                                  / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
1135                                  m_wilt(j,i)) / m_sat(j,i) )**(b_CH + 2.0_wp)
1136
1137!
1138!--             Calculate the hydraulic conductivity after Van Genuchten (1980)
1139                h_VG = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -          &
1140                           MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_VG(j,i)      &
1141                           / (n_VG(j,i)-1.0_wp)) - 1.0_wp                      &
1142                       )**(1.0_wp/n_VG(j,i)) / alpha_VG(j,i)
1143
1144                gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +               &
1145                                (alpha_VG(j,i)*h_VG)**n_VG(j,i))**(1.0_wp      &
1146                                -1.0_wp/n_VG(j,i)) - (alpha_VG(j,i)*h_VG       &
1147                                )**(n_VG(j,i)-1.0_wp))**2 )                    &
1148                                / ( (1.0_wp + (alpha_VG(j,i)*h_VG)**n_VG(j,i)  &
1149                                )**((1.0_wp - 1.0_wp/n_VG(j,i))*(l_VG(j,i)     &
1150                                + 2.0)) )
1151
1152             ENDDO
1153
1154
1155             IF ( humidity )  THEN
1156!
1157!--             Calculate soil diffusivity (lambda_w) at the _stag level
1158!--             using linear interpolation
1159                DO k = 0, soil_layers-2
1160                     
1161                   lambda_w(k,j,i) = lambda_temp(k) +                          &
1162                                     ( lambda_temp(k+1) - lambda_temp(k) )     &
1163                                     * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
1164                   gamma_w(k,j,i)  = gamma_temp(k) +                           &
1165                                     ( gamma_temp(k+1) - gamma_temp(k) )       &
1166                                     * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)                 
1167
1168                ENDDO
1169
1170!
1171!
1172!--             In case of a closed bottom (= water content is conserved), set
1173!--             hydraulic conductivity to zero to that no water will be lost
1174!--             in the bottom layer.
1175                IF ( conserve_water_content )  THEN
1176                   gamma_w(soil_layers-1,j,i) = 0.0_wp
1177                ELSE
1178                   gamma_w(soil_layers-1,j,i) = lambda_temp(soil_layers-1)
1179                ENDIF     
1180
1181!--             The root extraction (= root_extr * LE_veg / (rho_l * l_v))
1182!--             ensures the mass conservation for water. The transpiration of
1183!--             plants equals the cumulative withdrawals by the roots in the
1184!--             soil. The scheme takes into account the availability of water
1185!--             in the soil layers as well as the root fraction in the
1186!--             respective layer
1187
1188!
1189!--             Calculate the root extraction (ECMWF 7.69, with some
1190!--             modifications)
1191                m_total = 0.0_wp
1192                DO k = 0, soil_layers-1
1193                    m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) *       &
1194                              dz_soil_stag(k) 
1195
1196                ENDDO 
1197
1198!
1199!--             For conservation of mass, the sum of root_extr must be 1
1200                DO k = 0, soil_layers-1 
1201                   root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)               &
1202                                  * dz_soil_stag(k) / m_total
1203                ENDDO
1204
1205
1206!
1207!--             Prognostic equation for soil water content m_soil
1208                tend(:) = 0.0_wp
1209                tend(0) = ( lambda_w(0,j,i) * ( m_soil(1,j,i) - m_soil(0,j,i) )&
1210                            * ddz_soil(0) - gamma_w(0,j,i) - ( root_extr(0)    &
1211                            * LE_veg(j,i) + LE_soil(j,i) ) * drho_l_lv         &
1212                          ) * ddz_soil_stag(0)
1213
1214                DO  k = 1, soil_layers-2
1215                   tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
1216                               - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&
1217                               - lambda_w(k-1,j,i) * (m_soil(k,j,i) -          &
1218                               m_soil(k-1,j,i)) * ddz_soil(k-1)                &
1219                               + gamma_w(k-1,j,i) - (root_extr(k) * LE_veg(j,i)&
1220                               * drho_l_lv)                                    &
1221                             ) * ddz_soil_stag(k)
1222
1223                ENDDO
1224                tend(soil_layers-1) = ( - gamma_w(soil_layers-1,j,i)           &
1225                                        - lambda_w(soil_layers-2,j,i)          &
1226                                        * (m_soil(soil_layers-1,j,i)           &
1227                                        - m_soil(soil_layers-2,j,i))           &
1228                                        * ddz_soil(soil_layers-2)              &
1229                                        + gamma_w(soil_layers-2,j,i) - (       &
1230                                          root_extr(soil_layers-1)             &
1231                                        * LE_veg(j,i) * drho_l_lv      )       &
1232                                      ) * ddz_soil_stag(soil_layers-1)             
1233
1234                m_soil_p(0:soil_layers-1,j,i) = m_soil(0:soil_layers-1,j,i)    &
1235                                                + dt_3d * ( tsc(2) * tend(:)   &
1236                                                + tsc(3) * tm_soil_m(:,j,i) )   
1237
1238!
1239!--             Account for dry soils (find a better solution here!)
1240                m_soil_p(:,j,i) = MAX(m_soil_p(:,j,i),0.0_wp)
1241
1242!
1243!--             Calculate m_soil tendencies for the next Runge-Kutta step
1244                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1245                   IF ( intermediate_timestep_count == 1 )  THEN
1246                      DO  k = 0, soil_layers-1
1247                         tm_soil_m(k,j,i) = tend(k)
1248                      ENDDO
1249                   ELSEIF ( intermediate_timestep_count <                      &
1250                            intermediate_timestep_count_max )  THEN
1251                      DO  k = 0, soil_layers-1
1252                         tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
1253                                     * tm_soil_m(k,j,i)
1254                      ENDDO
1255                   ENDIF
1256                ENDIF
1257
1258             ENDIF
1259
1260          ENDDO
1261       ENDDO
1262
1263!
1264!--    Calculate surface specific humidity
1265       IF ( humidity )  THEN
1266          CALL calc_q0
1267       ENDIF
1268
1269
1270    END SUBROUTINE lsm_soil_model
1271
1272
1273!------------------------------------------------------------------------------!
1274! Description:
1275! ------------
1276!
1277!------------------------------------------------------------------------------!
1278    SUBROUTINE calc_q0
1279
1280       IMPLICIT NONE
1281
1282       INTEGER :: i              !: running index
1283       INTEGER :: j              !: running index
1284       INTEGER :: k              !: running index
1285       REAL(wp) :: resistance    !: aerodynamic and soil resistance term
1286
1287       DO i = nxlg, nxrg   
1288          DO j = nysg, nyng
1289             k = nzb_s_inner(j,i)
1290!
1291!--          Temporary solution as long as T_0 is prescribed
1292
1293             pt_p(k,j,i) = T_0(j,i) / exn
1294!
1295!--          Calculate water vapour pressure at saturation
1296             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) -         &
1297                                              273.16_wp ) /  ( T_0(j,i) -      &
1298                                              35.86_wp ) )
1299
1300!
1301!--          Calculate specific humidity at saturation
1302             q_s = 0.622_wp * e_s / surface_pressure
1303
1304
1305             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
1306
1307!
1308!--          Calculate specific humidity at surface
1309             q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance)             &
1310                          * q_p(k+1,j,i)
1311
1312          ENDDO
1313       ENDDO
1314
1315    END SUBROUTINE calc_q0
1316
1317
1318 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.