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