[3908] | 1 | ! radiation_cloud.F90 - Derived type to store cloud/precip properties |
---|
| 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 | ! 2019-01-14 R. Hogan Added inv_inhom_effective_size variable |
---|
| 17 | ! 2019-01-14 R. Hogan Added out_of_physical_bounds routine |
---|
[4584] | 18 | ! 2019-06-14 R. Hogan Added capability to store any number of cloud/precip types |
---|
[3908] | 19 | |
---|
| 20 | module radiation_cloud |
---|
| 21 | |
---|
| 22 | use parkind1, only : jprb |
---|
| 23 | |
---|
| 24 | implicit none |
---|
| 25 | public |
---|
| 26 | |
---|
| 27 | !--------------------------------------------------------------------- |
---|
| 28 | ! The intention is that all variables describing clouds and |
---|
| 29 | ! radiatively-active precipitation are contained in this derived |
---|
| 30 | ! type, and if cloud variables are to be added in future, they can |
---|
| 31 | ! be added to this type without requiring extra variables to be |
---|
| 32 | ! passed between subroutines elsewhere in the program. |
---|
| 33 | type cloud_type |
---|
| 34 | ! For maximum flexibility, an arbitrary number "ntype" of |
---|
[4584] | 35 | ! hydrometeor types can be stored, dimensioned (ncol,nlev,ntype) |
---|
| 36 | integer :: ntype = 0 |
---|
| 37 | real(jprb), allocatable, dimension(:,:,:) :: & |
---|
| 38 | & mixing_ratio, & ! mass mixing ratio (kg/kg) |
---|
| 39 | & effective_radius ! (m) |
---|
| 40 | |
---|
| 41 | ! For backwards compatibility, we also allow for the two |
---|
| 42 | ! traditional cloud types, liquid cloud droplets and ice cloud |
---|
| 43 | ! particles, dimensioned (ncol,nlev) |
---|
| 44 | real(jprb), pointer, dimension(:,:) :: & |
---|
[3908] | 45 | & q_liq, q_ice, & ! mass mixing ratio (kg/kg) |
---|
[4584] | 46 | & re_liq, re_ice ! effective radius (m) |
---|
[3908] | 47 | |
---|
[4584] | 48 | ! For the moment, the different types of hydrometeor are assumed |
---|
| 49 | ! to be mixed with each other, so there is just one cloud fraction |
---|
| 50 | ! variable varying from 0 to 1 |
---|
| 51 | real(jprb), allocatable, dimension(:,:) :: fraction |
---|
| 52 | |
---|
[3908] | 53 | ! The fractional standard deviation of cloud optical depth in the |
---|
| 54 | ! cloudy part of the gridbox. In the Tripleclouds representation |
---|
| 55 | ! of cloud inhomogeneity, this is implemented by splitting the |
---|
| 56 | ! cloudy part of the gridbox into two equal-area regions, one |
---|
| 57 | ! with the cloud optical depth scaled by 1+fractional_std and the |
---|
| 58 | ! other scaled by 1-fractional_std. This variable is dimensioned |
---|
| 59 | ! (ncol,nlev) |
---|
| 60 | real(jprb), allocatable, dimension(:,:) :: fractional_std |
---|
| 61 | |
---|
| 62 | ! The inverse of the effective horizontal size of the clouds in |
---|
| 63 | ! the gridbox, used to compute the cloud edge length per unit |
---|
| 64 | ! gridbox area for use in representing 3D effects. This variable |
---|
| 65 | ! is dimensioned (ncol,nlev). |
---|
| 66 | real(jprb), allocatable, dimension(:,:) :: inv_cloud_effective_size ! m-1 |
---|
| 67 | |
---|
| 68 | ! Similarly for the in-cloud heterogeneities, used to compute the |
---|
| 69 | ! edge length between the optically thin and thick cloudy regions |
---|
| 70 | ! of the gridbox. |
---|
| 71 | real(jprb), allocatable, dimension(:,:) :: inv_inhom_effective_size ! m-1 |
---|
| 72 | |
---|
| 73 | ! The following variable describes the overlap of cloud boundaries |
---|
| 74 | ! in adjacent layers, with dimensions (ncol,nlev-1): 1 corresponds |
---|
| 75 | ! to maximum overlap and 0 to random overlap. Depending on the |
---|
| 76 | ! ecRad configuration, it may be the "alpha" overlap parameter of |
---|
| 77 | ! Hogan and Illingworth (2000) or the "beta" overlap parameter of |
---|
| 78 | ! Shonk et al. (2010). |
---|
| 79 | real(jprb), allocatable, dimension(:,:) :: overlap_param |
---|
| 80 | |
---|
| 81 | contains |
---|
| 82 | procedure :: allocate => allocate_cloud_arrays |
---|
| 83 | procedure :: deallocate => deallocate_cloud_arrays |
---|
| 84 | procedure :: set_overlap_param |
---|
| 85 | procedure :: set_overlap_param_approx |
---|
| 86 | procedure :: create_fractional_std |
---|
| 87 | procedure :: create_inv_cloud_effective_size |
---|
| 88 | procedure :: create_inv_cloud_effective_size_eta |
---|
| 89 | procedure :: param_cloud_effective_separation_eta |
---|
| 90 | procedure :: crop_cloud_fraction |
---|
| 91 | procedure :: out_of_physical_bounds |
---|
| 92 | |
---|
| 93 | end type cloud_type |
---|
| 94 | |
---|
| 95 | contains |
---|
| 96 | |
---|
| 97 | !--------------------------------------------------------------------- |
---|
| 98 | ! Allocate arrays for describing clouds and precipitation, although |
---|
| 99 | ! in the offline code these are allocated when they are read from |
---|
| 100 | ! the NetCDF file |
---|
[4584] | 101 | subroutine allocate_cloud_arrays(this, ncol, nlev, ntype, use_inhom_effective_size) |
---|
[3908] | 102 | |
---|
| 103 | use yomhook, only : lhook, dr_hook |
---|
| 104 | |
---|
[4584] | 105 | class(cloud_type), intent(inout), target :: this |
---|
| 106 | integer, intent(in) :: ncol ! Number of columns |
---|
| 107 | integer, intent(in) :: nlev ! Number of levels |
---|
| 108 | ! Number of cloud/precip particle types. If not present then the |
---|
| 109 | ! older cloud behaviour is assumed: two types are present, (1) |
---|
| 110 | ! liquid and (2) ice, and they can be accessed via q_liq, q_ice, |
---|
| 111 | ! re_liq and re_ice. |
---|
| 112 | integer, intent(in), optional :: ntype |
---|
[3908] | 113 | logical, intent(in), optional :: use_inhom_effective_size |
---|
| 114 | |
---|
| 115 | real(jprb) :: hook_handle |
---|
| 116 | |
---|
| 117 | if (lhook) call dr_hook('radiation_cloud:allocate',0,hook_handle) |
---|
| 118 | |
---|
[4584] | 119 | if (present(ntype)) then |
---|
| 120 | this%ntype = ntype |
---|
| 121 | else |
---|
| 122 | this%ntype = 2 |
---|
| 123 | end if |
---|
| 124 | allocate(this%mixing_ratio(ncol,nlev,this%ntype)) |
---|
| 125 | allocate(this%effective_radius(ncol,nlev,this%ntype)) |
---|
| 126 | nullify(this%q_liq) |
---|
| 127 | nullify(this%q_ice) |
---|
| 128 | nullify(this%re_liq) |
---|
| 129 | nullify(this%re_ice) |
---|
| 130 | if (.not. present(ntype)) then |
---|
| 131 | ! Older interface in which only liquid and ice are supported |
---|
| 132 | this%q_liq => this%mixing_ratio(:,:,1) |
---|
| 133 | this%q_ice => this%mixing_ratio(:,:,2) |
---|
| 134 | this%re_liq => this%effective_radius(:,:,1) |
---|
| 135 | this%re_ice => this%effective_radius(:,:,2) |
---|
| 136 | end if |
---|
| 137 | |
---|
[3908] | 138 | allocate(this%fraction(ncol,nlev)) |
---|
| 139 | allocate(this%overlap_param(ncol,nlev-1)) |
---|
| 140 | allocate(this%fractional_std(ncol,nlev)) |
---|
| 141 | allocate(this%inv_cloud_effective_size(ncol,nlev)) |
---|
| 142 | |
---|
| 143 | if (present(use_inhom_effective_size)) then |
---|
| 144 | if (use_inhom_effective_size) then |
---|
| 145 | allocate(this%inv_inhom_effective_size(ncol,nlev)) |
---|
| 146 | end if |
---|
| 147 | end if |
---|
| 148 | |
---|
| 149 | if (lhook) call dr_hook('radiation_cloud:allocate',1,hook_handle) |
---|
| 150 | |
---|
| 151 | end subroutine allocate_cloud_arrays |
---|
| 152 | |
---|
| 153 | |
---|
| 154 | !--------------------------------------------------------------------- |
---|
| 155 | ! Deallocate arrays |
---|
| 156 | subroutine deallocate_cloud_arrays(this) |
---|
| 157 | |
---|
| 158 | use yomhook, only : lhook, dr_hook |
---|
| 159 | |
---|
| 160 | class(cloud_type), intent(inout) :: this |
---|
| 161 | |
---|
| 162 | real(jprb) :: hook_handle |
---|
| 163 | |
---|
| 164 | if (lhook) call dr_hook('radiation_cloud:deallocate',0,hook_handle) |
---|
| 165 | |
---|
[4584] | 166 | nullify(this%q_liq) |
---|
| 167 | nullify(this%q_ice) |
---|
| 168 | nullify(this%re_liq) |
---|
| 169 | nullify(this%re_ice) |
---|
| 170 | |
---|
| 171 | if (allocated(this%mixing_ratio)) deallocate(this%mixing_ratio) |
---|
| 172 | if (allocated(this%effective_radius)) deallocate(this%effective_radius) |
---|
| 173 | if (allocated(this%fraction)) deallocate(this%fraction) |
---|
| 174 | if (allocated(this%overlap_param)) deallocate(this%overlap_param) |
---|
| 175 | if (allocated(this%fractional_std)) deallocate(this%fractional_std) |
---|
[3908] | 176 | if (allocated(this%inv_cloud_effective_size)) & |
---|
| 177 | & deallocate(this%inv_cloud_effective_size) |
---|
| 178 | if (allocated(this%inv_inhom_effective_size)) & |
---|
| 179 | & deallocate(this%inv_inhom_effective_size) |
---|
| 180 | |
---|
| 181 | if (lhook) call dr_hook('radiation_cloud:deallocate',1,hook_handle) |
---|
| 182 | |
---|
| 183 | end subroutine deallocate_cloud_arrays |
---|
| 184 | |
---|
| 185 | |
---|
| 186 | !--------------------------------------------------------------------- |
---|
| 187 | ! Compute and store the overlap parameter from the provided overlap |
---|
| 188 | ! decorrelation length (in metres). If istartcol and/or iendcol are |
---|
| 189 | ! provided then only columns in this range are computed. If the |
---|
| 190 | ! overlap_param array has not been allocated then it will be |
---|
| 191 | ! allocated to be of the correct size relative to the pressure |
---|
| 192 | ! field. |
---|
| 193 | subroutine set_overlap_param(this, thermodynamics, decorrelation_length, & |
---|
| 194 | & istartcol, iendcol) |
---|
| 195 | |
---|
| 196 | use yomhook, only : lhook, dr_hook |
---|
| 197 | use radiation_thermodynamics, only : thermodynamics_type |
---|
| 198 | use radiation_constants, only : GasConstantDryAir, AccelDueToGravity |
---|
| 199 | |
---|
| 200 | class(cloud_type), intent(inout) :: this |
---|
| 201 | type(thermodynamics_type), intent(in) :: thermodynamics |
---|
| 202 | real(jprb), intent(in) :: decorrelation_length ! m |
---|
| 203 | integer, optional, intent(in) :: istartcol, iendcol |
---|
| 204 | |
---|
| 205 | ! Ratio of gas constant for dry air to acceleration due to gravity |
---|
| 206 | real(jprb), parameter :: R_over_g = GasConstantDryAir / AccelDueToGravity |
---|
| 207 | |
---|
| 208 | ! Process only columns i1 to i2, which will be istartcol to |
---|
| 209 | ! iendcol if they were provided |
---|
| 210 | integer :: i1, i2 |
---|
| 211 | |
---|
| 212 | integer :: ncol, nlev |
---|
| 213 | |
---|
[4584] | 214 | integer :: jcol, jlev |
---|
[3908] | 215 | |
---|
| 216 | real(jprb) :: hook_handle |
---|
| 217 | |
---|
| 218 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param',0,hook_handle) |
---|
| 219 | |
---|
| 220 | ! Pressure at half-levels, pressure_hl, is defined at nlev+1 |
---|
| 221 | ! points |
---|
| 222 | ncol = size(thermodynamics%pressure_hl,dim=1) |
---|
| 223 | nlev = size(thermodynamics%pressure_hl,dim=2)-1 |
---|
| 224 | |
---|
| 225 | if (present(istartcol)) then |
---|
| 226 | i1 = istartcol |
---|
| 227 | else |
---|
| 228 | i1 = 1 |
---|
| 229 | end if |
---|
| 230 | |
---|
| 231 | if (present(iendcol)) then |
---|
| 232 | i2 = iendcol |
---|
| 233 | else |
---|
| 234 | i2 = ncol |
---|
| 235 | end if |
---|
| 236 | |
---|
| 237 | if (.not. allocated(this%overlap_param)) then |
---|
| 238 | ! If pressure is of size (ncol,nlev+1) then overlap_param is of |
---|
| 239 | ! size (ncol,nlev-1), since overlap parameter is only defined here |
---|
| 240 | ! for interfaces between model layers, not for the interface to |
---|
| 241 | ! space or the surface |
---|
| 242 | allocate(this%overlap_param(ncol, nlev-1)) |
---|
| 243 | end if |
---|
| 244 | |
---|
| 245 | if (thermodynamics%pressure_hl(i1,2) > thermodynamics%pressure_hl(i1,1)) then |
---|
| 246 | ! Pressure is increasing with index (order of layers is |
---|
| 247 | ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we |
---|
| 248 | ! don't take the logarithm of the first pressure in each column. |
---|
[4584] | 249 | do jcol = i1,i2 |
---|
| 250 | this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length) & |
---|
| 251 | & * thermodynamics%temperature_hl(jcol,2) & |
---|
| 252 | & *log(thermodynamics%pressure_hl(jcol,3) & |
---|
| 253 | & /thermodynamics%pressure_hl(jcol,2))) |
---|
| 254 | end do |
---|
[3908] | 255 | |
---|
| 256 | do jlev = 2,nlev-1 |
---|
[4584] | 257 | do jcol = i1,i2 |
---|
| 258 | this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) & |
---|
| 259 | & * thermodynamics%temperature_hl(jcol,jlev+1) & |
---|
| 260 | & *log(thermodynamics%pressure_hl(jcol,jlev+2) & |
---|
| 261 | & /thermodynamics%pressure_hl(jcol,jlev))) |
---|
| 262 | end do |
---|
[3908] | 263 | end do |
---|
| 264 | |
---|
| 265 | else |
---|
| 266 | ! Pressure is decreasing with index (order of layers is surface |
---|
| 267 | ! to top-of-atmosphere). In case pressure_hl(:,nlev+1)=0, we |
---|
| 268 | ! don't take the logarithm of the last pressure in each column. |
---|
| 269 | do jlev = 1,nlev-2 |
---|
[4584] | 270 | do jcol = i1,i2 |
---|
| 271 | this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) & |
---|
| 272 | & * thermodynamics%temperature_hl(jcol,jlev+1) & |
---|
| 273 | & *log(thermodynamics%pressure_hl(jcol,jlev) & |
---|
| 274 | & /thermodynamics%pressure_hl(jcol,jlev+2))) |
---|
| 275 | end do |
---|
[3908] | 276 | end do |
---|
| 277 | |
---|
[4584] | 278 | do jcol = i1,i2 |
---|
| 279 | this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length) & |
---|
| 280 | & * thermodynamics%temperature_hl(jcol,nlev) & |
---|
| 281 | & *log(thermodynamics%pressure_hl(jcol,nlev-1) & |
---|
| 282 | & /thermodynamics%pressure_hl(jcol,nlev))) |
---|
| 283 | end do |
---|
[3908] | 284 | end if |
---|
| 285 | |
---|
| 286 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param',1,hook_handle) |
---|
| 287 | |
---|
| 288 | end subroutine set_overlap_param |
---|
| 289 | |
---|
| 290 | |
---|
| 291 | !--------------------------------------------------------------------- |
---|
| 292 | ! Compute and store the overlap parameter from the provided overlap |
---|
| 293 | ! decorrelation length (in metres). If istartcol and/or iendcol are |
---|
| 294 | ! provided then only columns in this range are computed. If the |
---|
| 295 | ! overlap_param array has not been allocated then it will be |
---|
| 296 | ! allocated to be of the correct size relative to the pressure |
---|
| 297 | ! field. This is the APPROXIMATE method as it assumes a fixed |
---|
| 298 | ! atmospheric scale height, which leads to differences particularly |
---|
| 299 | ! in low cloud. |
---|
| 300 | subroutine set_overlap_param_approx(this, thermodynamics, decorrelation_length, & |
---|
| 301 | & istartcol, iendcol) |
---|
| 302 | |
---|
| 303 | use yomhook, only : lhook, dr_hook |
---|
| 304 | use radiation_thermodynamics, only : thermodynamics_type |
---|
| 305 | |
---|
| 306 | class(cloud_type), intent(inout) :: this |
---|
| 307 | type(thermodynamics_type), intent(in) :: thermodynamics |
---|
| 308 | real(jprb), intent(in) :: decorrelation_length ! m |
---|
| 309 | integer, optional, intent(in) :: istartcol, iendcol |
---|
| 310 | |
---|
| 311 | ! To convert decorrelation length (m) to overlap parameter between |
---|
| 312 | ! layers, we need an estimate for the thickness of the layer. This |
---|
| 313 | ! is found using the pressure difference between the edges of the |
---|
| 314 | ! layer, along with the approximate scale height of the atmosphere |
---|
| 315 | ! (m) given here: |
---|
| 316 | real(jprb), parameter :: scale_height = 8000.0_jprb |
---|
| 317 | |
---|
| 318 | ! Process only columns i1 to i2, which will be istartcol to |
---|
| 319 | ! iendcol if they were provided |
---|
| 320 | integer :: i1, i2 |
---|
| 321 | |
---|
| 322 | integer :: ncol, nlev |
---|
| 323 | |
---|
| 324 | real(jprb) :: hook_handle |
---|
| 325 | |
---|
| 326 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param_approx',0,hook_handle) |
---|
| 327 | |
---|
| 328 | ! Pressure at half-levels, pressure_hl, is defined at nlev+1 |
---|
| 329 | ! points |
---|
| 330 | ncol = size(thermodynamics%pressure_hl,dim=1) |
---|
| 331 | nlev = size(thermodynamics%pressure_hl,dim=2)-1 |
---|
| 332 | |
---|
| 333 | if (present(istartcol)) then |
---|
| 334 | i1 = istartcol |
---|
| 335 | else |
---|
| 336 | i1 = 1 |
---|
| 337 | end if |
---|
| 338 | |
---|
| 339 | if (present(iendcol)) then |
---|
| 340 | i2 = iendcol |
---|
| 341 | else |
---|
| 342 | i2 = ncol |
---|
| 343 | end if |
---|
| 344 | |
---|
| 345 | if (.not. allocated(this%overlap_param)) then |
---|
| 346 | ! If pressure is of size (ncol,nlev+1) then overlap_param is of |
---|
| 347 | ! size (ncol,nlev-1), since overlap parameter is only defined here |
---|
| 348 | ! for interfaces between model layers, not for the interface to |
---|
| 349 | ! space or the surface |
---|
| 350 | allocate(this%overlap_param(ncol, nlev-1)) |
---|
| 351 | end if |
---|
| 352 | |
---|
| 353 | if (thermodynamics%pressure_hl(i1,2) > thermodynamics%pressure_hl(i1,1)) then |
---|
| 354 | ! Pressure is increasing with index (order of layers is |
---|
| 355 | ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we |
---|
| 356 | ! don't take the logarithm of the first pressure in each |
---|
| 357 | ! column. |
---|
| 358 | this%overlap_param(i1:i2,:) = exp(-(scale_height/decorrelation_length) & |
---|
| 359 | & * ( log(thermodynamics%pressure_hl(i1:i2,3:nlev+1) & |
---|
| 360 | & /thermodynamics%pressure_hl(i1:i2,2:nlev )) ) ) |
---|
| 361 | else |
---|
| 362 | ! Pressure is decreasing with index (order of layers is surface |
---|
| 363 | ! to top-of-atmosphere). In case pressure_hl(:,nlev+1)=0, we |
---|
| 364 | ! don't take the logarithm of the last pressure in each column. |
---|
| 365 | this%overlap_param(i1:i2,:) = exp(-(scale_height/decorrelation_length) & |
---|
| 366 | & * ( log(thermodynamics%pressure_hl(i1:i2,1:nlev-1) & |
---|
| 367 | & /thermodynamics%pressure_hl(i1:i2,2:nlev )) ) ) |
---|
| 368 | end if |
---|
| 369 | |
---|
| 370 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param_approx',1,hook_handle) |
---|
| 371 | |
---|
| 372 | end subroutine set_overlap_param_approx |
---|
| 373 | |
---|
| 374 | |
---|
| 375 | !--------------------------------------------------------------------- |
---|
| 376 | ! Create a matrix of constant fractional standard deviations |
---|
| 377 | ! (dimensionless) |
---|
| 378 | subroutine create_fractional_std(this, ncol, nlev, frac_std) |
---|
| 379 | |
---|
| 380 | use yomhook, only : lhook, dr_hook |
---|
| 381 | |
---|
| 382 | class(cloud_type), intent(inout) :: this |
---|
| 383 | integer, intent(in) :: ncol, nlev |
---|
| 384 | real(jprb), intent(in) :: frac_std |
---|
| 385 | |
---|
| 386 | real(jprb) :: hook_handle |
---|
| 387 | |
---|
| 388 | if (lhook) call dr_hook('radiation_cloud:create_fractional_std',0,hook_handle) |
---|
| 389 | |
---|
| 390 | if (allocated(this%fractional_std)) then |
---|
| 391 | deallocate(this%fractional_std) |
---|
| 392 | end if |
---|
| 393 | |
---|
| 394 | allocate(this%fractional_std(ncol, nlev)) |
---|
| 395 | |
---|
| 396 | this%fractional_std = frac_std |
---|
| 397 | |
---|
| 398 | if (lhook) call dr_hook('radiation_cloud:create_fractional_std',1,hook_handle) |
---|
| 399 | |
---|
| 400 | end subroutine create_fractional_std |
---|
| 401 | |
---|
| 402 | |
---|
| 403 | !--------------------------------------------------------------------- |
---|
| 404 | ! Create a matrix of constant inverse cloud effective size (m-1) |
---|
| 405 | subroutine create_inv_cloud_effective_size(this, ncol, nlev, inv_eff_size) |
---|
| 406 | |
---|
| 407 | use yomhook, only : lhook, dr_hook |
---|
| 408 | |
---|
| 409 | class(cloud_type), intent(inout) :: this |
---|
| 410 | integer, intent(in) :: ncol, nlev |
---|
| 411 | real(jprb), intent(in) :: inv_eff_size |
---|
| 412 | |
---|
| 413 | real(jprb) :: hook_handle |
---|
| 414 | |
---|
| 415 | if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size',0,hook_handle) |
---|
| 416 | |
---|
| 417 | if (allocated(this%inv_cloud_effective_size)) then |
---|
| 418 | deallocate(this%inv_cloud_effective_size) |
---|
| 419 | end if |
---|
| 420 | |
---|
| 421 | allocate(this%inv_cloud_effective_size(ncol, nlev)) |
---|
| 422 | |
---|
| 423 | this%inv_cloud_effective_size = inv_eff_size |
---|
| 424 | |
---|
| 425 | if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size',1,hook_handle) |
---|
| 426 | |
---|
| 427 | end subroutine create_inv_cloud_effective_size |
---|
| 428 | |
---|
| 429 | |
---|
| 430 | !--------------------------------------------------------------------- |
---|
| 431 | ! Create a matrix of inverse cloud effective size (m-1) according to |
---|
| 432 | ! the value of eta (=pressure divided by surface pressure) |
---|
| 433 | subroutine create_inv_cloud_effective_size_eta(this, ncol, nlev, & |
---|
| 434 | & pressure_hl, inv_eff_size_low, inv_eff_size_mid, inv_eff_size_high, & |
---|
| 435 | & eta_low_mid, eta_mid_high, istartcol, iendcol) |
---|
| 436 | |
---|
| 437 | use yomhook, only : lhook, dr_hook |
---|
| 438 | |
---|
| 439 | class(cloud_type), intent(inout) :: this |
---|
| 440 | integer, intent(in) :: ncol, nlev |
---|
| 441 | ! Pressure on half levels (Pa) |
---|
| 442 | real(jprb), intent(in) :: pressure_hl(:,:) |
---|
| 443 | ! Inverse effective size for low, mid and high cloud (m-1) |
---|
| 444 | real(jprb), intent(in) :: inv_eff_size_low |
---|
| 445 | real(jprb), intent(in) :: inv_eff_size_mid |
---|
| 446 | real(jprb), intent(in) :: inv_eff_size_high |
---|
| 447 | ! Eta values at low-mid and mid-high interfaces |
---|
| 448 | real(jprb), intent(in) :: eta_low_mid, eta_mid_high |
---|
| 449 | integer, optional, intent(in) :: istartcol, iendcol |
---|
| 450 | |
---|
| 451 | ! Ratio of layer midpoint pressure to surface pressure |
---|
| 452 | real(jprb) :: eta(nlev) |
---|
| 453 | |
---|
| 454 | ! Indices of column, level and surface half-level |
---|
| 455 | integer :: jcol, isurf |
---|
| 456 | |
---|
| 457 | ! Local values of istartcol, iendcol |
---|
| 458 | integer :: i1, i2 |
---|
| 459 | |
---|
| 460 | real(jprb) :: hook_handle |
---|
| 461 | |
---|
| 462 | if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size_eta',0,hook_handle) |
---|
| 463 | |
---|
| 464 | if (allocated(this%inv_cloud_effective_size)) then |
---|
| 465 | deallocate(this%inv_cloud_effective_size) |
---|
| 466 | end if |
---|
| 467 | |
---|
| 468 | allocate(this%inv_cloud_effective_size(ncol, nlev)) |
---|
| 469 | |
---|
| 470 | if (present(istartcol)) then |
---|
| 471 | i1 = istartcol |
---|
| 472 | else |
---|
| 473 | i1 = 1 |
---|
| 474 | end if |
---|
| 475 | |
---|
| 476 | if (present(iendcol)) then |
---|
| 477 | i2 = iendcol |
---|
| 478 | else |
---|
| 479 | i2 = ncol |
---|
| 480 | end if |
---|
| 481 | |
---|
| 482 | ! Locate the surface half-level |
---|
| 483 | if (pressure_hl(1,1) > pressure_hl(1,2)) then |
---|
| 484 | isurf = 1 |
---|
| 485 | else |
---|
| 486 | isurf = nlev+1 |
---|
| 487 | end if |
---|
| 488 | |
---|
| 489 | do jcol = i1,i2 |
---|
| 490 | eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) & |
---|
| 491 | & * (0.5_jprb / pressure_hl(jcol,isurf)) |
---|
| 492 | where (eta > eta_low_mid) |
---|
| 493 | this%inv_cloud_effective_size(jcol,:) = inv_eff_size_low |
---|
| 494 | elsewhere (eta > eta_mid_high) |
---|
| 495 | this%inv_cloud_effective_size(jcol,:) = inv_eff_size_mid |
---|
| 496 | elsewhere |
---|
| 497 | this%inv_cloud_effective_size(jcol,:) = inv_eff_size_high |
---|
| 498 | end where |
---|
| 499 | end do |
---|
| 500 | |
---|
| 501 | if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size_eta',1,hook_handle) |
---|
| 502 | |
---|
| 503 | end subroutine create_inv_cloud_effective_size_eta |
---|
| 504 | |
---|
| 505 | |
---|
| 506 | !--------------------------------------------------------------------- |
---|
| 507 | ! Create a matrix of inverse cloud and inhomogeneity effective size |
---|
| 508 | ! (m-1) parameterized according to the value of eta (=pressure |
---|
| 509 | ! divided by surface pressure): effective_separation = |
---|
| 510 | ! coeff_a + coeff_b*exp(-(eta**power)). |
---|
| 511 | subroutine param_cloud_effective_separation_eta(this, ncol, nlev, & |
---|
| 512 | & pressure_hl, separation_surf, separation_toa, power, & |
---|
| 513 | & inhom_separation_factor, istartcol, iendcol) |
---|
| 514 | |
---|
| 515 | use yomhook, only : lhook, dr_hook |
---|
| 516 | |
---|
| 517 | class(cloud_type), intent(inout) :: this |
---|
| 518 | integer, intent(in) :: ncol, nlev |
---|
| 519 | ! Pressure on half levels (Pa) |
---|
| 520 | real(jprb), intent(in) :: pressure_hl(:,:) |
---|
| 521 | ! Separation distances at surface and top-of-atmosphere, and power |
---|
| 522 | ! on eta |
---|
| 523 | real(jprb), intent(in) :: separation_surf ! m |
---|
| 524 | real(jprb), intent(in) :: separation_toa ! m |
---|
| 525 | real(jprb), intent(in) :: power |
---|
| 526 | real(jprb), optional, intent(in) :: inhom_separation_factor |
---|
| 527 | integer, optional, intent(in) :: istartcol, iendcol |
---|
| 528 | |
---|
| 529 | ! Ratio of layer midpoint pressure to surface pressure |
---|
| 530 | real(jprb) :: eta(nlev) |
---|
| 531 | |
---|
| 532 | ! Effective cloud separation (m) |
---|
| 533 | real(jprb) :: eff_separation(nlev) |
---|
| 534 | |
---|
| 535 | ! Coefficients used to compute effective separation distance |
---|
| 536 | real(jprb) :: coeff_e, coeff_a, coeff_b, inhom_sep_factor |
---|
| 537 | |
---|
| 538 | ! Indices of column, level and surface half-level |
---|
| 539 | integer :: jcol, isurf |
---|
| 540 | |
---|
| 541 | ! Local values of istartcol, iendcol |
---|
| 542 | integer :: i1, i2 |
---|
| 543 | |
---|
| 544 | real(jprb) :: hook_handle |
---|
| 545 | |
---|
| 546 | if (lhook) call dr_hook('radiation_cloud:param_cloud_effective_separation_eta',0,hook_handle) |
---|
| 547 | |
---|
| 548 | if (present(inhom_separation_factor)) then |
---|
| 549 | inhom_sep_factor = inhom_separation_factor |
---|
| 550 | else |
---|
| 551 | inhom_sep_factor = 1.0_jprb |
---|
| 552 | end if |
---|
| 553 | |
---|
| 554 | coeff_e = 1.0_jprb - exp(-1.0_jprb) |
---|
| 555 | coeff_b = (separation_toa - separation_surf) / coeff_e |
---|
| 556 | coeff_a = separation_toa - coeff_b |
---|
| 557 | |
---|
| 558 | if (allocated(this%inv_cloud_effective_size)) then |
---|
| 559 | deallocate(this%inv_cloud_effective_size) |
---|
| 560 | end if |
---|
| 561 | if (allocated(this%inv_inhom_effective_size)) then |
---|
| 562 | deallocate(this%inv_inhom_effective_size) |
---|
| 563 | end if |
---|
| 564 | |
---|
| 565 | allocate(this%inv_cloud_effective_size(ncol, nlev)) |
---|
| 566 | allocate(this%inv_inhom_effective_size(ncol, nlev)) |
---|
| 567 | |
---|
| 568 | if (present(istartcol)) then |
---|
| 569 | i1 = istartcol |
---|
| 570 | else |
---|
| 571 | i1 = 1 |
---|
| 572 | end if |
---|
| 573 | |
---|
| 574 | if (present(iendcol)) then |
---|
| 575 | i2 = iendcol |
---|
| 576 | else |
---|
| 577 | i2 = ncol |
---|
| 578 | end if |
---|
| 579 | |
---|
| 580 | ! Locate the surface half-level |
---|
| 581 | if (pressure_hl(1,1) > pressure_hl(1,2)) then |
---|
| 582 | isurf = 1 |
---|
| 583 | else |
---|
| 584 | isurf = nlev+1 |
---|
| 585 | end if |
---|
| 586 | |
---|
| 587 | do jcol = i1,i2 |
---|
| 588 | eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) & |
---|
| 589 | & * (0.5_jprb / pressure_hl(jcol,isurf)) |
---|
| 590 | eff_separation = coeff_a + coeff_b * exp(-eta**power) |
---|
| 591 | this%inv_cloud_effective_size(jcol,:) = 1.0_jprb / (eff_separation & |
---|
| 592 | & * sqrt(max(1.0e-5_jprb,this%fraction(jcol,:)*(1.0_jprb-this%fraction(jcol,:))))) |
---|
| 593 | this%inv_inhom_effective_size(jcol,:) = 1.0_jprb / (eff_separation * inhom_sep_factor & |
---|
| 594 | & * sqrt(max(1.0e-5_jprb,0.5_jprb*this%fraction(jcol,:)*(1.0_jprb-0.5_jprb*this%fraction(jcol,:))))) |
---|
| 595 | end do |
---|
| 596 | |
---|
| 597 | if (lhook) call dr_hook('radiation_cloud:param_cloud_effective_separation_eta',1,hook_handle) |
---|
| 598 | |
---|
| 599 | end subroutine param_cloud_effective_separation_eta |
---|
| 600 | |
---|
| 601 | |
---|
| 602 | !--------------------------------------------------------------------- |
---|
| 603 | ! Remove "ghost" clouds: those with a cloud fraction that is too |
---|
| 604 | ! small to treat sensibly (e.g. because it implies that the |
---|
| 605 | ! "in-cloud" water content is too high), or with a cloud water |
---|
| 606 | ! content that is too small. We do this in one place to ensure that |
---|
| 607 | ! all subsequent subroutines can assume that if cloud_fraction > 0.0 |
---|
| 608 | ! then cloud is really present and should be treated. |
---|
| 609 | subroutine crop_cloud_fraction(this, istartcol, iendcol, & |
---|
| 610 | & cloud_fraction_threshold, cloud_mixing_ratio_threshold) |
---|
| 611 | |
---|
| 612 | use yomhook, only : lhook, dr_hook |
---|
| 613 | |
---|
| 614 | class(cloud_type), intent(inout) :: this |
---|
| 615 | integer, intent(in) :: istartcol, iendcol |
---|
| 616 | |
---|
[4584] | 617 | integer :: nlev, ntype |
---|
| 618 | integer :: jcol, jlev, jh |
---|
[3908] | 619 | |
---|
| 620 | real(jprb) :: cloud_fraction_threshold, cloud_mixing_ratio_threshold |
---|
[4584] | 621 | real(jprb) :: sum_mixing_ratio(istartcol:iendcol) |
---|
[3908] | 622 | |
---|
| 623 | real(jprb) :: hook_handle |
---|
| 624 | |
---|
| 625 | if (lhook) call dr_hook('radiation_cloud:crop_cloud_fraction',0,hook_handle) |
---|
| 626 | |
---|
[4584] | 627 | nlev = size(this%fraction,2) |
---|
| 628 | ntype = size(this%mixing_ratio,3) |
---|
| 629 | |
---|
[3908] | 630 | do jlev = 1,nlev |
---|
| 631 | do jcol = istartcol,iendcol |
---|
[4584] | 632 | sum_mixing_ratio(jcol) = 0.0_jprb |
---|
| 633 | end do |
---|
| 634 | do jh = 1, ntype |
---|
| 635 | do jcol = istartcol,iendcol |
---|
| 636 | sum_mixing_ratio(jcol) = sum_mixing_ratio(jcol) + this%mixing_ratio(jcol,jlev,jh) |
---|
| 637 | end do |
---|
| 638 | end do |
---|
| 639 | do jcol = istartcol,iendcol |
---|
| 640 | if (this%fraction(jcol,jlev) < cloud_fraction_threshold & |
---|
| 641 | & .or. sum_mixing_ratio(jcol) < cloud_mixing_ratio_threshold) then |
---|
[3908] | 642 | this%fraction(jcol,jlev) = 0.0_jprb |
---|
| 643 | end if |
---|
| 644 | end do |
---|
| 645 | end do |
---|
| 646 | |
---|
| 647 | if (lhook) call dr_hook('radiation_cloud:crop_cloud_fraction',1,hook_handle) |
---|
| 648 | |
---|
| 649 | end subroutine crop_cloud_fraction |
---|
| 650 | |
---|
| 651 | |
---|
| 652 | !--------------------------------------------------------------------- |
---|
| 653 | ! Return .true. if variables are out of a physically sensible range, |
---|
| 654 | ! optionally only considering columns between istartcol and iendcol |
---|
| 655 | function out_of_physical_bounds(this, istartcol, iendcol, do_fix) result(is_bad) |
---|
| 656 | |
---|
| 657 | use yomhook, only : lhook, dr_hook |
---|
[4584] | 658 | use radiation_check, only : out_of_bounds_2d, out_of_bounds_3d |
---|
[3908] | 659 | |
---|
| 660 | class(cloud_type), intent(inout) :: this |
---|
| 661 | integer, optional,intent(in) :: istartcol, iendcol |
---|
| 662 | logical, optional,intent(in) :: do_fix |
---|
| 663 | logical :: is_bad |
---|
| 664 | |
---|
| 665 | logical :: do_fix_local |
---|
| 666 | |
---|
| 667 | real(jprb) :: hook_handle |
---|
| 668 | |
---|
| 669 | if (lhook) call dr_hook('radiation_cloud:out_of_physical_bounds',0,hook_handle) |
---|
| 670 | |
---|
| 671 | if (present(do_fix)) then |
---|
| 672 | do_fix_local = do_fix |
---|
| 673 | else |
---|
| 674 | do_fix_local = .false. |
---|
| 675 | end if |
---|
| 676 | |
---|
[4584] | 677 | is_bad = out_of_bounds_3d(this%mixing_ratio, 'cloud%mixing_ratio', 0.0_jprb, 1.0_jprb, & |
---|
[3908] | 678 | & do_fix_local, i1=istartcol, i2=iendcol) & |
---|
[4584] | 679 | & .or. out_of_bounds_3d(this%effective_radius, 'cloud%effective_radius', 0.0_jprb, 0.1_jprb, & |
---|
[3908] | 680 | & do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 681 | & .or. out_of_bounds_2d(this%fraction, 'cloud%fraction', 0.0_jprb, 1.0_jprb, & |
---|
| 682 | & do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 683 | & .or. out_of_bounds_2d(this%fractional_std, 'fractional_std', 0.0_jprb, 10.0_jprb, & |
---|
| 684 | & do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 685 | & .or. out_of_bounds_2d(this%inv_cloud_effective_size, 'inv_cloud_effective_size', & |
---|
| 686 | & 0.0_jprb, 1.0_jprb, do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 687 | & .or. out_of_bounds_2d(this%inv_inhom_effective_size, 'inv_inhom_effective_size', & |
---|
| 688 | & 0.0_jprb, 1.0_jprb, do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 689 | & .or. out_of_bounds_2d(this%overlap_param, 'overlap_param', -0.5_jprb, 1.0_jprb, & |
---|
| 690 | & do_fix_local, i1=istartcol, i2=iendcol) |
---|
| 691 | |
---|
| 692 | if (lhook) call dr_hook('radiation_cloud:out_of_physical_bounds',1,hook_handle) |
---|
| 693 | |
---|
| 694 | end function out_of_physical_bounds |
---|
| 695 | |
---|
| 696 | end module radiation_cloud |
---|