[3908] | 1 | ! radiation_cloud_generator.F90 - Generate water-content or optical-depth scalings for McICA |
---|
| 2 | ! |
---|
| 3 | ! (C) Copyright 2015- 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 | ! Generate clouds for McICA using a method modified from Raisanen et |
---|
| 16 | ! al. (2002) |
---|
| 17 | ! |
---|
| 18 | ! Modifications |
---|
| 19 | ! 2018-02-22 R. Hogan Call masked version of PDF sampler for speed |
---|
[4489] | 20 | ! 2020-03-31 R. Hogan More vectorizable version of Exp-Ran |
---|
[3908] | 21 | |
---|
| 22 | module radiation_cloud_generator |
---|
| 23 | |
---|
| 24 | public |
---|
| 25 | |
---|
| 26 | contains |
---|
| 27 | |
---|
| 28 | !--------------------------------------------------------------------- |
---|
| 29 | ! Generate scaling factors for the cloud optical depth to represent |
---|
| 30 | ! cloud overlap, the overlap of internal cloud inhomogeneities and |
---|
| 31 | ! the fractional standard deviation of these inhomogeneities, for |
---|
| 32 | ! use in a Monte Carlo Independent Column Approximation radiation |
---|
| 33 | ! scheme. All returned profiles contain cloud, and the total cloud |
---|
| 34 | ! cover is also returned, so the calling function can then do a |
---|
| 35 | ! weighted average of clear and cloudy skies; this is a way to |
---|
| 36 | ! reduce the Monte Carlo noise in profiles with low cloud cover. |
---|
| 37 | subroutine cloud_generator(ng, nlev, i_overlap_scheme, & |
---|
| 38 | & iseed, frac_threshold, & |
---|
| 39 | & frac, overlap_param, decorrelation_scaling, & |
---|
| 40 | & fractional_std, pdf_sampler, & |
---|
| 41 | & od_scaling, total_cloud_cover, & |
---|
[4489] | 42 | & use_beta_overlap, use_vectorizable_generator) |
---|
[3908] | 43 | |
---|
| 44 | use parkind1, only : jprb |
---|
| 45 | use yomhook, only : lhook, dr_hook |
---|
| 46 | use radiation_io, only : nulerr, radiation_abort |
---|
| 47 | use random_numbers_mix, only : randomnumberstream, & |
---|
| 48 | initialize_random_numbers, uniform_distribution |
---|
| 49 | use radiation_pdf_sampler, only : pdf_sampler_type |
---|
| 50 | use radiation_cloud_cover, only : cum_cloud_cover_exp_ran, & |
---|
| 51 | & cum_cloud_cover_max_ran, cum_cloud_cover_exp_exp, & |
---|
| 52 | & IOverlapMaximumRandom, IOverlapExponentialRandom, & |
---|
| 53 | & IOverlapExponential |
---|
| 54 | |
---|
| 55 | implicit none |
---|
| 56 | |
---|
| 57 | ! Inputs |
---|
| 58 | integer, intent(in) :: ng ! number of g points |
---|
| 59 | integer, intent(in) :: nlev ! number of model levels |
---|
| 60 | integer, intent(in) :: i_overlap_scheme |
---|
| 61 | integer, intent(in) :: iseed ! seed for random number generator |
---|
| 62 | |
---|
| 63 | ! Only cloud fractions above this threshold are considered to be |
---|
| 64 | ! clouds |
---|
| 65 | real(jprb), intent(in) :: frac_threshold |
---|
| 66 | |
---|
| 67 | ! Cloud fraction on full levels |
---|
| 68 | real(jprb), intent(in) :: frac(nlev) |
---|
| 69 | |
---|
| 70 | ! Cloud overlap parameter for interfaces between model layers, |
---|
| 71 | ! where 0 indicates random overlap and 1 indicates maximum-random |
---|
| 72 | ! overlap |
---|
| 73 | real(jprb), intent(in) :: overlap_param(nlev-1) |
---|
| 74 | |
---|
| 75 | ! Overlap parameter for internal inhomogeneities |
---|
| 76 | real(jprb), intent(in) :: decorrelation_scaling |
---|
| 77 | |
---|
| 78 | ! Fractional standard deviation at each layer |
---|
| 79 | real(jprb), intent(in) :: fractional_std(nlev) |
---|
| 80 | |
---|
| 81 | ! Object for sampling from a lognormal or gamma distribution |
---|
| 82 | type(pdf_sampler_type), intent(in) :: pdf_sampler |
---|
| 83 | |
---|
| 84 | ! This routine has been coded using the "alpha" overlap parameter |
---|
| 85 | ! of Hogan and Illingworth (2000). If the following logical is |
---|
| 86 | ! present and true then the input is interpretted to be the "beta" |
---|
| 87 | ! overlap parameter of Shonk et al. (2010), and needs to be |
---|
| 88 | ! converted to alpha. |
---|
[4489] | 89 | logical, intent(in), optional :: use_beta_overlap |
---|
[3908] | 90 | |
---|
[4489] | 91 | ! Do we use the more vectorizable cloud generator, at the expense |
---|
| 92 | ! of more random numbers being needed? |
---|
| 93 | logical, intent(in), optional :: use_vectorizable_generator |
---|
| 94 | |
---|
[3908] | 95 | ! Outputs |
---|
| 96 | |
---|
| 97 | ! Cloud optical depth scaling factor, with 0 indicating clear sky |
---|
| 98 | real(jprb), intent(out) :: od_scaling(ng,nlev) |
---|
| 99 | |
---|
| 100 | ! Total cloud cover using cloud fraction and overlap parameter |
---|
| 101 | real(jprb), intent(out) :: total_cloud_cover |
---|
| 102 | |
---|
| 103 | ! Local variables |
---|
| 104 | |
---|
| 105 | ! Cumulative cloud cover from TOA to the base of each layer |
---|
| 106 | real(jprb) :: cum_cloud_cover(nlev) |
---|
| 107 | |
---|
| 108 | ! Scaled random number for finding cloud |
---|
| 109 | real(jprb) :: trigger |
---|
| 110 | |
---|
| 111 | ! Uniform deviates between 0 and 1 |
---|
| 112 | real(jprb) :: rand_top(ng) |
---|
| 113 | |
---|
| 114 | ! Overlap parameter of inhomogeneities |
---|
| 115 | real(jprb) :: overlap_param_inhom(nlev-1) |
---|
| 116 | |
---|
| 117 | ! Seed for random number generator and stream for producing random |
---|
| 118 | ! numbers |
---|
| 119 | type(randomnumberstream) :: random_stream |
---|
| 120 | |
---|
| 121 | ! First and last cloudy layers |
---|
| 122 | integer :: ibegin, iend |
---|
| 123 | |
---|
| 124 | integer :: itrigger |
---|
| 125 | |
---|
| 126 | ! Loop index for model level and g-point |
---|
| 127 | integer :: jlev, jg |
---|
| 128 | |
---|
| 129 | ! Cloud cover of a pair of layers, and amount by which cloud at |
---|
| 130 | ! next level increases total cloud cover as seen from above |
---|
| 131 | real(jprb), dimension(nlev-1) :: pair_cloud_cover, overhang |
---|
| 132 | |
---|
[4489] | 133 | logical :: use_vec_gen |
---|
| 134 | |
---|
[3908] | 135 | real(jprb) :: hook_handle |
---|
| 136 | |
---|
| 137 | if (lhook) call dr_hook('radiation_cloud_generator:cloud_generator',0,hook_handle) |
---|
| 138 | |
---|
| 139 | if (i_overlap_scheme == IOverlapExponentialRandom) then |
---|
| 140 | call cum_cloud_cover_exp_ran(nlev, frac, overlap_param, & |
---|
[4489] | 141 | & cum_cloud_cover, pair_cloud_cover, use_beta_overlap) |
---|
[3908] | 142 | else if (i_overlap_scheme == IOverlapMaximumRandom) then |
---|
| 143 | call cum_cloud_cover_max_ran(nlev, frac, & |
---|
| 144 | & cum_cloud_cover, pair_cloud_cover) |
---|
| 145 | else if (i_overlap_scheme == IOverlapExponential) then |
---|
| 146 | call cum_cloud_cover_exp_exp(nlev, frac, overlap_param, & |
---|
[4489] | 147 | & cum_cloud_cover, pair_cloud_cover, use_beta_overlap) |
---|
[3908] | 148 | else |
---|
| 149 | write(nulerr,'(a)') '*** Error: cloud overlap scheme not recognised' |
---|
| 150 | call radiation_abort() |
---|
| 151 | end if |
---|
| 152 | |
---|
| 153 | total_cloud_cover = cum_cloud_cover(nlev); |
---|
| 154 | do jlev = 1,nlev-1 |
---|
| 155 | overhang(jlev) = cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev) |
---|
| 156 | end do |
---|
| 157 | |
---|
| 158 | if (total_cloud_cover < frac_threshold) then |
---|
| 159 | ! Treat column as clear sky: calling function therefore will not |
---|
| 160 | ! use od_scaling so we don't need to calculate it |
---|
| 161 | total_cloud_cover = 0.0_jprb |
---|
| 162 | |
---|
| 163 | else |
---|
| 164 | ! Cloud is present: need to calculate od_scaling |
---|
| 165 | |
---|
| 166 | ! Find range of cloudy layers |
---|
| 167 | jlev = 1 |
---|
| 168 | do while (frac(jlev) <= 0.0_jprb) |
---|
| 169 | jlev = jlev + 1 |
---|
| 170 | end do |
---|
| 171 | ibegin = jlev |
---|
| 172 | iend = jlev |
---|
| 173 | do jlev = jlev+1,nlev |
---|
| 174 | if (frac(jlev) > 0.0_jprb) then |
---|
| 175 | iend = jlev |
---|
| 176 | end if |
---|
| 177 | end do |
---|
| 178 | |
---|
| 179 | ! Set overlap parameter of inhomogeneities |
---|
| 180 | overlap_param_inhom = overlap_param |
---|
| 181 | |
---|
| 182 | do jlev = ibegin,iend-1 |
---|
| 183 | if (overlap_param(jlev) > 0.0_jprb) then |
---|
| 184 | overlap_param_inhom(jlev) & |
---|
| 185 | & = overlap_param(jlev)**(1.0_jprb/decorrelation_scaling) |
---|
| 186 | end if |
---|
| 187 | end do |
---|
| 188 | |
---|
| 189 | ! Reset optical depth scaling to clear skies |
---|
| 190 | od_scaling = 0.0_jprb |
---|
| 191 | |
---|
[4489] | 192 | if (present(use_vectorizable_generator)) then |
---|
| 193 | use_vec_gen = use_vectorizable_generator |
---|
| 194 | else |
---|
| 195 | use_vec_gen = .false. |
---|
| 196 | end if |
---|
[3908] | 197 | |
---|
[4489] | 198 | if (.not. use_vec_gen) then |
---|
| 199 | ! Original generator that minimizes the number of random |
---|
| 200 | ! numbers used, but is not vectorizable |
---|
[3908] | 201 | |
---|
[4489] | 202 | ! Expensive operation: initialize random number generator for |
---|
| 203 | ! this column |
---|
| 204 | call initialize_random_numbers(iseed, random_stream) |
---|
| 205 | |
---|
| 206 | ! Compute ng random numbers to use to locate cloud top |
---|
| 207 | call uniform_distribution(rand_top, random_stream) |
---|
| 208 | |
---|
| 209 | ! Loop over ng columns |
---|
| 210 | do jg = 1,ng |
---|
| 211 | ! Find the cloud top height corresponding to the current |
---|
| 212 | ! random number, and store in itrigger |
---|
| 213 | trigger = rand_top(jg) * total_cloud_cover |
---|
| 214 | jlev = ibegin |
---|
| 215 | do while (trigger > cum_cloud_cover(jlev) .and. jlev < iend) |
---|
| 216 | jlev = jlev + 1 |
---|
| 217 | end do |
---|
| 218 | itrigger = jlev |
---|
| 219 | |
---|
| 220 | if (i_overlap_scheme /= IOverlapExponential) then |
---|
| 221 | call generate_column_exp_ran(ng, nlev, jg, random_stream, pdf_sampler, & |
---|
| 222 | & frac, pair_cloud_cover, & |
---|
| 223 | & cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, & |
---|
| 224 | & itrigger, iend, od_scaling) |
---|
| 225 | else |
---|
| 226 | call generate_column_exp_exp(ng, nlev, jg, random_stream, pdf_sampler, & |
---|
| 227 | & frac, pair_cloud_cover, & |
---|
| 228 | & cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, & |
---|
| 229 | & itrigger, iend, od_scaling) |
---|
| 230 | end if |
---|
| 231 | |
---|
[3908] | 232 | end do |
---|
| 233 | |
---|
[4489] | 234 | else |
---|
| 235 | ! Alternative generator (only for Exp-Ran overlap so far) that |
---|
| 236 | ! should be vectorizable but generates more random numbers, |
---|
| 237 | ! some of which are not used |
---|
| 238 | |
---|
| 239 | if (i_overlap_scheme == IOverlapExponential) then |
---|
| 240 | write(nulerr,'(a)') '*** Error: vectorizable cloud generator is not available with Exp-Exp overlap' |
---|
| 241 | call radiation_abort() |
---|
[3908] | 242 | end if |
---|
| 243 | |
---|
[4489] | 244 | call generate_columns_exp_ran(ng, nlev, iseed, pdf_sampler, & |
---|
| 245 | & total_cloud_cover, frac_threshold, frac, pair_cloud_cover, & |
---|
| 246 | & cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, & |
---|
| 247 | & ibegin, iend, od_scaling) |
---|
| 248 | |
---|
| 249 | end if |
---|
| 250 | |
---|
[3908] | 251 | end if |
---|
| 252 | |
---|
| 253 | if (lhook) call dr_hook('radiation_cloud_generator:cloud_generator',1,hook_handle) |
---|
| 254 | |
---|
| 255 | end subroutine cloud_generator |
---|
| 256 | |
---|
| 257 | |
---|
| 258 | !--------------------------------------------------------------------- |
---|
| 259 | ! Generate a column of optical depth scalings using |
---|
| 260 | ! exponential-random overlap (which includes maximum-random overlap |
---|
| 261 | ! as a limiting case) |
---|
| 262 | subroutine generate_column_exp_ran(ng, nlev, ig, random_stream, pdf_sampler, & |
---|
| 263 | & frac, pair_cloud_cover, & |
---|
| 264 | & cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, & |
---|
| 265 | & itrigger, iend, od_scaling) |
---|
| 266 | |
---|
| 267 | use parkind1, only : jprb |
---|
| 268 | use radiation_pdf_sampler, only : pdf_sampler_type |
---|
| 269 | use random_numbers_mix, only : randomnumberstream, & |
---|
| 270 | initialize_random_numbers, uniform_distribution |
---|
| 271 | |
---|
| 272 | |
---|
| 273 | implicit none |
---|
| 274 | |
---|
| 275 | ! Number of g points / columns, and number of current column |
---|
| 276 | integer, intent(in) :: ng, ig |
---|
| 277 | |
---|
| 278 | ! Number of levels |
---|
| 279 | integer, intent(in) :: nlev |
---|
| 280 | |
---|
| 281 | ! Stream for producing random numbers |
---|
| 282 | type(randomnumberstream), intent(inout) :: random_stream |
---|
| 283 | |
---|
| 284 | ! Object for sampling from a lognormal or gamma distribution |
---|
| 285 | type(pdf_sampler_type), intent(in) :: pdf_sampler |
---|
| 286 | |
---|
| 287 | ! Cloud fraction, cumulative cloud cover and fractional standard |
---|
| 288 | ! deviation in each layer |
---|
| 289 | real(jprb), intent(in), dimension(nlev) :: frac, cum_cloud_cover, fractional_std |
---|
| 290 | |
---|
| 291 | ! Cloud cover of a pair of layers, and amount by which cloud at |
---|
| 292 | ! next level increases total cloud cover as seen from above |
---|
| 293 | real(jprb), intent(in), dimension(nlev-1) :: pair_cloud_cover, overhang |
---|
| 294 | |
---|
| 295 | ! Overlap parameter of inhomogeneities |
---|
| 296 | real(jprb), intent(in), dimension(nlev-1) :: overlap_param_inhom |
---|
| 297 | |
---|
| 298 | ! Top of highest cloudy layer (in this subcolumn) and base of |
---|
| 299 | ! lowest |
---|
| 300 | integer, intent(in) :: itrigger, iend |
---|
| 301 | |
---|
| 302 | ! Optical depth scaling to output |
---|
| 303 | real(jprb), intent(inout), dimension(ng,nlev) :: od_scaling |
---|
| 304 | |
---|
| 305 | ! Height indices |
---|
| 306 | integer :: jlev, jcloud |
---|
| 307 | |
---|
| 308 | ! Number of contiguous cloudy layers for which to compute optical |
---|
| 309 | ! depth scaling |
---|
| 310 | integer :: n_layers_to_scale |
---|
| 311 | |
---|
| 312 | integer :: iy |
---|
| 313 | |
---|
| 314 | ! Is it time to fill the od_scaling variable? |
---|
| 315 | logical :: do_fill_od_scaling |
---|
| 316 | |
---|
| 317 | real(jprb) :: rand_cloud(nlev) |
---|
| 318 | real(jprb) :: rand_inhom1(nlev), rand_inhom2(nlev) |
---|
| 319 | |
---|
| 320 | ! So far our vertically contiguous cloud contains only one layer |
---|
| 321 | n_layers_to_scale = 1 |
---|
| 322 | iy = 0 |
---|
| 323 | |
---|
| 324 | ! Locate the clouds below this layer: first generate some more |
---|
| 325 | ! random numbers |
---|
| 326 | call uniform_distribution(rand_cloud(1:(iend+1-itrigger)),random_stream) |
---|
| 327 | |
---|
| 328 | ! Loop from the layer below the local cloud top down to the |
---|
| 329 | ! bottom-most cloudy layer |
---|
| 330 | do jlev = itrigger+1,iend+1 |
---|
| 331 | do_fill_od_scaling = .false. |
---|
| 332 | if (jlev <= iend) then |
---|
| 333 | iy = iy+1 |
---|
| 334 | if (n_layers_to_scale > 0) then |
---|
| 335 | ! There is a cloud above, in which case the probability |
---|
| 336 | ! of cloud in the layer below is as follows |
---|
| 337 | if (rand_cloud(iy)*frac(jlev-1) & |
---|
| 338 | & < frac(jlev) + frac(jlev-1) - pair_cloud_cover(jlev-1)) then |
---|
| 339 | ! Add another cloudy layer |
---|
| 340 | n_layers_to_scale = n_layers_to_scale + 1 |
---|
| 341 | else |
---|
| 342 | ! Reached the end of a contiguous set of cloudy layers and |
---|
| 343 | ! will compute the optical depth scaling immediately. |
---|
| 344 | do_fill_od_scaling = .true. |
---|
| 345 | end if |
---|
| 346 | else |
---|
| 347 | ! There is clear-sky above, in which case the |
---|
| 348 | ! probability of cloud in the layer below is as follows |
---|
| 349 | if (rand_cloud(iy)*(cum_cloud_cover(jlev-1) - frac(jlev-1)) & |
---|
| 350 | & < pair_cloud_cover(jlev-1) - overhang(jlev-1) - frac(jlev-1)) then |
---|
| 351 | ! A new cloud top |
---|
| 352 | n_layers_to_scale = 1 |
---|
| 353 | end if |
---|
| 354 | end if |
---|
| 355 | else |
---|
| 356 | ! We are at the bottom of the cloudy layers in the model, |
---|
| 357 | ! so in a moment need to populate the od_scaling array |
---|
| 358 | do_fill_od_scaling = .true. |
---|
| 359 | end if |
---|
| 360 | |
---|
| 361 | if (do_fill_od_scaling) then |
---|
| 362 | ! We have a contiguous range of layers for which we |
---|
| 363 | ! compute the od_scaling elements using some random |
---|
| 364 | ! numbers |
---|
| 365 | call uniform_distribution(rand_inhom1(1:n_layers_to_scale),random_stream) |
---|
| 366 | call uniform_distribution(rand_inhom2(1:n_layers_to_scale),random_stream) |
---|
| 367 | |
---|
| 368 | ! Loop through the sequence of cloudy layers |
---|
| 369 | do jcloud = 2,n_layers_to_scale |
---|
| 370 | ! Use second random number, and inhomogeneity overlap |
---|
| 371 | ! parameter, to decide whether the first random number |
---|
| 372 | ! should be repeated (corresponding to maximum overlap) |
---|
| 373 | ! or not (corresponding to random overlap) |
---|
| 374 | if (rand_inhom2(jcloud) & |
---|
| 375 | & < overlap_param_inhom(jlev-n_layers_to_scale+jcloud-2)) then |
---|
| 376 | rand_inhom1(jcloud) = rand_inhom1(jcloud-1) |
---|
| 377 | end if |
---|
| 378 | end do |
---|
| 379 | |
---|
| 380 | ! Sample from a lognormal or gamma distribution to obtain |
---|
| 381 | ! the optical depth scalings |
---|
| 382 | call pdf_sampler%sample(fractional_std(jlev-n_layers_to_scale:jlev-1), & |
---|
| 383 | & rand_inhom1(1:n_layers_to_scale), od_scaling(ig,jlev-n_layers_to_scale:jlev-1)) |
---|
| 384 | |
---|
| 385 | n_layers_to_scale = 0 |
---|
| 386 | end if |
---|
| 387 | |
---|
| 388 | end do |
---|
| 389 | |
---|
| 390 | end subroutine generate_column_exp_ran |
---|
| 391 | |
---|
| 392 | |
---|
| 393 | |
---|
| 394 | !--------------------------------------------------------------------- |
---|
| 395 | ! Generate a column of optical depth scalings using |
---|
| 396 | ! exponential-exponential overlap |
---|
| 397 | subroutine generate_column_exp_exp(ng, nlev, ig, random_stream, pdf_sampler, & |
---|
| 398 | & frac, pair_cloud_cover, & |
---|
| 399 | & cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, & |
---|
| 400 | & itrigger, iend, od_scaling) |
---|
| 401 | |
---|
| 402 | use parkind1, only : jprb |
---|
| 403 | use radiation_pdf_sampler, only : pdf_sampler_type |
---|
| 404 | use random_numbers_mix, only : randomnumberstream, & |
---|
| 405 | initialize_random_numbers, uniform_distribution |
---|
| 406 | |
---|
| 407 | implicit none |
---|
| 408 | |
---|
| 409 | ! Number of g points / columns, and number of current column |
---|
| 410 | integer, intent(in) :: ng, ig |
---|
| 411 | |
---|
| 412 | ! Number of levels |
---|
| 413 | integer, intent(in) :: nlev |
---|
| 414 | |
---|
| 415 | ! Stream for producing random numbers |
---|
| 416 | type(randomnumberstream), intent(inout) :: random_stream |
---|
| 417 | |
---|
| 418 | ! Object for sampling from a lognormal or gamma distribution |
---|
| 419 | type(pdf_sampler_type), intent(in) :: pdf_sampler |
---|
| 420 | |
---|
| 421 | ! Cloud fraction, cumulative cloud cover and fractional standard |
---|
| 422 | ! deviation in each layer |
---|
| 423 | real(jprb), intent(in), dimension(nlev) :: frac, cum_cloud_cover, fractional_std |
---|
| 424 | |
---|
| 425 | ! Cloud cover of a pair of layers, and amount by which cloud at |
---|
| 426 | ! next level increases total cloud cover as seen from above |
---|
| 427 | real(jprb), intent(in), dimension(nlev-1) :: pair_cloud_cover, overhang |
---|
| 428 | |
---|
| 429 | ! Overlap parameter of inhomogeneities |
---|
| 430 | real(jprb), intent(in), dimension(nlev-1) :: overlap_param_inhom |
---|
| 431 | |
---|
| 432 | ! Top of highest cloudy layer (in this subcolumn) and base of |
---|
| 433 | ! lowest |
---|
| 434 | integer, intent(in) :: itrigger, iend |
---|
| 435 | |
---|
| 436 | ! Optical depth scaling to output |
---|
| 437 | real(jprb), intent(inout), dimension(ng,nlev) :: od_scaling |
---|
| 438 | |
---|
| 439 | ! Height indices |
---|
| 440 | integer :: jlev, jcloud |
---|
| 441 | |
---|
| 442 | integer :: iy |
---|
| 443 | |
---|
| 444 | real(jprb) :: rand_cloud(nlev) |
---|
| 445 | real(jprb) :: rand_inhom1(nlev), rand_inhom2(nlev) |
---|
| 446 | |
---|
| 447 | ! For each column analysed, this vector locates the clouds. It is |
---|
| 448 | ! only actually used for Exp-Exp overlap |
---|
| 449 | logical :: is_cloudy(nlev) |
---|
| 450 | |
---|
| 451 | ! Number of contiguous cloudy layers for which to compute optical |
---|
| 452 | ! depth scaling |
---|
| 453 | integer :: n_layers_to_scale |
---|
| 454 | |
---|
| 455 | iy = 0 |
---|
| 456 | |
---|
| 457 | is_cloudy = .false. |
---|
| 458 | is_cloudy(itrigger) = .true. |
---|
| 459 | |
---|
| 460 | ! Locate the clouds below this layer: first generate some more |
---|
| 461 | ! random numbers |
---|
| 462 | call uniform_distribution(rand_cloud(1:(iend+1-itrigger)),random_stream) |
---|
| 463 | |
---|
| 464 | ! Loop from the layer below the local cloud top down to the |
---|
| 465 | ! bottom-most cloudy layer |
---|
| 466 | do jlev = itrigger+1,iend |
---|
| 467 | iy = iy+1 |
---|
| 468 | if (is_cloudy(jlev-1)) then |
---|
| 469 | ! There is a cloud above, in which case the probability |
---|
| 470 | ! of cloud in the layer below is as follows |
---|
| 471 | if (rand_cloud(iy)*frac(jlev-1) & |
---|
| 472 | & < frac(jlev) + frac(jlev-1) - pair_cloud_cover(jlev-1)) then |
---|
| 473 | ! Add another cloudy layer |
---|
| 474 | is_cloudy(jlev) = .true. |
---|
| 475 | end if |
---|
| 476 | else |
---|
| 477 | ! There is clear-sky above, in which case the |
---|
| 478 | ! probability of cloud in the layer below is as follows |
---|
| 479 | if (rand_cloud(iy)*(cum_cloud_cover(jlev-1) - frac(jlev-1)) & |
---|
| 480 | & < pair_cloud_cover(jlev-1) - overhang(jlev-1) - frac(jlev-1)) then |
---|
| 481 | ! A new cloud top |
---|
| 482 | is_cloudy(jlev) = .true. |
---|
| 483 | end if |
---|
| 484 | end if |
---|
| 485 | end do |
---|
| 486 | |
---|
| 487 | ! We have a contiguous range of layers for which we compute the |
---|
| 488 | ! od_scaling elements using some random numbers |
---|
| 489 | |
---|
| 490 | ! In the Exp-Exp overlap scheme we do all layers at once |
---|
| 491 | n_layers_to_scale = iend+1 - itrigger |
---|
| 492 | |
---|
| 493 | call uniform_distribution(rand_inhom1(1:n_layers_to_scale),random_stream) |
---|
| 494 | call uniform_distribution(rand_inhom2(1:n_layers_to_scale),random_stream) |
---|
| 495 | |
---|
| 496 | ! Loop through the sequence of cloudy layers |
---|
| 497 | do jcloud = 2,n_layers_to_scale |
---|
| 498 | ! Use second random number, and inhomogeneity overlap |
---|
| 499 | ! parameter, to decide whether the first random number |
---|
| 500 | ! should be repeated (corresponding to maximum overlap) |
---|
| 501 | ! or not (corresponding to random overlap) |
---|
| 502 | if (rand_inhom2(jcloud) & |
---|
| 503 | & < overlap_param_inhom(iend-n_layers_to_scale+jcloud-1)) then |
---|
| 504 | rand_inhom1(jcloud) = rand_inhom1(jcloud-1) |
---|
| 505 | end if |
---|
| 506 | end do |
---|
| 507 | |
---|
| 508 | ! Sample from a lognormal or gamma distribution to obtain the |
---|
[4489] | 509 | ! optical depth scalings |
---|
| 510 | |
---|
| 511 | ! Masked version assuming values outside the range itrigger:iend |
---|
| 512 | ! are already zero: |
---|
[3908] | 513 | call pdf_sampler%masked_sample(n_layers_to_scale, & |
---|
| 514 | & fractional_std(itrigger:iend), & |
---|
| 515 | & rand_inhom1(1:n_layers_to_scale), od_scaling(ig,itrigger:iend), & |
---|
| 516 | & is_cloudy(itrigger:iend)) |
---|
| 517 | |
---|
[4489] | 518 | ! ! IFS version: |
---|
| 519 | ! !$omp simd |
---|
| 520 | ! do jlev=itrigger,iend |
---|
| 521 | ! if (.not. is_cloudy(jlev)) then |
---|
| 522 | ! od_scaling(ig,jlev) = 0.0_jprb |
---|
| 523 | ! else |
---|
| 524 | ! call sample_from_pdf_simd(& |
---|
| 525 | ! pdf_sampler,fractional_std(jlev),& |
---|
| 526 | ! rand_inhom1(jlev-itrigger+1), & |
---|
| 527 | ! od_scaling(ig,jlev)) |
---|
| 528 | ! end if |
---|
| 529 | ! end do |
---|
| 530 | |
---|
[3908] | 531 | end subroutine generate_column_exp_exp |
---|
| 532 | |
---|
[4489] | 533 | |
---|
| 534 | !--------------------------------------------------------------------- |
---|
| 535 | ! Extract the value of a lognormal distribution with fractional |
---|
| 536 | ! standard deviation "fsd" corresponding to the cumulative |
---|
| 537 | ! distribution function value "cdf", and return it in x. Since this |
---|
| 538 | ! is an elemental subroutine, fsd, cdf and x may be arrays. SIMD version. |
---|
| 539 | subroutine sample_from_pdf_simd(this, fsd, cdf, x) |
---|
| 540 | use parkind1, only : jprb |
---|
| 541 | use radiation_pdf_sampler, only : pdf_sampler_type |
---|
| 542 | implicit none |
---|
| 543 | type(pdf_sampler_type), intent(in) :: this |
---|
| 544 | |
---|
| 545 | ! Fractional standard deviation (0 to 4) and cumulative |
---|
| 546 | ! distribution function (0 to 1) |
---|
| 547 | real(jprb), intent(in) :: fsd, cdf |
---|
| 548 | |
---|
| 549 | ! Sample from distribution |
---|
| 550 | real(jprb), intent(out) :: x |
---|
| 551 | |
---|
| 552 | ! Index to look-up table |
---|
| 553 | integer :: ifsd, icdf |
---|
| 554 | |
---|
| 555 | ! Weights in bilinear interpolation |
---|
| 556 | real(jprb) :: wfsd, wcdf |
---|
| 557 | |
---|
| 558 | ! Bilinear interpolation with bounds |
---|
| 559 | wcdf = cdf * (this%ncdf-1) + 1.0_jprb |
---|
| 560 | icdf = max(1, min(int(wcdf), this%ncdf-1)) |
---|
| 561 | wcdf = max(0.0_jprb, min(wcdf - icdf, 1.0_jprb)) |
---|
| 562 | |
---|
| 563 | wfsd = (fsd-this%fsd1) * this%inv_fsd_interval + 1.0_jprb |
---|
| 564 | ifsd = max(1, min(int(wfsd), this%nfsd-1)) |
---|
| 565 | wfsd = max(0.0_jprb, min(wfsd - ifsd, 1.0_jprb)) |
---|
| 566 | |
---|
| 567 | x = (1.0_jprb-wcdf)*(1.0_jprb-wfsd) * this%val(icdf ,ifsd) & |
---|
| 568 | & + (1.0_jprb-wcdf)* wfsd * this%val(icdf ,ifsd+1) & |
---|
| 569 | & + wcdf *(1.0_jprb-wfsd) * this%val(icdf+1,ifsd) & |
---|
| 570 | & + wcdf * wfsd * this%val(icdf+1,ifsd+1) |
---|
| 571 | |
---|
| 572 | end subroutine sample_from_pdf_simd |
---|
| 573 | |
---|
| 574 | |
---|
| 575 | !--------------------------------------------------------------------- |
---|
| 576 | ! Generate columns of optical depth scalings using |
---|
| 577 | ! exponential-random overlap (which includes maximum-random overlap |
---|
| 578 | ! as a limiting case). This version is intended to work better on |
---|
| 579 | ! hardware with long vector lengths. As with all calculations in |
---|
| 580 | ! this file, we zoom into the fraction of the column with cloud at |
---|
| 581 | ! any height, so that all spectral intervals see a cloud somewhere. |
---|
| 582 | ! In the McICA solver, this is combined appropriately with the |
---|
| 583 | ! clear-sky calculation. |
---|
| 584 | subroutine generate_columns_exp_ran(ng, nlev, iseed, pdf_sampler, & |
---|
| 585 | & total_cloud_cover, frac_threshold, frac, pair_cloud_cover, & |
---|
| 586 | & cum_cloud_cover, overhang, fractional_std, overlap_param_inhom, & |
---|
| 587 | & ibegin, iend, od_scaling) |
---|
| 588 | |
---|
| 589 | use parkind1, only : jprb |
---|
| 590 | use radiation_pdf_sampler, only : pdf_sampler_type |
---|
| 591 | use radiation_random_numbers, only : rng_type, IRngMinstdVector, IRngNative |
---|
| 592 | |
---|
| 593 | implicit none |
---|
| 594 | |
---|
| 595 | ! Number of g points / columns |
---|
| 596 | integer, intent(in) :: ng |
---|
| 597 | |
---|
| 598 | ! Number of levels |
---|
| 599 | integer, intent(in) :: nlev |
---|
| 600 | |
---|
| 601 | integer, intent(in) :: iseed ! seed for random number generator |
---|
| 602 | |
---|
| 603 | ! Stream for producing random numbers |
---|
| 604 | !type(randomnumberstream) :: random_stream |
---|
| 605 | type(rng_type) :: random_number_generator |
---|
| 606 | |
---|
| 607 | ! Object for sampling from a lognormal or gamma distribution |
---|
| 608 | type(pdf_sampler_type), intent(in) :: pdf_sampler |
---|
| 609 | |
---|
| 610 | ! Total cloud cover using cloud fraction and overlap parameter |
---|
| 611 | real(jprb), intent(in) :: total_cloud_cover |
---|
| 612 | |
---|
| 613 | real(jprb), intent(in) :: frac_threshold |
---|
| 614 | |
---|
| 615 | ! Cloud fraction, cumulative cloud cover and fractional standard |
---|
| 616 | ! deviation in each layer |
---|
| 617 | real(jprb), intent(in), dimension(nlev) :: frac, cum_cloud_cover, fractional_std |
---|
| 618 | |
---|
| 619 | ! Cloud cover of a pair of layers, and amount by which cloud at |
---|
| 620 | ! next level increases total cloud cover as seen from above |
---|
| 621 | real(jprb), intent(in), dimension(nlev-1) :: pair_cloud_cover, overhang |
---|
| 622 | |
---|
| 623 | ! Overlap parameter of inhomogeneities |
---|
| 624 | real(jprb), intent(in), dimension(nlev-1) :: overlap_param_inhom |
---|
| 625 | |
---|
| 626 | ! Top of highest cloudy layer and base of lowest |
---|
| 627 | integer, intent(inout) :: ibegin, iend |
---|
| 628 | |
---|
| 629 | ! Optical depth scaling to output |
---|
| 630 | real(jprb), intent(inout), dimension(ng,nlev) :: od_scaling |
---|
| 631 | |
---|
| 632 | ! Loop indices |
---|
| 633 | integer :: jlev, jg |
---|
| 634 | |
---|
| 635 | real(jprb) :: rand_cloud(ng,ibegin:iend) |
---|
| 636 | real(jprb) :: rand_inhom(ng,ibegin-1:iend), rand_inhom2(ng,ibegin:iend) |
---|
| 637 | |
---|
| 638 | ! Is the cloud fraction above the minimum threshold at each level |
---|
| 639 | logical :: is_any_cloud(ibegin:iend) |
---|
| 640 | |
---|
| 641 | ! Scaled random number for finding cloud |
---|
| 642 | real(jprb) :: trigger(ng) |
---|
| 643 | |
---|
| 644 | logical :: is_cloud(ng) ! Is there cloud at this level and spectral interval? |
---|
| 645 | logical :: prev_cloud(ng) ! Was there cloud at level above? |
---|
| 646 | logical :: first_cloud(ng) ! At level of first cloud counting down from top? |
---|
| 647 | logical :: found_cloud(ng) ! Cloud found in this column counting down from top? |
---|
| 648 | |
---|
| 649 | is_any_cloud = (frac(ibegin:iend) >= frac_threshold) |
---|
| 650 | |
---|
| 651 | ! Initialize random number generator for this column, and state |
---|
| 652 | ! that random numbers will be requested in blocks of length the |
---|
| 653 | ! number of spectral intervals ng. |
---|
| 654 | call random_number_generator%initialize(IRngMinstdVector, iseed=iseed, & |
---|
| 655 | & nmaxstreams=ng) |
---|
| 656 | |
---|
| 657 | ! Random numbers to use to locate cloud top |
---|
| 658 | call random_number_generator%uniform_distribution(trigger) |
---|
| 659 | |
---|
| 660 | ! Random numbers to work out whether to transition vertically from |
---|
| 661 | ! clear to cloudy, cloudy to clear, clear to clear or cloudy to |
---|
| 662 | ! cloudy |
---|
| 663 | call random_number_generator%uniform_distribution(rand_cloud, is_any_cloud) |
---|
| 664 | |
---|
| 665 | ! Random numbers to generate sub-grid cloud structure |
---|
| 666 | call random_number_generator%uniform_distribution(rand_inhom) |
---|
| 667 | call random_number_generator%uniform_distribution(rand_inhom2, is_any_cloud) |
---|
| 668 | |
---|
| 669 | trigger = trigger * total_cloud_cover |
---|
| 670 | |
---|
| 671 | ! Initialize logicals for clear-sky above first cloudy layer |
---|
| 672 | found_cloud = .false. |
---|
| 673 | is_cloud = .false. |
---|
| 674 | first_cloud = .false. |
---|
| 675 | |
---|
| 676 | ! Loop down through layers starting at the first cloudy layer |
---|
| 677 | do jlev = ibegin,iend |
---|
| 678 | |
---|
| 679 | if (is_any_cloud(jlev)) then |
---|
| 680 | |
---|
| 681 | ! Added for DWD (2020) |
---|
| 682 | !NEC$ shortloop |
---|
| 683 | do jg = 1,ng |
---|
| 684 | ! The intention is that all these operations are vectorizable, |
---|
| 685 | ! since all are vector operations on vectors of length ng... |
---|
| 686 | |
---|
| 687 | ! Copy the cloud mask between levels |
---|
| 688 | prev_cloud(jg) = is_cloud(jg) |
---|
| 689 | |
---|
| 690 | ! For each spectral interval, has the first cloud appeared at this level? |
---|
| 691 | first_cloud(jg) = (trigger(jg) <= cum_cloud_cover(jlev) .and. .not. found_cloud(jg)) |
---|
| 692 | |
---|
| 693 | ! ...if so, add to found_cloud |
---|
| 694 | found_cloud(jg) = found_cloud(jg) .or. first_cloud(jg) |
---|
| 695 | |
---|
| 696 | ! There is cloud at this level either if a first cloud has |
---|
| 697 | ! appeared, or using separate probability calculations |
---|
| 698 | ! depending on whether there is a cloud above (given by |
---|
| 699 | ! prev_cloud) |
---|
| 700 | is_cloud(jg) = first_cloud(jg) & |
---|
| 701 | & .or. found_cloud(jg) .and. merge(rand_cloud(jg,jlev)*frac(jlev-1) & |
---|
| 702 | & < frac(jlev)+frac(jlev-1)-pair_cloud_cover(jlev-1), & |
---|
| 703 | & rand_cloud(jg,jlev)*(cum_cloud_cover(jlev-1) - frac(jlev-1)) & |
---|
| 704 | & < pair_cloud_cover(jlev-1) - overhang(jlev-1) - frac(jlev-1), & |
---|
| 705 | & prev_cloud(jg)) |
---|
| 706 | ! The random number determining cloud structure decorrelates |
---|
| 707 | ! with the one above it according to the overlap parameter, |
---|
| 708 | ! but always decorrelates if there is clear-sky above. If |
---|
| 709 | ! there is clear-sky in the present level, the random number |
---|
| 710 | ! is set to zero to ensure that the optical depth scaling is |
---|
| 711 | ! also zero. |
---|
| 712 | rand_inhom(jg,jlev) = merge(merge(rand_inhom(jg,jlev-1), rand_inhom(jg,jlev), & |
---|
| 713 | & rand_inhom2(jg,jlev) < overlap_param_inhom(jlev-1) & |
---|
| 714 | & .and. prev_cloud(jg)), & |
---|
| 715 | & 0.0_jprb, is_cloud(jg)) |
---|
| 716 | end do |
---|
| 717 | else |
---|
| 718 | ! No cloud at this level |
---|
| 719 | is_cloud = .false. |
---|
| 720 | end if |
---|
| 721 | end do |
---|
| 722 | |
---|
| 723 | ! Sample from a lognormal or gamma distribution to obtain the |
---|
| 724 | ! optical depth scalings, calling the faster masked version and |
---|
| 725 | ! assuming values outside the range ibegin:iend are already zero |
---|
| 726 | call pdf_sampler%masked_block_sample(iend-ibegin+1, ng, & |
---|
| 727 | & fractional_std(ibegin:iend), & |
---|
| 728 | & rand_inhom(:,ibegin:iend), od_scaling(:,ibegin:iend), & |
---|
| 729 | & is_any_cloud) |
---|
| 730 | |
---|
| 731 | end subroutine generate_columns_exp_ran |
---|
| 732 | |
---|
[3908] | 733 | end module radiation_cloud_generator |
---|