| 1 | ! radiation_cloud_cover.F90 - Compute cumulative cloud cover for McICA |
|---|
| 2 | ! |
|---|
| 3 | ! (C) Copyright 2016- 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 profiles of the cumulative cloud cover as seen from TOA, |
|---|
| 16 | ! used in the McICA cloud generator. |
|---|
| 17 | ! |
|---|
| 18 | ! Modifications |
|---|
| 19 | ! 2020-10-07 R. Hogan Ensure iobj1 initialized in case of alpha_obj==0 |
|---|
| 20 | |
|---|
| 21 | module radiation_cloud_cover |
|---|
| 22 | |
|---|
| 23 | use parkind1, only : jprb |
|---|
| 24 | |
|---|
| 25 | public |
|---|
| 26 | |
|---|
| 27 | ! Three overlap schemes. Note that "Exponential" means that |
|---|
| 28 | ! clear-sky regions have no special significance for computing the |
|---|
| 29 | ! cumulative cloud cover: non-contiguous clouds are exponentially |
|---|
| 30 | ! rather than randomly overlapped. This is the situaition in the |
|---|
| 31 | ! McRad radiation scheme at ECMWF. |
|---|
| 32 | enum, bind(c) |
|---|
| 33 | enumerator IOverlapMaximumRandom, IOverlapExponentialRandom, & |
|---|
| 34 | & IOverlapExponential |
|---|
| 35 | end enum |
|---|
| 36 | character(len=*), parameter :: OverlapName(0:2) = (/ 'Max-Ran', & |
|---|
| 37 | & 'Exp-Ran', & |
|---|
| 38 | & 'Exp-Exp' /) |
|---|
| 39 | |
|---|
| 40 | ! Maximum cloud fraction distinguishable from 1 |
|---|
| 41 | real(jprb), parameter :: MaxCloudFrac = 1.0_jprb-epsilon(1.0_jprb)*10.0_jprb |
|---|
| 42 | |
|---|
| 43 | |
|---|
| 44 | contains |
|---|
| 45 | |
|---|
| 46 | !--------------------------------------------------------------------- |
|---|
| 47 | ! Convert "beta" overlap parameter of Shonk et al. (2010) to "alpha" |
|---|
| 48 | ! overlap parameter of Hogan and Illingworth (2000) |
|---|
| 49 | elemental function beta2alpha(beta, frac1, frac2) |
|---|
| 50 | |
|---|
| 51 | implicit none |
|---|
| 52 | |
|---|
| 53 | ! Beta overlap parameter and the cloud fractions in the upper and |
|---|
| 54 | ! lower layers |
|---|
| 55 | real(jprb), intent(in) :: beta, frac1, frac2 |
|---|
| 56 | |
|---|
| 57 | real(jprb) :: beta2alpha |
|---|
| 58 | |
|---|
| 59 | ! Absolute difference in cloud fraction |
|---|
| 60 | real(jprb) :: frac_diff |
|---|
| 61 | |
|---|
| 62 | if (beta < 1.0_jprb) then |
|---|
| 63 | frac_diff = abs(frac1-frac2) |
|---|
| 64 | beta2alpha = beta & |
|---|
| 65 | & + (1.0_jprb-beta)*frac_diff & |
|---|
| 66 | & / (frac_diff + 1.0_jprb/beta - 1.0_jprb) |
|---|
| 67 | else |
|---|
| 68 | beta2alpha = 1.0_jprb |
|---|
| 69 | end if |
|---|
| 70 | |
|---|
| 71 | end function beta2alpha |
|---|
| 72 | |
|---|
| 73 | |
|---|
| 74 | !--------------------------------------------------------------------- |
|---|
| 75 | ! Compute total cloud cover according to the specified overlap |
|---|
| 76 | ! rule. This can be used to compute the high, mid and low cloud |
|---|
| 77 | ! cover by passing in subsets of the cloud fraction array |
|---|
| 78 | function cloud_cover(nlev, i_overlap_scheme, frac, overlap_param, & |
|---|
| 79 | & is_beta_overlap) |
|---|
| 80 | |
|---|
| 81 | implicit none |
|---|
| 82 | |
|---|
| 83 | ! Number of levels and the overlap scheme to be applied |
|---|
| 84 | integer, intent(in) :: nlev, i_overlap_scheme |
|---|
| 85 | |
|---|
| 86 | ! Cloud fraction and the overlap parameter between adjacent pairs |
|---|
| 87 | ! of levels |
|---|
| 88 | real(jprb), intent(in) :: frac(nlev), overlap_param(nlev-1) |
|---|
| 89 | |
|---|
| 90 | ! Do we use the "beta" overlap scheme of Shonk et al. (2010)? |
|---|
| 91 | ! Default is false. |
|---|
| 92 | logical, intent(in), optional :: is_beta_overlap |
|---|
| 93 | |
|---|
| 94 | ! Return cloud cover |
|---|
| 95 | real(jprb) :: cloud_cover |
|---|
| 96 | |
|---|
| 97 | ! Cumulative cloud cover from TOA to the base of each layer |
|---|
| 98 | real(jprb) :: cum_cloud_cover(nlev) |
|---|
| 99 | |
|---|
| 100 | ! Cloud cover of a pair of layers |
|---|
| 101 | real(jprb) :: pair_cloud_cover(nlev-1) |
|---|
| 102 | |
|---|
| 103 | if (i_overlap_scheme == IOverlapExponentialRandom) then |
|---|
| 104 | call cum_cloud_cover_exp_ran(nlev, frac, overlap_param, & |
|---|
| 105 | & cum_cloud_cover, pair_cloud_cover, is_beta_overlap) |
|---|
| 106 | else if (i_overlap_scheme == IOverlapExponential) then |
|---|
| 107 | call cum_cloud_cover_exp_exp(nlev, frac, overlap_param, & |
|---|
| 108 | & cum_cloud_cover, pair_cloud_cover, is_beta_overlap) |
|---|
| 109 | else |
|---|
| 110 | call cum_cloud_cover_max_ran(nlev, frac, cum_cloud_cover, & |
|---|
| 111 | & pair_cloud_cover) |
|---|
| 112 | end if |
|---|
| 113 | |
|---|
| 114 | cloud_cover = cum_cloud_cover(nlev) |
|---|
| 115 | |
|---|
| 116 | end function cloud_cover |
|---|
| 117 | |
|---|
| 118 | |
|---|
| 119 | !--------------------------------------------------------------------- |
|---|
| 120 | ! Maximum-random overlap: Geleyn & Hollingsworth formula |
|---|
| 121 | subroutine cum_cloud_cover_max_ran(nlev, frac, & |
|---|
| 122 | & cum_cloud_cover, pair_cloud_cover) |
|---|
| 123 | |
|---|
| 124 | use yomhook, only : lhook, dr_hook |
|---|
| 125 | |
|---|
| 126 | implicit none |
|---|
| 127 | |
|---|
| 128 | ! Inputs |
|---|
| 129 | integer, intent(in) :: nlev ! number of model levels |
|---|
| 130 | |
|---|
| 131 | ! Cloud fraction on full levels |
|---|
| 132 | real(jprb), intent(in) :: frac(nlev) |
|---|
| 133 | |
|---|
| 134 | ! Outputs |
|---|
| 135 | |
|---|
| 136 | ! Cumulative cloud cover from TOA to the base of each layer |
|---|
| 137 | real(jprb), intent(out) :: cum_cloud_cover(nlev) |
|---|
| 138 | |
|---|
| 139 | ! Cloud cover of a pair of layers |
|---|
| 140 | real(jprb), intent(out) :: pair_cloud_cover(nlev-1) |
|---|
| 141 | |
|---|
| 142 | ! Local variables |
|---|
| 143 | |
|---|
| 144 | ! Cumulative product needed in computation of total_cloud_cover |
|---|
| 145 | real(jprb) :: cum_product |
|---|
| 146 | |
|---|
| 147 | ! Loop index for model level |
|---|
| 148 | integer :: jlev |
|---|
| 149 | |
|---|
| 150 | real(jprb) :: hook_handle |
|---|
| 151 | |
|---|
| 152 | if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_max_ran',0,hook_handle) |
|---|
| 153 | |
|---|
| 154 | ! Loop to compute total cloud cover and the cumulative cloud cover |
|---|
| 155 | |
|---|
| 156 | ! down to the base of each layer |
|---|
| 157 | cum_product = 1.0_jprb - frac(1) |
|---|
| 158 | cum_cloud_cover(1) = frac(1) |
|---|
| 159 | do jlev = 1,nlev-1 |
|---|
| 160 | ! Compute the combined cloud cover of layers jlev and jlev+1 |
|---|
| 161 | pair_cloud_cover(jlev) = max(frac(jlev),frac(jlev+1)) |
|---|
| 162 | |
|---|
| 163 | if (frac(jlev) >= MaxCloudFrac) then |
|---|
| 164 | ! Cloud cover has reached one |
|---|
| 165 | cum_product = 0.0_jprb |
|---|
| 166 | else |
|---|
| 167 | cum_product = cum_product * (1.0_jprb - pair_cloud_cover(jlev)) & |
|---|
| 168 | & / (1.0_jprb - frac(jlev)) |
|---|
| 169 | end if |
|---|
| 170 | cum_cloud_cover(jlev+1) = 1.0_jprb - cum_product |
|---|
| 171 | end do |
|---|
| 172 | |
|---|
| 173 | if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_max_ran',1,hook_handle) |
|---|
| 174 | |
|---|
| 175 | end subroutine cum_cloud_cover_max_ran |
|---|
| 176 | |
|---|
| 177 | |
|---|
| 178 | !--------------------------------------------------------------------- |
|---|
| 179 | ! Exponential-random overlap: exponential overlap for contiguous |
|---|
| 180 | ! clouds, random overlap for non-contiguous clouds |
|---|
| 181 | subroutine cum_cloud_cover_exp_ran(nlev, frac, overlap_param, & |
|---|
| 182 | & cum_cloud_cover, pair_cloud_cover, is_beta_overlap) |
|---|
| 183 | |
|---|
| 184 | use yomhook, only : lhook, dr_hook |
|---|
| 185 | |
|---|
| 186 | implicit none |
|---|
| 187 | |
|---|
| 188 | ! Inputs |
|---|
| 189 | integer, intent(in) :: nlev ! number of model levels |
|---|
| 190 | |
|---|
| 191 | ! Cloud fraction on full levels |
|---|
| 192 | real(jprb), intent(in) :: frac(nlev) |
|---|
| 193 | |
|---|
| 194 | ! Cloud overlap parameter for interfaces between model layers, |
|---|
| 195 | ! where 0 indicates random overlap and 1 indicates maximum-random |
|---|
| 196 | ! overlap |
|---|
| 197 | real(jprb), intent(in) :: overlap_param(nlev-1) |
|---|
| 198 | |
|---|
| 199 | ! This routine has been coded using the "alpha" overlap parameter |
|---|
| 200 | ! of Hogan and Illingworth (2000). If the following logical is |
|---|
| 201 | ! present and true then the input is interpretted to be the "beta" |
|---|
| 202 | ! overlap parameter of Shonk et al. (2010), and needs to be |
|---|
| 203 | ! converted to alpha. |
|---|
| 204 | logical, intent(in), optional :: is_beta_overlap |
|---|
| 205 | |
|---|
| 206 | ! Outputs |
|---|
| 207 | |
|---|
| 208 | ! Cumulative cloud cover from TOA to the base of each layer |
|---|
| 209 | real(jprb), intent(out) :: cum_cloud_cover(nlev) |
|---|
| 210 | |
|---|
| 211 | ! Cloud cover of a pair of layers |
|---|
| 212 | real(jprb), intent(out) :: pair_cloud_cover(nlev-1) |
|---|
| 213 | |
|---|
| 214 | ! Local variables |
|---|
| 215 | |
|---|
| 216 | ! Cumulative product needed in computation of total_cloud_cover |
|---|
| 217 | real(jprb) :: cum_product |
|---|
| 218 | |
|---|
| 219 | ! "Alpha" overlap parameter |
|---|
| 220 | real(jprb) :: overlap_alpha |
|---|
| 221 | logical :: do_overlap_conversion |
|---|
| 222 | |
|---|
| 223 | ! Loop index for model level |
|---|
| 224 | integer :: jlev |
|---|
| 225 | |
|---|
| 226 | real(jprb) :: hook_handle |
|---|
| 227 | |
|---|
| 228 | if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_exp_ran',0,hook_handle) |
|---|
| 229 | |
|---|
| 230 | |
|---|
| 231 | if (present(is_beta_overlap)) then |
|---|
| 232 | do_overlap_conversion = is_beta_overlap |
|---|
| 233 | else |
|---|
| 234 | do_overlap_conversion = .false. |
|---|
| 235 | end if |
|---|
| 236 | |
|---|
| 237 | ! Loop to compute total cloud cover and the cumulative cloud cover |
|---|
| 238 | |
|---|
| 239 | ! down to the base of each layer |
|---|
| 240 | cum_product = 1.0_jprb - frac(1) |
|---|
| 241 | cum_cloud_cover(1) = frac(1) |
|---|
| 242 | do jlev = 1,nlev-1 |
|---|
| 243 | ! Convert to "alpha" overlap parameter if necessary |
|---|
| 244 | if (do_overlap_conversion) then |
|---|
| 245 | overlap_alpha = beta2alpha(overlap_param(jlev), & |
|---|
| 246 | & frac(jlev), frac(jlev+1)) |
|---|
| 247 | else |
|---|
| 248 | overlap_alpha = overlap_param(jlev) |
|---|
| 249 | end if |
|---|
| 250 | |
|---|
| 251 | ! Compute the combined cloud cover of layers jlev and jlev+1 |
|---|
| 252 | pair_cloud_cover(jlev) = overlap_alpha*max(frac(jlev),frac(jlev+1)) & |
|---|
| 253 | & + (1.0_jprb - overlap_alpha) & |
|---|
| 254 | & * (frac(jlev)+frac(jlev+1)-frac(jlev)*frac(jlev+1)) |
|---|
| 255 | |
|---|
| 256 | if (frac(jlev) >= MaxCloudFrac) then |
|---|
| 257 | ! Cloud cover has reached one |
|---|
| 258 | cum_product = 0.0_jprb |
|---|
| 259 | else |
|---|
| 260 | cum_product = cum_product * (1.0_jprb - pair_cloud_cover(jlev)) & |
|---|
| 261 | & / (1.0_jprb - frac(jlev)) |
|---|
| 262 | end if |
|---|
| 263 | cum_cloud_cover(jlev+1) = 1.0_jprb - cum_product |
|---|
| 264 | end do |
|---|
| 265 | |
|---|
| 266 | |
|---|
| 267 | if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_exp_ran',1,hook_handle) |
|---|
| 268 | |
|---|
| 269 | end subroutine cum_cloud_cover_exp_ran |
|---|
| 270 | |
|---|
| 271 | |
|---|
| 272 | |
|---|
| 273 | !--------------------------------------------------------------------- |
|---|
| 274 | ! Exponential-exponential overlap: exponential overlap for both |
|---|
| 275 | ! contiguous and non-contiguous clouds. This is the result of the |
|---|
| 276 | ! simple Raisanen cloud generator, but unfortunately it has no |
|---|
| 277 | ! (known) analytic formula for the total cloud cover, or the |
|---|
| 278 | ! cumulative cloud cover. In partially cloudy columns, The McICA |
|---|
| 279 | ! scheme needs this info in order to devote all the cloudy g-points |
|---|
| 280 | ! to columns containing cloud, which reduces McICA noise. The |
|---|
| 281 | ! following routine provides an approximate estimate of cumulative |
|---|
| 282 | ! cloud cover consistent with the exponential-exponential scheme. |
|---|
| 283 | subroutine cum_cloud_cover_exp_exp(nlev, frac, overlap_param, & |
|---|
| 284 | & cum_cloud_cover, pair_cloud_cover, is_beta_overlap) |
|---|
| 285 | |
|---|
| 286 | use yomhook, only : lhook, dr_hook |
|---|
| 287 | |
|---|
| 288 | implicit none |
|---|
| 289 | |
|---|
| 290 | ! Inputs |
|---|
| 291 | integer, intent(in) :: nlev ! number of model levels |
|---|
| 292 | |
|---|
| 293 | ! Cloud fraction on full levels |
|---|
| 294 | real(jprb), intent(in) :: frac(nlev) |
|---|
| 295 | |
|---|
| 296 | ! Cloud overlap parameter for interfaces between model layers, |
|---|
| 297 | ! where 0 indicates random overlap and 1 indicates maximum-random |
|---|
| 298 | ! overlap |
|---|
| 299 | real(jprb), intent(in) :: overlap_param(nlev-1) |
|---|
| 300 | |
|---|
| 301 | ! This routine has been coded using the "alpha" overlap parameter |
|---|
| 302 | ! of Hogan and Illingworth (2000). If the following logical is |
|---|
| 303 | ! present and true then the input is interpretted to be the "beta" |
|---|
| 304 | ! overlap parameter of Shonk et al. (2010), and needs to be |
|---|
| 305 | ! converted to alpha. |
|---|
| 306 | logical, intent(in), optional :: is_beta_overlap |
|---|
| 307 | |
|---|
| 308 | ! Outputs |
|---|
| 309 | |
|---|
| 310 | ! Cumulative cloud cover from TOA to the base of each layer |
|---|
| 311 | real(jprb), intent(out) :: cum_cloud_cover(nlev) |
|---|
| 312 | |
|---|
| 313 | ! Cloud cover of a pair of layers |
|---|
| 314 | real(jprb), intent(out) :: pair_cloud_cover(nlev-1) |
|---|
| 315 | |
|---|
| 316 | ! Local variables |
|---|
| 317 | |
|---|
| 318 | ! If this routine is called from the radiation_interface tree then |
|---|
| 319 | ! very low cloud fractions have already been set to zero, but if |
|---|
| 320 | ! it is called as a cloud cover diagnostic then this can't be |
|---|
| 321 | ! guaranteed so a small non-zero numbers is required |
|---|
| 322 | real(jprb), parameter :: min_frac = 1.0e-6_jprb |
|---|
| 323 | |
|---|
| 324 | ! "Alpha" overlap parameter |
|---|
| 325 | real(jprb) :: overlap_alpha(nlev-1) |
|---|
| 326 | logical :: do_overlap_conversion |
|---|
| 327 | |
|---|
| 328 | ! Variables describing "concave cloud objects", i.e. contiguous |
|---|
| 329 | ! layers where the cloud fraction monotonically increases then |
|---|
| 330 | ! monotonically decreases |
|---|
| 331 | |
|---|
| 332 | ! Number of objects |
|---|
| 333 | integer :: nobj |
|---|
| 334 | |
|---|
| 335 | ! Indices to the location of the top, maximum cloud fraction, and |
|---|
| 336 | ! base, of each concave cloud object |
|---|
| 337 | integer :: i_top_obj(nlev) |
|---|
| 338 | integer :: i_max_obj(nlev) |
|---|
| 339 | integer :: i_base_obj(nlev) |
|---|
| 340 | ! Poor-man's linked list to allow for deletion of objects: this |
|---|
| 341 | ! variable points to the index of the next active object |
|---|
| 342 | integer :: i_next_obj(nlev) |
|---|
| 343 | |
|---|
| 344 | ! Cloud cover of object |
|---|
| 345 | real(jprb) :: cc_obj(nlev) |
|---|
| 346 | |
|---|
| 347 | ! Overlap parameter between objects |
|---|
| 348 | real(jprb) :: alpha_obj(nlev) |
|---|
| 349 | |
|---|
| 350 | ! Do (while) loop index for model level |
|---|
| 351 | integer :: jlev |
|---|
| 352 | |
|---|
| 353 | ! Do loop index for object |
|---|
| 354 | integer :: jobj |
|---|
| 355 | |
|---|
| 356 | ! Maximum correlation between adjacent objects |
|---|
| 357 | real(jprb) :: alpha_max |
|---|
| 358 | |
|---|
| 359 | ! Indices to pair of objects |
|---|
| 360 | integer :: iobj1, iobj2 |
|---|
| 361 | |
|---|
| 362 | ! Combined cloud cover of pair of objects, and scaling to modify |
|---|
| 363 | ! cumulative cloud cover of lower layer |
|---|
| 364 | real(jprb) :: cc_pair, scaling |
|---|
| 365 | |
|---|
| 366 | real(jprb) :: hook_handle |
|---|
| 367 | |
|---|
| 368 | if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_exp_exp',0,hook_handle) |
|---|
| 369 | |
|---|
| 370 | |
|---|
| 371 | if (present(is_beta_overlap)) then |
|---|
| 372 | do_overlap_conversion = is_beta_overlap |
|---|
| 373 | else |
|---|
| 374 | do_overlap_conversion = .false. |
|---|
| 375 | end if |
|---|
| 376 | |
|---|
| 377 | ! Loop down through atmosphere to locate objects and compute their |
|---|
| 378 | ! basic properties |
|---|
| 379 | jlev = 1 |
|---|
| 380 | nobj = 0 |
|---|
| 381 | do while (jlev <= nlev) |
|---|
| 382 | if (frac(jlev) > min_frac) then |
|---|
| 383 | ! Starting a new object: store its top |
|---|
| 384 | nobj = nobj + 1 |
|---|
| 385 | i_top_obj(nobj) = jlev; |
|---|
| 386 | ! Find its maximum cloud fraction |
|---|
| 387 | jlev = jlev + 1 |
|---|
| 388 | do while (jlev <= nlev) |
|---|
| 389 | if (frac(jlev) < frac(jlev-1)) then |
|---|
| 390 | exit |
|---|
| 391 | end if |
|---|
| 392 | jlev = jlev + 1 |
|---|
| 393 | end do |
|---|
| 394 | i_max_obj(nobj) = jlev - 1 |
|---|
| 395 | ! Find its base |
|---|
| 396 | do while (jlev <= nlev) |
|---|
| 397 | if (frac(jlev) > frac(jlev-1) .or. frac(jlev) <= min_frac) then |
|---|
| 398 | exit |
|---|
| 399 | end if |
|---|
| 400 | jlev = jlev + 1 |
|---|
| 401 | end do |
|---|
| 402 | ! In the case of cloud fraction profile starting from the top |
|---|
| 403 | ! like this: 0.1 0.2 0.1 0.2 0.1, we may want the object grouping |
|---|
| 404 | ! to be (0.1 0.2) (0.1 0.2 0.1), not (0.1 0.2 0.1) (0.2 0.1), |
|---|
| 405 | ! in which case the following should be uncommented |
|---|
| 406 | ! if (jlev < nlev) then |
|---|
| 407 | ! if (frac(jlev) > frac(jlev-1)) then |
|---|
| 408 | ! jlev = jlev - 1 |
|---|
| 409 | ! end if |
|---|
| 410 | ! end if |
|---|
| 411 | i_base_obj(nobj) = jlev - 1 |
|---|
| 412 | ! Index to the next object |
|---|
| 413 | i_next_obj(nobj) = nobj+1 |
|---|
| 414 | else |
|---|
| 415 | jlev = jlev + 1 |
|---|
| 416 | end if |
|---|
| 417 | end do |
|---|
| 418 | |
|---|
| 419 | ! Array assignments |
|---|
| 420 | cum_cloud_cover = 0.0_jprb |
|---|
| 421 | pair_cloud_cover = 0.0_jprb |
|---|
| 422 | |
|---|
| 423 | if (nobj > 0) then |
|---|
| 424 | ! Only do any more work if there is cloud present |
|---|
| 425 | |
|---|
| 426 | ! To minimize the potential calls to beta2alpha, we do all the |
|---|
| 427 | ! computations related to overlap parameter here |
|---|
| 428 | if (.not. do_overlap_conversion) then |
|---|
| 429 | |
|---|
| 430 | ! Compute the combined cloud cover of pairs of layers |
|---|
| 431 | do jlev = 1,nlev-1 |
|---|
| 432 | pair_cloud_cover(jlev) & |
|---|
| 433 | & = overlap_param(jlev)*max(frac(jlev),frac(jlev+1)) & |
|---|
| 434 | & + (1.0_jprb - overlap_param(jlev)) & |
|---|
| 435 | & * (frac(jlev)+frac(jlev+1)-frac(jlev)*frac(jlev+1)) |
|---|
| 436 | end do |
|---|
| 437 | ! Estimate the effective overlap parameter "alpha_obj" between |
|---|
| 438 | ! adjacent objects as the product of the layerwise overlap |
|---|
| 439 | ! parameters between their layers of maximum cloud fraction |
|---|
| 440 | do jobj = 1,nobj-1 |
|---|
| 441 | alpha_obj(jobj) & |
|---|
| 442 | & = product(overlap_param(i_max_obj(jobj):i_max_obj(jobj+1)-1)) |
|---|
| 443 | end do |
|---|
| 444 | |
|---|
| 445 | else |
|---|
| 446 | |
|---|
| 447 | ! Convert Shonk et al overlap parameter to Hogan and |
|---|
| 448 | ! Illingworth definition |
|---|
| 449 | overlap_alpha = beta2alpha(overlap_param, & |
|---|
| 450 | & frac(1:nlev-1), frac(2:nlev)) |
|---|
| 451 | ! Compute the combined cloud cover of pairs of layers |
|---|
| 452 | do jlev = 1,nlev-1 |
|---|
| 453 | pair_cloud_cover(jlev) & |
|---|
| 454 | & = overlap_alpha(jlev)*max(frac(jlev),frac(jlev+1)) & |
|---|
| 455 | & + (1.0_jprb - overlap_alpha(jlev)) & |
|---|
| 456 | & * (frac(jlev)+frac(jlev+1)-frac(jlev)*frac(jlev+1)) |
|---|
| 457 | end do |
|---|
| 458 | ! Estimate the effective overlap parameter "alpha_obj" between |
|---|
| 459 | ! adjacent objects as the product of the layerwise overlap |
|---|
| 460 | ! parameters between their layers of maximum cloud fraction |
|---|
| 461 | do jobj = 1,nobj-1 |
|---|
| 462 | alpha_obj(jobj) & |
|---|
| 463 | & = product(overlap_alpha(i_max_obj(jobj):i_max_obj(jobj+1)-1)) |
|---|
| 464 | end do |
|---|
| 465 | |
|---|
| 466 | end if |
|---|
| 467 | |
|---|
| 468 | |
|---|
| 469 | ! Compute the cumulative cloud cover working down from the top |
|---|
| 470 | ! of each object: this will later be converted to the cumulative |
|---|
| 471 | ! cloud cover working down from TOA |
|---|
| 472 | do jobj = 1,nobj |
|---|
| 473 | cum_cloud_cover(i_top_obj(jobj)) = frac(i_top_obj(jobj)) |
|---|
| 474 | do jlev = i_top_obj(jobj), i_base_obj(jobj)-1 |
|---|
| 475 | if (frac(jlev) >= MaxCloudFrac) then |
|---|
| 476 | ! Cloud cover has reached one |
|---|
| 477 | cum_cloud_cover(jlev+1) = 1.0_jprb |
|---|
| 478 | else |
|---|
| 479 | cum_cloud_cover(jlev+1) = 1.0_jprb & |
|---|
| 480 | & - (1.0_jprb - cum_cloud_cover(jlev)) & |
|---|
| 481 | & * (1.0_jprb - pair_cloud_cover(jlev)) & |
|---|
| 482 | & / (1.0_jprb - frac(jlev)) |
|---|
| 483 | end if |
|---|
| 484 | end do |
|---|
| 485 | cc_obj(jobj) = cum_cloud_cover(i_base_obj(jobj)) |
|---|
| 486 | end do |
|---|
| 487 | |
|---|
| 488 | iobj1 = 1 |
|---|
| 489 | |
|---|
| 490 | ! Sequentially combine objects until there is only one left |
|---|
| 491 | ! covering the full vertical extent of clouds in the profile |
|---|
| 492 | do while (nobj > 1) |
|---|
| 493 | ! Find the most correlated adjacent pair of objects |
|---|
| 494 | alpha_max = 0.0_jprb |
|---|
| 495 | |
|---|
| 496 | ! Need to re-initialize iobj1 here in case alpha_obj(:)==0.0, |
|---|
| 497 | ! which would mean that the "if" statement in the following |
|---|
| 498 | ! loop would never get triggered |
|---|
| 499 | iobj1 = 1 |
|---|
| 500 | |
|---|
| 501 | jobj = 1 |
|---|
| 502 | do while (jobj < nobj) |
|---|
| 503 | if (alpha_obj(jobj) > alpha_max) then |
|---|
| 504 | alpha_max = alpha_obj(jobj) |
|---|
| 505 | iobj1 = jobj |
|---|
| 506 | end if |
|---|
| 507 | jobj = i_next_obj(jobj) |
|---|
| 508 | end do |
|---|
| 509 | |
|---|
| 510 | ! iobj1 is the index to the first object in the pair, set |
|---|
| 511 | ! iobj2 to the second |
|---|
| 512 | iobj2 = i_next_obj(iobj1) |
|---|
| 513 | |
|---|
| 514 | ! Set the cumulative cloud cover in the clear-sky gap between |
|---|
| 515 | ! the objects to the value at the base of the upper object |
|---|
| 516 | cum_cloud_cover(i_base_obj(iobj1)+1:i_top_obj(iobj2)-1) & |
|---|
| 517 | & = cum_cloud_cover(i_base_obj(iobj1)) |
|---|
| 518 | |
|---|
| 519 | ! Calculate the combined cloud cover of the pair of objects |
|---|
| 520 | cc_pair = alpha_obj(iobj1)*max(cc_obj(iobj1), cc_obj(iobj2)) & |
|---|
| 521 | & + (1.0_jprb - alpha_obj(iobj1)) & |
|---|
| 522 | & * (cc_obj(iobj1) + cc_obj(iobj2) - cc_obj(iobj1)*cc_obj(iobj2)) |
|---|
| 523 | scaling = min(max((cc_pair-cc_obj(iobj1)) / max(min_frac, cc_obj(iobj2)), & |
|---|
| 524 | & 0.0_jprb), & |
|---|
| 525 | & 1.0_jprb) |
|---|
| 526 | |
|---|
| 527 | ! Scale the combined cloud cover of the lower object to |
|---|
| 528 | ! account for its overlap with the upper object |
|---|
| 529 | do jlev = i_top_obj(iobj2),i_base_obj(iobj2) |
|---|
| 530 | cum_cloud_cover(jlev) = cum_cloud_cover(i_base_obj(iobj1)) & |
|---|
| 531 | + cum_cloud_cover(jlev) * scaling |
|---|
| 532 | end do |
|---|
| 533 | |
|---|
| 534 | ! Merge the objects by setting the properties of the upper |
|---|
| 535 | ! object to the combined properties of both. Note that |
|---|
| 536 | ! i_max_obj is not modified because it is no longer needed. |
|---|
| 537 | cc_obj(iobj1) = cc_pair |
|---|
| 538 | i_base_obj(iobj1) = i_base_obj(iobj2) |
|---|
| 539 | i_next_obj(iobj1) = i_next_obj(iobj2) |
|---|
| 540 | alpha_obj(iobj1) = alpha_obj(iobj2) |
|---|
| 541 | nobj = nobj - 1 |
|---|
| 542 | end do |
|---|
| 543 | |
|---|
| 544 | ! Finish off the total cloud cover below cloud |
|---|
| 545 | cum_cloud_cover(i_base_obj(iobj1)+1:nlev) & |
|---|
| 546 | & = cum_cloud_cover(i_base_obj(iobj1)) |
|---|
| 547 | |
|---|
| 548 | ! Ensure that the combined cloud cover of pairs of layers is |
|---|
| 549 | ! consistent with the overhang |
|---|
| 550 | do jlev = 1,nlev-1 |
|---|
| 551 | pair_cloud_cover(jlev) = max(pair_cloud_cover(jlev), & |
|---|
| 552 | & frac(jlev)+cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev)) |
|---|
| 553 | end do |
|---|
| 554 | |
|---|
| 555 | ! Sometimes round-off error can lead to cloud cover just above |
|---|
| 556 | ! one, which in turn can lead to direct shortwave fluxes just |
|---|
| 557 | ! below zero |
|---|
| 558 | cum_cloud_cover = min(cum_cloud_cover, 1.0_jprb) |
|---|
| 559 | |
|---|
| 560 | end if ! cloud is present in profile |
|---|
| 561 | |
|---|
| 562 | if (lhook) call dr_hook('radiation_cloud_cover:cum_cloud_cover_exp_exp',1,hook_handle) |
|---|
| 563 | |
|---|
| 564 | end subroutine cum_cloud_cover_exp_exp |
|---|
| 565 | |
|---|
| 566 | end module radiation_cloud_cover |
|---|
