| 1 | ! radiation_pdf_sampler.F90 - Get samples from a PDF for McICA |
|---|
| 2 | ! |
|---|
| 3 | ! (C) Copyright 2015- ECMWF. |
|---|
| 4 | ! |
|---|
| 5 | ! This software is licensed under the terms of the Apache Licence Version 2.0 |
|---|
| 6 | ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
|---|
| 7 | ! |
|---|
| 8 | ! In applying this licence, ECMWF does not waive the privileges and immunities |
|---|
| 9 | ! granted to it by virtue of its status as an intergovernmental organisation |
|---|
| 10 | ! nor does it submit to any jurisdiction. |
|---|
| 11 | ! |
|---|
| 12 | ! Author: Robin Hogan |
|---|
| 13 | ! Email: r.j.hogan@ecmwf.int |
|---|
| 14 | ! |
|---|
| 15 | |
|---|
| 16 | module radiation_pdf_sampler |
|---|
| 17 | |
|---|
| 18 | use parkind1, only : jprb |
|---|
| 19 | |
|---|
| 20 | implicit none |
|---|
| 21 | public |
|---|
| 22 | |
|---|
| 23 | !--------------------------------------------------------------------- |
|---|
| 24 | ! Derived type for sampling from a lognormal or gamma distribution, |
|---|
| 25 | ! or other PDF, used to generate water content or optical depth |
|---|
| 26 | ! scalings for use in the Monte Carlo Independent Column |
|---|
| 27 | ! Approximation (McICA) |
|---|
| 28 | type pdf_sampler_type |
|---|
| 29 | ! Number of points in look-up table for cumulative distribution |
|---|
| 30 | ! function (CDF) and fractional standard deviation (FSD) |
|---|
| 31 | ! dimensions |
|---|
| 32 | integer :: ncdf, nfsd |
|---|
| 33 | |
|---|
| 34 | ! First value of FSD and the reciprocal of the interval between |
|---|
| 35 | ! FSD values (which are assumed to be uniformly distributed) |
|---|
| 36 | real(jprb) :: fsd1, inv_fsd_interval |
|---|
| 37 | |
|---|
| 38 | ! Value of the distribution for each CDF and FSD bin |
|---|
| 39 | real(jprb), allocatable, dimension(:,:) :: val |
|---|
| 40 | |
|---|
| 41 | contains |
|---|
| 42 | |
|---|
| 43 | procedure :: setup => setup_pdf_sampler |
|---|
| 44 | procedure :: sample => sample_from_pdf |
|---|
| 45 | procedure :: masked_sample => sample_from_pdf_masked |
|---|
| 46 | procedure :: block_sample => sample_from_pdf_block |
|---|
| 47 | procedure :: masked_block_sample => sample_from_pdf_masked_block |
|---|
| 48 | procedure :: deallocate => deallocate_pdf_sampler |
|---|
| 49 | |
|---|
| 50 | end type pdf_sampler_type |
|---|
| 51 | |
|---|
| 52 | contains |
|---|
| 53 | |
|---|
| 54 | !--------------------------------------------------------------------- |
|---|
| 55 | ! Load look-up table from a file |
|---|
| 56 | subroutine setup_pdf_sampler(this, file_name, iverbose) |
|---|
| 57 | |
|---|
| 58 | use yomhook, only : lhook, dr_hook, jphook |
|---|
| 59 | use easy_netcdf, only : netcdf_file |
|---|
| 60 | |
|---|
| 61 | class(pdf_sampler_type), intent(inout) :: this |
|---|
| 62 | character(len=*), intent(in) :: file_name |
|---|
| 63 | integer, optional, intent(in) :: iverbose |
|---|
| 64 | |
|---|
| 65 | type(netcdf_file) :: file |
|---|
| 66 | integer :: iverb |
|---|
| 67 | real(jprb), allocatable :: fsd(:) |
|---|
| 68 | |
|---|
| 69 | real(jphook) :: hook_handle |
|---|
| 70 | |
|---|
| 71 | if (lhook) call dr_hook('radiation_pdf_sampler:setup',0,hook_handle) |
|---|
| 72 | |
|---|
| 73 | if (present(iverbose)) then |
|---|
| 74 | iverb = iverbose |
|---|
| 75 | else |
|---|
| 76 | iverb = 2 |
|---|
| 77 | end if |
|---|
| 78 | |
|---|
| 79 | if (allocated(this%val)) then |
|---|
| 80 | deallocate(this%val) |
|---|
| 81 | end if |
|---|
| 82 | |
|---|
| 83 | call file%open(trim(file_name), iverbose=iverb) |
|---|
| 84 | |
|---|
| 85 | call file%get('fsd',fsd) |
|---|
| 86 | call file%get('x', this%val) |
|---|
| 87 | |
|---|
| 88 | call file%close() |
|---|
| 89 | |
|---|
| 90 | this%ncdf = size(this%val,1) |
|---|
| 91 | this%nfsd = size(this%val,2) |
|---|
| 92 | this%fsd1 = fsd(1) |
|---|
| 93 | this%inv_fsd_interval = 1.0_jprb / (fsd(2)-fsd(1)) |
|---|
| 94 | |
|---|
| 95 | deallocate(fsd) |
|---|
| 96 | |
|---|
| 97 | if (lhook) call dr_hook('radiation_pdf_sampler:setup',1,hook_handle) |
|---|
| 98 | |
|---|
| 99 | end subroutine setup_pdf_sampler |
|---|
| 100 | |
|---|
| 101 | !--------------------------------------------------------------------- |
|---|
| 102 | ! Deallocate data in pdf_sampler_type derived type |
|---|
| 103 | subroutine deallocate_pdf_sampler(this) |
|---|
| 104 | |
|---|
| 105 | use yomhook, only : lhook, dr_hook, jphook |
|---|
| 106 | |
|---|
| 107 | class(pdf_sampler_type), intent(inout) :: this |
|---|
| 108 | real(jphook) :: hook_handle |
|---|
| 109 | |
|---|
| 110 | if (lhook) call dr_hook('radiation_pdf_sampler:deallocate',0,hook_handle) |
|---|
| 111 | |
|---|
| 112 | if (allocated(this%val)) then |
|---|
| 113 | deallocate(this%val) |
|---|
| 114 | end if |
|---|
| 115 | |
|---|
| 116 | if (lhook) call dr_hook('radiation_pdf_sampler:deallocate',1,hook_handle) |
|---|
| 117 | |
|---|
| 118 | end subroutine deallocate_pdf_sampler |
|---|
| 119 | |
|---|
| 120 | |
|---|
| 121 | !--------------------------------------------------------------------- |
|---|
| 122 | ! Extract the value from a PDF with fractional standard deviation |
|---|
| 123 | ! "fsd" corresponding to the cumulative distribution function value |
|---|
| 124 | ! "cdf", and return it in val. Since this is an elemental |
|---|
| 125 | ! subroutine, fsd, cdf and val may be arrays. |
|---|
| 126 | elemental subroutine sample_from_pdf(this, fsd, cdf, val) |
|---|
| 127 | |
|---|
| 128 | class(pdf_sampler_type), intent(in) :: this |
|---|
| 129 | |
|---|
| 130 | ! Fractional standard deviation (0 to 4) and cumulative |
|---|
| 131 | ! distribution function (0 to 1) |
|---|
| 132 | real(jprb), intent(in) :: fsd, cdf |
|---|
| 133 | |
|---|
| 134 | ! Sample from distribution |
|---|
| 135 | real(jprb), intent(out) :: val |
|---|
| 136 | |
|---|
| 137 | ! Index to look-up table |
|---|
| 138 | integer :: ifsd, icdf |
|---|
| 139 | |
|---|
| 140 | ! Weights in bilinear interpolation |
|---|
| 141 | real(jprb) :: wfsd, wcdf |
|---|
| 142 | |
|---|
| 143 | ! Bilinear interpolation with bounds |
|---|
| 144 | wcdf = cdf * (this%ncdf-1) + 1.0_jprb |
|---|
| 145 | icdf = max(1, min(int(wcdf), this%ncdf-1)) |
|---|
| 146 | wcdf = max(0.0_jprb, min(wcdf - icdf, 1.0_jprb)) |
|---|
| 147 | |
|---|
| 148 | wfsd = (fsd-this%fsd1) * this%inv_fsd_interval + 1.0_jprb |
|---|
| 149 | ifsd = max(1, min(int(wfsd), this%nfsd-1)) |
|---|
| 150 | wfsd = max(0.0_jprb, min(wfsd - ifsd, 1.0_jprb)) |
|---|
| 151 | |
|---|
| 152 | val = (1.0_jprb-wcdf)*(1.0_jprb-wfsd) * this%val(icdf ,ifsd) & |
|---|
| 153 | & + (1.0_jprb-wcdf)* wfsd * this%val(icdf ,ifsd+1) & |
|---|
| 154 | & + wcdf *(1.0_jprb-wfsd) * this%val(icdf+1,ifsd) & |
|---|
| 155 | & + wcdf * wfsd * this%val(icdf+1,ifsd+1) |
|---|
| 156 | |
|---|
| 157 | end subroutine sample_from_pdf |
|---|
| 158 | |
|---|
| 159 | |
|---|
| 160 | !--------------------------------------------------------------------- |
|---|
| 161 | ! For true elements of mask, extract the values of a PDF with |
|---|
| 162 | ! fractional standard deviation "fsd" corresponding to the |
|---|
| 163 | ! cumulative distribution function values "cdf", and return in |
|---|
| 164 | ! val. For false elements of mask, return zero in val. |
|---|
| 165 | subroutine sample_from_pdf_masked(this, nsamp, fsd, cdf, val, mask) |
|---|
| 166 | |
|---|
| 167 | class(pdf_sampler_type), intent(in) :: this |
|---|
| 168 | |
|---|
| 169 | ! Number of samples |
|---|
| 170 | integer, intent(in) :: nsamp |
|---|
| 171 | |
|---|
| 172 | ! Fractional standard deviation (0 to 4) and cumulative |
|---|
| 173 | ! distribution function (0 to 1) |
|---|
| 174 | real(jprb), intent(in) :: fsd(nsamp), cdf(nsamp) |
|---|
| 175 | |
|---|
| 176 | ! Sample from distribution |
|---|
| 177 | real(jprb), intent(out) :: val(:) |
|---|
| 178 | |
|---|
| 179 | ! Mask |
|---|
| 180 | logical, intent(in) :: mask(nsamp) |
|---|
| 181 | |
|---|
| 182 | ! Loop index |
|---|
| 183 | integer :: jsamp |
|---|
| 184 | |
|---|
| 185 | ! Index to look-up table |
|---|
| 186 | integer :: ifsd, icdf |
|---|
| 187 | |
|---|
| 188 | ! Weights in bilinear interpolation |
|---|
| 189 | real(jprb) :: wfsd, wcdf |
|---|
| 190 | |
|---|
| 191 | do jsamp = 1,nsamp |
|---|
| 192 | if (mask(jsamp)) then |
|---|
| 193 | ! Bilinear interpolation with bounds |
|---|
| 194 | wcdf = cdf(jsamp) * (this%ncdf-1) + 1.0_jprb |
|---|
| 195 | icdf = max(1, min(int(wcdf), this%ncdf-1)) |
|---|
| 196 | wcdf = max(0.0_jprb, min(wcdf - icdf, 1.0_jprb)) |
|---|
| 197 | |
|---|
| 198 | wfsd = (fsd(jsamp)-this%fsd1) * this%inv_fsd_interval + 1.0_jprb |
|---|
| 199 | ifsd = max(1, min(int(wfsd), this%nfsd-1)) |
|---|
| 200 | wfsd = max(0.0_jprb, min(wfsd - ifsd, 1.0_jprb)) |
|---|
| 201 | |
|---|
| 202 | val(jsamp)=(1.0_jprb-wcdf)*(1.0_jprb-wfsd) * this%val(icdf ,ifsd) & |
|---|
| 203 | & +(1.0_jprb-wcdf)* wfsd * this%val(icdf ,ifsd+1) & |
|---|
| 204 | & + wcdf *(1.0_jprb-wfsd) * this%val(icdf+1,ifsd) & |
|---|
| 205 | & + wcdf * wfsd * this%val(icdf+1,ifsd+1) |
|---|
| 206 | else |
|---|
| 207 | val(jsamp) = 0.0_jprb |
|---|
| 208 | end if |
|---|
| 209 | end do |
|---|
| 210 | end subroutine sample_from_pdf_masked |
|---|
| 211 | |
|---|
| 212 | !--------------------------------------------------------------------- |
|---|
| 213 | ! Extract the values of a PDF with fractional standard deviation |
|---|
| 214 | ! "fsd" corresponding to the cumulative distribution function values |
|---|
| 215 | ! "cdf", and return in val. This version works on 2D blocks of data. |
|---|
| 216 | subroutine sample_from_pdf_block(this, nz, ng, fsd, cdf, val) |
|---|
| 217 | |
|---|
| 218 | class(pdf_sampler_type), intent(in) :: this |
|---|
| 219 | |
|---|
| 220 | ! Number of samples |
|---|
| 221 | integer, intent(in) :: nz, ng |
|---|
| 222 | |
|---|
| 223 | ! Fractional standard deviation (0 to 4) and cumulative |
|---|
| 224 | ! distribution function (0 to 1) |
|---|
| 225 | real(jprb), intent(in) :: fsd(nz), cdf(ng, nz) |
|---|
| 226 | |
|---|
| 227 | ! Sample from distribution |
|---|
| 228 | real(jprb), intent(out) :: val(:,:) |
|---|
| 229 | |
|---|
| 230 | ! Loop index |
|---|
| 231 | integer :: jz, jg |
|---|
| 232 | |
|---|
| 233 | ! Index to look-up table |
|---|
| 234 | integer :: ifsd, icdf |
|---|
| 235 | |
|---|
| 236 | ! Weights in bilinear interpolation |
|---|
| 237 | real(jprb) :: wfsd, wcdf |
|---|
| 238 | |
|---|
| 239 | do jz = 1,nz |
|---|
| 240 | do jg = 1,ng |
|---|
| 241 | if (cdf(jg, jz) > 0.0_jprb) then |
|---|
| 242 | ! Bilinear interpolation with bounds |
|---|
| 243 | wcdf = cdf(jg,jz) * (this%ncdf-1) + 1.0_jprb |
|---|
| 244 | icdf = max(1, min(int(wcdf), this%ncdf-1)) |
|---|
| 245 | wcdf = max(0.0_jprb, min(wcdf - icdf, 1.0_jprb)) |
|---|
| 246 | |
|---|
| 247 | wfsd = (fsd(jz)-this%fsd1) * this%inv_fsd_interval + 1.0_jprb |
|---|
| 248 | ifsd = max(1, min(int(wfsd), this%nfsd-1)) |
|---|
| 249 | wfsd = max(0.0_jprb, min(wfsd - ifsd, 1.0_jprb)) |
|---|
| 250 | |
|---|
| 251 | val(jg,jz)=(1.0_jprb-wcdf)*(1.0_jprb-wfsd) * this%val(icdf ,ifsd) & |
|---|
| 252 | & +(1.0_jprb-wcdf)* wfsd * this%val(icdf ,ifsd+1) & |
|---|
| 253 | & + wcdf *(1.0_jprb-wfsd) * this%val(icdf+1,ifsd) & |
|---|
| 254 | & + wcdf * wfsd * this%val(icdf+1,ifsd+1) |
|---|
| 255 | else |
|---|
| 256 | val(jg,jz) = 0.0_jprb |
|---|
| 257 | end if |
|---|
| 258 | end do |
|---|
| 259 | end do |
|---|
| 260 | |
|---|
| 261 | end subroutine sample_from_pdf_block |
|---|
| 262 | |
|---|
| 263 | !--------------------------------------------------------------------- |
|---|
| 264 | ! Extract the values of a PDF with fractional standard deviation |
|---|
| 265 | ! "fsd" corresponding to the cumulative distribution function values |
|---|
| 266 | ! "cdf", and return in val. This version works on 2D blocks of data. |
|---|
| 267 | subroutine sample_from_pdf_masked_block(this, nz, ng, fsd, cdf, val, mask) |
|---|
| 268 | |
|---|
| 269 | class(pdf_sampler_type), intent(in) :: this |
|---|
| 270 | |
|---|
| 271 | ! Number of samples |
|---|
| 272 | integer, intent(in) :: nz, ng |
|---|
| 273 | |
|---|
| 274 | ! Fractional standard deviation (0 to 4) and cumulative |
|---|
| 275 | ! distribution function (0 to 1) |
|---|
| 276 | real(jprb), intent(in) :: fsd(nz), cdf(ng, nz) |
|---|
| 277 | |
|---|
| 278 | ! Sample from distribution |
|---|
| 279 | real(jprb), intent(out) :: val(:,:) |
|---|
| 280 | |
|---|
| 281 | ! Mask |
|---|
| 282 | logical, intent(in), optional :: mask(nz) |
|---|
| 283 | |
|---|
| 284 | ! Loop index |
|---|
| 285 | integer :: jz, jg |
|---|
| 286 | |
|---|
| 287 | ! Index to look-up table |
|---|
| 288 | integer :: ifsd, icdf |
|---|
| 289 | |
|---|
| 290 | ! Weights in bilinear interpolation |
|---|
| 291 | real(jprb) :: wfsd, wcdf |
|---|
| 292 | |
|---|
| 293 | do jz = 1,nz |
|---|
| 294 | |
|---|
| 295 | if (mask(jz)) then |
|---|
| 296 | |
|---|
| 297 | do jg = 1,ng |
|---|
| 298 | if (cdf(jg, jz) > 0.0_jprb) then |
|---|
| 299 | ! Bilinear interpolation with bounds |
|---|
| 300 | wcdf = cdf(jg,jz) * (this%ncdf-1) + 1.0_jprb |
|---|
| 301 | icdf = max(1, min(int(wcdf), this%ncdf-1)) |
|---|
| 302 | wcdf = max(0.0_jprb, min(wcdf - icdf, 1.0_jprb)) |
|---|
| 303 | |
|---|
| 304 | wfsd = (fsd(jz)-this%fsd1) * this%inv_fsd_interval + 1.0_jprb |
|---|
| 305 | ifsd = max(1, min(int(wfsd), this%nfsd-1)) |
|---|
| 306 | wfsd = max(0.0_jprb, min(wfsd - ifsd, 1.0_jprb)) |
|---|
| 307 | |
|---|
| 308 | val(jg,jz)=(1.0_jprb-wcdf)*(1.0_jprb-wfsd) * this%val(icdf ,ifsd) & |
|---|
| 309 | & +(1.0_jprb-wcdf)* wfsd * this%val(icdf ,ifsd+1) & |
|---|
| 310 | & + wcdf *(1.0_jprb-wfsd) * this%val(icdf+1,ifsd) & |
|---|
| 311 | & + wcdf * wfsd * this%val(icdf+1,ifsd+1) |
|---|
| 312 | else |
|---|
| 313 | val(jg,jz) = 0.0_jprb |
|---|
| 314 | end if |
|---|
| 315 | end do |
|---|
| 316 | |
|---|
| 317 | end if |
|---|
| 318 | |
|---|
| 319 | end do |
|---|
| 320 | |
|---|
| 321 | end subroutine sample_from_pdf_masked_block |
|---|
| 322 | |
|---|
| 323 | end module radiation_pdf_sampler |
|---|
