[4444] | 1 | ! radiation_general_cloud_optics_data.F90 - Type to store generalized cloud optical properties |
---|
| 2 | ! |
---|
| 3 | ! (C) Copyright 2019- 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 | ! License: see the COPYING file for details |
---|
| 15 | ! |
---|
| 16 | |
---|
| 17 | module radiation_general_cloud_optics_data |
---|
| 18 | |
---|
| 19 | use parkind1, only : jprb |
---|
| 20 | |
---|
| 21 | implicit none |
---|
| 22 | |
---|
| 23 | public |
---|
| 24 | |
---|
| 25 | !--------------------------------------------------------------------- |
---|
| 26 | ! This type holds the configuration information to compute optical |
---|
| 27 | ! properties for a particular type of cloud or hydrometeor in one of |
---|
| 28 | ! the shortwave or longwave |
---|
| 29 | type general_cloud_optics_type |
---|
| 30 | ! Band-specific (or g-point-specific) values as a look-up table |
---|
| 31 | ! versus effective radius dimensioned (nband,n_effective_radius) |
---|
| 32 | |
---|
| 33 | ! Extinction coefficient per unit mass (m2 kg-1) |
---|
| 34 | real(jprb), allocatable, dimension(:,:) :: & |
---|
| 35 | & mass_ext |
---|
| 36 | |
---|
| 37 | ! Single-scattering albedo and asymmetry factor (dimensionless) |
---|
| 38 | real(jprb), allocatable, dimension(:,:) :: & |
---|
| 39 | & ssa, asymmetry |
---|
| 40 | |
---|
| 41 | ! Number of effective radius coefficients, start value and |
---|
| 42 | ! interval in look-up table |
---|
| 43 | integer :: n_effective_radius = 0 |
---|
| 44 | real(jprb) :: effective_radius_0, d_effective_radius |
---|
| 45 | |
---|
| 46 | ! Name of cloud/precip type (e.g. "liquid", "ice", "rain", "snow") |
---|
| 47 | ! and the name of the optics scheme. These two are used to |
---|
| 48 | ! generate the name of the data file from which the coefficients |
---|
| 49 | ! are read. |
---|
| 50 | character(len=511) :: type_name, scheme_name |
---|
| 51 | |
---|
| 52 | ! Do we use bands or g-points? |
---|
| 53 | logical :: use_bands = .false. |
---|
| 54 | |
---|
| 55 | contains |
---|
| 56 | procedure :: setup => setup_general_cloud_optics |
---|
| 57 | procedure :: add_optical_properties |
---|
| 58 | |
---|
| 59 | end type general_cloud_optics_type |
---|
| 60 | |
---|
| 61 | contains |
---|
| 62 | |
---|
| 63 | ! Provides elemental function "delta_eddington" |
---|
| 64 | #include "radiation_delta_eddington.h" |
---|
| 65 | |
---|
| 66 | !--------------------------------------------------------------------- |
---|
| 67 | ! Setup cloud optics coefficients by reading them from a file |
---|
| 68 | subroutine setup_general_cloud_optics(this, file_name, specdef, & |
---|
| 69 | & use_bands, use_thick_averaging, & |
---|
| 70 | & weighting_temperature, & |
---|
| 71 | & iverbose) |
---|
| 72 | |
---|
| 73 | use yomhook, only : lhook, dr_hook |
---|
| 74 | use easy_netcdf, only : netcdf_file |
---|
| 75 | use radiation_spectral_definition, only : spectral_definition_type |
---|
| 76 | use radiation_io, only : nulout, nulerr, radiation_abort |
---|
| 77 | |
---|
| 78 | class(general_cloud_optics_type), intent(inout) :: this |
---|
| 79 | character(len=*), intent(in) :: file_name |
---|
| 80 | type(spectral_definition_type), intent(in) :: specdef |
---|
| 81 | logical, intent(in), optional :: use_bands, use_thick_averaging |
---|
| 82 | real(jprb), intent(in), optional :: weighting_temperature ! K |
---|
| 83 | integer, intent(in), optional :: iverbose |
---|
| 84 | |
---|
| 85 | ! Spectral properties read from file, dimensioned (wavenumber, |
---|
| 86 | ! n_effective_radius) |
---|
| 87 | real(jprb), dimension(:,:), allocatable :: mass_ext, & ! m2 kg-1 |
---|
| 88 | & ssa, asymmetry |
---|
| 89 | |
---|
| 90 | ! Reflectance of an infinitely thick cloud, needed for thick |
---|
| 91 | ! averaging |
---|
| 92 | real(jprb), dimension(:,:), allocatable :: ref_inf |
---|
| 93 | |
---|
| 94 | ! Coordinate variables from file |
---|
| 95 | real(jprb), dimension(:), allocatable :: wavenumber ! cm-1 |
---|
| 96 | real(jprb), dimension(:), allocatable :: effective_radius ! m |
---|
| 97 | |
---|
| 98 | ! Matrix mapping optical properties in the file to values per |
---|
| 99 | ! g-point or band, such that in the thin-averaging case, |
---|
| 100 | ! this%mass_ext=matmul(mapping,file%mass_ext), so mapping is |
---|
| 101 | ! dimensioned (ngpoint,nwav) |
---|
| 102 | real(jprb), dimension(:,:), allocatable :: mapping |
---|
| 103 | |
---|
| 104 | ! The NetCDF file containing the coefficients |
---|
| 105 | type(netcdf_file) :: file |
---|
| 106 | |
---|
| 107 | real(jprb) :: diff_spread |
---|
| 108 | integer :: iverb |
---|
| 109 | integer :: nre ! Number of effective radii |
---|
| 110 | integer :: nwav ! Number of wavenumbers describing cloud |
---|
| 111 | |
---|
| 112 | logical :: use_bands_local, use_thick_averaging_local |
---|
| 113 | |
---|
| 114 | real(jprb) :: hook_handle |
---|
| 115 | |
---|
| 116 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:setup',0,hook_handle) |
---|
| 117 | |
---|
| 118 | ! Set local values of optional inputs |
---|
| 119 | if (present(iverbose)) then |
---|
| 120 | iverb = iverbose |
---|
| 121 | else |
---|
| 122 | iverb = 2 |
---|
| 123 | end if |
---|
| 124 | |
---|
| 125 | if (present(use_bands)) then |
---|
| 126 | use_bands_local = use_bands |
---|
| 127 | else |
---|
| 128 | use_bands_local = .false. |
---|
| 129 | end if |
---|
| 130 | |
---|
| 131 | if (present(use_thick_averaging)) then |
---|
| 132 | use_thick_averaging_local = use_thick_averaging |
---|
| 133 | else |
---|
| 134 | use_thick_averaging_local = .false. |
---|
| 135 | end if |
---|
| 136 | |
---|
| 137 | ! Open the scattering file and configure the way it is read |
---|
| 138 | call file%open(trim(file_name), iverbose=iverb) |
---|
| 139 | !call file%transpose_matrices() |
---|
| 140 | |
---|
| 141 | ! Read coordinate variables |
---|
| 142 | call file%get('wavenumber', wavenumber) |
---|
| 143 | call file%get('effective_radius', effective_radius) |
---|
| 144 | |
---|
| 145 | ! Read the band-specific coefficients |
---|
| 146 | call file%get('mass_extinction_coefficient', mass_ext) |
---|
| 147 | call file%get('single_scattering_albedo', ssa) |
---|
| 148 | call file%get('asymmetry_factor', asymmetry) |
---|
| 149 | |
---|
| 150 | ! Close scattering file |
---|
| 151 | call file%close() |
---|
| 152 | |
---|
| 153 | ! Check effective radius is evenly spaced |
---|
| 154 | nre = size(effective_radius) |
---|
| 155 | ! Fractional range of differences, should be near zero for evenly |
---|
| 156 | ! spaced data |
---|
| 157 | diff_spread = (maxval(effective_radius(2:nre)-effective_radius(1:nre-1)) & |
---|
| 158 | & -minval(effective_radius(2:nre)-effective_radius(1:nre-1))) & |
---|
| 159 | & / minval(abs(effective_radius(2:nre)-effective_radius(1:nre-1))) |
---|
| 160 | if (diff_spread > 0.01_jprb) then |
---|
| 161 | write(nulerr, '(a,a,a)') '*** Error: effective_radius in ', & |
---|
| 162 | & trim(file_name), ', is not evenly spaced to 1%' |
---|
| 163 | call radiation_abort('Radiation configuration error') |
---|
| 164 | end if |
---|
| 165 | |
---|
| 166 | ! Set up effective radius coordinate variable |
---|
| 167 | this%n_effective_radius = nre |
---|
| 168 | this%effective_radius_0 = effective_radius(1) |
---|
| 169 | this%d_effective_radius = effective_radius(2) - effective_radius(1) |
---|
| 170 | |
---|
| 171 | ! Set up weighting |
---|
| 172 | if (.not. present(weighting_temperature)) then |
---|
| 173 | write(nulerr, '(a)') '*** Error: weighting_temperature not provided' |
---|
| 174 | call radiation_abort('Radiation configuration error') |
---|
| 175 | end if |
---|
| 176 | |
---|
| 177 | nwav = size(wavenumber) |
---|
| 178 | |
---|
| 179 | ! Define the mapping matrix |
---|
| 180 | call specdef%calc_mapping(weighting_temperature, & |
---|
| 181 | & wavenumber, mapping, use_bands=use_bands) |
---|
| 182 | |
---|
| 183 | ! Thick averaging should be performed on delta-Eddington scaled |
---|
| 184 | ! quantities (it makes no difference to thin averaging) |
---|
| 185 | call delta_eddington(mass_ext, ssa, asymmetry) |
---|
| 186 | |
---|
| 187 | ! Thin averaging |
---|
| 188 | this%mass_ext = matmul(mapping, mass_ext) |
---|
| 189 | this%ssa = matmul(mapping, mass_ext*ssa) / this%mass_ext |
---|
| 190 | this%asymmetry = matmul(mapping, mass_ext*ssa*asymmetry) / (this%mass_ext*this%ssa) |
---|
| 191 | |
---|
| 192 | if (use_thick_averaging_local) then |
---|
| 193 | ! Thick averaging as described by Edwards and Slingo (1996), |
---|
| 194 | ! modifying only the single-scattering albedo |
---|
| 195 | allocate(ref_inf(nwav, nre)) |
---|
| 196 | |
---|
| 197 | ! Eqs. 18 and 17 of Edwards & Slingo (1996) |
---|
| 198 | ref_inf = sqrt((1.0_jprb - ssa) / (1.0_jprb - ssa*asymmetry)) |
---|
| 199 | ref_inf = (1.0_jprb - ref_inf) / (1.0_jprb + ref_inf) |
---|
| 200 | ! Here the left-hand side is actually the averaged ref_inf |
---|
| 201 | this%ssa = matmul(mapping, ref_inf) |
---|
| 202 | ! Eq. 19 of Edwards and Slingo (1996) |
---|
| 203 | this%ssa = 4.0_jprb * this%ssa / ((1.0_jprb + this%ssa)**2 & |
---|
| 204 | & - this%asymmetry * (1.0_jprb - this%ssa)**2) |
---|
| 205 | |
---|
| 206 | deallocate(ref_inf) |
---|
| 207 | end if |
---|
| 208 | |
---|
| 209 | deallocate(mapping) |
---|
| 210 | |
---|
| 211 | ! Revert back to unscaled quantities |
---|
| 212 | call revert_delta_eddington(this%mass_ext, this%ssa, this%asymmetry) |
---|
| 213 | |
---|
| 214 | if (iverb >= 2) then |
---|
| 215 | write(nulout,'(a,a)') ' File: ', trim(file_name) |
---|
| 216 | write(nulout,'(a,f7.1,a)') ' Weighting temperature: ', weighting_temperature, ' K' |
---|
| 217 | if (use_thick_averaging_local) then |
---|
| 218 | write(nulout,'(a)') ' SSA averaging: optically thick limit' |
---|
| 219 | else |
---|
| 220 | write(nulout,'(a)') ' SSA averaging: optically thin limit' |
---|
| 221 | end if |
---|
| 222 | if (use_bands_local) then |
---|
| 223 | write(nulout,'(a,i0,a)') ' Spectral discretization: ', specdef%nband, ' bands' |
---|
| 224 | else |
---|
| 225 | write(nulout,'(a,i0,a)') ' Spectral discretization: ', specdef%ng, ' g-points' |
---|
| 226 | end if |
---|
| 227 | write(nulout,'(a,i0,a,f6.1,a,f6.1,a)') ' Effective radius look-up: ', nre, ' points in range ', & |
---|
| 228 | & effective_radius(1)*1.0e6_jprb, '-', effective_radius(nre)*1.0e6_jprb, ' um' |
---|
| 229 | write(nulout,'(a,i0,a,i0,a)') ' Wavenumber range: ', int(specdef%min_wavenumber()), '-', & |
---|
| 230 | & int(specdef%max_wavenumber()), ' cm-1' |
---|
| 231 | end if |
---|
| 232 | |
---|
| 233 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:setup',1,hook_handle) |
---|
| 234 | |
---|
| 235 | end subroutine setup_general_cloud_optics |
---|
| 236 | |
---|
| 237 | |
---|
| 238 | !--------------------------------------------------------------------- |
---|
| 239 | ! Add the optical properties of a particular cloud type to the |
---|
| 240 | ! accumulated optical properties of all cloud types |
---|
| 241 | subroutine add_optical_properties(this, ng, nlev, ncol, & |
---|
| 242 | & cloud_fraction, & |
---|
| 243 | & water_path, effective_radius, & |
---|
| 244 | & od, scat_od, scat_asymmetry) |
---|
| 245 | |
---|
| 246 | use yomhook, only : lhook, dr_hook |
---|
| 247 | |
---|
| 248 | class(general_cloud_optics_type), intent(in) :: this |
---|
| 249 | |
---|
| 250 | ! Number of g points, levels and columns |
---|
| 251 | integer, intent(in) :: ng, nlev, ncol |
---|
| 252 | |
---|
| 253 | ! Properties of present cloud type, dimensioned (ncol,nlev) |
---|
| 254 | real(jprb), intent(in) :: cloud_fraction(:,:) |
---|
| 255 | real(jprb), intent(in) :: water_path(:,:) ! kg m-2 |
---|
| 256 | real(jprb), intent(in) :: effective_radius(:,:) ! m |
---|
| 257 | |
---|
| 258 | ! Optical properties which are additive per cloud type, |
---|
| 259 | ! dimensioned (ng,nlev,ncol) |
---|
| 260 | real(jprb), intent(inout), dimension(ng,nlev,ncol) & |
---|
| 261 | & :: od ! Optical depth of layer |
---|
| 262 | real(jprb), intent(inout), dimension(ng,nlev,ncol), optional & |
---|
| 263 | & :: scat_od, & ! Scattering optical depth of layer |
---|
| 264 | & scat_asymmetry ! Scattering optical depth x asymmetry factor |
---|
| 265 | |
---|
| 266 | real(jprb) :: od_local(ng) |
---|
| 267 | |
---|
| 268 | real(jprb) :: re_index, weight1, weight2 |
---|
| 269 | integer :: ire |
---|
| 270 | |
---|
| 271 | integer :: jcol, jlev |
---|
| 272 | |
---|
| 273 | real(jprb) :: hook_handle |
---|
| 274 | |
---|
| 275 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:add_optical_properties',0,hook_handle) |
---|
| 276 | |
---|
| 277 | if (present(scat_od)) then |
---|
| 278 | do jcol = 1,ncol |
---|
| 279 | do jlev = 1,nlev |
---|
| 280 | if (cloud_fraction(jcol, jlev) > 0.0_jprb) then |
---|
| 281 | re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) & |
---|
| 282 | & / this%d_effective_radius, this%n_effective_radius-0.0001_jprb)) |
---|
| 283 | ire = int(re_index) |
---|
| 284 | weight2 = re_index - ire |
---|
| 285 | weight1 = 1.0_jprb - weight2 |
---|
| 286 | od_local = water_path(jcol, jlev) * (weight1*this%mass_ext(:,ire) & |
---|
| 287 | & +weight2*this%mass_ext(:,ire+1)) |
---|
| 288 | od(:,jlev,jcol) = od(:,jlev,jcol) + od_local |
---|
| 289 | od_local = od_local * (weight1*this%ssa(:,ire) & |
---|
| 290 | & +weight2*this%ssa(:,ire+1)) |
---|
| 291 | scat_od(:,jlev,jcol) = scat_od(:,jlev,jcol) + od_local |
---|
| 292 | scat_asymmetry(:,jlev,jcol) = scat_asymmetry(:,jlev,jcol) & |
---|
| 293 | & + od_local * (weight1*this%asymmetry(:,ire) & |
---|
| 294 | & +weight2*this%asymmetry(:,ire+1)) |
---|
| 295 | end if |
---|
| 296 | end do |
---|
| 297 | end do |
---|
| 298 | else |
---|
| 299 | ! No scattering: return the absorption optical depth |
---|
| 300 | do jcol = 1,ncol |
---|
| 301 | do jlev = 1,nlev |
---|
| 302 | if (water_path(jcol, jlev) > 0.0_jprb) then |
---|
| 303 | re_index = max(1.0, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) & |
---|
| 304 | & / this%d_effective_radius, this%n_effective_radius-0.0001_jprb)) |
---|
| 305 | ire = int(re_index) |
---|
| 306 | weight2 = re_index - ire |
---|
| 307 | weight1 = 1.0_jprb - weight2 |
---|
| 308 | od(:,jlev,jcol) = od(:,jlev,jcol) & |
---|
| 309 | & + water_path(jcol, jlev) * (weight1*this%mass_ext(:,ire) & |
---|
| 310 | & +weight2*this%mass_ext(:,ire+1)) & |
---|
| 311 | & * (1.0_jprb - (weight1*this%ssa(:,ire)+weight2*this%ssa(:,ire+1))) |
---|
| 312 | end if |
---|
| 313 | end do |
---|
| 314 | end do |
---|
| 315 | end if |
---|
| 316 | |
---|
| 317 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:add_optical_properties',1,hook_handle) |
---|
| 318 | |
---|
| 319 | end subroutine add_optical_properties |
---|
| 320 | |
---|
| 321 | |
---|
| 322 | !--------------------------------------------------------------------- |
---|
| 323 | ! Return the Planck function (in W m-2 (cm-1)-1) for a given |
---|
| 324 | ! wavenumber (cm-1) and temperature (K), ensuring double precision |
---|
| 325 | ! for internal calculation |
---|
| 326 | elemental function calc_planck_function_wavenumber(wavenumber, temperature) |
---|
| 327 | |
---|
| 328 | use parkind1, only : jprb, jprd |
---|
| 329 | use radiation_constants, only : SpeedOfLight, BoltzmannConstant, PlanckConstant |
---|
| 330 | |
---|
| 331 | real(jprb), intent(in) :: wavenumber ! cm-1 |
---|
| 332 | real(jprb), intent(in) :: temperature ! K |
---|
| 333 | real(jprb) :: calc_planck_function_wavenumber |
---|
| 334 | |
---|
| 335 | real(jprd) :: freq ! Hz |
---|
| 336 | real(jprd) :: planck_fn_freq ! W m-2 Hz-1 |
---|
| 337 | |
---|
| 338 | freq = 100.0_jprd * real(SpeedOfLight,jprd) * real(wavenumber,jprd) |
---|
| 339 | planck_fn_freq = 2.0_jprd * real(PlanckConstant,jprd) * freq**3 & |
---|
| 340 | & / (real(SpeedOfLight,jprd)**2 * (exp(real(PlanckConstant,jprd)*freq & |
---|
| 341 | & /(real(BoltzmannConstant,jprd)*real(temperature,jprd))) - 1.0_jprd)) |
---|
| 342 | calc_planck_function_wavenumber = real(planck_fn_freq * 100.0_jprd * real(SpeedOfLight,jprd), jprb) |
---|
| 343 | |
---|
| 344 | end function calc_planck_function_wavenumber |
---|
| 345 | |
---|
| 346 | end module radiation_general_cloud_optics_data |
---|