[4773] | 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 |
---|
| 18 | ! 2019-06-14 R. Hogan Added capability to store any number of cloud/precip types |
---|
| 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 |
---|
| 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(:,:) :: & |
---|
| 45 | & q_liq, q_ice, & ! mass mixing ratio (kg/kg) |
---|
| 46 | & re_liq, re_ice ! effective radius (m) |
---|
| 47 | |
---|
| 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 | |
---|
| 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). |
---|
[4853] | 66 | real(jprb), allocatable, dimension(:,:) :: inv_cloud_effective_size ! m-1 |
---|
[4773] | 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_fix |
---|
| 85 | procedure :: set_overlap_param_var |
---|
| 86 | generic :: set_overlap_param => set_overlap_param_fix, set_overlap_param_var |
---|
| 87 | procedure :: set_overlap_param_approx |
---|
| 88 | procedure :: create_fractional_std |
---|
| 89 | procedure :: create_inv_cloud_effective_size |
---|
| 90 | procedure :: create_inv_cloud_effective_size_eta |
---|
| 91 | procedure :: param_cloud_effective_separation_eta |
---|
| 92 | procedure :: crop_cloud_fraction |
---|
| 93 | procedure :: out_of_physical_bounds |
---|
| 94 | |
---|
| 95 | end type cloud_type |
---|
| 96 | |
---|
| 97 | contains |
---|
| 98 | |
---|
| 99 | !--------------------------------------------------------------------- |
---|
| 100 | ! Allocate arrays for describing clouds and precipitation, although |
---|
| 101 | ! in the offline code these are allocated when they are read from |
---|
| 102 | ! the NetCDF file |
---|
| 103 | subroutine allocate_cloud_arrays(this, ncol, nlev, ntype, use_inhom_effective_size) |
---|
| 104 | |
---|
| 105 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 106 | |
---|
| 107 | class(cloud_type), intent(inout), target :: this |
---|
| 108 | integer, intent(in) :: ncol ! Number of columns |
---|
| 109 | integer, intent(in) :: nlev ! Number of levels |
---|
| 110 | ! Number of cloud/precip particle types. If not present then the |
---|
| 111 | ! older cloud behaviour is assumed: two types are present, (1) |
---|
| 112 | ! liquid and (2) ice, and they can be accessed via q_liq, q_ice, |
---|
| 113 | ! re_liq and re_ice. |
---|
| 114 | integer, intent(in), optional :: ntype |
---|
| 115 | logical, intent(in), optional :: use_inhom_effective_size |
---|
| 116 | |
---|
| 117 | real(jphook) :: hook_handle |
---|
| 118 | |
---|
| 119 | if (lhook) call dr_hook('radiation_cloud:allocate',0,hook_handle) |
---|
| 120 | |
---|
| 121 | if (present(ntype)) then |
---|
| 122 | this%ntype = ntype |
---|
| 123 | else |
---|
| 124 | this%ntype = 2 |
---|
| 125 | end if |
---|
| 126 | allocate(this%mixing_ratio(ncol,nlev,this%ntype)) |
---|
| 127 | allocate(this%effective_radius(ncol,nlev,this%ntype)) |
---|
| 128 | nullify(this%q_liq) |
---|
| 129 | nullify(this%q_ice) |
---|
| 130 | nullify(this%re_liq) |
---|
| 131 | nullify(this%re_ice) |
---|
| 132 | if (.not. present(ntype)) then |
---|
| 133 | ! Older interface in which only liquid and ice are supported |
---|
| 134 | this%q_liq => this%mixing_ratio(:,:,1) |
---|
| 135 | this%q_ice => this%mixing_ratio(:,:,2) |
---|
| 136 | this%re_liq => this%effective_radius(:,:,1) |
---|
| 137 | this%re_ice => this%effective_radius(:,:,2) |
---|
| 138 | end if |
---|
| 139 | |
---|
| 140 | allocate(this%fraction(ncol,nlev)) |
---|
| 141 | allocate(this%overlap_param(ncol,nlev-1)) |
---|
| 142 | allocate(this%fractional_std(ncol,nlev)) |
---|
| 143 | allocate(this%inv_cloud_effective_size(ncol,nlev)) |
---|
| 144 | |
---|
| 145 | if (present(use_inhom_effective_size)) then |
---|
| 146 | if (use_inhom_effective_size) then |
---|
| 147 | allocate(this%inv_inhom_effective_size(ncol,nlev)) |
---|
| 148 | end if |
---|
| 149 | end if |
---|
| 150 | |
---|
| 151 | if (lhook) call dr_hook('radiation_cloud:allocate',1,hook_handle) |
---|
| 152 | |
---|
| 153 | end subroutine allocate_cloud_arrays |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | !--------------------------------------------------------------------- |
---|
| 157 | ! Deallocate arrays |
---|
| 158 | subroutine deallocate_cloud_arrays(this) |
---|
| 159 | |
---|
| 160 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 161 | |
---|
| 162 | class(cloud_type), intent(inout) :: this |
---|
| 163 | |
---|
| 164 | real(jphook) :: hook_handle |
---|
| 165 | |
---|
| 166 | if (lhook) call dr_hook('radiation_cloud:deallocate',0,hook_handle) |
---|
| 167 | |
---|
| 168 | nullify(this%q_liq) |
---|
| 169 | nullify(this%q_ice) |
---|
| 170 | nullify(this%re_liq) |
---|
| 171 | nullify(this%re_ice) |
---|
| 172 | |
---|
| 173 | if (allocated(this%mixing_ratio)) deallocate(this%mixing_ratio) |
---|
| 174 | if (allocated(this%effective_radius)) deallocate(this%effective_radius) |
---|
| 175 | if (allocated(this%fraction)) deallocate(this%fraction) |
---|
| 176 | if (allocated(this%overlap_param)) deallocate(this%overlap_param) |
---|
| 177 | if (allocated(this%fractional_std)) deallocate(this%fractional_std) |
---|
| 178 | if (allocated(this%inv_cloud_effective_size)) & |
---|
| 179 | & deallocate(this%inv_cloud_effective_size) |
---|
| 180 | if (allocated(this%inv_inhom_effective_size)) & |
---|
| 181 | & deallocate(this%inv_inhom_effective_size) |
---|
| 182 | |
---|
| 183 | if (lhook) call dr_hook('radiation_cloud:deallocate',1,hook_handle) |
---|
| 184 | |
---|
| 185 | end subroutine deallocate_cloud_arrays |
---|
| 186 | |
---|
| 187 | |
---|
| 188 | !--------------------------------------------------------------------- |
---|
| 189 | ! Compute and store the overlap parameter from the provided overlap |
---|
| 190 | ! decorrelation length (in metres). If istartcol and/or iendcol are |
---|
| 191 | ! provided then only columns in this range are computed. If the |
---|
| 192 | ! overlap_param array has not been allocated then it will be |
---|
| 193 | ! allocated to be of the correct size relative to the pressure |
---|
| 194 | ! field. This version assumes a fixed decorrelation_length for all |
---|
| 195 | ! columns. |
---|
| 196 | subroutine set_overlap_param_fix(this, thermodynamics, decorrelation_length, & |
---|
| 197 | & istartcol, iendcol) |
---|
| 198 | |
---|
| 199 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 200 | use radiation_thermodynamics, only : thermodynamics_type |
---|
| 201 | use radiation_constants, only : GasConstantDryAir, AccelDueToGravity |
---|
| 202 | |
---|
| 203 | class(cloud_type), intent(inout) :: this |
---|
| 204 | type(thermodynamics_type), intent(in) :: thermodynamics |
---|
| 205 | real(jprb), intent(in) :: decorrelation_length ! m |
---|
| 206 | integer, optional, intent(in) :: istartcol, iendcol |
---|
| 207 | |
---|
| 208 | ! Ratio of gas constant for dry air to acceleration due to gravity |
---|
| 209 | real(jprb), parameter :: R_over_g = GasConstantDryAir / AccelDueToGravity |
---|
| 210 | |
---|
| 211 | ! Process only columns i1 to i2, which will be istartcol to |
---|
| 212 | ! iendcol if they were provided |
---|
| 213 | integer :: i1, i2 |
---|
| 214 | |
---|
| 215 | integer :: ncol, nlev |
---|
| 216 | |
---|
| 217 | integer :: jcol, jlev |
---|
| 218 | |
---|
| 219 | real(jphook) :: hook_handle |
---|
| 220 | |
---|
| 221 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param_fix',0,hook_handle) |
---|
| 222 | |
---|
| 223 | ! Pressure at half-levels, pressure_hl, is defined at nlev+1 |
---|
| 224 | ! points |
---|
| 225 | ncol = size(thermodynamics%pressure_hl,dim=1) |
---|
| 226 | nlev = size(thermodynamics%pressure_hl,dim=2)-1 |
---|
| 227 | |
---|
| 228 | if (present(istartcol)) then |
---|
| 229 | i1 = istartcol |
---|
| 230 | else |
---|
| 231 | i1 = 1 |
---|
| 232 | end if |
---|
| 233 | |
---|
| 234 | if (present(iendcol)) then |
---|
| 235 | i2 = iendcol |
---|
| 236 | else |
---|
| 237 | i2 = ncol |
---|
| 238 | end if |
---|
| 239 | |
---|
| 240 | if (.not. allocated(this%overlap_param)) then |
---|
| 241 | ! If pressure is of size (ncol,nlev+1) then overlap_param is of |
---|
| 242 | ! size (ncol,nlev-1), since overlap parameter is only defined here |
---|
| 243 | ! for interfaces between model layers, not for the interface to |
---|
| 244 | ! space or the surface |
---|
| 245 | allocate(this%overlap_param(ncol, nlev-1)) |
---|
| 246 | end if |
---|
| 247 | |
---|
| 248 | if (thermodynamics%pressure_hl(i1,2) > thermodynamics%pressure_hl(i1,1)) then |
---|
| 249 | ! Pressure is increasing with index (order of layers is |
---|
| 250 | ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we |
---|
| 251 | ! don't take the logarithm of the first pressure in each column. |
---|
| 252 | do jcol = i1,i2 |
---|
| 253 | this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length) & |
---|
| 254 | & * thermodynamics%temperature_hl(jcol,2) & |
---|
| 255 | & *log(thermodynamics%pressure_hl(jcol,3) & |
---|
| 256 | & /thermodynamics%pressure_hl(jcol,2))) |
---|
| 257 | end do |
---|
| 258 | |
---|
| 259 | do jlev = 2,nlev-1 |
---|
| 260 | do jcol = i1,i2 |
---|
| 261 | this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) & |
---|
| 262 | & * thermodynamics%temperature_hl(jcol,jlev+1) & |
---|
| 263 | & *log(thermodynamics%pressure_hl(jcol,jlev+2) & |
---|
| 264 | & /thermodynamics%pressure_hl(jcol,jlev))) |
---|
| 265 | end do |
---|
| 266 | end do |
---|
| 267 | |
---|
| 268 | else |
---|
| 269 | ! Pressure is decreasing with index (order of layers is surface |
---|
| 270 | ! to top-of-atmosphere). In case pressure_hl(:,nlev+1)=0, we |
---|
| 271 | ! don't take the logarithm of the last pressure in each column. |
---|
| 272 | do jlev = 1,nlev-2 |
---|
| 273 | do jcol = i1,i2 |
---|
| 274 | this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) & |
---|
| 275 | & * thermodynamics%temperature_hl(jcol,jlev+1) & |
---|
| 276 | & *log(thermodynamics%pressure_hl(jcol,jlev) & |
---|
| 277 | & /thermodynamics%pressure_hl(jcol,jlev+2))) |
---|
| 278 | end do |
---|
| 279 | end do |
---|
| 280 | |
---|
| 281 | do jcol = i1,i2 |
---|
| 282 | this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length) & |
---|
| 283 | & * thermodynamics%temperature_hl(jcol,nlev) & |
---|
| 284 | & *log(thermodynamics%pressure_hl(jcol,nlev-1) & |
---|
| 285 | & /thermodynamics%pressure_hl(jcol,nlev))) |
---|
| 286 | end do |
---|
| 287 | end if |
---|
| 288 | |
---|
| 289 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param_fix',1,hook_handle) |
---|
| 290 | |
---|
| 291 | end subroutine set_overlap_param_fix |
---|
| 292 | |
---|
| 293 | |
---|
| 294 | !--------------------------------------------------------------------- |
---|
| 295 | ! Compute and store the overlap parameter from the provided overlap |
---|
| 296 | ! decorrelation length (in metres), which may vary with column. Only |
---|
| 297 | ! columns from istartcol to iendcol are computed. If the |
---|
| 298 | ! overlap_param array has not been allocated then it will be |
---|
| 299 | ! allocated to be of the correct size relative to the pressure |
---|
| 300 | ! field. |
---|
| 301 | subroutine set_overlap_param_var(this, thermodynamics, decorrelation_length, & |
---|
| 302 | & istartcol, iendcol) |
---|
| 303 | |
---|
| 304 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 305 | use radiation_thermodynamics, only : thermodynamics_type |
---|
| 306 | use radiation_constants, only : GasConstantDryAir, AccelDueToGravity |
---|
| 307 | |
---|
| 308 | class(cloud_type), intent(inout) :: this |
---|
| 309 | type(thermodynamics_type), intent(in) :: thermodynamics |
---|
| 310 | integer, intent(in) :: istartcol, iendcol |
---|
| 311 | real(jprb), intent(in) :: decorrelation_length(istartcol:iendcol) ! m |
---|
| 312 | |
---|
| 313 | ! Ratio of gas constant for dry air to acceleration due to gravity |
---|
| 314 | real(jprb), parameter :: R_over_g = GasConstantDryAir / AccelDueToGravity |
---|
| 315 | |
---|
| 316 | integer :: ncol, nlev |
---|
| 317 | |
---|
| 318 | integer :: jcol, jlev |
---|
| 319 | |
---|
| 320 | real(jphook) :: hook_handle |
---|
| 321 | |
---|
| 322 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param_var',0,hook_handle) |
---|
| 323 | |
---|
| 324 | ! Pressure at half-levels, pressure_hl, is defined at nlev+1 |
---|
| 325 | ! points |
---|
| 326 | ncol = size(thermodynamics%pressure_hl,dim=1) |
---|
| 327 | nlev = size(thermodynamics%pressure_hl,dim=2)-1 |
---|
| 328 | |
---|
| 329 | if (.not. allocated(this%overlap_param)) then |
---|
| 330 | ! If pressure is of size (ncol,nlev+1) then overlap_param is of |
---|
| 331 | ! size (ncol,nlev-1), since overlap parameter is only defined here |
---|
| 332 | ! for interfaces between model layers, not for the interface to |
---|
| 333 | ! space or the surface |
---|
| 334 | allocate(this%overlap_param(ncol, nlev-1)) |
---|
| 335 | end if |
---|
| 336 | |
---|
| 337 | if (thermodynamics%pressure_hl(istartcol,2) > thermodynamics%pressure_hl(istartcol,1)) then |
---|
| 338 | ! Pressure is increasing with index (order of layers is |
---|
| 339 | ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we |
---|
| 340 | ! don't take the logarithm of the first pressure in each column. |
---|
| 341 | do jcol = istartcol,iendcol |
---|
| 342 | this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length(jcol)) & |
---|
| 343 | & * thermodynamics%temperature_hl(jcol,2) & |
---|
| 344 | & *log(thermodynamics%pressure_hl(jcol,3) & |
---|
| 345 | & /thermodynamics%pressure_hl(jcol,2))) |
---|
| 346 | end do |
---|
| 347 | |
---|
| 348 | do jlev = 2,nlev-1 |
---|
| 349 | do jcol = istartcol,iendcol |
---|
| 350 | this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol)) & |
---|
| 351 | & * thermodynamics%temperature_hl(jcol,jlev+1) & |
---|
| 352 | & *log(thermodynamics%pressure_hl(jcol,jlev+2) & |
---|
| 353 | & /thermodynamics%pressure_hl(jcol,jlev))) |
---|
| 354 | end do |
---|
| 355 | end do |
---|
| 356 | |
---|
| 357 | else |
---|
| 358 | ! Pressure is decreasing with index (order of layers is surface |
---|
| 359 | ! to top-of-atmosphere). In case pressure_hl(:,nlev+1)=0, we |
---|
| 360 | ! don't take the logarithm of the last pressure in each column. |
---|
| 361 | do jlev = 1,nlev-2 |
---|
| 362 | do jcol = istartcol,iendcol |
---|
| 363 | this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol)) & |
---|
| 364 | & * thermodynamics%temperature_hl(jcol,jlev+1) & |
---|
| 365 | & *log(thermodynamics%pressure_hl(jcol,jlev) & |
---|
| 366 | & /thermodynamics%pressure_hl(jcol,jlev+2))) |
---|
| 367 | end do |
---|
| 368 | end do |
---|
| 369 | |
---|
| 370 | do jcol = istartcol,iendcol |
---|
| 371 | this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length(jcol)) & |
---|
| 372 | & * thermodynamics%temperature_hl(jcol,nlev) & |
---|
| 373 | & *log(thermodynamics%pressure_hl(jcol,nlev-1) & |
---|
| 374 | & /thermodynamics%pressure_hl(jcol,nlev))) |
---|
| 375 | end do |
---|
| 376 | end if |
---|
| 377 | |
---|
| 378 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param_var',1,hook_handle) |
---|
| 379 | |
---|
| 380 | end subroutine set_overlap_param_var |
---|
| 381 | |
---|
| 382 | |
---|
| 383 | !--------------------------------------------------------------------- |
---|
| 384 | ! Compute and store the overlap parameter from the provided overlap |
---|
| 385 | ! decorrelation length (in metres). If istartcol and/or iendcol are |
---|
| 386 | ! provided then only columns in this range are computed. If the |
---|
| 387 | ! overlap_param array has not been allocated then it will be |
---|
| 388 | ! allocated to be of the correct size relative to the pressure |
---|
| 389 | ! field. This is the APPROXIMATE method as it assumes a fixed |
---|
| 390 | ! atmospheric scale height, which leads to differences particularly |
---|
| 391 | ! in low cloud. |
---|
| 392 | subroutine set_overlap_param_approx(this, thermodynamics, decorrelation_length, & |
---|
| 393 | & istartcol, iendcol) |
---|
| 394 | |
---|
| 395 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 396 | use radiation_thermodynamics, only : thermodynamics_type |
---|
| 397 | |
---|
| 398 | class(cloud_type), intent(inout) :: this |
---|
| 399 | type(thermodynamics_type), intent(in) :: thermodynamics |
---|
| 400 | real(jprb), intent(in) :: decorrelation_length ! m |
---|
| 401 | integer, optional, intent(in) :: istartcol, iendcol |
---|
| 402 | |
---|
| 403 | ! To convert decorrelation length (m) to overlap parameter between |
---|
| 404 | ! layers, we need an estimate for the thickness of the layer. This |
---|
| 405 | ! is found using the pressure difference between the edges of the |
---|
| 406 | ! layer, along with the approximate scale height of the atmosphere |
---|
| 407 | ! (m) given here: |
---|
| 408 | real(jprb), parameter :: scale_height = 8000.0_jprb |
---|
| 409 | |
---|
| 410 | ! Process only columns i1 to i2, which will be istartcol to |
---|
| 411 | ! iendcol if they were provided |
---|
| 412 | integer :: i1, i2 |
---|
| 413 | |
---|
| 414 | integer :: ncol, nlev |
---|
| 415 | |
---|
| 416 | real(jphook) :: hook_handle |
---|
| 417 | |
---|
| 418 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param_approx',0,hook_handle) |
---|
| 419 | |
---|
| 420 | ! Pressure at half-levels, pressure_hl, is defined at nlev+1 |
---|
| 421 | ! points |
---|
| 422 | ncol = size(thermodynamics%pressure_hl,dim=1) |
---|
| 423 | nlev = size(thermodynamics%pressure_hl,dim=2)-1 |
---|
| 424 | |
---|
| 425 | if (present(istartcol)) then |
---|
| 426 | i1 = istartcol |
---|
| 427 | else |
---|
| 428 | i1 = 1 |
---|
| 429 | end if |
---|
| 430 | |
---|
| 431 | if (present(iendcol)) then |
---|
| 432 | i2 = iendcol |
---|
| 433 | else |
---|
| 434 | i2 = ncol |
---|
| 435 | end if |
---|
| 436 | |
---|
| 437 | if (.not. allocated(this%overlap_param)) then |
---|
| 438 | ! If pressure is of size (ncol,nlev+1) then overlap_param is of |
---|
| 439 | ! size (ncol,nlev-1), since overlap parameter is only defined here |
---|
| 440 | ! for interfaces between model layers, not for the interface to |
---|
| 441 | ! space or the surface |
---|
| 442 | allocate(this%overlap_param(ncol, nlev-1)) |
---|
| 443 | end if |
---|
| 444 | |
---|
| 445 | if (thermodynamics%pressure_hl(i1,2) > thermodynamics%pressure_hl(i1,1)) then |
---|
| 446 | ! Pressure is increasing with index (order of layers is |
---|
| 447 | ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we |
---|
| 448 | ! don't take the logarithm of the first pressure in each |
---|
| 449 | ! column. |
---|
| 450 | this%overlap_param(i1:i2,:) = exp(-(scale_height/decorrelation_length) & |
---|
| 451 | & * ( log(thermodynamics%pressure_hl(i1:i2,3:nlev+1) & |
---|
| 452 | & /thermodynamics%pressure_hl(i1:i2,2:nlev )) ) ) |
---|
| 453 | else |
---|
| 454 | ! Pressure is decreasing with index (order of layers is surface |
---|
| 455 | ! to top-of-atmosphere). In case pressure_hl(:,nlev+1)=0, we |
---|
| 456 | ! don't take the logarithm of the last pressure in each column. |
---|
| 457 | this%overlap_param(i1:i2,:) = exp(-(scale_height/decorrelation_length) & |
---|
| 458 | & * ( log(thermodynamics%pressure_hl(i1:i2,1:nlev-1) & |
---|
| 459 | & /thermodynamics%pressure_hl(i1:i2,2:nlev )) ) ) |
---|
| 460 | end if |
---|
| 461 | |
---|
| 462 | if (lhook) call dr_hook('radiation_cloud:set_overlap_param_approx',1,hook_handle) |
---|
| 463 | |
---|
| 464 | end subroutine set_overlap_param_approx |
---|
| 465 | |
---|
| 466 | |
---|
| 467 | !--------------------------------------------------------------------- |
---|
| 468 | ! Create a matrix of constant fractional standard deviations |
---|
| 469 | ! (dimensionless) |
---|
| 470 | subroutine create_fractional_std(this, ncol, nlev, frac_std) |
---|
| 471 | |
---|
| 472 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 473 | |
---|
| 474 | class(cloud_type), intent(inout) :: this |
---|
| 475 | integer, intent(in) :: ncol, nlev |
---|
| 476 | real(jprb), intent(in) :: frac_std |
---|
| 477 | |
---|
| 478 | real(jphook) :: hook_handle |
---|
| 479 | |
---|
| 480 | if (lhook) call dr_hook('radiation_cloud:create_fractional_std',0,hook_handle) |
---|
| 481 | |
---|
| 482 | if (allocated(this%fractional_std)) then |
---|
| 483 | deallocate(this%fractional_std) |
---|
| 484 | end if |
---|
| 485 | |
---|
| 486 | allocate(this%fractional_std(ncol, nlev)) |
---|
| 487 | |
---|
| 488 | this%fractional_std = frac_std |
---|
| 489 | |
---|
| 490 | if (lhook) call dr_hook('radiation_cloud:create_fractional_std',1,hook_handle) |
---|
| 491 | |
---|
| 492 | end subroutine create_fractional_std |
---|
| 493 | |
---|
| 494 | |
---|
| 495 | !--------------------------------------------------------------------- |
---|
| 496 | ! Create a matrix of constant inverse cloud effective size (m-1) |
---|
| 497 | subroutine create_inv_cloud_effective_size(this, ncol, nlev, inv_eff_size) |
---|
| 498 | |
---|
| 499 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 500 | |
---|
| 501 | class(cloud_type), intent(inout) :: this |
---|
| 502 | integer, intent(in) :: ncol, nlev |
---|
| 503 | real(jprb), intent(in) :: inv_eff_size |
---|
| 504 | |
---|
| 505 | real(jphook) :: hook_handle |
---|
| 506 | |
---|
| 507 | if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size',0,hook_handle) |
---|
| 508 | |
---|
| 509 | if (allocated(this%inv_cloud_effective_size)) then |
---|
| 510 | deallocate(this%inv_cloud_effective_size) |
---|
| 511 | end if |
---|
| 512 | |
---|
| 513 | allocate(this%inv_cloud_effective_size(ncol, nlev)) |
---|
| 514 | |
---|
| 515 | this%inv_cloud_effective_size = inv_eff_size |
---|
| 516 | |
---|
| 517 | if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size',1,hook_handle) |
---|
| 518 | |
---|
| 519 | end subroutine create_inv_cloud_effective_size |
---|
| 520 | |
---|
| 521 | |
---|
| 522 | !--------------------------------------------------------------------- |
---|
| 523 | ! Create a matrix of inverse cloud effective size (m-1) according to |
---|
| 524 | ! the value of eta (=pressure divided by surface pressure) |
---|
| 525 | subroutine create_inv_cloud_effective_size_eta(this, ncol, nlev, & |
---|
| 526 | & pressure_hl, inv_eff_size_low, inv_eff_size_mid, inv_eff_size_high, & |
---|
| 527 | & eta_low_mid, eta_mid_high, istartcol, iendcol) |
---|
| 528 | |
---|
| 529 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 530 | |
---|
| 531 | class(cloud_type), intent(inout) :: this |
---|
| 532 | integer, intent(in) :: ncol, nlev |
---|
| 533 | ! Pressure on half levels (Pa) |
---|
| 534 | real(jprb), intent(in) :: pressure_hl(:,:) |
---|
| 535 | ! Inverse effective size for low, mid and high cloud (m-1) |
---|
| 536 | real(jprb), intent(in) :: inv_eff_size_low |
---|
| 537 | real(jprb), intent(in) :: inv_eff_size_mid |
---|
| 538 | real(jprb), intent(in) :: inv_eff_size_high |
---|
| 539 | ! Eta values at low-mid and mid-high interfaces |
---|
| 540 | real(jprb), intent(in) :: eta_low_mid, eta_mid_high |
---|
| 541 | integer, optional, intent(in) :: istartcol, iendcol |
---|
| 542 | |
---|
| 543 | ! Ratio of layer midpoint pressure to surface pressure |
---|
| 544 | real(jprb) :: eta(nlev) |
---|
| 545 | |
---|
| 546 | ! Indices of column, level and surface half-level |
---|
| 547 | integer :: jcol, isurf |
---|
| 548 | |
---|
| 549 | ! Local values of istartcol, iendcol |
---|
| 550 | integer :: i1, i2 |
---|
| 551 | |
---|
| 552 | real(jphook) :: hook_handle |
---|
| 553 | |
---|
| 554 | if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size_eta',0,hook_handle) |
---|
| 555 | |
---|
| 556 | if (allocated(this%inv_cloud_effective_size)) then |
---|
| 557 | deallocate(this%inv_cloud_effective_size) |
---|
| 558 | end if |
---|
| 559 | |
---|
| 560 | allocate(this%inv_cloud_effective_size(ncol, nlev)) |
---|
| 561 | |
---|
| 562 | if (present(istartcol)) then |
---|
| 563 | i1 = istartcol |
---|
| 564 | else |
---|
| 565 | i1 = 1 |
---|
| 566 | end if |
---|
| 567 | |
---|
| 568 | if (present(iendcol)) then |
---|
| 569 | i2 = iendcol |
---|
| 570 | else |
---|
| 571 | i2 = ncol |
---|
| 572 | end if |
---|
| 573 | |
---|
| 574 | ! Locate the surface half-level |
---|
| 575 | if (pressure_hl(1,1) > pressure_hl(1,2)) then |
---|
| 576 | isurf = 1 |
---|
| 577 | else |
---|
| 578 | isurf = nlev+1 |
---|
| 579 | end if |
---|
| 580 | |
---|
| 581 | do jcol = i1,i2 |
---|
| 582 | eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) & |
---|
| 583 | & * (0.5_jprb / pressure_hl(jcol,isurf)) |
---|
| 584 | where (eta > eta_low_mid) |
---|
| 585 | this%inv_cloud_effective_size(jcol,:) = inv_eff_size_low |
---|
| 586 | elsewhere (eta > eta_mid_high) |
---|
| 587 | this%inv_cloud_effective_size(jcol,:) = inv_eff_size_mid |
---|
| 588 | elsewhere |
---|
| 589 | this%inv_cloud_effective_size(jcol,:) = inv_eff_size_high |
---|
| 590 | end where |
---|
| 591 | end do |
---|
| 592 | |
---|
| 593 | if (lhook) call dr_hook('radiation_cloud:create_inv_cloud_effective_size_eta',1,hook_handle) |
---|
| 594 | |
---|
| 595 | end subroutine create_inv_cloud_effective_size_eta |
---|
| 596 | |
---|
| 597 | |
---|
| 598 | !--------------------------------------------------------------------- |
---|
| 599 | ! Create a matrix of inverse cloud and inhomogeneity effective size |
---|
| 600 | ! (m-1) parameterized according to the value of eta (=pressure |
---|
| 601 | ! divided by surface pressure): effective_separation = |
---|
| 602 | ! coeff_a + coeff_b*exp(-(eta**power)). |
---|
| 603 | subroutine param_cloud_effective_separation_eta(this, ncol, nlev, & |
---|
| 604 | & pressure_hl, separation_surf, separation_toa, power, & |
---|
| 605 | & inhom_separation_factor, istartcol, iendcol) |
---|
| 606 | |
---|
| 607 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 608 | |
---|
| 609 | class(cloud_type), intent(inout) :: this |
---|
| 610 | integer, intent(in) :: ncol, nlev |
---|
| 611 | ! Pressure on half levels (Pa) |
---|
| 612 | real(jprb), intent(in) :: pressure_hl(:,:) |
---|
| 613 | ! Separation distances at surface and top-of-atmosphere, and power |
---|
| 614 | ! on eta |
---|
| 615 | real(jprb), intent(in) :: separation_surf ! m |
---|
| 616 | real(jprb), intent(in) :: separation_toa ! m |
---|
| 617 | real(jprb), intent(in) :: power |
---|
| 618 | real(jprb), optional, intent(in) :: inhom_separation_factor |
---|
| 619 | integer, optional, intent(in) :: istartcol, iendcol |
---|
| 620 | |
---|
| 621 | ! Ratio of layer midpoint pressure to surface pressure |
---|
| 622 | real(jprb) :: eta(nlev) |
---|
| 623 | |
---|
| 624 | ! Effective cloud separation (m) |
---|
| 625 | real(jprb) :: eff_separation(nlev) |
---|
| 626 | |
---|
| 627 | ! Coefficients used to compute effective separation distance |
---|
| 628 | real(jprb) :: coeff_e, coeff_a, coeff_b, inhom_sep_factor |
---|
| 629 | |
---|
| 630 | ! Indices of column, level and surface half-level |
---|
| 631 | integer :: jcol, isurf |
---|
| 632 | |
---|
| 633 | ! Local values of istartcol, iendcol |
---|
| 634 | integer :: i1, i2 |
---|
| 635 | |
---|
| 636 | real(jphook) :: hook_handle |
---|
| 637 | |
---|
| 638 | if (lhook) call dr_hook('radiation_cloud:param_cloud_effective_separation_eta',0,hook_handle) |
---|
| 639 | |
---|
| 640 | if (present(inhom_separation_factor)) then |
---|
| 641 | inhom_sep_factor = inhom_separation_factor |
---|
| 642 | else |
---|
| 643 | inhom_sep_factor = 1.0_jprb |
---|
| 644 | end if |
---|
| 645 | |
---|
| 646 | coeff_e = 1.0_jprb - exp(-1.0_jprb) |
---|
| 647 | coeff_b = (separation_toa - separation_surf) / coeff_e |
---|
| 648 | coeff_a = separation_toa - coeff_b |
---|
| 649 | |
---|
| 650 | if (allocated(this%inv_cloud_effective_size)) then |
---|
| 651 | deallocate(this%inv_cloud_effective_size) |
---|
| 652 | end if |
---|
| 653 | if (allocated(this%inv_inhom_effective_size)) then |
---|
| 654 | deallocate(this%inv_inhom_effective_size) |
---|
| 655 | end if |
---|
| 656 | |
---|
| 657 | allocate(this%inv_cloud_effective_size(ncol, nlev)) |
---|
| 658 | allocate(this%inv_inhom_effective_size(ncol, nlev)) |
---|
| 659 | |
---|
| 660 | if (present(istartcol)) then |
---|
| 661 | i1 = istartcol |
---|
| 662 | else |
---|
| 663 | i1 = 1 |
---|
| 664 | end if |
---|
| 665 | |
---|
| 666 | if (present(iendcol)) then |
---|
| 667 | i2 = iendcol |
---|
| 668 | else |
---|
| 669 | i2 = ncol |
---|
| 670 | end if |
---|
| 671 | |
---|
| 672 | ! Locate the surface half-level |
---|
| 673 | if (pressure_hl(1,1) > pressure_hl(1,2)) then |
---|
| 674 | isurf = 1 |
---|
| 675 | else |
---|
| 676 | isurf = nlev+1 |
---|
| 677 | end if |
---|
| 678 | |
---|
| 679 | do jcol = i1,i2 |
---|
| 680 | eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) & |
---|
| 681 | & * (0.5_jprb / pressure_hl(jcol,isurf)) |
---|
| 682 | eff_separation = coeff_a + coeff_b * exp(-eta**power) |
---|
| 683 | this%inv_cloud_effective_size(jcol,:) = 1.0_jprb / (eff_separation & |
---|
| 684 | & * sqrt(max(1.0e-5_jprb,this%fraction(jcol,:)*(1.0_jprb-this%fraction(jcol,:))))) |
---|
| 685 | this%inv_inhom_effective_size(jcol,:) = 1.0_jprb / (eff_separation * inhom_sep_factor & |
---|
| 686 | & * sqrt(max(1.0e-5_jprb,0.5_jprb*this%fraction(jcol,:)*(1.0_jprb-0.5_jprb*this%fraction(jcol,:))))) |
---|
| 687 | end do |
---|
| 688 | |
---|
| 689 | if (lhook) call dr_hook('radiation_cloud:param_cloud_effective_separation_eta',1,hook_handle) |
---|
| 690 | |
---|
| 691 | end subroutine param_cloud_effective_separation_eta |
---|
| 692 | |
---|
| 693 | |
---|
| 694 | !--------------------------------------------------------------------- |
---|
| 695 | ! Remove "ghost" clouds: those with a cloud fraction that is too |
---|
| 696 | ! small to treat sensibly (e.g. because it implies that the |
---|
| 697 | ! "in-cloud" water content is too high), or with a cloud water |
---|
| 698 | ! content that is too small. We do this in one place to ensure that |
---|
| 699 | ! all subsequent subroutines can assume that if cloud_fraction > 0.0 |
---|
| 700 | ! then cloud is really present and should be treated. |
---|
| 701 | subroutine crop_cloud_fraction(this, istartcol, iendcol, & |
---|
| 702 | & cloud_fraction_threshold, cloud_mixing_ratio_threshold) |
---|
| 703 | |
---|
| 704 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 705 | |
---|
| 706 | class(cloud_type), intent(inout) :: this |
---|
| 707 | integer, intent(in) :: istartcol, iendcol |
---|
| 708 | |
---|
| 709 | integer :: nlev, ntype |
---|
| 710 | integer :: jcol, jlev, jh |
---|
| 711 | |
---|
| 712 | real(jprb) :: cloud_fraction_threshold, cloud_mixing_ratio_threshold |
---|
| 713 | real(jprb) :: sum_mixing_ratio(istartcol:iendcol) |
---|
| 714 | |
---|
| 715 | real(jphook) :: hook_handle |
---|
| 716 | |
---|
| 717 | if (lhook) call dr_hook('radiation_cloud:crop_cloud_fraction',0,hook_handle) |
---|
| 718 | |
---|
| 719 | nlev = size(this%fraction,2) |
---|
| 720 | ntype = size(this%mixing_ratio,3) |
---|
| 721 | |
---|
| 722 | do jlev = 1,nlev |
---|
| 723 | do jcol = istartcol,iendcol |
---|
| 724 | sum_mixing_ratio(jcol) = 0.0_jprb |
---|
| 725 | end do |
---|
| 726 | do jh = 1, ntype |
---|
| 727 | do jcol = istartcol,iendcol |
---|
| 728 | sum_mixing_ratio(jcol) = sum_mixing_ratio(jcol) + this%mixing_ratio(jcol,jlev,jh) |
---|
| 729 | end do |
---|
| 730 | end do |
---|
| 731 | do jcol = istartcol,iendcol |
---|
| 732 | if (this%fraction(jcol,jlev) < cloud_fraction_threshold & |
---|
| 733 | & .or. sum_mixing_ratio(jcol) < cloud_mixing_ratio_threshold) then |
---|
| 734 | this%fraction(jcol,jlev) = 0.0_jprb |
---|
| 735 | end if |
---|
| 736 | end do |
---|
| 737 | end do |
---|
| 738 | |
---|
| 739 | if (lhook) call dr_hook('radiation_cloud:crop_cloud_fraction',1,hook_handle) |
---|
| 740 | |
---|
| 741 | end subroutine crop_cloud_fraction |
---|
| 742 | |
---|
| 743 | |
---|
| 744 | !--------------------------------------------------------------------- |
---|
| 745 | ! Return .true. if variables are out of a physically sensible range, |
---|
| 746 | ! optionally only considering columns between istartcol and iendcol |
---|
| 747 | function out_of_physical_bounds(this, istartcol, iendcol, do_fix) result(is_bad) |
---|
| 748 | |
---|
| 749 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 750 | use radiation_check, only : out_of_bounds_2d, out_of_bounds_3d |
---|
| 751 | |
---|
| 752 | class(cloud_type), intent(inout) :: this |
---|
| 753 | integer, optional,intent(in) :: istartcol, iendcol |
---|
| 754 | logical, optional,intent(in) :: do_fix |
---|
| 755 | logical :: is_bad |
---|
| 756 | |
---|
| 757 | logical :: do_fix_local |
---|
| 758 | |
---|
| 759 | real(jphook) :: hook_handle |
---|
| 760 | |
---|
| 761 | if (lhook) call dr_hook('radiation_cloud:out_of_physical_bounds',0,hook_handle) |
---|
| 762 | |
---|
| 763 | if (present(do_fix)) then |
---|
| 764 | do_fix_local = do_fix |
---|
| 765 | else |
---|
| 766 | do_fix_local = .false. |
---|
| 767 | end if |
---|
| 768 | |
---|
| 769 | is_bad = out_of_bounds_3d(this%mixing_ratio, 'cloud%mixing_ratio', 0.0_jprb, 1.0_jprb, & |
---|
| 770 | & do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 771 | & .or. out_of_bounds_3d(this%effective_radius, 'cloud%effective_radius', 0.0_jprb, 0.1_jprb, & |
---|
| 772 | & do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 773 | & .or. out_of_bounds_2d(this%fraction, 'cloud%fraction', 0.0_jprb, 1.0_jprb, & |
---|
| 774 | & do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 775 | & .or. out_of_bounds_2d(this%fractional_std, 'fractional_std', 0.0_jprb, 10.0_jprb, & |
---|
| 776 | & do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 777 | & .or. out_of_bounds_2d(this%inv_cloud_effective_size, 'inv_cloud_effective_size', & |
---|
| 778 | & 0.0_jprb, 1.0_jprb, do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 779 | & .or. out_of_bounds_2d(this%inv_inhom_effective_size, 'inv_inhom_effective_size', & |
---|
| 780 | & 0.0_jprb, 1.0_jprb, do_fix_local, i1=istartcol, i2=iendcol) & |
---|
| 781 | & .or. out_of_bounds_2d(this%overlap_param, 'overlap_param', -0.5_jprb, 1.0_jprb, & |
---|
| 782 | & do_fix_local, i1=istartcol, i2=iendcol) |
---|
| 783 | |
---|
| 784 | if (lhook) call dr_hook('radiation_cloud:out_of_physical_bounds',1,hook_handle) |
---|
| 785 | |
---|
| 786 | end function out_of_physical_bounds |
---|
| 787 | |
---|
| 788 | end module radiation_cloud |
---|