[1682] | 1 | !> @file land_surface_model.f90 |
---|
[1496] | 2 | !--------------------------------------------------------------------------------! |
---|
| 3 | ! This file is part of PALM. |
---|
| 4 | ! |
---|
| 5 | ! PALM is free software: you can redistribute it and/or modify it under the terms |
---|
| 6 | ! of the GNU General Public License as published by the Free Software Foundation, |
---|
| 7 | ! either version 3 of the License, or (at your option) any later version. |
---|
| 8 | ! |
---|
| 9 | ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
| 10 | ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
---|
| 11 | ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
---|
| 12 | ! |
---|
| 13 | ! You should have received a copy of the GNU General Public License along with |
---|
| 14 | ! PALM. If not, see <http://www.gnu.org/licenses/>. |
---|
| 15 | ! |
---|
[1691] | 16 | ! Copyright 1997-2015 Leibniz Universitaet Hannover |
---|
[1496] | 17 | !--------------------------------------------------------------------------------! |
---|
| 18 | ! |
---|
| 19 | ! Current revisions: |
---|
| 20 | ! ----------------- |
---|
[1788] | 21 | ! Bugfix: calculate lambda_surface based on temperature gradient between skin |
---|
| 22 | ! layer and soil layer instead of Obukhov length |
---|
| 23 | ! Changed: moved calculation of surface specific humidity to energy balance solver |
---|
| 24 | ! New: water surfaces are available by using a fixed sea surface temperature. |
---|
| 25 | ! The roughness lengths are calculated dynamically using the Charnock |
---|
| 26 | ! parameterization. This involves the new roughness length for moisture z0q. |
---|
| 27 | ! New: modified solution of the energy balance solver and soil model for |
---|
| 28 | ! paved surfaces (i.e. asphalt concrete). |
---|
| 29 | ! Syntax layout improved. |
---|
| 30 | ! Changed: parameter dewfall removed. |
---|
[1758] | 31 | ! |
---|
| 32 | ! Former revisions: |
---|
| 33 | ! ----------------- |
---|
| 34 | ! $Id: land_surface_model.f90 1788 2016-03-10 11:01:04Z maronga $ |
---|
| 35 | ! |
---|
[1784] | 36 | ! 1783 2016-03-06 18:36:17Z raasch |
---|
| 37 | ! netcdf variables moved to netcdf module |
---|
| 38 | ! |
---|
[1758] | 39 | ! 1757 2016-02-22 15:49:32Z maronga |
---|
[1757] | 40 | ! Bugfix: set tm_soil_m to zero after allocation. Added parameter |
---|
| 41 | ! unscheduled_radiation_calls to control calls of the radiation model based on |
---|
| 42 | ! the skin temperature change during one time step (preliminary version). Set |
---|
| 43 | ! qsws_soil_eb to zero at model start (previously set to qsws_eb). Removed MAX |
---|
| 44 | ! function as it cannot be vectorized. |
---|
[1710] | 45 | ! |
---|
| 46 | ! 1709 2015-11-04 14:47:01Z maronga |
---|
[1709] | 47 | ! Renamed pt_1 and qv_1 to pt1 and qv1. |
---|
| 48 | ! Bugfix: set initial values for t_surface_p in case of restart runs |
---|
| 49 | ! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky' |
---|
| 50 | ! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky' |
---|
| 51 | ! Added todo action |
---|
[1698] | 52 | ! |
---|
| 53 | ! 1697 2015-10-28 17:14:10Z raasch |
---|
| 54 | ! bugfix: misplaced cpp-directive |
---|
| 55 | ! |
---|
| 56 | ! 1695 2015-10-27 10:03:11Z maronga |
---|
[1692] | 57 | ! Bugfix: REAL constants provided with KIND-attribute in call of |
---|
[1696] | 58 | ! Replaced rif with ol |
---|
| 59 | ! |
---|
[1692] | 60 | ! 1691 2015-10-26 16:17:44Z maronga |
---|
[1691] | 61 | ! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes: |
---|
| 62 | ! Soil temperatures are now defined at the edges of the layers, calculation of |
---|
| 63 | ! shb_eb corrected, prognostic equation for skin temperature corrected. Surface |
---|
| 64 | ! fluxes are now directly transfered to atmosphere |
---|
[1552] | 65 | ! |
---|
[1683] | 66 | ! 1682 2015-10-07 23:56:08Z knoop |
---|
| 67 | ! Code annotations made doxygen readable |
---|
| 68 | ! |
---|
[1591] | 69 | ! 1590 2015-05-08 13:56:27Z maronga |
---|
| 70 | ! Bugfix: definition of character strings requires same length for all elements |
---|
| 71 | ! |
---|
[1586] | 72 | ! 1585 2015-04-30 07:05:52Z maronga |
---|
| 73 | ! Modifications for RRTMG. Changed tables to PARAMETER type. |
---|
| 74 | ! |
---|
[1572] | 75 | ! 1571 2015-03-12 16:12:49Z maronga |
---|
| 76 | ! Removed upper-case variable names. Corrected distribution of precipitation to |
---|
| 77 | ! the liquid water reservoir and the bare soil fractions. |
---|
| 78 | ! |
---|
[1556] | 79 | ! 1555 2015-03-04 17:44:27Z maronga |
---|
| 80 | ! Added output of r_a and r_s |
---|
| 81 | ! |
---|
[1554] | 82 | ! 1553 2015-03-03 17:33:54Z maronga |
---|
| 83 | ! Improved better treatment of roughness lengths. Added default soil temperature |
---|
| 84 | ! profile |
---|
| 85 | ! |
---|
[1552] | 86 | ! 1551 2015-03-03 14:18:16Z maronga |
---|
[1551] | 87 | ! Flux calculation is now done in prandtl_fluxes. Added support for data output. |
---|
| 88 | ! Vertical indices have been replaced. Restart runs are now possible. Some |
---|
| 89 | ! variables have beem renamed. Bugfix in the prognostic equation for the surface |
---|
| 90 | ! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of |
---|
| 91 | ! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for |
---|
| 92 | ! the hydraulic conductivity. Bugfix for root fraction and extraction |
---|
| 93 | ! calculation |
---|
[1496] | 94 | ! |
---|
[1514] | 95 | ! intrinsic function MAX and MIN |
---|
[1496] | 96 | ! |
---|
[1501] | 97 | ! 1500 2014-12-03 17:42:41Z maronga |
---|
| 98 | ! Corrected calculation of aerodynamic resistance (r_a). |
---|
| 99 | ! Precipitation is now added to liquid water reservoir using LE_liq. |
---|
| 100 | ! Added support for dry runs. |
---|
| 101 | ! |
---|
[1497] | 102 | ! 1496 2014-12-02 17:25:50Z maronga |
---|
| 103 | ! Initial revision |
---|
| 104 | ! |
---|
[1496] | 105 | ! |
---|
| 106 | ! Description: |
---|
| 107 | ! ------------ |
---|
[1682] | 108 | !> Land surface model, consisting of a solver for the energy balance at the |
---|
| 109 | !> surface and a four layer soil scheme. The scheme is similar to the TESSEL |
---|
| 110 | !> scheme implemented in the ECMWF IFS model, with modifications according to |
---|
| 111 | !> H-TESSEL. The implementation is based on the formulation implemented in the |
---|
| 112 | !> DALES and UCLA-LES models. |
---|
| 113 | !> |
---|
[1691] | 114 | !> @todo Consider partial absorption of the net shortwave radiation by the |
---|
[1709] | 115 | !> skin layer. |
---|
[1788] | 116 | !> @todo Improve surface water parameterization |
---|
[1691] | 117 | !> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, |
---|
| 118 | !> nzt_soil=3)). |
---|
| 119 | !> @todo Implement surface runoff model (required when performing long-term LES |
---|
| 120 | !> with considerable precipitation. |
---|
[1709] | 121 | !> @todo Fix crashes with radiation_scheme == 'constant' |
---|
[1682] | 122 | !> |
---|
[1709] | 123 | !> @note No time step criterion is required as long as the soil layers do not |
---|
| 124 | !> become too thin. |
---|
[1496] | 125 | !------------------------------------------------------------------------------! |
---|
[1682] | 126 | MODULE land_surface_model_mod |
---|
| 127 | |
---|
[1691] | 128 | USE arrays_3d, & |
---|
[1788] | 129 | ONLY: hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h, & |
---|
| 130 | z0q |
---|
[1496] | 131 | |
---|
[1691] | 132 | USE cloud_parameters, & |
---|
| 133 | ONLY: cp, hyrho, l_d_cp, l_d_r, l_v, prr, pt_d_t, rho_l, r_d, r_v |
---|
[1496] | 134 | |
---|
[1691] | 135 | USE control_parameters, & |
---|
| 136 | ONLY: cloud_physics, dt_3d, humidity, intermediate_timestep_count, & |
---|
| 137 | initializing_actions, intermediate_timestep_count_max, & |
---|
| 138 | max_masks, precipitation, pt_surface, rho_surface, & |
---|
| 139 | roughness_length, surface_pressure, timestep_scheme, tsc, & |
---|
| 140 | z0h_factor, time_since_reference_point |
---|
[1496] | 141 | |
---|
[1691] | 142 | USE indices, & |
---|
| 143 | ONLY: nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner |
---|
[1496] | 144 | |
---|
[1691] | 145 | USE kinds |
---|
[1496] | 146 | |
---|
[1691] | 147 | USE pegrid |
---|
[1496] | 148 | |
---|
[1691] | 149 | USE radiation_model_mod, & |
---|
| 150 | ONLY: force_radiation_call, radiation_scheme, rad_net, rad_sw_in, & |
---|
[1757] | 151 | rad_lw_out, rad_lw_out_change_0, sigma_sb, & |
---|
| 152 | unscheduled_radiation_calls |
---|
[1691] | 153 | |
---|
| 154 | #if defined ( __rrtmg ) |
---|
| 155 | USE radiation_model_mod, & |
---|
[1709] | 156 | ONLY: rrtm_idrv |
---|
[1691] | 157 | #endif |
---|
[1496] | 158 | |
---|
[1691] | 159 | USE statistics, & |
---|
| 160 | ONLY: hom, statistic_regions |
---|
| 161 | |
---|
[1496] | 162 | IMPLICIT NONE |
---|
| 163 | |
---|
| 164 | ! |
---|
| 165 | !-- LSM model constants |
---|
[1682] | 166 | INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !< bottom of the soil model (to be switched) |
---|
| 167 | nzt_soil = 3, & !< top of the soil model (to be switched) |
---|
| 168 | nzs = 4 !< number of soil layers (fixed for now) |
---|
[1496] | 169 | |
---|
| 170 | REAL(wp), PARAMETER :: & |
---|
[1551] | 171 | b_ch = 6.04_wp, & ! Clapp & Hornberger exponent |
---|
| 172 | lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil |
---|
[1496] | 173 | lambda_h_sm = 3.44_wp, & ! heat conductivity of the soil matrix |
---|
| 174 | lambda_h_water = 0.57_wp, & ! heat conductivity of water |
---|
| 175 | psi_sat = -0.388_wp, & ! soil matrix potential at saturation |
---|
[1551] | 176 | rho_c_soil = 2.19E6_wp, & ! volumetric heat capacity of soil |
---|
| 177 | rho_c_water = 4.20E6_wp, & ! volumetric heat capacity of water |
---|
[1496] | 178 | m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m) |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | ! |
---|
| 182 | !-- LSM variables |
---|
[1788] | 183 | INTEGER(iwp) :: veg_type = 2, & !< NAMELIST veg_type_2d |
---|
| 184 | soil_type = 3 !< NAMELIST soil_type_2d |
---|
[1496] | 185 | |
---|
[1788] | 186 | INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: soil_type_2d, & !< soil type, 0: user-defined, 1-7: generic (see list) |
---|
| 187 | veg_type_2d !< vegetation type, 0: user-defined, 1-19: generic (see list) |
---|
| 188 | |
---|
| 189 | LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface, & !< flag parameter for water surfaces (classes 14+15) |
---|
| 190 | pave_surface, & !< flag parameter for pavements (asphalt etc.) (class 20) |
---|
| 191 | building_surface !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet) |
---|
| 192 | |
---|
[1691] | 193 | LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model |
---|
| 194 | force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls |
---|
| 195 | land_surface = .FALSE. !< flag parameter indicating wheather the lsm is used |
---|
[1496] | 196 | |
---|
| 197 | ! value 9999999.9_wp -> generic available or user-defined value must be set |
---|
| 198 | ! otherwise -> no generic variable and user setting is optional |
---|
[1682] | 199 | REAL(wp) :: alpha_vangenuchten = 9999999.9_wp, & !< NAMELIST alpha_vg |
---|
| 200 | canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d |
---|
| 201 | c_surface = 20000.0_wp, & !< Surface (skin) heat capacity |
---|
| 202 | drho_l_lv, & !< (rho_l * l_v)**-1 |
---|
| 203 | exn, & !< value of the Exner function |
---|
| 204 | e_s = 0.0_wp, & !< saturation water vapour pressure |
---|
| 205 | field_capacity = 9999999.9_wp, & !< NAMELIST m_fc |
---|
| 206 | f_shortwave_incoming = 9999999.9_wp, & !< NAMELIST f_sw_in |
---|
| 207 | hydraulic_conductivity = 9999999.9_wp, & !< NAMELIST gamma_w_sat |
---|
| 208 | ke = 0.0_wp, & !< Kersten number |
---|
[1691] | 209 | lambda_h_sat = 0.0_wp, & !< heat conductivity for saturated soil |
---|
[1682] | 210 | lambda_surface_stable = 9999999.9_wp, & !< NAMELIST lambda_surface_s |
---|
| 211 | lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u |
---|
| 212 | leaf_area_index = 9999999.9_wp, & !< NAMELIST lai |
---|
| 213 | l_vangenuchten = 9999999.9_wp, & !< NAMELIST l_vg |
---|
| 214 | min_canopy_resistance = 9999999.9_wp, & !< NAMELIST r_canopy_min |
---|
| 215 | min_soil_resistance = 50.0_wp, & !< NAMELIST r_soil_min |
---|
| 216 | m_total = 0.0_wp, & !< weighted total water content of the soil (m3/m3) |
---|
| 217 | n_vangenuchten = 9999999.9_wp, & !< NAMELIST n_vg |
---|
[1788] | 218 | pave_depth = 9999999.9_wp, & !< depth of the pavement |
---|
| 219 | pave_heat_capacity = 1.94E6_wp, & !< volumetric heat capacity of pavement (e.g. roads) |
---|
| 220 | pave_heat_conductivity = 1.00_wp, & !< heat conductivity for pavements (e.g. roads) |
---|
[1682] | 221 | q_s = 0.0_wp, & !< saturation specific humidity |
---|
| 222 | residual_moisture = 9999999.9_wp, & !< NAMELIST m_res |
---|
| 223 | rho_cp, & !< rho_surface * cp |
---|
| 224 | rho_lv, & !< rho * l_v |
---|
| 225 | rd_d_rv, & !< r_d / r_v |
---|
| 226 | saturation_moisture = 9999999.9_wp, & !< NAMELIST m_sat |
---|
[1691] | 227 | skip_time_do_lsm = 0.0_wp, & !< LSM is not called before this time |
---|
[1682] | 228 | vegetation_coverage = 9999999.9_wp, & !< NAMELIST c_veg |
---|
| 229 | wilting_point = 9999999.9_wp, & !< NAMELIST m_wilt |
---|
| 230 | z0_eb = 9999999.9_wp, & !< NAMELIST z0 (lsm_par) |
---|
[1788] | 231 | z0h_eb = 9999999.9_wp, & !< NAMELIST z0h (lsm_par) |
---|
| 232 | z0q_eb = 9999999.9_wp !< NAMELIST z0q (lsm_par) |
---|
[1496] | 233 | |
---|
[1551] | 234 | REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: & |
---|
[1691] | 235 | ddz_soil_stag, & !< 1/dz_soil_stag |
---|
| 236 | dz_soil_stag, & !< soil grid spacing (center-center) |
---|
| 237 | root_extr = 0.0_wp, & !< root extraction |
---|
[1551] | 238 | root_fraction = (/9999999.9_wp, 9999999.9_wp, & |
---|
[1682] | 239 | 9999999.9_wp, 9999999.9_wp /), & !< distribution of root surface area to the individual soil layers |
---|
| 240 | zs = (/0.07_wp, 0.28_wp, 1.00_wp, 2.89_wp/), & !< soil layer depths (m) |
---|
| 241 | soil_moisture = 0.0_wp !< soil moisture content (m3/m3) |
---|
[1496] | 242 | |
---|
[1551] | 243 | REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) :: & |
---|
[1691] | 244 | soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp, 283.0_wp, & !< soil temperature (K) |
---|
| 245 | 283.0_wp /), & |
---|
| 246 | ddz_soil, & !< 1/dz_soil |
---|
| 247 | dz_soil !< soil grid spacing (edge-edge) |
---|
[1496] | 248 | |
---|
| 249 | #if defined( __nopointer ) |
---|
[1682] | 250 | REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface, & !< surface temperature (K) |
---|
| 251 | t_surface_p, & !< progn. surface temperature (K) |
---|
| 252 | m_liq_eb, & !< liquid water reservoir (m) |
---|
| 253 | m_liq_eb_av, & !< liquid water reservoir (m) |
---|
| 254 | m_liq_eb_p !< progn. liquid water reservoir (m) |
---|
[1496] | 255 | #else |
---|
[1551] | 256 | REAL(wp), DIMENSION(:,:), POINTER :: t_surface, & |
---|
| 257 | t_surface_p, & |
---|
| 258 | m_liq_eb, & |
---|
| 259 | m_liq_eb_p |
---|
[1496] | 260 | |
---|
[1551] | 261 | REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, & |
---|
| 262 | m_liq_eb_av, & |
---|
| 263 | m_liq_eb_1, m_liq_eb_2 |
---|
[1496] | 264 | #endif |
---|
| 265 | |
---|
| 266 | ! |
---|
| 267 | !-- Temporal tendencies for time stepping |
---|
[1682] | 268 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m, & !< surface temperature tendency (K) |
---|
| 269 | tm_liq_eb_m !< liquid water reservoir tendency (m) |
---|
[1496] | 270 | |
---|
| 271 | ! |
---|
| 272 | !-- Energy balance variables |
---|
[1691] | 273 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & |
---|
[1682] | 274 | alpha_vg, & !< coef. of Van Genuchten |
---|
| 275 | c_liq, & !< liquid water coverage (of vegetated area) |
---|
| 276 | c_liq_av, & !< average of c_liq |
---|
| 277 | c_soil_av, & !< average of c_soil |
---|
| 278 | c_veg, & !< vegetation coverage |
---|
| 279 | c_veg_av, & !< average of c_veg |
---|
| 280 | f_sw_in, & !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet) |
---|
| 281 | ghf_eb, & !< ground heat flux |
---|
| 282 | ghf_eb_av, & !< average of ghf_eb |
---|
| 283 | gamma_w_sat, & !< hydraulic conductivity at saturation |
---|
| 284 | g_d, & !< coefficient for dependence of r_canopy on water vapour pressure deficit |
---|
| 285 | lai, & !< leaf area index |
---|
| 286 | lai_av, & !< average of lai |
---|
[1691] | 287 | lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type) |
---|
| 288 | lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type) |
---|
[1682] | 289 | l_vg, & !< coef. of Van Genuchten |
---|
| 290 | m_fc, & !< soil moisture at field capacity (m3/m3) |
---|
| 291 | m_res, & !< residual soil moisture |
---|
| 292 | m_sat, & !< saturation soil moisture (m3/m3) |
---|
| 293 | m_wilt, & !< soil moisture at permanent wilting point (m3/m3) |
---|
| 294 | n_vg, & !< coef. Van Genuchten |
---|
| 295 | qsws_eb, & !< surface flux of latent heat (total) |
---|
| 296 | qsws_eb_av, & !< average of qsws_eb |
---|
| 297 | qsws_liq_eb, & !< surface flux of latent heat (liquid water portion) |
---|
| 298 | qsws_liq_eb_av, & !< average of qsws_liq_eb |
---|
| 299 | qsws_soil_eb, & !< surface flux of latent heat (soil portion) |
---|
| 300 | qsws_soil_eb_av, & !< average of qsws_soil_eb |
---|
| 301 | qsws_veg_eb, & !< surface flux of latent heat (vegetation portion) |
---|
| 302 | qsws_veg_eb_av, & !< average of qsws_veg_eb |
---|
| 303 | rad_net_l, & !< local copy of rad_net (net radiation at surface) |
---|
| 304 | r_a, & !< aerodynamic resistance |
---|
[1709] | 305 | r_a_av, & !< average of r_a |
---|
[1682] | 306 | r_canopy, & !< canopy resistance |
---|
[1691] | 307 | r_soil, & !< soil resistance |
---|
[1682] | 308 | r_soil_min, & !< minimum soil resistance |
---|
| 309 | r_s, & !< total surface resistance (combination of r_soil and r_canopy) |
---|
[1691] | 310 | r_s_av, & !< average of r_s |
---|
[1682] | 311 | r_canopy_min, & !< minimum canopy (stomatal) resistance |
---|
| 312 | shf_eb, & !< surface flux of sensible heat |
---|
| 313 | shf_eb_av !< average of shf_eb |
---|
[1496] | 314 | |
---|
[1691] | 315 | |
---|
[1496] | 316 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & |
---|
[1788] | 317 | lambda_h, & !< heat conductivity of soil (W/m/K) |
---|
[1682] | 318 | lambda_w, & !< hydraulic diffusivity of soil (?) |
---|
[1788] | 319 | gamma_w, & !< hydraulic conductivity of soil (W/m/K) |
---|
[1682] | 320 | rho_c_total !< volumetric heat capacity of the actual soil matrix (?) |
---|
[1496] | 321 | |
---|
| 322 | #if defined( __nopointer ) |
---|
| 323 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & |
---|
[1682] | 324 | t_soil, & !< Soil temperature (K) |
---|
| 325 | t_soil_av, & !< Average of t_soil |
---|
| 326 | t_soil_p, & !< Prog. soil temperature (K) |
---|
| 327 | m_soil, & !< Soil moisture (m3/m3) |
---|
| 328 | m_soil_av, & !< Average of m_soil |
---|
| 329 | m_soil_p !< Prog. soil moisture (m3/m3) |
---|
[1496] | 330 | #else |
---|
| 331 | REAL(wp), DIMENSION(:,:,:), POINTER :: & |
---|
[1551] | 332 | t_soil, t_soil_p, & |
---|
[1496] | 333 | m_soil, m_soil_p |
---|
| 334 | |
---|
| 335 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & |
---|
[1551] | 336 | t_soil_av, t_soil_1, t_soil_2, & |
---|
| 337 | m_soil_av, m_soil_1, m_soil_2 |
---|
[1496] | 338 | #endif |
---|
| 339 | |
---|
| 340 | |
---|
| 341 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & |
---|
[1757] | 342 | tt_soil_m, & !< t_soil storage array |
---|
| 343 | tm_soil_m, & !< m_soil storage array |
---|
[1682] | 344 | root_fr !< root fraction (sum=1) |
---|
[1496] | 345 | |
---|
[1551] | 346 | |
---|
[1496] | 347 | ! |
---|
[1551] | 348 | !-- Predefined Land surface classes (veg_type) |
---|
[1788] | 349 | CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ & |
---|
[1590] | 350 | 'user defined ', & ! 0 |
---|
| 351 | 'crops, mixed farming ', & ! 1 |
---|
| 352 | 'short grass ', & ! 2 |
---|
[1585] | 353 | 'evergreen needleleaf trees', & ! 3 |
---|
| 354 | 'deciduous needleleaf trees', & ! 4 |
---|
[1590] | 355 | 'evergreen broadleaf trees ', & ! 5 |
---|
| 356 | 'deciduous broadleaf trees ', & ! 6 |
---|
| 357 | 'tall grass ', & ! 7 |
---|
| 358 | 'desert ', & ! 8 |
---|
| 359 | 'tundra ', & ! 9 |
---|
| 360 | 'irrigated crops ', & ! 10 |
---|
| 361 | 'semidesert ', & ! 11 |
---|
| 362 | 'ice caps and glaciers ', & ! 12 |
---|
| 363 | 'bogs and marshes ', & ! 13 |
---|
| 364 | 'inland water ', & ! 14 |
---|
| 365 | 'ocean ', & ! 15 |
---|
| 366 | 'evergreen shrubs ', & ! 16 |
---|
| 367 | 'deciduous shrubs ', & ! 17 |
---|
| 368 | 'mixed forest/woodland ', & ! 18 |
---|
[1788] | 369 | 'interrupted forest ', & ! 19 |
---|
| 370 | 'pavements/roads ' & ! 20 |
---|
[1585] | 371 | /) |
---|
[1496] | 372 | |
---|
| 373 | ! |
---|
[1551] | 374 | !-- Soil model classes (soil_type) |
---|
[1585] | 375 | CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ & |
---|
| 376 | 'user defined', & ! 0 |
---|
[1590] | 377 | 'coarse ', & ! 1 |
---|
| 378 | 'medium ', & ! 2 |
---|
| 379 | 'medium-fine ', & ! 3 |
---|
| 380 | 'fine ', & ! 4 |
---|
| 381 | 'very fine ', & ! 5 |
---|
| 382 | 'organic ', & ! 6 |
---|
| 383 | 'loamy (CH) ' & ! 7 |
---|
[1585] | 384 | /) |
---|
[1551] | 385 | ! |
---|
| 386 | !-- Land surface parameters according to the respective classes (veg_type) |
---|
| 387 | |
---|
| 388 | ! |
---|
| 389 | !-- Land surface parameters I |
---|
| 390 | !-- r_canopy_min, lai, c_veg, g_d |
---|
[1788] | 391 | REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ & |
---|
[1585] | 392 | 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 1 |
---|
| 393 | 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & ! 2 |
---|
| 394 | 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 3 |
---|
| 395 | 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 4 |
---|
| 396 | 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 5 |
---|
| 397 | 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & ! 6 |
---|
| 398 | 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & ! 7 |
---|
| 399 | 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, & ! 8 |
---|
| 400 | 80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & ! 9 |
---|
| 401 | 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10 |
---|
| 402 | 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, & ! 11 |
---|
| 403 | 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12 |
---|
| 404 | 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, & ! 13 |
---|
| 405 | 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14 |
---|
| 406 | 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15 |
---|
| 407 | 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, & ! 16 |
---|
| 408 | 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17 |
---|
| 409 | 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18 |
---|
[1788] | 410 | 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp, & ! 19 |
---|
| 411 | 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp & ! 20 |
---|
| 412 | /), (/ 4, 20 /) ) |
---|
[1496] | 413 | |
---|
| 414 | ! |
---|
[1788] | 415 | !-- Land surface parameters II z0, z0h, z0q |
---|
| 416 | REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ & |
---|
| 417 | 0.25_wp, 0.25E-2_wp, 0.25E-2_wp, & ! 1 |
---|
| 418 | 0.20_wp, 0.20E-2_wp, 0.20E-2_wp, & ! 2 |
---|
| 419 | 2.00_wp, 2.00_wp, 2.00_wp, & ! 3 |
---|
| 420 | 2.00_wp, 2.00_wp, 2.00_wp, & ! 4 |
---|
| 421 | 2.00_wp, 2.00_wp, 2.00_wp, & ! 5 |
---|
| 422 | 2.00_wp, 2.00_wp, 2.00_wp, & ! 6 |
---|
| 423 | 0.47_wp, 0.47E-2_wp, 0.47E-2_wp, & ! 7 |
---|
| 424 | 0.013_wp, 0.013E-2_wp, 0.013E-2_wp, & ! 8 |
---|
| 425 | 0.034_wp, 0.034E-2_wp, 0.034E-2_wp, & ! 9 |
---|
| 426 | 0.5_wp, 0.50E-2_wp, 0.50E-2_wp, & ! 10 |
---|
| 427 | 0.17_wp, 0.17E-2_wp, 0.17E-2_wp, & ! 11 |
---|
| 428 | 1.3E-3_wp, 1.3E-4_wp, 1.3E-4_wp, & ! 12 |
---|
| 429 | 0.83_wp, 0.83E-2_wp, 0.83E-2_wp, & ! 13 |
---|
| 430 | 0.00_wp, 0.00_wp, 0.00_wp, & ! 14 |
---|
| 431 | 0.00_wp, 0.00_wp, 0.00_wp, & ! 15 |
---|
| 432 | 0.10_wp, 0.10E-2_wp, 0.10E-2_wp, & ! 16 |
---|
| 433 | 0.25_wp, 0.25E-2_wp, 0.25E-2_wp, & ! 17 |
---|
| 434 | 2.00_wp, 2.00E-2_wp, 2.00E-2_wp, & ! 18 |
---|
| 435 | 1.10_wp, 1.10E-2_wp, 1.10E-2_wp, & ! 19 |
---|
| 436 | 1.0E-4_wp, 1.0E-5_wp, 1.0E-5_wp & ! 20 |
---|
| 437 | /), (/ 3, 20 /) ) |
---|
[1496] | 438 | |
---|
| 439 | ! |
---|
[1551] | 440 | !-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in |
---|
[1788] | 441 | REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ & |
---|
[1585] | 442 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 1 |
---|
| 443 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 2 |
---|
| 444 | 20.0_wp, 15.0_wp, 0.03_wp, & ! 3 |
---|
| 445 | 20.0_wp, 15.0_wp, 0.03_wp, & ! 4 |
---|
| 446 | 20.0_wp, 15.0_wp, 0.03_wp, & ! 5 |
---|
| 447 | 20.0_wp, 15.0_wp, 0.03_wp, & ! 6 |
---|
| 448 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 7 |
---|
| 449 | 15.0_wp, 15.0_wp, 0.00_wp, & ! 8 |
---|
| 450 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 9 |
---|
| 451 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 10 |
---|
| 452 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 11 |
---|
| 453 | 58.0_wp, 58.0_wp, 0.00_wp, & ! 12 |
---|
| 454 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 13 |
---|
[1788] | 455 | 1.0E10_wp, 1.0E10_wp, 0.00_wp, & ! 14 |
---|
| 456 | 1.0E10_wp, 1.0E10_wp, 0.00_wp, & ! 15 |
---|
[1585] | 457 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 16 |
---|
| 458 | 10.0_wp, 10.0_wp, 0.05_wp, & ! 17 |
---|
| 459 | 20.0_wp, 15.0_wp, 0.03_wp, & ! 18 |
---|
[1788] | 460 | 20.0_wp, 15.0_wp, 0.03_wp, & ! 19 |
---|
| 461 | 0.0_wp, 0.0_wp, 0.00_wp & ! 20 |
---|
| 462 | /), (/ 3, 20 /) ) |
---|
[1496] | 463 | |
---|
| 464 | ! |
---|
| 465 | !-- Root distribution (sum = 1) level 1, level 2, level 3, level 4, |
---|
[1788] | 466 | REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ & |
---|
[1585] | 467 | 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 1 |
---|
| 468 | 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 2 |
---|
| 469 | 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & ! 3 |
---|
| 470 | 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & ! 4 |
---|
| 471 | 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & ! 5 |
---|
| 472 | 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & ! 6 |
---|
| 473 | 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & ! 7 |
---|
| 474 | 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 8 |
---|
| 475 | 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 9 |
---|
| 476 | 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 10 |
---|
| 477 | 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 11 |
---|
| 478 | 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12 |
---|
| 479 | 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 13 |
---|
| 480 | 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14 |
---|
| 481 | 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15 |
---|
| 482 | 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16 |
---|
| 483 | 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17 |
---|
| 484 | 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18 |
---|
[1788] | 485 | 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 19 |
---|
| 486 | 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp & ! 20 |
---|
| 487 | /), (/ 4, 20 /) ) |
---|
[1496] | 488 | |
---|
| 489 | ! |
---|
| 490 | !-- Soil parameters according to the following porosity classes (soil_type) |
---|
[1551] | 491 | |
---|
[1496] | 492 | ! |
---|
[1551] | 493 | !-- Soil parameters I alpha_vg, l_vg, n_vg, gamma_w_sat |
---|
[1585] | 494 | REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/ & |
---|
[1496] | 495 | 3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, & ! 1 |
---|
| 496 | 3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, & ! 2 |
---|
| 497 | 0.83_wp, -0.588_wp, 1.25_wp, 0.26E-6_wp, & ! 3 |
---|
| 498 | 3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, & ! 4 |
---|
| 499 | 2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, & ! 5 |
---|
[1551] | 500 | 1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp, & ! 6 |
---|
| 501 | 0.00_wp, 0.00_wp, 0.00_wp, 0.57E-6_wp & ! 7 |
---|
| 502 | /), (/ 4, 7 /) ) |
---|
[1496] | 503 | |
---|
| 504 | ! |
---|
| 505 | !-- Soil parameters II m_sat, m_fc, m_wilt, m_res |
---|
[1585] | 506 | REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ & |
---|
[1496] | 507 | 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1 |
---|
| 508 | 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2 |
---|
| 509 | 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3 |
---|
| 510 | 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4 |
---|
| 511 | 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5 |
---|
[1551] | 512 | 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6 |
---|
| 513 | 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp & ! 7 |
---|
| 514 | /), (/ 4, 7 /) ) |
---|
[1496] | 515 | |
---|
| 516 | |
---|
| 517 | SAVE |
---|
| 518 | |
---|
| 519 | |
---|
| 520 | PRIVATE |
---|
| 521 | |
---|
| 522 | |
---|
[1551] | 523 | ! |
---|
| 524 | !-- Public parameters, constants and initial values |
---|
| 525 | PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient, & |
---|
[1788] | 526 | conserve_water_content, field_capacity, & |
---|
[1551] | 527 | f_shortwave_incoming, hydraulic_conductivity, init_lsm, & |
---|
| 528 | init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable, & |
---|
| 529 | land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model, & |
---|
| 530 | lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance, & |
---|
[1788] | 531 | min_soil_resistance, n_vangenuchten, pave_heat_capacity, & |
---|
| 532 | pave_depth, pave_heat_conductivity, residual_moisture, rho_cp, & |
---|
[1691] | 533 | rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm, & |
---|
| 534 | soil_moisture, soil_temperature, soil_type, soil_type_name, & |
---|
| 535 | vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, & |
---|
[1788] | 536 | z0h_eb, z0q_eb |
---|
[1496] | 537 | |
---|
[1551] | 538 | ! |
---|
[1783] | 539 | !-- Public grid variables |
---|
| 540 | PUBLIC nzb_soil, nzs, nzt_soil, zs |
---|
[1496] | 541 | |
---|
[1551] | 542 | ! |
---|
| 543 | !-- Public 2D output variables |
---|
| 544 | PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av, & |
---|
| 545 | lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av, & |
---|
| 546 | qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av, & |
---|
[1555] | 547 | r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av |
---|
[1551] | 548 | |
---|
| 549 | ! |
---|
| 550 | !-- Public prognostic variables |
---|
| 551 | PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av |
---|
| 552 | |
---|
[1496] | 553 | INTERFACE init_lsm |
---|
| 554 | MODULE PROCEDURE init_lsm |
---|
| 555 | END INTERFACE init_lsm |
---|
| 556 | |
---|
| 557 | INTERFACE lsm_energy_balance |
---|
| 558 | MODULE PROCEDURE lsm_energy_balance |
---|
| 559 | END INTERFACE lsm_energy_balance |
---|
| 560 | |
---|
| 561 | INTERFACE lsm_soil_model |
---|
| 562 | MODULE PROCEDURE lsm_soil_model |
---|
| 563 | END INTERFACE lsm_soil_model |
---|
| 564 | |
---|
[1551] | 565 | INTERFACE lsm_swap_timelevel |
---|
| 566 | MODULE PROCEDURE lsm_swap_timelevel |
---|
| 567 | END INTERFACE lsm_swap_timelevel |
---|
[1496] | 568 | |
---|
| 569 | CONTAINS |
---|
| 570 | |
---|
| 571 | |
---|
| 572 | !------------------------------------------------------------------------------! |
---|
| 573 | ! Description: |
---|
| 574 | ! ------------ |
---|
[1682] | 575 | !> Allocate land surface model arrays and define pointers |
---|
[1496] | 576 | !------------------------------------------------------------------------------! |
---|
[1551] | 577 | SUBROUTINE init_lsm_arrays |
---|
[1496] | 578 | |
---|
| 579 | |
---|
| 580 | IMPLICIT NONE |
---|
| 581 | |
---|
| 582 | ! |
---|
[1551] | 583 | !-- Allocate surface and soil temperature / humidity |
---|
[1496] | 584 | #if defined( __nopointer ) |
---|
[1551] | 585 | ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) ) |
---|
| 586 | ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) ) |
---|
| 587 | ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
| 588 | ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
| 589 | ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) ) |
---|
| 590 | ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) ) |
---|
| 591 | ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 592 | ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 593 | #else |
---|
[1551] | 594 | ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) ) |
---|
| 595 | ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) ) |
---|
| 596 | ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
| 597 | ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
| 598 | ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) ) |
---|
| 599 | ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) ) |
---|
| 600 | ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 601 | ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 602 | #endif |
---|
| 603 | |
---|
| 604 | ! |
---|
[1551] | 605 | !-- Allocate intermediate timestep arrays |
---|
| 606 | ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) ) |
---|
| 607 | ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
| 608 | ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) ) |
---|
| 609 | ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 610 | |
---|
| 611 | ! |
---|
| 612 | !-- Allocate 2D vegetation model arrays |
---|
[1551] | 613 | ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) ) |
---|
[1788] | 614 | ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 615 | ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) ) |
---|
| 616 | ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) ) |
---|
[1551] | 617 | ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) ) |
---|
| 618 | ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 619 | ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) ) |
---|
[1551] | 620 | ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) ) |
---|
| 621 | ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) ) |
---|
| 622 | ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) ) |
---|
| 623 | ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) ) |
---|
| 624 | ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 625 | ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) ) |
---|
| 626 | ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) ) |
---|
| 627 | ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) ) |
---|
| 628 | ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) ) |
---|
[1551] | 629 | ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) ) |
---|
[1788] | 630 | ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) ) |
---|
[1551] | 631 | ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) ) |
---|
| 632 | ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) ) |
---|
| 633 | ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) ) |
---|
| 634 | ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) ) |
---|
[1585] | 635 | ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 636 | ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) ) |
---|
| 637 | ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) ) |
---|
| 638 | ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) ) |
---|
| 639 | ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) ) |
---|
| 640 | ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) ) |
---|
[1551] | 641 | ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) ) |
---|
| 642 | ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) ) |
---|
[1788] | 643 | ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) ) |
---|
| 644 | ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) ) |
---|
| 645 | ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 646 | |
---|
[1551] | 647 | #if ! defined( __nopointer ) |
---|
[1496] | 648 | ! |
---|
[1585] | 649 | !-- Initial assignment of the pointers |
---|
[1551] | 650 | t_soil => t_soil_1; t_soil_p => t_soil_2 |
---|
| 651 | t_surface => t_surface_1; t_surface_p => t_surface_2 |
---|
| 652 | m_soil => m_soil_1; m_soil_p => m_soil_2 |
---|
| 653 | m_liq_eb => m_liq_eb_1; m_liq_eb_p => m_liq_eb_2 |
---|
| 654 | #endif |
---|
[1496] | 655 | |
---|
| 656 | |
---|
[1551] | 657 | END SUBROUTINE init_lsm_arrays |
---|
[1500] | 658 | |
---|
[1551] | 659 | !------------------------------------------------------------------------------! |
---|
| 660 | ! Description: |
---|
| 661 | ! ------------ |
---|
[1682] | 662 | !> Initialization of the land surface model |
---|
[1551] | 663 | !------------------------------------------------------------------------------! |
---|
| 664 | SUBROUTINE init_lsm |
---|
| 665 | |
---|
| 666 | |
---|
| 667 | IMPLICIT NONE |
---|
| 668 | |
---|
[1682] | 669 | INTEGER(iwp) :: i !< running index |
---|
| 670 | INTEGER(iwp) :: j !< running index |
---|
| 671 | INTEGER(iwp) :: k !< running index |
---|
[1551] | 672 | |
---|
[1709] | 673 | REAL(wp) :: pt1 !< potential temperature at first grid level |
---|
[1551] | 674 | |
---|
[1691] | 675 | |
---|
[1551] | 676 | ! |
---|
| 677 | !-- Calculate Exner function |
---|
[1691] | 678 | exn = ( surface_pressure / 1000.0_wp )**0.286_wp |
---|
[1551] | 679 | |
---|
| 680 | |
---|
| 681 | ! |
---|
| 682 | !-- If no cloud physics is used, rho_surface has not been calculated before |
---|
[1788] | 683 | IF ( .NOT. cloud_physics ) THEN |
---|
[1551] | 684 | rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn ) |
---|
| 685 | ENDIF |
---|
| 686 | |
---|
| 687 | ! |
---|
| 688 | !-- Calculate frequently used parameters |
---|
| 689 | rho_cp = cp * rho_surface |
---|
| 690 | rd_d_rv = r_d / r_v |
---|
| 691 | rho_lv = rho_surface * l_v |
---|
| 692 | drho_l_lv = 1.0_wp / (rho_l * l_v) |
---|
| 693 | |
---|
| 694 | ! |
---|
| 695 | !-- Set inital values for prognostic quantities |
---|
| 696 | tt_surface_m = 0.0_wp |
---|
| 697 | tt_soil_m = 0.0_wp |
---|
[1757] | 698 | tm_soil_m = 0.0_wp |
---|
[1551] | 699 | tm_liq_eb_m = 0.0_wp |
---|
| 700 | c_liq = 0.0_wp |
---|
| 701 | |
---|
| 702 | ghf_eb = 0.0_wp |
---|
| 703 | shf_eb = rho_cp * shf |
---|
| 704 | |
---|
[1500] | 705 | IF ( humidity ) THEN |
---|
[1551] | 706 | qsws_eb = rho_l * l_v * qsws |
---|
[1500] | 707 | ELSE |
---|
[1551] | 708 | qsws_eb = 0.0_wp |
---|
[1500] | 709 | ENDIF |
---|
| 710 | |
---|
[1551] | 711 | qsws_liq_eb = 0.0_wp |
---|
[1757] | 712 | qsws_soil_eb = 0.0_wp |
---|
[1551] | 713 | qsws_veg_eb = 0.0_wp |
---|
[1496] | 714 | |
---|
| 715 | r_a = 50.0_wp |
---|
[1551] | 716 | r_s = 50.0_wp |
---|
[1496] | 717 | r_canopy = 0.0_wp |
---|
| 718 | r_soil = 0.0_wp |
---|
| 719 | |
---|
| 720 | ! |
---|
| 721 | !-- Allocate 3D soil model arrays |
---|
[1551] | 722 | ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
| 723 | ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
| 724 | ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 725 | |
---|
| 726 | lambda_h = 0.0_wp |
---|
| 727 | ! |
---|
| 728 | !-- If required, allocate humidity-related variables for the soil model |
---|
| 729 | IF ( humidity ) THEN |
---|
[1551] | 730 | ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
| 731 | ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) |
---|
[1496] | 732 | |
---|
| 733 | lambda_w = 0.0_wp |
---|
| 734 | ENDIF |
---|
| 735 | |
---|
| 736 | ! |
---|
| 737 | !-- Calculate grid spacings. Temperature and moisture are defined at |
---|
[1691] | 738 | !-- the edges of the soil layers (_stag), whereas gradients/fluxes are defined |
---|
| 739 | !-- at the centers |
---|
| 740 | dz_soil(nzb_soil) = zs(nzb_soil) |
---|
[1496] | 741 | |
---|
[1691] | 742 | DO k = nzb_soil+1, nzt_soil |
---|
| 743 | dz_soil(k) = zs(k) - zs(k-1) |
---|
[1496] | 744 | ENDDO |
---|
[1691] | 745 | dz_soil(nzt_soil+1) = dz_soil(nzt_soil) |
---|
[1496] | 746 | |
---|
[1691] | 747 | DO k = nzb_soil, nzt_soil-1 |
---|
| 748 | dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k)) |
---|
[1496] | 749 | ENDDO |
---|
[1691] | 750 | dz_soil_stag(nzt_soil) = dz_soil(nzt_soil) |
---|
[1496] | 751 | |
---|
[1551] | 752 | ddz_soil = 1.0_wp / dz_soil |
---|
| 753 | ddz_soil_stag = 1.0_wp / dz_soil_stag |
---|
| 754 | |
---|
[1496] | 755 | ! |
---|
[1551] | 756 | !-- Initialize standard soil types. It is possible to overwrite each |
---|
| 757 | !-- parameter by setting the respecticy NAMELIST variable to a |
---|
| 758 | !-- value /= 9999999.9. |
---|
| 759 | IF ( soil_type /= 0 ) THEN |
---|
| 760 | |
---|
| 761 | IF ( alpha_vangenuchten == 9999999.9_wp ) THEN |
---|
| 762 | alpha_vangenuchten = soil_pars(0,soil_type) |
---|
| 763 | ENDIF |
---|
| 764 | |
---|
| 765 | IF ( l_vangenuchten == 9999999.9_wp ) THEN |
---|
| 766 | l_vangenuchten = soil_pars(1,soil_type) |
---|
| 767 | ENDIF |
---|
| 768 | |
---|
| 769 | IF ( n_vangenuchten == 9999999.9_wp ) THEN |
---|
| 770 | n_vangenuchten = soil_pars(2,soil_type) |
---|
| 771 | ENDIF |
---|
| 772 | |
---|
| 773 | IF ( hydraulic_conductivity == 9999999.9_wp ) THEN |
---|
| 774 | hydraulic_conductivity = soil_pars(3,soil_type) |
---|
| 775 | ENDIF |
---|
| 776 | |
---|
| 777 | IF ( saturation_moisture == 9999999.9_wp ) THEN |
---|
| 778 | saturation_moisture = m_soil_pars(0,soil_type) |
---|
| 779 | ENDIF |
---|
| 780 | |
---|
| 781 | IF ( field_capacity == 9999999.9_wp ) THEN |
---|
| 782 | field_capacity = m_soil_pars(1,soil_type) |
---|
| 783 | ENDIF |
---|
| 784 | |
---|
| 785 | IF ( wilting_point == 9999999.9_wp ) THEN |
---|
| 786 | wilting_point = m_soil_pars(2,soil_type) |
---|
| 787 | ENDIF |
---|
| 788 | |
---|
| 789 | IF ( residual_moisture == 9999999.9_wp ) THEN |
---|
| 790 | residual_moisture = m_soil_pars(3,soil_type) |
---|
| 791 | ENDIF |
---|
| 792 | |
---|
[1496] | 793 | ENDIF |
---|
| 794 | |
---|
[1788] | 795 | ! |
---|
| 796 | !-- Map values to the respective 2D arrays |
---|
[1551] | 797 | alpha_vg = alpha_vangenuchten |
---|
| 798 | l_vg = l_vangenuchten |
---|
| 799 | n_vg = n_vangenuchten |
---|
| 800 | gamma_w_sat = hydraulic_conductivity |
---|
| 801 | m_sat = saturation_moisture |
---|
| 802 | m_fc = field_capacity |
---|
| 803 | m_wilt = wilting_point |
---|
| 804 | m_res = residual_moisture |
---|
| 805 | r_soil_min = min_soil_resistance |
---|
| 806 | |
---|
[1496] | 807 | ! |
---|
[1551] | 808 | !-- Initial run actions |
---|
| 809 | IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN |
---|
[1496] | 810 | |
---|
[1551] | 811 | t_soil = 0.0_wp |
---|
| 812 | m_liq_eb = 0.0_wp |
---|
| 813 | m_soil = 0.0_wp |
---|
[1496] | 814 | |
---|
[1551] | 815 | ! |
---|
| 816 | !-- Map user settings of T and q for each soil layer |
---|
| 817 | !-- (make sure that the soil moisture does not drop below the permanent |
---|
| 818 | !-- wilting point) -> problems with devision by zero) |
---|
[1691] | 819 | DO k = nzb_soil, nzt_soil |
---|
[1551] | 820 | t_soil(k,:,:) = soil_temperature(k) |
---|
| 821 | m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:)) |
---|
| 822 | soil_moisture(k) = MAX(soil_moisture(k),wilting_point) |
---|
| 823 | ENDDO |
---|
| 824 | t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1) |
---|
[1496] | 825 | |
---|
[1551] | 826 | ! |
---|
| 827 | !-- Calculate surface temperature |
---|
[1691] | 828 | t_surface = pt_surface * exn |
---|
[1496] | 829 | |
---|
| 830 | ! |
---|
[1788] | 831 | !-- Set artifical values for ts and us so that r_a has its initial value |
---|
| 832 | !-- for the first time step |
---|
[1551] | 833 | DO i = nxlg, nxrg |
---|
| 834 | DO j = nysg, nyng |
---|
| 835 | k = nzb_s_inner(j,i) |
---|
| 836 | |
---|
[1691] | 837 | IF ( cloud_physics ) THEN |
---|
[1709] | 838 | pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) |
---|
[1691] | 839 | ELSE |
---|
[1709] | 840 | pt1 = pt(k+1,j,i) |
---|
[1691] | 841 | ENDIF |
---|
| 842 | |
---|
[1551] | 843 | ! |
---|
| 844 | !-- Assure that r_a cannot be zero at model start |
---|
[1709] | 845 | IF ( pt1 == pt(k,j,i) ) pt1 = pt1 + 1.0E-10_wp |
---|
[1551] | 846 | |
---|
[1691] | 847 | us(j,i) = 0.1_wp |
---|
[1709] | 848 | ts(j,i) = (pt1 - pt(k,j,i)) / r_a(j,i) |
---|
[1551] | 849 | shf(j,i) = - us(j,i) * ts(j,i) |
---|
| 850 | ENDDO |
---|
| 851 | ENDDO |
---|
| 852 | |
---|
| 853 | ! |
---|
| 854 | !-- Actions for restart runs |
---|
| 855 | ELSE |
---|
| 856 | |
---|
| 857 | DO i = nxlg, nxrg |
---|
| 858 | DO j = nysg, nyng |
---|
| 859 | k = nzb_s_inner(j,i) |
---|
| 860 | t_surface(j,i) = pt(k,j,i) * exn |
---|
| 861 | ENDDO |
---|
| 862 | ENDDO |
---|
| 863 | |
---|
| 864 | ENDIF |
---|
| 865 | |
---|
[1691] | 866 | DO k = nzb_soil, nzt_soil |
---|
[1551] | 867 | root_fr(k,:,:) = root_fraction(k) |
---|
| 868 | ENDDO |
---|
[1496] | 869 | |
---|
[1551] | 870 | IF ( veg_type /= 0 ) THEN |
---|
| 871 | IF ( min_canopy_resistance == 9999999.9_wp ) THEN |
---|
| 872 | min_canopy_resistance = veg_pars(0,veg_type) |
---|
| 873 | ENDIF |
---|
| 874 | IF ( leaf_area_index == 9999999.9_wp ) THEN |
---|
| 875 | leaf_area_index = veg_pars(1,veg_type) |
---|
| 876 | ENDIF |
---|
| 877 | IF ( vegetation_coverage == 9999999.9_wp ) THEN |
---|
| 878 | vegetation_coverage = veg_pars(2,veg_type) |
---|
| 879 | ENDIF |
---|
| 880 | IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN |
---|
| 881 | canopy_resistance_coefficient= veg_pars(3,veg_type) |
---|
| 882 | ENDIF |
---|
| 883 | IF ( lambda_surface_stable == 9999999.9_wp ) THEN |
---|
| 884 | lambda_surface_stable = surface_pars(0,veg_type) |
---|
| 885 | ENDIF |
---|
| 886 | IF ( lambda_surface_unstable == 9999999.9_wp ) THEN |
---|
| 887 | lambda_surface_unstable = surface_pars(1,veg_type) |
---|
| 888 | ENDIF |
---|
| 889 | IF ( f_shortwave_incoming == 9999999.9_wp ) THEN |
---|
| 890 | f_shortwave_incoming = surface_pars(2,veg_type) |
---|
| 891 | ENDIF |
---|
| 892 | IF ( z0_eb == 9999999.9_wp ) THEN |
---|
| 893 | roughness_length = roughness_par(0,veg_type) |
---|
| 894 | z0_eb = roughness_par(0,veg_type) |
---|
| 895 | ENDIF |
---|
| 896 | IF ( z0h_eb == 9999999.9_wp ) THEN |
---|
| 897 | z0h_eb = roughness_par(1,veg_type) |
---|
| 898 | ENDIF |
---|
[1788] | 899 | IF ( z0q_eb == 9999999.9_wp ) THEN |
---|
| 900 | z0q_eb = roughness_par(2,veg_type) |
---|
| 901 | ENDIF |
---|
| 902 | z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp ) |
---|
[1496] | 903 | |
---|
[1551] | 904 | IF ( ANY( root_fraction == 9999999.9_wp ) ) THEN |
---|
[1691] | 905 | DO k = nzb_soil, nzt_soil |
---|
[1551] | 906 | root_fr(k,:,:) = root_distribution(k,veg_type) |
---|
| 907 | root_fraction(k) = root_distribution(k,veg_type) |
---|
| 908 | ENDDO |
---|
| 909 | ENDIF |
---|
[1496] | 910 | |
---|
[1553] | 911 | ELSE |
---|
| 912 | |
---|
| 913 | IF ( z0_eb == 9999999.9_wp ) THEN |
---|
| 914 | z0_eb = roughness_length |
---|
| 915 | ENDIF |
---|
| 916 | IF ( z0h_eb == 9999999.9_wp ) THEN |
---|
| 917 | z0h_eb = z0_eb * z0h_factor |
---|
| 918 | ENDIF |
---|
[1788] | 919 | IF ( z0q_eb == 9999999.9_wp ) THEN |
---|
| 920 | z0q_eb = z0_eb * z0h_factor |
---|
| 921 | ENDIF |
---|
[1553] | 922 | |
---|
[1496] | 923 | ENDIF |
---|
| 924 | |
---|
| 925 | ! |
---|
[1788] | 926 | !-- For surfaces covered with pavement, set depth of the pavement (with dry |
---|
| 927 | !-- soil below). The depth must be greater than the first soil layer depth |
---|
| 928 | IF ( veg_type == 20 ) THEN |
---|
| 929 | IF ( pave_depth == 9999999.9_wp ) THEN |
---|
| 930 | pave_depth = zs(nzb_soil) |
---|
| 931 | ELSE |
---|
| 932 | pave_depth = MAX( zs(nzb_soil), pave_depth ) |
---|
| 933 | ENDIF |
---|
| 934 | ENDIF |
---|
| 935 | |
---|
| 936 | ! |
---|
| 937 | !-- Map vegetation and soil types to 2D array to allow for heterogeneous |
---|
| 938 | !-- surfaces via user interface see below |
---|
| 939 | veg_type_2d = veg_type |
---|
| 940 | soil_type_2d = soil_type |
---|
| 941 | |
---|
| 942 | ! |
---|
| 943 | !-- Map vegetation parameters to the respective 2D arrays |
---|
[1551] | 944 | r_canopy_min = min_canopy_resistance |
---|
| 945 | lai = leaf_area_index |
---|
| 946 | c_veg = vegetation_coverage |
---|
| 947 | g_d = canopy_resistance_coefficient |
---|
| 948 | lambda_surface_s = lambda_surface_stable |
---|
| 949 | lambda_surface_u = lambda_surface_unstable |
---|
| 950 | f_sw_in = f_shortwave_incoming |
---|
| 951 | z0 = z0_eb |
---|
| 952 | z0h = z0h_eb |
---|
[1788] | 953 | z0q = z0q_eb |
---|
[1551] | 954 | |
---|
| 955 | ! |
---|
[1496] | 956 | !-- Possibly do user-defined actions (e.g. define heterogeneous land surface) |
---|
| 957 | CALL user_init_land_surface |
---|
| 958 | |
---|
[1788] | 959 | ! |
---|
| 960 | !-- Set flag parameter if vegetation type was set to a water surface. Also |
---|
| 961 | !-- set temperature to a constant value in all "soil" layers. |
---|
| 962 | DO i = nxlg, nxrg |
---|
| 963 | DO j = nysg, nyng |
---|
| 964 | IF ( veg_type_2d(j,i) == 14 .OR. veg_type_2d(j,i) == 15 ) THEN |
---|
| 965 | water_surface(j,i) = .TRUE. |
---|
| 966 | ELSEIF ( veg_type_2d(j,i) == 20 ) THEN |
---|
| 967 | pave_surface(j,i) = .TRUE. |
---|
| 968 | m_soil(:,j,i) = 0.0_wp |
---|
| 969 | ENDIF |
---|
[1500] | 970 | |
---|
[1788] | 971 | ENDDO |
---|
| 972 | ENDDO |
---|
| 973 | |
---|
| 974 | ! |
---|
| 975 | !-- Calculate new roughness lengths (for water surfaces only) |
---|
| 976 | CALL calc_z0_water_surface |
---|
| 977 | |
---|
[1551] | 978 | t_soil_p = t_soil |
---|
| 979 | m_soil_p = m_soil |
---|
| 980 | m_liq_eb_p = m_liq_eb |
---|
[1709] | 981 | t_surface_p = t_surface |
---|
[1496] | 982 | |
---|
[1709] | 983 | |
---|
| 984 | |
---|
[1551] | 985 | !-- Store initial profiles of t_soil and m_soil (assuming they are |
---|
| 986 | !-- horizontally homogeneous on this PE) |
---|
| 987 | hom(nzb_soil:nzt_soil,1,90,:) = SPREAD( t_soil(nzb_soil:nzt_soil, & |
---|
| 988 | nysg,nxlg), 2, & |
---|
| 989 | statistic_regions+1 ) |
---|
| 990 | hom(nzb_soil:nzt_soil,1,92,:) = SPREAD( m_soil(nzb_soil:nzt_soil, & |
---|
| 991 | nysg,nxlg), 2, & |
---|
| 992 | statistic_regions+1 ) |
---|
| 993 | |
---|
[1496] | 994 | END SUBROUTINE init_lsm |
---|
| 995 | |
---|
| 996 | |
---|
| 997 | |
---|
| 998 | !------------------------------------------------------------------------------! |
---|
| 999 | ! Description: |
---|
| 1000 | ! ------------ |
---|
[1682] | 1001 | !> Solver for the energy balance at the surface. |
---|
[1496] | 1002 | !------------------------------------------------------------------------------! |
---|
| 1003 | SUBROUTINE lsm_energy_balance |
---|
| 1004 | |
---|
| 1005 | |
---|
| 1006 | IMPLICIT NONE |
---|
| 1007 | |
---|
[1682] | 1008 | INTEGER(iwp) :: i !< running index |
---|
| 1009 | INTEGER(iwp) :: j !< running index |
---|
| 1010 | INTEGER(iwp) :: k, ks !< running index |
---|
[1496] | 1011 | |
---|
[1788] | 1012 | REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface |
---|
| 1013 | f1, & !< resistance correction term 1 |
---|
[1682] | 1014 | f2, & !< resistance correction term 2 |
---|
| 1015 | f3, & !< resistance correction term 3 |
---|
| 1016 | m_min, & !< minimum soil moisture |
---|
| 1017 | e, & !< water vapour pressure |
---|
| 1018 | e_s, & !< water vapour saturation pressure |
---|
| 1019 | e_s_dt, & !< derivate of e_s with respect to T |
---|
| 1020 | tend, & !< tendency |
---|
| 1021 | dq_s_dt, & !< derivate of q_s with respect to T |
---|
| 1022 | coef_1, & !< coef. for prognostic equation |
---|
| 1023 | coef_2, & !< coef. for prognostic equation |
---|
| 1024 | f_qsws, & !< factor for qsws_eb |
---|
| 1025 | f_qsws_veg, & !< factor for qsws_veg_eb |
---|
| 1026 | f_qsws_soil, & !< factor for qsws_soil_eb |
---|
| 1027 | f_qsws_liq, & !< factor for qsws_liq_eb |
---|
| 1028 | f_shf, & !< factor for shf_eb |
---|
| 1029 | lambda_surface, & !< Current value of lambda_surface |
---|
[1691] | 1030 | m_liq_eb_max, & !< maxmimum value of the liq. water reservoir |
---|
[1709] | 1031 | pt1, & !< potential temperature at first grid level |
---|
| 1032 | qv1 !< specific humidity at first grid level |
---|
[1496] | 1033 | |
---|
| 1034 | ! |
---|
| 1035 | !-- Calculate the exner function for the current time step |
---|
| 1036 | exn = ( surface_pressure / 1000.0_wp )**0.286_wp |
---|
| 1037 | |
---|
[1691] | 1038 | DO i = nxlg, nxrg |
---|
| 1039 | DO j = nysg, nyng |
---|
[1500] | 1040 | k = nzb_s_inner(j,i) |
---|
[1496] | 1041 | |
---|
| 1042 | ! |
---|
[1788] | 1043 | !-- Set lambda_surface according to stratification between skin layer and soil |
---|
| 1044 | IF ( .NOT. pave_surface(j,i) ) THEN |
---|
| 1045 | |
---|
| 1046 | c_surface_tmp = c_surface |
---|
| 1047 | |
---|
| 1048 | IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i)) THEN |
---|
| 1049 | lambda_surface = lambda_surface_s(j,i) |
---|
| 1050 | ELSE |
---|
| 1051 | lambda_surface = lambda_surface_u(j,i) |
---|
| 1052 | ENDIF |
---|
[1496] | 1053 | ELSE |
---|
[1788] | 1054 | |
---|
| 1055 | c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp |
---|
| 1056 | lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil) |
---|
| 1057 | |
---|
[1496] | 1058 | ENDIF |
---|
[1500] | 1059 | |
---|
[1496] | 1060 | ! |
---|
[1500] | 1061 | !-- First step: calculate aerodyamic resistance. As pt, us, ts |
---|
| 1062 | !-- are not available for the prognostic time step, data from the last |
---|
| 1063 | !-- time step is used here. Note that this formulation is the |
---|
| 1064 | !-- equivalent to the ECMWF formulation using drag coefficients |
---|
[1691] | 1065 | IF ( cloud_physics ) THEN |
---|
[1709] | 1066 | pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) |
---|
| 1067 | qv1 = q(k+1,j,i) - ql(k+1,j,i) |
---|
[1691] | 1068 | ELSE |
---|
[1709] | 1069 | pt1 = pt(k+1,j,i) |
---|
| 1070 | qv1 = q(k+1,j,i) |
---|
[1691] | 1071 | ENDIF |
---|
[1496] | 1072 | |
---|
[1709] | 1073 | r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp) |
---|
[1691] | 1074 | |
---|
[1496] | 1075 | ! |
---|
[1709] | 1076 | !-- Make sure that the resistance does not drop to zero |
---|
| 1077 | IF ( ABS(r_a(j,i)) < 1.0E-10_wp ) r_a(j,i) = 1.0E-10_wp |
---|
| 1078 | |
---|
| 1079 | ! |
---|
[1496] | 1080 | !-- Second step: calculate canopy resistance r_canopy |
---|
| 1081 | !-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation |
---|
| 1082 | |
---|
[1551] | 1083 | !-- f1: correction for incoming shortwave radiation (stomata close at |
---|
| 1084 | !-- night) |
---|
[1585] | 1085 | IF ( radiation_scheme /= 'constant' ) THEN |
---|
[1709] | 1086 | f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) / & |
---|
| 1087 | (0.81_wp * (0.004_wp * rad_sw_in(k,j,i) & |
---|
[1585] | 1088 | + 1.0_wp)) ) |
---|
[1551] | 1089 | ELSE |
---|
| 1090 | f1 = 1.0_wp |
---|
| 1091 | ENDIF |
---|
[1496] | 1092 | |
---|
[1709] | 1093 | |
---|
[1496] | 1094 | ! |
---|
[1551] | 1095 | !-- f2: correction for soil moisture availability to plants (the |
---|
| 1096 | !-- integrated soil moisture must thus be considered here) |
---|
| 1097 | !-- f2 = 0 for very dry soils |
---|
[1496] | 1098 | m_total = 0.0_wp |
---|
[1691] | 1099 | DO ks = nzb_soil, nzt_soil |
---|
[1551] | 1100 | m_total = m_total + root_fr(ks,j,i) & |
---|
| 1101 | * MAX(m_soil(ks,j,i),m_wilt(j,i)) |
---|
[1496] | 1102 | ENDDO |
---|
| 1103 | |
---|
[1788] | 1104 | IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) ) THEN |
---|
[1496] | 1105 | f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) |
---|
[1691] | 1106 | ELSEIF ( m_total >= m_fc(j,i) ) THEN |
---|
[1551] | 1107 | f2 = 1.0_wp |
---|
[1496] | 1108 | ELSE |
---|
| 1109 | f2 = 1.0E-20_wp |
---|
| 1110 | ENDIF |
---|
| 1111 | |
---|
| 1112 | ! |
---|
| 1113 | !-- Calculate water vapour pressure at saturation |
---|
[1551] | 1114 | e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) & |
---|
| 1115 | - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) ) |
---|
[1496] | 1116 | |
---|
| 1117 | ! |
---|
| 1118 | !-- f3: correction for vapour pressure deficit |
---|
[1551] | 1119 | IF ( g_d(j,i) /= 0.0_wp ) THEN |
---|
[1496] | 1120 | ! |
---|
| 1121 | !-- Calculate vapour pressure |
---|
[1709] | 1122 | e = qv1 * surface_pressure / 0.622_wp |
---|
[1551] | 1123 | f3 = EXP ( -g_d(j,i) * (e_s - e) ) |
---|
[1496] | 1124 | ELSE |
---|
| 1125 | f3 = 1.0_wp |
---|
| 1126 | ENDIF |
---|
| 1127 | |
---|
| 1128 | ! |
---|
[1551] | 1129 | !-- Calculate canopy resistance. In case that c_veg is 0 (bare soils), |
---|
| 1130 | !-- this calculation is obsolete, as r_canopy is not used below. |
---|
[1496] | 1131 | !-- To do: check for very dry soil -> r_canopy goes to infinity |
---|
[1551] | 1132 | r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3 & |
---|
| 1133 | + 1.0E-20_wp) |
---|
[1496] | 1134 | |
---|
| 1135 | ! |
---|
[1551] | 1136 | !-- Third step: calculate bare soil resistance r_soil. The Clapp & |
---|
| 1137 | !-- Hornberger parametrization does not consider c_veg. |
---|
[1788] | 1138 | IF ( soil_type_2d(j,i) /= 7 ) THEN |
---|
[1551] | 1139 | m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & |
---|
| 1140 | m_res(j,i) |
---|
| 1141 | ELSE |
---|
| 1142 | m_min = m_wilt(j,i) |
---|
| 1143 | ENDIF |
---|
[1496] | 1144 | |
---|
[1551] | 1145 | f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min ) |
---|
[1513] | 1146 | f2 = MAX(f2,1.0E-20_wp) |
---|
[1551] | 1147 | f2 = MIN(f2,1.0_wp) |
---|
[1496] | 1148 | |
---|
| 1149 | r_soil(j,i) = r_soil_min(j,i) / f2 |
---|
| 1150 | |
---|
| 1151 | ! |
---|
[1788] | 1152 | !-- Calculate the maximum possible liquid water amount on plants and |
---|
| 1153 | !-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is |
---|
| 1154 | !-- assumed, while paved surfaces might hold up 1 mm of water. The |
---|
| 1155 | !-- liquid water fraction for paved surfaces is calculated after |
---|
| 1156 | !-- Noilhan & Planton (1989), while the ECMWF formulation is used for |
---|
| 1157 | !-- vegetated surfaces and bare soils. |
---|
| 1158 | IF ( pave_surface(j,i) ) THEN |
---|
| 1159 | m_liq_eb_max = m_max_depth * 5.0_wp |
---|
| 1160 | c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 ) |
---|
| 1161 | ELSE |
---|
| 1162 | m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i) & |
---|
| 1163 | + (1.0_wp - c_veg(j,i)) ) |
---|
| 1164 | c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max ) |
---|
| 1165 | ENDIF |
---|
[1496] | 1166 | |
---|
[1551] | 1167 | ! |
---|
| 1168 | !-- Calculate saturation specific humidity |
---|
[1496] | 1169 | q_s = 0.622_wp * e_s / surface_pressure |
---|
[1500] | 1170 | |
---|
| 1171 | ! |
---|
[1551] | 1172 | !-- In case of dewfall, set evapotranspiration to zero |
---|
| 1173 | !-- All super-saturated water is then removed from the air |
---|
[1788] | 1174 | IF ( humidity .AND. q_s <= qv1 ) THEN |
---|
[1551] | 1175 | r_canopy(j,i) = 0.0_wp |
---|
| 1176 | r_soil(j,i) = 0.0_wp |
---|
[1691] | 1177 | ENDIF |
---|
[1496] | 1178 | |
---|
| 1179 | ! |
---|
| 1180 | !-- Calculate coefficients for the total evapotranspiration |
---|
[1788] | 1181 | !-- In case of water surface, set vegetation and soil fluxes to zero. |
---|
| 1182 | !-- For pavements, only evaporation of liquid water is possible. |
---|
| 1183 | IF ( water_surface(j,i) ) THEN |
---|
| 1184 | f_qsws_veg = 0.0_wp |
---|
| 1185 | f_qsws_soil = 0.0_wp |
---|
| 1186 | f_qsws_liq = rho_lv / r_a(j,i) |
---|
| 1187 | ELSEIF ( pave_surface (j,i) ) THEN |
---|
| 1188 | f_qsws_veg = 0.0_wp |
---|
| 1189 | f_qsws_soil = 0.0_wp |
---|
| 1190 | f_qsws_liq = rho_lv * c_liq(j,i) / r_a(j,i) |
---|
| 1191 | ELSE |
---|
| 1192 | f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & |
---|
| 1193 | (r_a(j,i) + r_canopy(j,i)) |
---|
| 1194 | f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & |
---|
[1551] | 1195 | r_soil(j,i)) |
---|
[1788] | 1196 | f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) |
---|
| 1197 | ENDIF |
---|
[1500] | 1198 | ! |
---|
| 1199 | !-- If soil moisture is below wilting point, plants do no longer |
---|
| 1200 | !-- transpirate. |
---|
[1691] | 1201 | ! IF ( m_soil(k,j,i) < m_wilt(j,i) ) THEN |
---|
[1551] | 1202 | ! f_qsws_veg = 0.0_wp |
---|
| 1203 | ! ENDIF |
---|
[1496] | 1204 | |
---|
[1551] | 1205 | f_shf = rho_cp / r_a(j,i) |
---|
| 1206 | f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq |
---|
| 1207 | |
---|
| 1208 | ! |
---|
[1496] | 1209 | !-- Calculate derivative of q_s for Taylor series expansion |
---|
[1571] | 1210 | e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) - & |
---|
[1551] | 1211 | 17.269_wp*(t_surface(j,i) - 273.16_wp) & |
---|
| 1212 | / (t_surface(j,i) - 35.86_wp)**2 ) |
---|
[1496] | 1213 | |
---|
[1571] | 1214 | dq_s_dt = 0.622_wp * e_s_dt / surface_pressure |
---|
[1496] | 1215 | |
---|
| 1216 | ! |
---|
| 1217 | !-- Add LW up so that it can be removed in prognostic equation |
---|
[1691] | 1218 | rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) |
---|
[1496] | 1219 | |
---|
[1788] | 1220 | ! |
---|
| 1221 | !-- Calculate new skin temperature |
---|
[1500] | 1222 | IF ( humidity ) THEN |
---|
[1691] | 1223 | #if defined ( __rrtmg ) |
---|
[1496] | 1224 | ! |
---|
[1500] | 1225 | !-- Numerator of the prognostic equation |
---|
[1709] | 1226 | coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & |
---|
[1691] | 1227 | * t_surface(j,i) - rad_lw_out(nzb,j,i) & |
---|
[1709] | 1228 | + f_shf * pt1 + f_qsws * ( qv1 - q_s & |
---|
[1571] | 1229 | + dq_s_dt * t_surface(j,i) ) + lambda_surface & |
---|
[1551] | 1230 | * t_soil(nzb_soil,j,i) |
---|
[1496] | 1231 | |
---|
| 1232 | ! |
---|
[1500] | 1233 | !-- Denominator of the prognostic equation |
---|
[1709] | 1234 | coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt & |
---|
[1691] | 1235 | + lambda_surface + f_shf / exn |
---|
| 1236 | #else |
---|
[1496] | 1237 | |
---|
[1691] | 1238 | ! |
---|
| 1239 | !-- Numerator of the prognostic equation |
---|
| 1240 | coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & |
---|
| 1241 | * t_surface(j,i) ** 4 & |
---|
[1709] | 1242 | + f_shf * pt1 + f_qsws * ( qv1 & |
---|
[1691] | 1243 | - q_s + dq_s_dt * t_surface(j,i) ) & |
---|
| 1244 | + lambda_surface * t_soil(nzb_soil,j,i) |
---|
| 1245 | |
---|
| 1246 | ! |
---|
[1709] | 1247 | !-- Denominator of the prognostic equation |
---|
| 1248 | coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws & |
---|
| 1249 | * dq_s_dt + lambda_surface + f_shf / exn |
---|
[1691] | 1250 | |
---|
| 1251 | #endif |
---|
[1500] | 1252 | ELSE |
---|
| 1253 | |
---|
[1691] | 1254 | #if defined ( __rrtmg ) |
---|
[1500] | 1255 | ! |
---|
| 1256 | !-- Numerator of the prognostic equation |
---|
[1709] | 1257 | coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & |
---|
[1691] | 1258 | * t_surface(j,i) - rad_lw_out(nzb,j,i) & |
---|
[1709] | 1259 | + f_shf * pt1 + lambda_surface & |
---|
[1551] | 1260 | * t_soil(nzb_soil,j,i) |
---|
[1500] | 1261 | |
---|
| 1262 | ! |
---|
| 1263 | !-- Denominator of the prognostic equation |
---|
[1709] | 1264 | coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn |
---|
[1691] | 1265 | #else |
---|
| 1266 | |
---|
| 1267 | ! |
---|
| 1268 | !-- Numerator of the prognostic equation |
---|
| 1269 | coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & |
---|
[1709] | 1270 | * t_surface(j,i) ** 4 + f_shf * pt1 & |
---|
[1691] | 1271 | + lambda_surface * t_soil(nzb_soil,j,i) |
---|
| 1272 | |
---|
| 1273 | ! |
---|
| 1274 | !-- Denominator of the prognostic equation |
---|
[1571] | 1275 | coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 & |
---|
[1551] | 1276 | + lambda_surface + f_shf / exn |
---|
[1691] | 1277 | #endif |
---|
[1500] | 1278 | ENDIF |
---|
| 1279 | |
---|
[1496] | 1280 | tend = 0.0_wp |
---|
| 1281 | |
---|
| 1282 | ! |
---|
[1551] | 1283 | !-- Implicit solution when the surface layer has no heat capacity, |
---|
[1496] | 1284 | !-- otherwise use RK3 scheme. |
---|
[1788] | 1285 | t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * & |
---|
| 1286 | t_surface(j,i) ) / ( c_surface_tmp + coef_2 & |
---|
| 1287 | * dt_3d * tsc(2) ) |
---|
| 1288 | |
---|
[1496] | 1289 | ! |
---|
| 1290 | !-- Add RK3 term |
---|
[1788] | 1291 | IF ( c_surface_tmp /= 0.0_wp ) THEN |
---|
[1496] | 1292 | |
---|
[1691] | 1293 | t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & |
---|
| 1294 | * tt_surface_m(j,i) |
---|
| 1295 | |
---|
[1496] | 1296 | ! |
---|
[1691] | 1297 | !-- Calculate true tendency |
---|
| 1298 | tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) & |
---|
| 1299 | * tt_surface_m(j,i)) / (dt_3d * tsc(2)) |
---|
[1496] | 1300 | ! |
---|
[1691] | 1301 | !-- Calculate t_surface tendencies for the next Runge-Kutta step |
---|
| 1302 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1303 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 1304 | tt_surface_m(j,i) = tend |
---|
| 1305 | ELSEIF ( intermediate_timestep_count < & |
---|
| 1306 | intermediate_timestep_count_max ) THEN |
---|
| 1307 | tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp & |
---|
| 1308 | * tt_surface_m(j,i) |
---|
| 1309 | ENDIF |
---|
[1496] | 1310 | ENDIF |
---|
| 1311 | ENDIF |
---|
| 1312 | |
---|
| 1313 | ! |
---|
[1757] | 1314 | !-- In case of fast changes in the skin temperature, it is possible to |
---|
| 1315 | !-- update the radiative fluxes independently from the prescribed |
---|
| 1316 | !-- radiation call frequency. This effectively prevents oscillations, |
---|
| 1317 | !-- especially when setting skip_time_do_radiation /= 0. The threshold |
---|
| 1318 | !-- value of 0.2 used here is just a first guess. This method should be |
---|
| 1319 | !-- revised in the future as tests have shown that the threshold is |
---|
| 1320 | !-- often reached, when no oscillations would occur (causes immense |
---|
| 1321 | !-- computing time for the radiation code). |
---|
[1788] | 1322 | IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND. & |
---|
[1757] | 1323 | unscheduled_radiation_calls ) THEN |
---|
[1691] | 1324 | force_radiation_call_l = .TRUE. |
---|
| 1325 | ENDIF |
---|
| 1326 | |
---|
| 1327 | pt(k,j,i) = t_surface_p(j,i) / exn |
---|
| 1328 | |
---|
| 1329 | ! |
---|
[1496] | 1330 | !-- Calculate fluxes |
---|
[1691] | 1331 | #if defined ( __rrtmg ) |
---|
[1709] | 1332 | rad_net_l(j,i) = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & |
---|
[1691] | 1333 | * t_surface(j,i) - rad_lw_out(nzb,j,i) & |
---|
[1709] | 1334 | - rad_lw_out_change_0(j,i) * t_surface_p(j,i) |
---|
[1691] | 1335 | |
---|
| 1336 | IF ( rrtm_idrv == 1 ) THEN |
---|
| 1337 | rad_net(j,i) = rad_net_l(j,i) |
---|
| 1338 | rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) & |
---|
[1709] | 1339 | + rad_lw_out_change_0(j,i) & |
---|
[1691] | 1340 | * ( t_surface_p(j,i) - t_surface(j,i) ) |
---|
| 1341 | ENDIF |
---|
| 1342 | #else |
---|
| 1343 | rad_net_l(j,i) = rad_net_l(j,i) + 3.0_wp * sigma_sb & |
---|
| 1344 | * t_surface(j,i)**4 - 4.0_wp * sigma_sb & |
---|
| 1345 | * t_surface(j,i)**3 * t_surface_p(j,i) |
---|
| 1346 | #endif |
---|
| 1347 | |
---|
[1551] | 1348 | ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & |
---|
| 1349 | - t_soil(nzb_soil,j,i)) |
---|
[1496] | 1350 | |
---|
[1709] | 1351 | shf_eb(j,i) = - f_shf * ( pt1 - pt(k,j,i) ) |
---|
[1691] | 1352 | |
---|
| 1353 | shf(j,i) = shf_eb(j,i) / rho_cp |
---|
| 1354 | |
---|
[1500] | 1355 | IF ( humidity ) THEN |
---|
[1709] | 1356 | qsws_eb(j,i) = - f_qsws * ( qv1 - q_s + dq_s_dt & |
---|
[1571] | 1357 | * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) |
---|
[1496] | 1358 | |
---|
[1691] | 1359 | qsws(j,i) = qsws_eb(j,i) / rho_lv |
---|
| 1360 | |
---|
[1709] | 1361 | qsws_veg_eb(j,i) = - f_qsws_veg * ( qv1 - q_s & |
---|
[1571] | 1362 | + dq_s_dt * t_surface(j,i) - dq_s_dt & |
---|
[1551] | 1363 | * t_surface_p(j,i) ) |
---|
| 1364 | |
---|
[1709] | 1365 | qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s & |
---|
[1571] | 1366 | + dq_s_dt * t_surface(j,i) - dq_s_dt & |
---|
[1551] | 1367 | * t_surface_p(j,i) ) |
---|
| 1368 | |
---|
[1709] | 1369 | qsws_liq_eb(j,i) = - f_qsws_liq * ( qv1 - q_s & |
---|
[1571] | 1370 | + dq_s_dt * t_surface(j,i) - dq_s_dt & |
---|
[1551] | 1371 | * t_surface_p(j,i) ) |
---|
[1500] | 1372 | ENDIF |
---|
[1788] | 1373 | |
---|
[1500] | 1374 | ! |
---|
| 1375 | !-- Calculate the true surface resistance |
---|
[1691] | 1376 | IF ( qsws_eb(j,i) == 0.0_wp ) THEN |
---|
[1551] | 1377 | r_s(j,i) = 1.0E10_wp |
---|
[1496] | 1378 | ELSE |
---|
[1788] | 1379 | r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt & |
---|
[1571] | 1380 | * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & |
---|
[1551] | 1381 | / qsws_eb(j,i) - r_a(j,i) |
---|
[1496] | 1382 | ENDIF |
---|
| 1383 | |
---|
| 1384 | ! |
---|
| 1385 | !-- Calculate change in liquid water reservoir due to dew fall or |
---|
[1500] | 1386 | !-- evaporation of liquid water |
---|
| 1387 | IF ( humidity ) THEN |
---|
[1496] | 1388 | ! |
---|
[1571] | 1389 | !-- If precipitation is activated, add rain water to qsws_liq_eb |
---|
| 1390 | !-- and qsws_soil_eb according the the vegetation coverage. |
---|
[1500] | 1391 | !-- precipitation_rate is given in mm. |
---|
| 1392 | IF ( precipitation ) THEN |
---|
[1571] | 1393 | |
---|
| 1394 | ! |
---|
| 1395 | !-- Add precipitation to liquid water reservoir, if possible. |
---|
[1788] | 1396 | !-- Otherwise, add the water to soil. In case of |
---|
| 1397 | !-- pavements, the exceeding water amount is implicitely removed |
---|
| 1398 | !-- as runoff as qsws_soil_eb is then not used in the soil model |
---|
[1691] | 1399 | IF ( m_liq_eb(j,i) /= m_liq_eb_max ) THEN |
---|
[1571] | 1400 | qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & |
---|
[1691] | 1401 | + c_veg(j,i) * prr(k,j,i) * hyrho(k) & |
---|
[1571] | 1402 | * 0.001_wp * rho_l * l_v |
---|
| 1403 | ELSE |
---|
| 1404 | qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & |
---|
[1691] | 1405 | + c_veg(j,i) * prr(k,j,i) * hyrho(k) & |
---|
[1571] | 1406 | * 0.001_wp * rho_l * l_v |
---|
| 1407 | ENDIF |
---|
| 1408 | |
---|
| 1409 | !-- Add precipitation to bare soil according to the bare soil |
---|
| 1410 | !-- coverage. |
---|
| 1411 | qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp & |
---|
[1691] | 1412 | - c_veg(j,i)) * prr(k,j,i) * hyrho(k) & |
---|
[1571] | 1413 | * 0.001_wp * rho_l * l_v |
---|
[1496] | 1414 | ENDIF |
---|
[1691] | 1415 | |
---|
[1500] | 1416 | ! |
---|
| 1417 | !-- If the air is saturated, check the reservoir water level |
---|
[1691] | 1418 | IF ( qsws_eb(j,i) < 0.0_wp ) THEN |
---|
| 1419 | |
---|
[1500] | 1420 | ! |
---|
[1551] | 1421 | !-- Check if reservoir is full (avoid values > m_liq_eb_max) |
---|
| 1422 | !-- In that case, qsws_liq_eb goes to qsws_soil_eb. In this |
---|
| 1423 | !-- case qsws_veg_eb is zero anyway (because c_liq = 1), |
---|
| 1424 | !-- so that tend is zero and no further check is needed |
---|
[1691] | 1425 | IF ( m_liq_eb(j,i) == m_liq_eb_max ) THEN |
---|
[1551] | 1426 | qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & |
---|
| 1427 | + qsws_liq_eb(j,i) |
---|
| 1428 | qsws_liq_eb(j,i) = 0.0_wp |
---|
[1500] | 1429 | ENDIF |
---|
[1496] | 1430 | |
---|
| 1431 | ! |
---|
[1551] | 1432 | !-- In case qsws_veg_eb becomes negative (unphysical behavior), |
---|
| 1433 | !-- let the water enter the liquid water reservoir as dew on the |
---|
[1500] | 1434 | !-- plant |
---|
[1691] | 1435 | IF ( qsws_veg_eb(j,i) < 0.0_wp ) THEN |
---|
[1551] | 1436 | qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i) |
---|
| 1437 | qsws_veg_eb(j,i) = 0.0_wp |
---|
[1500] | 1438 | ENDIF |
---|
| 1439 | ENDIF |
---|
[1496] | 1440 | |
---|
[1551] | 1441 | tend = - qsws_liq_eb(j,i) * drho_l_lv |
---|
[1496] | 1442 | |
---|
[1551] | 1443 | m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend & |
---|
| 1444 | + tsc(3) * tm_liq_eb_m(j,i) ) |
---|
[1496] | 1445 | |
---|
| 1446 | ! |
---|
[1500] | 1447 | !-- Check if reservoir is overfull -> reduce to maximum |
---|
| 1448 | !-- (conservation of water is violated here) |
---|
[1551] | 1449 | m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max) |
---|
[1496] | 1450 | |
---|
| 1451 | ! |
---|
[1500] | 1452 | !-- Check if reservoir is empty (avoid values < 0.0) |
---|
| 1453 | !-- (conservation of water is violated here) |
---|
[1551] | 1454 | m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp) |
---|
[1496] | 1455 | |
---|
| 1456 | |
---|
| 1457 | ! |
---|
[1551] | 1458 | !-- Calculate m_liq_eb tendencies for the next Runge-Kutta step |
---|
[1500] | 1459 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1460 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
[1551] | 1461 | tm_liq_eb_m(j,i) = tend |
---|
[1500] | 1462 | ELSEIF ( intermediate_timestep_count < & |
---|
| 1463 | intermediate_timestep_count_max ) THEN |
---|
[1551] | 1464 | tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp & |
---|
| 1465 | * tm_liq_eb_m(j,i) |
---|
[1500] | 1466 | ENDIF |
---|
[1496] | 1467 | ENDIF |
---|
| 1468 | |
---|
[1500] | 1469 | ENDIF |
---|
| 1470 | |
---|
[1496] | 1471 | ENDDO |
---|
[1500] | 1472 | ENDDO |
---|
[1496] | 1473 | |
---|
[1691] | 1474 | ! |
---|
| 1475 | !-- Make a logical OR for all processes. Force radiation call if at |
---|
[1757] | 1476 | !-- least one processor reached the threshold change in skin temperature |
---|
[1788] | 1477 | IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count & |
---|
[1757] | 1478 | == intermediate_timestep_count_max-1 ) THEN |
---|
[1697] | 1479 | #if defined( __parallel ) |
---|
[1691] | 1480 | IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) |
---|
| 1481 | CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, & |
---|
| 1482 | 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr ) |
---|
| 1483 | #else |
---|
| 1484 | force_radiation_call = force_radiation_call_l |
---|
| 1485 | #endif |
---|
| 1486 | force_radiation_call_l = .FALSE. |
---|
| 1487 | ENDIF |
---|
| 1488 | |
---|
[1788] | 1489 | ! |
---|
| 1490 | !-- Calculate surface specific humidity |
---|
| 1491 | IF ( humidity ) THEN |
---|
| 1492 | CALL calc_q_surface |
---|
| 1493 | ENDIF |
---|
[1691] | 1494 | |
---|
[1788] | 1495 | ! |
---|
| 1496 | !-- Calculate new roughness lengths (for water surfaces only) |
---|
| 1497 | CALL calc_z0_water_surface |
---|
| 1498 | |
---|
| 1499 | |
---|
[1496] | 1500 | END SUBROUTINE lsm_energy_balance |
---|
| 1501 | |
---|
| 1502 | |
---|
| 1503 | !------------------------------------------------------------------------------! |
---|
| 1504 | ! Description: |
---|
| 1505 | ! ------------ |
---|
[1682] | 1506 | !> Soil model as part of the land surface model. The model predicts soil |
---|
| 1507 | !> temperature and water content. |
---|
[1496] | 1508 | !------------------------------------------------------------------------------! |
---|
| 1509 | SUBROUTINE lsm_soil_model |
---|
| 1510 | |
---|
| 1511 | |
---|
| 1512 | IMPLICIT NONE |
---|
| 1513 | |
---|
[1682] | 1514 | INTEGER(iwp) :: i !< running index |
---|
| 1515 | INTEGER(iwp) :: j !< running index |
---|
| 1516 | INTEGER(iwp) :: k !< running index |
---|
[1496] | 1517 | |
---|
[1682] | 1518 | REAL(wp) :: h_vg !< Van Genuchten coef. h |
---|
[1496] | 1519 | |
---|
[1682] | 1520 | REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp, & !< temp. gamma |
---|
[1691] | 1521 | lambda_temp, & !< temp. lambda |
---|
| 1522 | tend !< tendency |
---|
[1496] | 1523 | |
---|
[1691] | 1524 | DO i = nxlg, nxrg |
---|
| 1525 | DO j = nysg, nyng |
---|
[1788] | 1526 | |
---|
| 1527 | IF ( pave_surface(j,i) ) THEN |
---|
| 1528 | rho_c_total(nzb_soil,j,i) = pave_heat_capacity |
---|
| 1529 | lambda_temp(nzb_soil) = pave_heat_conductivity |
---|
| 1530 | ENDIF |
---|
| 1531 | |
---|
| 1532 | IF ( .NOT. water_surface(j,i) ) THEN |
---|
| 1533 | DO k = nzb_soil, nzt_soil |
---|
| 1534 | |
---|
| 1535 | |
---|
| 1536 | IF ( pave_surface(j,i) .AND. zs(k) <= pave_depth ) THEN |
---|
| 1537 | |
---|
| 1538 | rho_c_total(k,j,i) = pave_heat_capacity |
---|
| 1539 | lambda_temp(k) = pave_heat_conductivity |
---|
| 1540 | |
---|
| 1541 | ELSE |
---|
[1496] | 1542 | ! |
---|
[1788] | 1543 | !-- Calculate volumetric heat capacity of the soil, taking |
---|
| 1544 | !-- into account water content |
---|
| 1545 | rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) & |
---|
| 1546 | + rho_c_water * m_soil(k,j,i)) |
---|
[1496] | 1547 | |
---|
| 1548 | ! |
---|
[1788] | 1549 | !-- Calculate soil heat conductivity at the center of the soil |
---|
| 1550 | !-- layers |
---|
| 1551 | lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) * & |
---|
| 1552 | lambda_h_water ** m_soil(k,j,i) |
---|
[1691] | 1553 | |
---|
[1788] | 1554 | ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) |
---|
[1691] | 1555 | |
---|
[1788] | 1556 | lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + & |
---|
| 1557 | lambda_h_dry |
---|
| 1558 | ENDIF |
---|
[1496] | 1559 | |
---|
[1788] | 1560 | ENDDO |
---|
[1496] | 1561 | |
---|
| 1562 | ! |
---|
[1788] | 1563 | !-- Calculate soil heat conductivity (lambda_h) at the _stag level |
---|
| 1564 | !-- using linear interpolation. For pavement surface, the |
---|
| 1565 | !-- true pavement depth is considered |
---|
| 1566 | DO k = nzb_soil, nzt_soil-1 |
---|
| 1567 | IF ( pave_surface(j,i) .AND. zs(k) < pave_depth & |
---|
| 1568 | .AND. zs(k+1) > pave_depth ) THEN |
---|
| 1569 | lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1) & |
---|
| 1570 | * lambda_temp(k) & |
---|
| 1571 | + ( 1.0_wp - ( pave_depth - zs(k) ) & |
---|
| 1572 | / dz_soil(k+1) ) * lambda_temp(k+1) |
---|
| 1573 | ELSE |
---|
| 1574 | lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & |
---|
| 1575 | * 0.5_wp |
---|
| 1576 | ENDIF |
---|
| 1577 | ENDDO |
---|
| 1578 | lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil) |
---|
[1496] | 1579 | |
---|
[1788] | 1580 | |
---|
| 1581 | |
---|
| 1582 | |
---|
[1496] | 1583 | ! |
---|
[1788] | 1584 | !-- Prognostic equation for soil temperature t_soil |
---|
| 1585 | tend(:) = 0.0_wp |
---|
[1496] | 1586 | |
---|
[1788] | 1587 | tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) * & |
---|
| 1588 | ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) & |
---|
| 1589 | - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1) & |
---|
| 1590 | + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil) |
---|
[1691] | 1591 | |
---|
[1788] | 1592 | DO k = nzb_soil+1, nzt_soil |
---|
| 1593 | tend(k) = (1.0_wp/rho_c_total(k,j,i)) & |
---|
| 1594 | * ( lambda_h(k,j,i) & |
---|
| 1595 | * ( t_soil(k+1,j,i) - t_soil(k,j,i) ) & |
---|
| 1596 | * ddz_soil(k+1) & |
---|
| 1597 | - lambda_h(k-1,j,i) & |
---|
| 1598 | * ( t_soil(k,j,i) - t_soil(k-1,j,i) ) & |
---|
| 1599 | * ddz_soil(k) & |
---|
| 1600 | ) * ddz_soil_stag(k) |
---|
[1691] | 1601 | |
---|
[1788] | 1602 | ENDDO |
---|
[1691] | 1603 | |
---|
[1788] | 1604 | t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)& |
---|
| 1605 | + dt_3d * ( tsc(2) & |
---|
| 1606 | * tend(nzb_soil:nzt_soil) & |
---|
| 1607 | + tsc(3) & |
---|
| 1608 | * tt_soil_m(:,j,i) ) |
---|
[1496] | 1609 | |
---|
| 1610 | ! |
---|
[1788] | 1611 | !-- Calculate t_soil tendencies for the next Runge-Kutta step |
---|
| 1612 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1613 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 1614 | DO k = nzb_soil, nzt_soil |
---|
| 1615 | tt_soil_m(k,j,i) = tend(k) |
---|
| 1616 | ENDDO |
---|
| 1617 | ELSEIF ( intermediate_timestep_count < & |
---|
| 1618 | intermediate_timestep_count_max ) THEN |
---|
| 1619 | DO k = nzb_soil, nzt_soil |
---|
| 1620 | tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & |
---|
[1551] | 1621 | * tt_soil_m(k,j,i) |
---|
[1788] | 1622 | ENDDO |
---|
| 1623 | ENDIF |
---|
[1496] | 1624 | ENDIF |
---|
| 1625 | |
---|
| 1626 | |
---|
[1788] | 1627 | DO k = nzb_soil, nzt_soil |
---|
[1551] | 1628 | |
---|
[1496] | 1629 | ! |
---|
[1788] | 1630 | !-- Calculate soil diffusivity at the center of the soil layers |
---|
| 1631 | lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat & |
---|
| 1632 | / m_sat(j,i) ) * ( MAX( m_soil(k,j,i), & |
---|
| 1633 | m_wilt(j,i) ) / m_sat(j,i) )**( & |
---|
| 1634 | b_ch + 2.0_wp ) |
---|
[1496] | 1635 | |
---|
| 1636 | ! |
---|
[1788] | 1637 | !-- Parametrization of Van Genuchten |
---|
| 1638 | IF ( soil_type /= 7 ) THEN |
---|
[1551] | 1639 | ! |
---|
[1788] | 1640 | !-- Calculate the hydraulic conductivity after Van Genuchten |
---|
| 1641 | !-- (1980) |
---|
| 1642 | h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - & |
---|
| 1643 | MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**( & |
---|
| 1644 | n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp & |
---|
| 1645 | )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i) |
---|
[1496] | 1646 | |
---|
[1691] | 1647 | |
---|
[1788] | 1648 | gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & |
---|
| 1649 | ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**( & |
---|
| 1650 | 1.0_wp - 1.0_wp / n_vg(j,i) ) - ( & |
---|
| 1651 | alpha_vg(j,i) * h_vg )**( n_vg(j,i) & |
---|
| 1652 | - 1.0_wp) )**2 ) & |
---|
| 1653 | / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg & |
---|
| 1654 | )**n_vg(j,i) )**( ( 1.0_wp - 1.0_wp & |
---|
| 1655 | / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) ) |
---|
[1496] | 1656 | |
---|
[1551] | 1657 | ! |
---|
[1788] | 1658 | !-- Parametrization of Clapp & Hornberger |
---|
| 1659 | ELSE |
---|
| 1660 | gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i) & |
---|
| 1661 | / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp) |
---|
| 1662 | ENDIF |
---|
[1551] | 1663 | |
---|
[1788] | 1664 | ENDDO |
---|
[1496] | 1665 | |
---|
| 1666 | ! |
---|
[1788] | 1667 | !-- Prognostic equation for soil moisture content. Only performed, |
---|
| 1668 | !-- when humidity is enabled in the atmosphere and the surface type |
---|
| 1669 | !-- is not pavement (implies dry soil below). |
---|
| 1670 | IF ( humidity .AND. .NOT. pave_surface(j,i) ) THEN |
---|
| 1671 | ! |
---|
| 1672 | !-- Calculate soil diffusivity (lambda_w) at the _stag level |
---|
| 1673 | !-- using linear interpolation. To do: replace this with |
---|
| 1674 | !-- ECMWF-IFS Eq. 8.81 |
---|
| 1675 | DO k = nzb_soil, nzt_soil-1 |
---|
[1496] | 1676 | |
---|
[1788] | 1677 | lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & |
---|
| 1678 | * 0.5_wp |
---|
| 1679 | gamma_w(k,j,i) = ( gamma_temp(k+1) + gamma_temp(k) ) & |
---|
| 1680 | * 0.5_wp |
---|
[1496] | 1681 | |
---|
[1788] | 1682 | ENDDO |
---|
[1496] | 1683 | |
---|
| 1684 | ! |
---|
| 1685 | ! |
---|
[1788] | 1686 | !-- In case of a closed bottom (= water content is conserved), set |
---|
| 1687 | !-- hydraulic conductivity to zero to that no water will be lost |
---|
| 1688 | !-- in the bottom layer. |
---|
| 1689 | IF ( conserve_water_content ) THEN |
---|
| 1690 | gamma_w(nzt_soil,j,i) = 0.0_wp |
---|
| 1691 | ELSE |
---|
| 1692 | gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil) |
---|
| 1693 | ENDIF |
---|
[1496] | 1694 | |
---|
[1788] | 1695 | !-- The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v)) |
---|
| 1696 | !-- ensures the mass conservation for water. The transpiration of |
---|
| 1697 | !-- plants equals the cumulative withdrawals by the roots in the |
---|
| 1698 | !-- soil. The scheme takes into account the availability of water |
---|
| 1699 | !-- in the soil layers as well as the root fraction in the |
---|
| 1700 | !-- respective layer. Layer with moisture below wilting point will |
---|
| 1701 | !-- not contribute, which reflects the preference of plants to |
---|
| 1702 | !-- take water from moister layers. |
---|
[1496] | 1703 | |
---|
| 1704 | ! |
---|
[1788] | 1705 | !-- Calculate the root extraction (ECMWF 7.69, the sum of root_extr |
---|
| 1706 | !-- = 1). The energy balance solver guarantees a positive |
---|
| 1707 | !-- transpiration, so that there is no need for an additional check. |
---|
[1691] | 1708 | DO k = nzb_soil, nzt_soil |
---|
[1788] | 1709 | IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN |
---|
| 1710 | m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) |
---|
| 1711 | ENDIF |
---|
| 1712 | ENDDO |
---|
[1496] | 1713 | |
---|
[1788] | 1714 | IF ( m_total > 0.0_wp ) THEN |
---|
| 1715 | DO k = nzb_soil, nzt_soil |
---|
| 1716 | IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN |
---|
| 1717 | root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total |
---|
| 1718 | ELSE |
---|
| 1719 | root_extr(k) = 0.0_wp |
---|
| 1720 | ENDIF |
---|
| 1721 | ENDDO |
---|
| 1722 | ENDIF |
---|
| 1723 | |
---|
[1496] | 1724 | ! |
---|
[1788] | 1725 | !-- Prognostic equation for soil water content m_soil. |
---|
| 1726 | tend(:) = 0.0_wp |
---|
| 1727 | |
---|
| 1728 | tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * ( & |
---|
[1551] | 1729 | m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) ) & |
---|
[1691] | 1730 | * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( & |
---|
[1788] | 1731 | root_extr(nzb_soil) * qsws_veg_eb(j,i) & |
---|
| 1732 | + qsws_soil_eb(j,i) ) * drho_l_lv ) & |
---|
| 1733 | * ddz_soil_stag(nzb_soil) |
---|
[1496] | 1734 | |
---|
[1788] | 1735 | DO k = nzb_soil+1, nzt_soil-1 |
---|
| 1736 | tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i) & |
---|
| 1737 | - m_soil(k,j,i) ) * ddz_soil(k+1) & |
---|
| 1738 | - gamma_w(k,j,i) & |
---|
| 1739 | - lambda_w(k-1,j,i) * (m_soil(k,j,i) - & |
---|
| 1740 | m_soil(k-1,j,i)) * ddz_soil(k) & |
---|
| 1741 | + gamma_w(k-1,j,i) - (root_extr(k) & |
---|
| 1742 | * qsws_veg_eb(j,i) * drho_l_lv) & |
---|
| 1743 | ) * ddz_soil_stag(k) |
---|
[1496] | 1744 | |
---|
[1788] | 1745 | ENDDO |
---|
| 1746 | tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i) & |
---|
| 1747 | - lambda_w(nzt_soil-1,j,i) & |
---|
| 1748 | * (m_soil(nzt_soil,j,i) & |
---|
| 1749 | - m_soil(nzt_soil-1,j,i)) & |
---|
| 1750 | * ddz_soil(nzt_soil) & |
---|
| 1751 | + gamma_w(nzt_soil-1,j,i) - ( & |
---|
| 1752 | root_extr(nzt_soil) & |
---|
| 1753 | * qsws_veg_eb(j,i) * drho_l_lv ) & |
---|
| 1754 | ) * ddz_soil_stag(nzt_soil) |
---|
[1496] | 1755 | |
---|
[1788] | 1756 | m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)& |
---|
| 1757 | + dt_3d * ( tsc(2) * tend(:) & |
---|
| 1758 | + tsc(3) * tm_soil_m(:,j,i) ) |
---|
| 1759 | |
---|
[1496] | 1760 | ! |
---|
[1788] | 1761 | !-- Account for dry soils (find a better solution here!) |
---|
| 1762 | DO k = nzb_soil, nzt_soil |
---|
| 1763 | IF ( m_soil_p(k,j,i) < 0.0_wp ) m_soil_p(k,j,i) = 0.0_wp |
---|
| 1764 | ENDDO |
---|
[1496] | 1765 | |
---|
| 1766 | ! |
---|
[1788] | 1767 | !-- Calculate m_soil tendencies for the next Runge-Kutta step |
---|
| 1768 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1769 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 1770 | DO k = nzb_soil, nzt_soil |
---|
| 1771 | tm_soil_m(k,j,i) = tend(k) |
---|
| 1772 | ENDDO |
---|
| 1773 | ELSEIF ( intermediate_timestep_count < & |
---|
| 1774 | intermediate_timestep_count_max ) THEN |
---|
| 1775 | DO k = nzb_soil, nzt_soil |
---|
| 1776 | tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp& |
---|
[1496] | 1777 | * tm_soil_m(k,j,i) |
---|
[1788] | 1778 | ENDDO |
---|
| 1779 | ENDIF |
---|
[1496] | 1780 | ENDIF |
---|
[1788] | 1781 | |
---|
[1496] | 1782 | ENDIF |
---|
| 1783 | |
---|
| 1784 | ENDIF |
---|
| 1785 | |
---|
| 1786 | ENDDO |
---|
| 1787 | ENDDO |
---|
| 1788 | |
---|
[1788] | 1789 | END SUBROUTINE lsm_soil_model |
---|
| 1790 | |
---|
| 1791 | |
---|
| 1792 | !------------------------------------------------------------------------------! |
---|
| 1793 | ! Description: |
---|
| 1794 | ! ------------ |
---|
| 1795 | !> Calculation of roughness length for open water (lakes, ocean). The |
---|
| 1796 | !> parameterization follows Charnock (1955). Two different implementations |
---|
| 1797 | !> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al. |
---|
| 1798 | !> 2012) |
---|
| 1799 | !------------------------------------------------------------------------------! |
---|
| 1800 | SUBROUTINE calc_z0_water_surface |
---|
| 1801 | |
---|
| 1802 | USE control_parameters, & |
---|
| 1803 | ONLY: g, kappa, molecular_viscosity |
---|
| 1804 | |
---|
| 1805 | IMPLICIT NONE |
---|
| 1806 | |
---|
| 1807 | INTEGER :: i !< running index |
---|
| 1808 | INTEGER :: j !< running index |
---|
| 1809 | |
---|
| 1810 | REAL(wp), PARAMETER :: alpha_ch = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF |
---|
| 1811 | ! REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number) |
---|
| 1812 | ! REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization |
---|
| 1813 | ! REAL(wp) :: re_0 !< near-surface roughness Reynolds number |
---|
| 1814 | |
---|
| 1815 | |
---|
| 1816 | DO i = nxlg, nxrg |
---|
| 1817 | DO j = nysg, nyng |
---|
| 1818 | IF ( water_surface(j,i) ) THEN |
---|
| 1819 | |
---|
[1496] | 1820 | ! |
---|
[1788] | 1821 | !-- Disabled: FLake parameterization. Ideally, the Charnock |
---|
| 1822 | !-- coefficient should depend on the water depth and the fetch |
---|
| 1823 | !-- length |
---|
| 1824 | ! re_0 = z0(j,i) * us(j,i) / molecular_viscosity |
---|
| 1825 | ! |
---|
| 1826 | ! z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i), & |
---|
| 1827 | ! alpha_ch * us(j,i) / g ) |
---|
| 1828 | ! |
---|
| 1829 | ! z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) ) |
---|
| 1830 | ! z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) ) |
---|
[1496] | 1831 | |
---|
[1788] | 1832 | ! |
---|
| 1833 | !-- Set minimum roughness length for u* > 0.2 |
---|
| 1834 | ! IF ( us(j,i) > 0.2_wp ) THEN |
---|
| 1835 | ! z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) ) |
---|
| 1836 | ! z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) ) |
---|
| 1837 | ! ENDIF |
---|
[1496] | 1838 | |
---|
[1788] | 1839 | ! |
---|
| 1840 | !-- ECMWF IFS model parameterization after Beljaars (1994). At low |
---|
| 1841 | !-- wind speed, the sea surface becomes aerodynamically smooth and |
---|
| 1842 | !-- the roughness scales with the viscosity. At high wind speed, the |
---|
| 1843 | !-- Charnock relation is used. |
---|
| 1844 | z0(j,i) = ( 0.11_wp * molecular_viscosity / us(j,i) ) & |
---|
| 1845 | + ( alpha_ch * us(j,i)**2 / g ) |
---|
[1496] | 1846 | |
---|
[1788] | 1847 | z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i) |
---|
| 1848 | z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i) |
---|
[1496] | 1849 | |
---|
[1788] | 1850 | ENDIF |
---|
| 1851 | ENDDO |
---|
| 1852 | ENDDO |
---|
| 1853 | |
---|
| 1854 | END SUBROUTINE calc_z0_water_surface |
---|
| 1855 | |
---|
| 1856 | |
---|
[1496] | 1857 | !------------------------------------------------------------------------------! |
---|
| 1858 | ! Description: |
---|
| 1859 | ! ------------ |
---|
[1788] | 1860 | !> Calculation of specific humidity of the skin layer (surface). It is assumend |
---|
| 1861 | !> that the skin is always saturated. |
---|
[1496] | 1862 | !------------------------------------------------------------------------------! |
---|
[1551] | 1863 | SUBROUTINE calc_q_surface |
---|
[1496] | 1864 | |
---|
| 1865 | IMPLICIT NONE |
---|
| 1866 | |
---|
[1682] | 1867 | INTEGER :: i !< running index |
---|
| 1868 | INTEGER :: j !< running index |
---|
| 1869 | INTEGER :: k !< running index |
---|
[1788] | 1870 | |
---|
[1682] | 1871 | REAL(wp) :: resistance !< aerodynamic and soil resistance term |
---|
[1496] | 1872 | |
---|
[1691] | 1873 | DO i = nxlg, nxrg |
---|
| 1874 | DO j = nysg, nyng |
---|
[1496] | 1875 | k = nzb_s_inner(j,i) |
---|
| 1876 | |
---|
| 1877 | ! |
---|
| 1878 | !-- Calculate water vapour pressure at saturation |
---|
[1691] | 1879 | e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i) & |
---|
| 1880 | - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) ) |
---|
[1496] | 1881 | |
---|
| 1882 | ! |
---|
| 1883 | !-- Calculate specific humidity at saturation |
---|
| 1884 | q_s = 0.622_wp * e_s / surface_pressure |
---|
| 1885 | |
---|
| 1886 | resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i)) |
---|
| 1887 | |
---|
| 1888 | ! |
---|
| 1889 | !-- Calculate specific humidity at surface |
---|
[1691] | 1890 | IF ( cloud_physics ) THEN |
---|
| 1891 | q(k,j,i) = resistance * q_s + (1.0_wp - resistance) & |
---|
| 1892 | * ( q(k+1,j,i) - ql(k+1,j,i) ) |
---|
| 1893 | ELSE |
---|
| 1894 | q(k,j,i) = resistance * q_s + (1.0_wp - resistance) & |
---|
| 1895 | * q(k+1,j,i) |
---|
| 1896 | ENDIF |
---|
[1496] | 1897 | |
---|
[1691] | 1898 | ! |
---|
| 1899 | !-- Update virtual potential temperature |
---|
| 1900 | vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) ) |
---|
| 1901 | |
---|
[1496] | 1902 | ENDDO |
---|
| 1903 | ENDDO |
---|
| 1904 | |
---|
[1551] | 1905 | END SUBROUTINE calc_q_surface |
---|
[1496] | 1906 | |
---|
[1788] | 1907 | |
---|
[1551] | 1908 | !------------------------------------------------------------------------------! |
---|
| 1909 | ! Description: |
---|
| 1910 | ! ------------ |
---|
[1682] | 1911 | !> Swapping of timelevels |
---|
[1551] | 1912 | !------------------------------------------------------------------------------! |
---|
| 1913 | SUBROUTINE lsm_swap_timelevel ( mod_count ) |
---|
[1496] | 1914 | |
---|
[1551] | 1915 | IMPLICIT NONE |
---|
| 1916 | |
---|
| 1917 | INTEGER, INTENT(IN) :: mod_count |
---|
| 1918 | |
---|
| 1919 | #if defined( __nopointer ) |
---|
| 1920 | |
---|
| 1921 | t_surface = t_surface_p |
---|
| 1922 | t_soil = t_soil_p |
---|
| 1923 | IF ( humidity ) THEN |
---|
| 1924 | m_soil = m_soil_p |
---|
| 1925 | m_liq_eb = m_liq_eb_p |
---|
| 1926 | ENDIF |
---|
| 1927 | |
---|
| 1928 | #else |
---|
| 1929 | |
---|
| 1930 | SELECT CASE ( mod_count ) |
---|
| 1931 | |
---|
| 1932 | CASE ( 0 ) |
---|
| 1933 | |
---|
[1585] | 1934 | t_surface => t_surface_1; t_surface_p => t_surface_2 |
---|
| 1935 | t_soil => t_soil_1; t_soil_p => t_soil_2 |
---|
[1551] | 1936 | IF ( humidity ) THEN |
---|
[1585] | 1937 | m_soil => m_soil_1; m_soil_p => m_soil_2 |
---|
| 1938 | m_liq_eb => m_liq_eb_1; m_liq_eb_p => m_liq_eb_2 |
---|
[1551] | 1939 | ENDIF |
---|
| 1940 | |
---|
| 1941 | |
---|
| 1942 | CASE ( 1 ) |
---|
| 1943 | |
---|
[1585] | 1944 | t_surface => t_surface_2; t_surface_p => t_surface_1 |
---|
| 1945 | t_soil => t_soil_2; t_soil_p => t_soil_1 |
---|
[1551] | 1946 | IF ( humidity ) THEN |
---|
[1585] | 1947 | m_soil => m_soil_2; m_soil_p => m_soil_1 |
---|
| 1948 | m_liq_eb => m_liq_eb_2; m_liq_eb_p => m_liq_eb_1 |
---|
[1551] | 1949 | ENDIF |
---|
| 1950 | |
---|
| 1951 | END SELECT |
---|
| 1952 | #endif |
---|
| 1953 | |
---|
| 1954 | END SUBROUTINE lsm_swap_timelevel |
---|
| 1955 | |
---|
| 1956 | |
---|
[1496] | 1957 | END MODULE land_surface_model_mod |
---|