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

Last change on this file since 1501 was 1501, checked in by maronga, 9 years ago

last commit documented

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