[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 |
---|