[3880] | 1 | ! radiation_two_stream.F90 - Compute two-stream coefficients |
---|
| 2 | ! |
---|
| 3 | ! (C) Copyright 2014- ECMWF. |
---|
| 4 | ! |
---|
| 5 | ! This software is licensed under the terms of the Apache Licence Version 2.0 |
---|
| 6 | ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
| 7 | ! |
---|
| 8 | ! In applying this licence, ECMWF does not waive the privileges and immunities |
---|
| 9 | ! granted to it by virtue of its status as an intergovernmental organisation |
---|
| 10 | ! nor does it submit to any jurisdiction. |
---|
| 11 | ! |
---|
| 12 | ! Author: Robin Hogan |
---|
| 13 | ! Email: r.j.hogan@ecmwf.int |
---|
| 14 | ! |
---|
| 15 | ! Modifications |
---|
| 16 | ! 2017-05-04 P Dueben/R Hogan Use JPRD where double precision essential |
---|
| 17 | ! 2017-07-12 R Hogan Optimized LW coeffs in low optical depth case |
---|
| 18 | ! 2017-07-26 R Hogan Added calc_frac_scattered_diffuse_sw routine |
---|
| 19 | ! 2017-10-23 R Hogan Renamed single-character variables |
---|
| 20 | |
---|
| 21 | module radiation_two_stream |
---|
| 22 | |
---|
| 23 | use parkind1, only : jprb, jprd |
---|
| 24 | |
---|
| 25 | implicit none |
---|
| 26 | public |
---|
| 27 | |
---|
| 28 | ! Elsasser's factor: the effective factor by which the zenith |
---|
| 29 | ! optical depth needs to be multiplied to account for longwave |
---|
| 30 | ! transmission at all angles through the atmosphere. Alternatively |
---|
| 31 | ! think of acos(1/lw_diffusivity) to be the effective zenith angle |
---|
| 32 | ! of longwave radiation. |
---|
| 33 | real(jprd), parameter :: LwDiffusivity = 1.66_jprd |
---|
| 34 | |
---|
| 35 | ! Shortwave diffusivity factor assumes hemispheric isotropy, assumed |
---|
| 36 | ! by Zdunkowski's scheme and most others; note that for efficiency |
---|
| 37 | ! this parameter is not used in the calculation of the gamma values, |
---|
| 38 | ! but is used in the SPARTACUS solver. |
---|
| 39 | real(jprb), parameter :: SwDiffusivity = 2.00_jprb |
---|
| 40 | |
---|
| 41 | ! The routines in this module can be called millions of times, so |
---|
| 42 | !calling Dr Hook for each one may be a significant overhead. |
---|
| 43 | !Uncomment the following to turn Dr Hook on. |
---|
| 44 | !#define DO_DR_HOOK_TWO_STREAM |
---|
| 45 | |
---|
| 46 | contains |
---|
| 47 | |
---|
| 48 | #ifdef FAST_EXPONENTIAL |
---|
| 49 | !--------------------------------------------------------------------- |
---|
| 50 | ! Fast exponential for negative arguments: a Pade approximant that |
---|
| 51 | ! doesn't go negative for negative arguments, applied to arg/8, and |
---|
| 52 | ! the result is then squared three times |
---|
| 53 | elemental function exp_fast(arg) result(ex) |
---|
| 54 | real(jprd) :: arg, ex |
---|
| 55 | ex = 1.0_jprd / (1.0_jprd + arg*(-0.125_jprd & |
---|
| 56 | + arg*(0.0078125_jprd - 0.000325520833333333_jprd * arg))) |
---|
| 57 | ex = ex*ex |
---|
| 58 | ex = ex*ex |
---|
| 59 | ex = ex*ex |
---|
| 60 | end function exp_fast |
---|
| 61 | #else |
---|
| 62 | #define exp_fast exp |
---|
| 63 | #endif |
---|
| 64 | |
---|
| 65 | !--------------------------------------------------------------------- |
---|
| 66 | ! Calculate the two-stream coefficients gamma1 and gamma2 for the |
---|
| 67 | ! longwave |
---|
| 68 | subroutine calc_two_stream_gammas_lw(ng, ssa, g, & |
---|
| 69 | & gamma1, gamma2) |
---|
| 70 | |
---|
| 71 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 72 | use yomhook, only : lhook, dr_hook |
---|
| 73 | #endif |
---|
| 74 | |
---|
| 75 | integer, intent(in) :: ng |
---|
| 76 | ! Sngle scattering albedo and asymmetry factor: |
---|
| 77 | real(jprb), intent(in), dimension(:) :: ssa, g |
---|
| 78 | real(jprb), intent(out), dimension(:) :: gamma1, gamma2 |
---|
| 79 | |
---|
| 80 | real(jprb) :: factor |
---|
| 81 | |
---|
| 82 | integer :: jg |
---|
| 83 | |
---|
| 84 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 85 | real(jprb) :: hook_handle |
---|
| 86 | |
---|
| 87 | if (lhook) call dr_hook('radiation_two_stream:calc_two_stream_gammas_lw',0,hook_handle) |
---|
| 88 | #endif |
---|
| 89 | |
---|
| 90 | do jg = 1, ng |
---|
| 91 | ! Fu et al. (1997), Eq 2.9 and 2.10: |
---|
| 92 | ! gamma1(jg) = LwDiffusivity * (1.0_jprb - 0.5_jprb*ssa(jg) & |
---|
| 93 | ! & * (1.0_jprb + g(jg))) |
---|
| 94 | ! gamma2(jg) = LwDiffusivity * 0.5_jprb * ssa(jg) & |
---|
| 95 | ! & * (1.0_jprb - g(jg)) |
---|
| 96 | ! Reduce number of multiplications |
---|
| 97 | factor = (LwDiffusivity * 0.5_jprb) * ssa(jg) |
---|
| 98 | gamma1(jg) = LwDiffusivity - factor*(1.0_jprb + g(jg)) |
---|
| 99 | gamma2(jg) = factor * (1.0_jprb - g(jg)) |
---|
| 100 | end do |
---|
| 101 | |
---|
| 102 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 103 | if (lhook) call dr_hook('radiation_two_stream:calc_two_stream_gammas_lw',1,hook_handle) |
---|
| 104 | #endif |
---|
| 105 | |
---|
| 106 | end subroutine calc_two_stream_gammas_lw |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | !--------------------------------------------------------------------- |
---|
| 110 | ! Calculate the two-stream coefficients gamma1-gamma4 in the |
---|
| 111 | ! shortwave |
---|
| 112 | subroutine calc_two_stream_gammas_sw(ng, mu0, ssa, g, & |
---|
| 113 | & gamma1, gamma2, gamma3) |
---|
| 114 | |
---|
| 115 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 116 | use yomhook, only : lhook, dr_hook |
---|
| 117 | #endif |
---|
| 118 | |
---|
| 119 | integer, intent(in) :: ng |
---|
| 120 | ! Cosine of solar zenith angle, single scattering albedo and |
---|
| 121 | ! asymmetry factor: |
---|
| 122 | real(jprb), intent(in) :: mu0 |
---|
| 123 | real(jprb), intent(in), dimension(:) :: ssa, g |
---|
| 124 | real(jprb), intent(out), dimension(:) :: gamma1, gamma2, gamma3 |
---|
| 125 | |
---|
| 126 | real(jprb) :: factor |
---|
| 127 | |
---|
| 128 | integer :: jg |
---|
| 129 | |
---|
| 130 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 131 | real(jprb) :: hook_handle |
---|
| 132 | |
---|
| 133 | if (lhook) call dr_hook('radiation_two_stream:calc_two_stream_gammas_sw',0,hook_handle) |
---|
| 134 | #endif |
---|
| 135 | |
---|
| 136 | ! Zdunkowski "PIFM" (Zdunkowski et al., 1980; Contributions to |
---|
| 137 | ! Atmospheric Physics 53, 147-66) |
---|
| 138 | do jg = 1, ng |
---|
| 139 | ! gamma1(jg) = 2.0_jprb - ssa(jg) * (1.25_jprb + 0.75_jprb*g(jg)) |
---|
| 140 | ! gamma2(jg) = 0.75_jprb *(ssa(jg) * (1.0_jprb - g(jg))) |
---|
| 141 | ! gamma3(jg) = 0.5_jprb - (0.75_jprb*mu0)*g(jg) |
---|
| 142 | ! Optimized version: |
---|
| 143 | factor = 0.75_jprb*g(jg) |
---|
| 144 | gamma1(jg) = 2.0_jprb - ssa(jg) * (1.25_jprb + factor) |
---|
| 145 | gamma2(jg) = ssa(jg) * (0.75_jprb - factor) |
---|
| 146 | gamma3(jg) = 0.5_jprb - mu0*factor |
---|
| 147 | end do |
---|
| 148 | |
---|
| 149 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 150 | if (lhook) call dr_hook('radiation_two_stream:calc_two_stream_gammas_sw',1,hook_handle) |
---|
| 151 | #endif |
---|
| 152 | |
---|
| 153 | end subroutine calc_two_stream_gammas_sw |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | !--------------------------------------------------------------------- |
---|
| 157 | ! Compute the longwave reflectance and transmittance to diffuse |
---|
| 158 | ! radiation using the Meador & Weaver formulas, as well as the |
---|
| 159 | ! upward flux at the top and the downward flux at the base of the |
---|
| 160 | ! layer due to emission from within the layer assuming a linear |
---|
| 161 | ! variation of Planck function within the layer. |
---|
| 162 | subroutine calc_reflectance_transmittance_lw(ng, & |
---|
| 163 | & od, gamma1, gamma2, planck_top, planck_bot, & |
---|
| 164 | & reflectance, transmittance, source_up, source_dn) |
---|
| 165 | |
---|
| 166 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 167 | use yomhook, only : lhook, dr_hook |
---|
| 168 | #endif |
---|
| 169 | |
---|
| 170 | integer, intent(in) :: ng |
---|
| 171 | |
---|
| 172 | ! Optical depth and single scattering albedo |
---|
| 173 | real(jprb), intent(in), dimension(ng) :: od |
---|
| 174 | |
---|
| 175 | ! The two transfer coefficients from the two-stream |
---|
| 176 | ! differentiatial equations (computed by |
---|
| 177 | ! calc_two_stream_gammas_lw) |
---|
| 178 | real(jprb), intent(in), dimension(ng) :: gamma1, gamma2 |
---|
| 179 | |
---|
| 180 | ! The Planck terms (functions of temperature) at the top and |
---|
| 181 | ! bottom of the layer |
---|
| 182 | real(jprb), intent(in), dimension(ng) :: planck_top, planck_bot |
---|
| 183 | |
---|
| 184 | ! The diffuse reflectance and transmittance, i.e. the fraction of |
---|
| 185 | ! diffuse radiation incident on a layer from either top or bottom |
---|
| 186 | ! that is reflected back or transmitted through |
---|
| 187 | real(jprb), intent(out), dimension(ng) :: reflectance, transmittance |
---|
| 188 | |
---|
| 189 | ! The upward emission at the top of the layer and the downward |
---|
| 190 | ! emission at its base, due to emission from within the layer |
---|
| 191 | real(jprb), intent(out), dimension(ng) :: source_up, source_dn |
---|
| 192 | |
---|
| 193 | real(jprd) :: k_exponent, reftrans_factor |
---|
| 194 | real(jprd) :: exponential ! = exp(-k_exponent*od) |
---|
| 195 | real(jprd) :: exponential2 ! = exp(-2*k_exponent*od) |
---|
| 196 | |
---|
| 197 | real(jprd) :: coeff, coeff_up_top, coeff_up_bot, coeff_dn_top, coeff_dn_bot |
---|
| 198 | |
---|
| 199 | integer :: jg |
---|
| 200 | |
---|
| 201 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 202 | real(jprb) :: hook_handle |
---|
| 203 | |
---|
| 204 | if (lhook) call dr_hook('radiation_two_stream:calc_reflectance_transmittance_lw',0,hook_handle) |
---|
| 205 | #endif |
---|
| 206 | |
---|
| 207 | do jg = 1, ng |
---|
| 208 | if (od(jg) > 1.0e-3_jprd) then |
---|
| 209 | k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & |
---|
| 210 | 1.E-12_jprd)) ! Eq 18 of Meador & Weaver (1980) |
---|
| 211 | exponential = exp_fast(-k_exponent*od(jg)) |
---|
| 212 | exponential2 = exponential*exponential |
---|
| 213 | reftrans_factor = 1.0 / (k_exponent + gamma1(jg) + (k_exponent - gamma1(jg))*exponential2) |
---|
| 214 | ! Meador & Weaver (1980) Eq. 25 |
---|
| 215 | reflectance(jg) = gamma2(jg) * (1.0_jprd - exponential2) * reftrans_factor |
---|
| 216 | ! Meador & Weaver (1980) Eq. 26 |
---|
| 217 | transmittance(jg) = 2.0_jprd * k_exponent * exponential * reftrans_factor |
---|
| 218 | |
---|
| 219 | ! Compute upward and downward emission assuming the Planck |
---|
| 220 | ! function to vary linearly with optical depth within the layer |
---|
| 221 | ! (e.g. Wiscombe , JQSRT 1976). |
---|
| 222 | |
---|
| 223 | ! Stackhouse and Stephens (JAS 1991) Eqs 5 & 12 |
---|
| 224 | coeff = (planck_bot(jg)-planck_top(jg)) / (od(jg)*(gamma1(jg)+gamma2(jg))) |
---|
| 225 | coeff_up_top = coeff + planck_top(jg) |
---|
| 226 | coeff_up_bot = coeff + planck_bot(jg) |
---|
| 227 | coeff_dn_top = -coeff + planck_top(jg) |
---|
| 228 | coeff_dn_bot = -coeff + planck_bot(jg) |
---|
| 229 | source_up(jg) = coeff_up_top - reflectance(jg) * coeff_dn_top - transmittance(jg) * coeff_up_bot |
---|
| 230 | source_dn(jg) = coeff_dn_bot - reflectance(jg) * coeff_up_bot - transmittance(jg) * coeff_dn_top |
---|
| 231 | else |
---|
| 232 | k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & |
---|
| 233 | 1.E-12_jprd)) ! Eq 18 of Meador & Weaver (1980) |
---|
| 234 | reflectance(jg) = gamma2(jg) * od(jg) |
---|
| 235 | transmittance(jg) = (1.0_jprb - k_exponent*od(jg)) / (1.0_jprb + od(jg)*(gamma1(jg)-k_exponent)) |
---|
| 236 | source_up(jg) = (1.0_jprb - reflectance(jg) - transmittance(jg)) & |
---|
| 237 | & * 0.5 * (planck_top(jg) + planck_bot(jg)) |
---|
| 238 | source_dn(jg) = source_up(jg) |
---|
| 239 | end if |
---|
| 240 | end do |
---|
| 241 | |
---|
| 242 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 243 | if (lhook) call dr_hook('radiation_two_stream:calc_reflectance_transmittance_lw',1,hook_handle) |
---|
| 244 | #endif |
---|
| 245 | |
---|
| 246 | end subroutine calc_reflectance_transmittance_lw |
---|
| 247 | |
---|
| 248 | |
---|
| 249 | |
---|
| 250 | !--------------------------------------------------------------------- |
---|
| 251 | ! As calc_reflectance_transmittance_lw but for an isothermal layer |
---|
| 252 | subroutine calc_reflectance_transmittance_isothermal_lw(ng, & |
---|
| 253 | & od, gamma1, gamma2, planck, & |
---|
| 254 | & reflectance, transmittance, source) |
---|
| 255 | |
---|
| 256 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 257 | use yomhook, only : lhook, dr_hook |
---|
| 258 | #endif |
---|
| 259 | |
---|
| 260 | integer, intent(in) :: ng |
---|
| 261 | |
---|
| 262 | ! Optical depth and single scattering albedo |
---|
| 263 | real(jprb), intent(in), dimension(ng) :: od |
---|
| 264 | |
---|
| 265 | ! The two transfer coefficients from the two-stream |
---|
| 266 | ! differentiatial equations (computed by |
---|
| 267 | ! calc_two_stream_gammas_lw) |
---|
| 268 | real(jprb), intent(in), dimension(ng) :: gamma1, gamma2 |
---|
| 269 | |
---|
| 270 | ! The Planck terms (functions of temperature) constant through the |
---|
| 271 | ! layer |
---|
| 272 | real(jprb), intent(in), dimension(ng) :: planck |
---|
| 273 | |
---|
| 274 | ! The diffuse reflectance and transmittance, i.e. the fraction of |
---|
| 275 | ! diffuse radiation incident on a layer from either top or bottom |
---|
| 276 | ! that is reflected back or transmitted through |
---|
| 277 | real(jprb), intent(out), dimension(ng) :: reflectance, transmittance |
---|
| 278 | |
---|
| 279 | ! The upward emission at the top of the layer and the downward |
---|
| 280 | ! emission at its base, due to emission from within the layer |
---|
| 281 | real(jprb), intent(out), dimension(ng) :: source |
---|
| 282 | |
---|
| 283 | real(jprd) :: k_exponent, reftrans_factor |
---|
| 284 | real(jprd) :: exponential ! = exp(-k_exponent*od) |
---|
| 285 | real(jprd) :: exponential2 ! = exp(-2*k_exponent*od) |
---|
| 286 | |
---|
| 287 | integer :: jg |
---|
| 288 | |
---|
| 289 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 290 | real(jprb) :: hook_handle |
---|
| 291 | |
---|
| 292 | if (lhook) call dr_hook('radiation_two_stream:calc_reflectance_transmittance_isothermal_lw',0,hook_handle) |
---|
| 293 | #endif |
---|
| 294 | |
---|
| 295 | do jg = 1, ng |
---|
| 296 | k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & |
---|
| 297 | 1.E-12_jprd)) ! Eq 18 of Meador & Weaver (1980) |
---|
| 298 | exponential = exp_fast(-k_exponent*od(jg)) |
---|
| 299 | exponential2 = exponential*exponential |
---|
| 300 | reftrans_factor = 1.0 / (k_exponent + gamma1(jg) + (k_exponent - gamma1(jg))*exponential2) |
---|
| 301 | ! Meador & Weaver (1980) Eq. 25 |
---|
| 302 | reflectance(jg) = gamma2(jg) * (1.0_jprd - exponential2) * reftrans_factor |
---|
| 303 | ! Meador & Weaver (1980) Eq. 26 |
---|
| 304 | transmittance(jg) = 2.0_jprd * k_exponent * exponential * reftrans_factor |
---|
| 305 | |
---|
| 306 | ! Emissivity of layer is one minus reflectance minus |
---|
| 307 | ! transmittance, multiply by Planck function to get emitted |
---|
| 308 | ! ousrce |
---|
| 309 | source(jg) = planck(jg) * (1.0_jprd - reflectance(jg) - transmittance(jg)) |
---|
| 310 | end do |
---|
| 311 | |
---|
| 312 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 313 | if (lhook) call dr_hook('radiation_two_stream:calc_reflectance_transmittance_isothermal_lw',1,hook_handle) |
---|
| 314 | #endif |
---|
| 315 | |
---|
| 316 | end subroutine calc_reflectance_transmittance_isothermal_lw |
---|
| 317 | |
---|
| 318 | |
---|
| 319 | |
---|
| 320 | !--------------------------------------------------------------------- |
---|
| 321 | ! Compute the longwave transmittance to diffuse radiation in the |
---|
| 322 | ! no-scattering case, as well as the upward flux at the top and the |
---|
| 323 | ! downward flux at the base of the layer due to emission from within |
---|
| 324 | ! the layer assuming a linear variation of Planck function within |
---|
| 325 | ! the layer. |
---|
| 326 | subroutine calc_no_scattering_transmittance_lw(ng, & |
---|
| 327 | & od, planck_top, planck_bot, transmittance, source_up, source_dn) |
---|
| 328 | |
---|
| 329 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 330 | use yomhook, only : lhook, dr_hook |
---|
| 331 | #endif |
---|
| 332 | |
---|
| 333 | integer, intent(in) :: ng |
---|
| 334 | |
---|
| 335 | ! Optical depth and single scattering albedo |
---|
| 336 | real(jprb), intent(in), dimension(ng) :: od |
---|
| 337 | |
---|
| 338 | ! The Planck terms (functions of temperature) at the top and |
---|
| 339 | ! bottom of the layer |
---|
| 340 | real(jprb), intent(in), dimension(ng) :: planck_top, planck_bot |
---|
| 341 | |
---|
| 342 | ! The diffuse transmittance, i.e. the fraction of diffuse |
---|
| 343 | ! radiation incident on a layer from either top or bottom that is |
---|
| 344 | ! reflected back or transmitted through |
---|
| 345 | real(jprb), intent(out), dimension(ng) :: transmittance |
---|
| 346 | |
---|
| 347 | ! The upward emission at the top of the layer and the downward |
---|
| 348 | ! emission at its base, due to emission from within the layer |
---|
| 349 | real(jprb), intent(out), dimension(ng) :: source_up, source_dn |
---|
| 350 | |
---|
| 351 | real(jprd) :: coeff, coeff_up_top, coeff_up_bot, coeff_dn_top, coeff_dn_bot !, planck_mean |
---|
| 352 | |
---|
| 353 | integer :: jg |
---|
| 354 | |
---|
| 355 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 356 | real(jprb) :: hook_handle |
---|
| 357 | |
---|
| 358 | if (lhook) call dr_hook('radiation_two_stream:calc_no_scattering_transmittance_lw',0,hook_handle) |
---|
| 359 | #endif |
---|
| 360 | |
---|
| 361 | do jg = 1, ng |
---|
| 362 | ! Compute upward and downward emission assuming the Planck |
---|
| 363 | ! function to vary linearly with optical depth within the layer |
---|
| 364 | ! (e.g. Wiscombe , JQSRT 1976). |
---|
| 365 | if (od(jg) > 1.0e-3) then |
---|
| 366 | ! Simplified from calc_reflectance_transmittance_lw above |
---|
| 367 | coeff = LwDiffusivity*od(jg) |
---|
| 368 | transmittance(jg) = exp_fast(-coeff) |
---|
| 369 | coeff = (planck_bot(jg)-planck_top(jg)) / coeff |
---|
| 370 | coeff_up_top = coeff + planck_top(jg) |
---|
| 371 | coeff_up_bot = coeff + planck_bot(jg) |
---|
| 372 | coeff_dn_top = -coeff + planck_top(jg) |
---|
| 373 | coeff_dn_bot = -coeff + planck_bot(jg) |
---|
| 374 | source_up(jg) = coeff_up_top - transmittance(jg) * coeff_up_bot |
---|
| 375 | source_dn(jg) = coeff_dn_bot - transmittance(jg) * coeff_dn_top |
---|
| 376 | else |
---|
| 377 | ! Linear limit at low optical depth |
---|
| 378 | coeff = LwDiffusivity*od(jg) |
---|
| 379 | transmittance(jg) = 1.0_jprb - coeff |
---|
| 380 | source_up(jg) = coeff * 0.5_jprb * (planck_top(jg)+planck_bot(jg)) |
---|
| 381 | source_dn(jg) = source_up(jg) |
---|
| 382 | end if |
---|
| 383 | end do |
---|
| 384 | |
---|
| 385 | ! Method in the older IFS radiation scheme |
---|
| 386 | ! do j = 1, n |
---|
| 387 | ! coeff = od(jg) / (3.59712_jprd + od(jg)) |
---|
| 388 | ! planck_mean = 0.5_jprd * (planck_top(jg) + planck_bot(jg)) |
---|
| 389 | ! |
---|
| 390 | ! source_up(jg) = (1.0_jprd-transmittance(jg)) * (planck_mean + (planck_top(jg) - planck_mean) * coeff) |
---|
| 391 | ! source_dn(jg) = (1.0_jprd-transmittance(jg)) * (planck_mean + (planck_bot(jg) - planck_mean) * coeff) |
---|
| 392 | ! end do |
---|
| 393 | |
---|
| 394 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 395 | if (lhook) call dr_hook('radiation_two_stream:calc_no_scattering_transmittance_lw',1,hook_handle) |
---|
| 396 | #endif |
---|
| 397 | |
---|
| 398 | end subroutine calc_no_scattering_transmittance_lw |
---|
| 399 | |
---|
| 400 | |
---|
| 401 | !--------------------------------------------------------------------- |
---|
| 402 | ! Compute the shortwave reflectance and transmittance to diffuse |
---|
| 403 | ! radiation using the Meador & Weaver formulas, as well as the |
---|
| 404 | ! "direct" reflection and transmission, which really means the rate |
---|
| 405 | ! of transfer of direct solar radiation (into a plane perpendicular |
---|
| 406 | ! to the direct beam) into diffuse upward and downward streams at |
---|
| 407 | ! the top and bottom of the layer, respectively. Finally, |
---|
| 408 | ! trans_dir_dir is the transmittance of the atmosphere to direct |
---|
| 409 | ! radiation with no scattering. |
---|
| 410 | subroutine calc_reflectance_transmittance_sw(ng, mu0, od, ssa, & |
---|
| 411 | & gamma1, gamma2, gamma3, ref_diff, trans_diff, & |
---|
| 412 | & ref_dir, trans_dir_diff, trans_dir_dir) |
---|
| 413 | |
---|
| 414 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 415 | use yomhook, only : lhook, dr_hook |
---|
| 416 | #endif |
---|
| 417 | |
---|
| 418 | integer, intent(in) :: ng |
---|
| 419 | |
---|
| 420 | ! Cosine of solar zenith angle |
---|
| 421 | real(jprb), intent(in) :: mu0 |
---|
| 422 | |
---|
| 423 | ! Optical depth and single scattering albedo |
---|
| 424 | real(jprb), intent(in), dimension(ng) :: od, ssa |
---|
| 425 | |
---|
| 426 | ! The three transfer coefficients from the two-stream |
---|
| 427 | ! differentiatial equations (computed by calc_two_stream_gammas) |
---|
| 428 | real(jprb), intent(in), dimension(ng) :: gamma1, gamma2, gamma3 |
---|
| 429 | |
---|
| 430 | ! The direct reflectance and transmittance, i.e. the fraction of |
---|
| 431 | ! incoming direct solar radiation incident at the top of a layer |
---|
| 432 | ! that is either reflected back (ref_dir) or scattered but |
---|
| 433 | ! transmitted through the layer to the base (trans_dir_diff) |
---|
| 434 | real(jprb), intent(out), dimension(ng) :: ref_dir, trans_dir_diff |
---|
| 435 | |
---|
| 436 | ! The diffuse reflectance and transmittance, i.e. the fraction of |
---|
| 437 | ! diffuse radiation incident on a layer from either top or bottom |
---|
| 438 | ! that is reflected back or transmitted through |
---|
| 439 | real(jprb), intent(out), dimension(ng) :: ref_diff, trans_diff |
---|
| 440 | |
---|
| 441 | ! Transmittance of the direct been with no scattering |
---|
| 442 | real(jprb), intent(out), dimension(ng) :: trans_dir_dir |
---|
| 443 | |
---|
| 444 | real(jprd) :: gamma4, alpha1, alpha2, k_exponent, reftrans_factor |
---|
| 445 | real(jprd) :: exponential0 ! = exp(-od/mu0) |
---|
| 446 | real(jprd) :: exponential ! = exp(-k_exponent*od) |
---|
| 447 | real(jprd) :: exponential2 ! = exp(-2*k_exponent*od) |
---|
| 448 | real(jprd) :: k_mu0, k_gamma3, k_gamma4 |
---|
| 449 | real(jprd) :: k_2_exponential, od_over_mu0 |
---|
| 450 | integer :: jg |
---|
| 451 | |
---|
| 452 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 453 | real(jprb) :: hook_handle |
---|
| 454 | |
---|
| 455 | if (lhook) call dr_hook('radiation_two_stream:calc_reflectance_transmittance_sw',0,hook_handle) |
---|
| 456 | #endif |
---|
| 457 | |
---|
| 458 | do jg = 1, ng |
---|
| 459 | od_over_mu0 = max(od(jg) / mu0, 0.0_jprd) |
---|
| 460 | ! In the IFS this appears to be faster without testing the value |
---|
| 461 | ! of od_over_mu0: |
---|
| 462 | if (.true.) then |
---|
| 463 | ! if (od_over_mu0 > 1.0e-6_jprd) then |
---|
| 464 | gamma4 = 1.0_jprd - gamma3(jg) |
---|
| 465 | alpha1 = gamma1(jg)*gamma4 + gamma2(jg)*gamma3(jg) ! Eq. 16 |
---|
| 466 | alpha2 = gamma1(jg)*gamma3(jg) + gamma2(jg)*gamma4 ! Eq. 17 |
---|
| 467 | |
---|
| 468 | ! Note that if the minimum value is reduced (e.g. to 1.0e-24) |
---|
| 469 | ! then noise starts to appear as a function of solar zenith |
---|
| 470 | ! angle |
---|
| 471 | k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & |
---|
| 472 | & 1.0e-12_jprd)) ! Eq 18 |
---|
| 473 | k_mu0 = k_exponent*mu0 |
---|
| 474 | k_gamma3 = k_exponent*gamma3(jg) |
---|
| 475 | k_gamma4 = k_exponent*gamma4 |
---|
| 476 | ! Check for mu0 <= 0! |
---|
| 477 | exponential0 = exp_fast(-od_over_mu0) |
---|
| 478 | trans_dir_dir(jg) = exponential0 |
---|
| 479 | exponential = exp_fast(-k_exponent*od(jg)) |
---|
| 480 | |
---|
| 481 | exponential2 = exponential*exponential |
---|
| 482 | k_2_exponential = 2.0_jprd * k_exponent * exponential |
---|
| 483 | |
---|
| 484 | if (k_mu0 == 1.0_jprd) then |
---|
| 485 | k_mu0 = 1.0_jprd - 10.0_jprd*epsilon(1.0_jprd) |
---|
| 486 | end if |
---|
| 487 | |
---|
| 488 | reftrans_factor = 1.0_jprd / (k_exponent + gamma1(jg) + (k_exponent - gamma1(jg))*exponential2) |
---|
| 489 | |
---|
| 490 | ! Meador & Weaver (1980) Eq. 25 |
---|
| 491 | ref_diff(jg) = gamma2(jg) * (1.0_jprd - exponential2) * reftrans_factor |
---|
| 492 | |
---|
| 493 | ! Meador & Weaver (1980) Eq. 26 |
---|
| 494 | trans_diff(jg) = k_2_exponential * reftrans_factor |
---|
| 495 | |
---|
| 496 | ! Here we need mu0 even though it wasn't in Meador and Weaver |
---|
| 497 | ! because we are assuming the incoming direct flux is defined |
---|
| 498 | ! to be the flux into a plane perpendicular to the direction of |
---|
| 499 | ! the sun, not into a horizontal plane |
---|
| 500 | reftrans_factor = mu0 * ssa(jg) * reftrans_factor / (1.0_jprd - k_mu0*k_mu0) |
---|
| 501 | |
---|
| 502 | ! Meador & Weaver (1980) Eq. 14, multiplying top & bottom by |
---|
| 503 | ! exp(-k_exponent*od) in case of very high optical depths |
---|
| 504 | ref_dir(jg) = reftrans_factor & |
---|
| 505 | & * ( (1.0_jprd - k_mu0) * (alpha2 + k_gamma3) & |
---|
| 506 | & -(1.0_jprd + k_mu0) * (alpha2 - k_gamma3)*exponential2 & |
---|
| 507 | & -k_2_exponential*(gamma3(jg) - alpha2*mu0)*exponential0) |
---|
| 508 | |
---|
| 509 | ! Meador & Weaver (1980) Eq. 15, multiplying top & bottom by |
---|
| 510 | ! exp(-k_exponent*od), minus the 1*exp(-od/mu0) term representing direct |
---|
| 511 | ! unscattered transmittance. |
---|
| 512 | trans_dir_diff(jg) = reftrans_factor * ( k_2_exponential*(gamma4 + alpha1*mu0) & |
---|
| 513 | & - exponential0 & |
---|
| 514 | & * ( (1.0_jprd + k_mu0) * (alpha1 + k_gamma4) & |
---|
| 515 | & -(1.0_jprd - k_mu0) * (alpha1 - k_gamma4) * exponential2) ) |
---|
| 516 | |
---|
| 517 | else |
---|
| 518 | ! Low optical-depth limit; see equations 19, 20 and 27 from |
---|
| 519 | ! Meador & Weaver (1980) |
---|
| 520 | trans_diff(jg) = 1.0_jprb - gamma1(jg) * od(jg) |
---|
| 521 | ref_diff(jg) = gamma2(jg) * od(jg) |
---|
| 522 | trans_dir_diff(jg) = (1.0_jprb - gamma3(jg)) * ssa(jg) * od(jg) |
---|
| 523 | ref_dir(jg) = gamma3(jg) * ssa(jg) * od(jg) |
---|
| 524 | trans_dir_dir(jg) = 1.0_jprd - od_over_mu0 |
---|
| 525 | end if |
---|
| 526 | end do |
---|
| 527 | |
---|
| 528 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 529 | if (lhook) call dr_hook('radiation_two_stream:calc_reflectance_transmittance_sw',1,hook_handle) |
---|
| 530 | #endif |
---|
| 531 | |
---|
| 532 | end subroutine calc_reflectance_transmittance_sw |
---|
| 533 | |
---|
| 534 | !--------------------------------------------------------------------- |
---|
| 535 | ! As above but with height as a vertical coordinate rather than |
---|
| 536 | ! optical depth |
---|
| 537 | subroutine calc_reflectance_transmittance_z_sw(ng, mu0, depth, & |
---|
| 538 | & gamma0, gamma1, gamma2, gamma3, gamma4, & |
---|
| 539 | & ref_diff, trans_diff, ref_dir, trans_dir_diff, trans_dir_dir) |
---|
| 540 | |
---|
| 541 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 542 | use yomhook, only : lhook, dr_hook |
---|
| 543 | #endif |
---|
| 544 | |
---|
| 545 | integer, intent(in) :: ng |
---|
| 546 | |
---|
| 547 | ! Cosine of solar zenith angle |
---|
| 548 | real(jprb), intent(in) :: mu0 |
---|
| 549 | |
---|
| 550 | ! Layer depth |
---|
| 551 | real(jprb), intent(in) :: depth |
---|
| 552 | |
---|
| 553 | ! The four transfer coefficients from the two-stream |
---|
| 554 | ! differentiatial equations |
---|
| 555 | real(jprb), intent(in), dimension(ng) :: gamma1, gamma2, gamma3, gamma4 |
---|
| 556 | |
---|
| 557 | ! An additional coefficient for direct unscattered flux "Fdir" |
---|
| 558 | ! such that dFdir/dz = -gamma0*Fdir |
---|
| 559 | real(jprb), intent(in), dimension(ng) :: gamma0 |
---|
| 560 | |
---|
| 561 | ! The direct reflectance and transmittance, i.e. the fraction of |
---|
| 562 | ! incoming direct solar radiation incident at the top of a layer |
---|
| 563 | ! that is either reflected back (ref_dir) or scattered but |
---|
| 564 | ! transmitted through the layer to the base (trans_dir_diff) |
---|
| 565 | real(jprb), intent(out), dimension(ng) :: ref_dir, trans_dir_diff |
---|
| 566 | |
---|
| 567 | ! The diffuse reflectance and transmittance, i.e. the fraction of |
---|
| 568 | ! diffuse radiation incident on a layer from either top or bottom |
---|
| 569 | ! that is reflected back or transmitted through |
---|
| 570 | real(jprb), intent(out), dimension(ng) :: ref_diff, trans_diff |
---|
| 571 | |
---|
| 572 | ! Transmittance of the direct been with no scattering |
---|
| 573 | real(jprb), intent(out), dimension(ng) :: trans_dir_dir |
---|
| 574 | |
---|
| 575 | real(jprd) :: alpha1, alpha2, k_exponent, reftrans_factor |
---|
| 576 | real(jprd) :: exponential0 ! = exp(-od/mu0) |
---|
| 577 | real(jprd) :: exponential ! = exp(-k_exponent*od) |
---|
| 578 | real(jprd) :: exponential2 ! = exp(-2*k_exponent*od) |
---|
| 579 | real(jprd) :: k_mu0, k_gamma3, k_gamma4 |
---|
| 580 | real(jprd) :: k_2_exponential, od_over_mu0 |
---|
| 581 | integer :: jg |
---|
| 582 | |
---|
| 583 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 584 | real(jprb) :: hook_handle |
---|
| 585 | |
---|
| 586 | if (lhook) call dr_hook('radiation_two_stream:calc_reflectance_transmittance_z_sw',0,hook_handle) |
---|
| 587 | #endif |
---|
| 588 | |
---|
| 589 | do jg = 1, ng |
---|
| 590 | od_over_mu0 = max(gamma0(jg) * depth, 0.0_jprd) |
---|
| 591 | ! In the IFS this appears to be faster without testing the value |
---|
| 592 | ! of od_over_mu0: |
---|
| 593 | if (.true.) then |
---|
| 594 | ! if (od_over_mu0 > 1.0e-6_jprd) then |
---|
| 595 | alpha1 = gamma1(jg)*gamma4(jg) + gamma2(jg)*gamma3(jg) ! Eq. 16 |
---|
| 596 | alpha2 = gamma1(jg)*gamma3(jg) + gamma2(jg)*gamma4(jg) ! Eq. 17 |
---|
| 597 | |
---|
| 598 | ! Note that if the minimum value is reduced (e.g. to 1.0e-24) |
---|
| 599 | ! then noise starts to appear as a function of solar zenith |
---|
| 600 | ! angle |
---|
| 601 | k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & |
---|
| 602 | & 1.0e-12_jprd)) ! Eq 18 |
---|
| 603 | k_mu0 = k_exponent*mu0 |
---|
| 604 | k_gamma3 = k_exponent*gamma3(jg) |
---|
| 605 | k_gamma4 = k_exponent*gamma4(jg) |
---|
| 606 | ! Check for mu0 <= 0! |
---|
| 607 | exponential0 = exp_fast(-od_over_mu0) |
---|
| 608 | trans_dir_dir(jg) = exponential0 |
---|
| 609 | exponential = exp_fast(-k_exponent*depth) |
---|
| 610 | |
---|
| 611 | exponential2 = exponential*exponential |
---|
| 612 | k_2_exponential = 2.0_jprd * k_exponent * exponential |
---|
| 613 | |
---|
| 614 | if (k_mu0 == 1.0_jprd) then |
---|
| 615 | k_mu0 = 1.0_jprd - 10.0_jprd*epsilon(1.0_jprd) |
---|
| 616 | end if |
---|
| 617 | |
---|
| 618 | reftrans_factor = 1.0_jprd / (k_exponent + gamma1(jg) + (k_exponent - gamma1(jg))*exponential2) |
---|
| 619 | |
---|
| 620 | ! Meador & Weaver (1980) Eq. 25 |
---|
| 621 | ref_diff(jg) = gamma2(jg) * (1.0_jprd - exponential2) * reftrans_factor |
---|
| 622 | |
---|
| 623 | ! Meador & Weaver (1980) Eq. 26 |
---|
| 624 | trans_diff(jg) = k_2_exponential * reftrans_factor |
---|
| 625 | |
---|
| 626 | ! Here we need mu0 even though it wasn't in Meador and Weaver |
---|
| 627 | ! because we are assuming the incoming direct flux is defined |
---|
| 628 | ! to be the flux into a plane perpendicular to the direction of |
---|
| 629 | ! the sun, not into a horizontal plane |
---|
| 630 | reftrans_factor = mu0 * reftrans_factor / (1.0_jprd - k_mu0*k_mu0) |
---|
| 631 | |
---|
| 632 | ! Meador & Weaver (1980) Eq. 14, multiplying top & bottom by |
---|
| 633 | ! exp(-k_exponent*od) in case of very high optical depths |
---|
| 634 | ref_dir(jg) = reftrans_factor & |
---|
| 635 | & * ( (1.0_jprd - k_mu0) * (alpha2 + k_gamma3) & |
---|
| 636 | & -(1.0_jprd + k_mu0) * (alpha2 - k_gamma3)*exponential2 & |
---|
| 637 | & -k_2_exponential*(gamma3(jg) - alpha2*mu0)*exponential0) |
---|
| 638 | |
---|
| 639 | ! Meador & Weaver (1980) Eq. 15, multiplying top & bottom by |
---|
| 640 | ! exp(-k_exponent*od), minus the 1*exp(-od/mu0) term representing direct |
---|
| 641 | ! unscattered transmittance. |
---|
| 642 | trans_dir_diff(jg) = reftrans_factor * ( k_2_exponential*(gamma4(jg) + alpha1*mu0) & |
---|
| 643 | & - exponential0 & |
---|
| 644 | & * ( (1.0_jprd + k_mu0) * (alpha1 + k_gamma4) & |
---|
| 645 | & -(1.0_jprd - k_mu0) * (alpha1 - k_gamma4) * exponential2) ) |
---|
| 646 | |
---|
| 647 | else |
---|
| 648 | ! Low optical-depth limit; see equations 19, 20 and 27 from |
---|
| 649 | ! Meador & Weaver (1980) |
---|
| 650 | trans_diff(jg) = 1.0_jprb - gamma1(jg) * depth |
---|
| 651 | ref_diff(jg) = gamma2(jg) * depth |
---|
| 652 | trans_dir_diff(jg) = (1.0_jprb - gamma3(jg)) * depth |
---|
| 653 | ref_dir(jg) = gamma3(jg) * depth |
---|
| 654 | trans_dir_dir(jg) = 1.0_jprd - od_over_mu0 |
---|
| 655 | end if |
---|
| 656 | end do |
---|
| 657 | |
---|
| 658 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 659 | if (lhook) call dr_hook('radiation_two_stream:calc_reflectance_transmittance_z_sw',1,hook_handle) |
---|
| 660 | #endif |
---|
| 661 | |
---|
| 662 | end subroutine calc_reflectance_transmittance_z_sw |
---|
| 663 | |
---|
| 664 | |
---|
| 665 | !--------------------------------------------------------------------- |
---|
| 666 | ! Compute the fraction of shortwave transmitted diffuse radiation |
---|
| 667 | ! that is scattered during its transmission, used to compute |
---|
| 668 | ! entrapment in SPARTACUS |
---|
| 669 | subroutine calc_frac_scattered_diffuse_sw(ng, od, & |
---|
| 670 | & gamma1, gamma2, frac_scat_diffuse) |
---|
| 671 | |
---|
| 672 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 673 | use yomhook, only : lhook, dr_hook |
---|
| 674 | #endif |
---|
| 675 | |
---|
| 676 | integer, intent(in) :: ng |
---|
| 677 | |
---|
| 678 | ! Optical depth |
---|
| 679 | real(jprb), intent(in), dimension(ng) :: od |
---|
| 680 | |
---|
| 681 | ! The first two transfer coefficients from the two-stream |
---|
| 682 | ! differentiatial equations (computed by calc_two_stream_gammas) |
---|
| 683 | real(jprb), intent(in), dimension(ng) :: gamma1, gamma2 |
---|
| 684 | |
---|
| 685 | ! The fraction of shortwave transmitted diffuse radiation that is |
---|
| 686 | ! scattered during its transmission |
---|
| 687 | real(jprb), intent(out), dimension(ng) :: frac_scat_diffuse |
---|
| 688 | |
---|
| 689 | real(jprd) :: k_exponent, reftrans_factor |
---|
| 690 | real(jprd) :: exponential ! = exp(-k_exponent*od) |
---|
| 691 | real(jprd) :: exponential2 ! = exp(-2*k_exponent*od) |
---|
| 692 | real(jprd) :: k_2_exponential |
---|
| 693 | integer :: jg |
---|
| 694 | |
---|
| 695 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 696 | real(jprb) :: hook_handle |
---|
| 697 | |
---|
| 698 | if (lhook) call dr_hook('radiation_two_stream:calc_frac_scattered_diffuse_sw',0,hook_handle) |
---|
| 699 | #endif |
---|
| 700 | |
---|
| 701 | do jg = 1, ng |
---|
| 702 | ! Note that if the minimum value is reduced (e.g. to 1.0e-24) |
---|
| 703 | ! then noise starts to appear as a function of solar zenith |
---|
| 704 | ! angle |
---|
| 705 | k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & |
---|
| 706 | & 1.0e-12_jprd)) ! Eq 18 |
---|
| 707 | exponential = exp_fast(-k_exponent*od(jg)) |
---|
| 708 | exponential2 = exponential*exponential |
---|
| 709 | k_2_exponential = 2.0_jprd * k_exponent * exponential |
---|
| 710 | |
---|
| 711 | reftrans_factor = 1.0_jprd / (k_exponent + gamma1(jg) + (k_exponent - gamma1(jg))*exponential2) |
---|
| 712 | |
---|
| 713 | ! Meador & Weaver (1980) Eq. 26. |
---|
| 714 | ! Until 1.1.8, used LwDiffusivity instead of 2.0, although the |
---|
| 715 | ! effect is very small |
---|
| 716 | ! frac_scat_diffuse(jg) = 1.0_jprb - min(1.0_jprb,exp_fast(-LwDiffusivity*od(jg)) & |
---|
| 717 | ! & / max(1.0e-8_jprb, k_2_exponential * reftrans_factor)) |
---|
| 718 | frac_scat_diffuse(jg) = 1.0_jprb & |
---|
| 719 | & - min(1.0_jprb,exp_fast(-2.0_jprb*od(jg)) & |
---|
| 720 | & / max(1.0e-8_jprb, k_2_exponential * reftrans_factor)) |
---|
| 721 | end do |
---|
| 722 | |
---|
| 723 | #ifdef DO_DR_HOOK_TWO_STREAM |
---|
| 724 | if (lhook) call dr_hook('radiation_two_stream:calc_frac_scattered_diffuse_sw',1,hook_handle) |
---|
| 725 | #endif |
---|
| 726 | |
---|
| 727 | end subroutine calc_frac_scattered_diffuse_sw |
---|
| 728 | |
---|
| 729 | end module radiation_two_stream |
---|