[4773] | 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 and scattering model |
---|
| 47 | ! (e.g. "mie_droplet", "fu-muskatel_ice"). These are used to |
---|
| 48 | ! generate the name of the data file from which the coefficients |
---|
| 49 | ! are read. |
---|
| 50 | character(len=511) :: type_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 | procedure :: save => save_general_cloud_optics_data |
---|
| 59 | |
---|
| 60 | end type general_cloud_optics_type |
---|
| 61 | |
---|
| 62 | contains |
---|
| 63 | |
---|
| 64 | ! Provides elemental function "delta_eddington" |
---|
| 65 | #include "radiation_delta_eddington.h" |
---|
| 66 | |
---|
| 67 | !--------------------------------------------------------------------- |
---|
| 68 | ! Setup cloud optics coefficients by reading them from a file |
---|
| 69 | subroutine setup_general_cloud_optics(this, file_name, specdef, & |
---|
| 70 | & use_bands, use_thick_averaging, & |
---|
| 71 | & weighting_temperature, & |
---|
| 72 | & iverbose) |
---|
| 73 | |
---|
| 74 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 75 | use easy_netcdf, only : netcdf_file |
---|
| 76 | use radiation_spectral_definition, only : spectral_definition_type |
---|
| 77 | use radiation_io, only : nulout, nulerr, radiation_abort |
---|
| 78 | |
---|
| 79 | class(general_cloud_optics_type), intent(inout) :: this |
---|
| 80 | character(len=*), intent(in) :: file_name |
---|
| 81 | type(spectral_definition_type), intent(in) :: specdef |
---|
| 82 | logical, intent(in), optional :: use_bands, use_thick_averaging |
---|
| 83 | real(jprb), intent(in), optional :: weighting_temperature ! K |
---|
| 84 | integer, intent(in), optional :: iverbose |
---|
| 85 | |
---|
| 86 | ! Spectral properties read from file, dimensioned (wavenumber, |
---|
| 87 | ! n_effective_radius) |
---|
| 88 | real(jprb), dimension(:,:), allocatable :: mass_ext, & ! m2 kg-1 |
---|
| 89 | & ssa, asymmetry |
---|
| 90 | |
---|
| 91 | ! Reflectance of an infinitely thick cloud, needed for thick |
---|
| 92 | ! averaging |
---|
| 93 | real(jprb), dimension(:,:), allocatable :: ref_inf |
---|
| 94 | |
---|
| 95 | ! Coordinate variables from file |
---|
| 96 | real(jprb), dimension(:), allocatable :: wavenumber ! cm-1 |
---|
| 97 | real(jprb), dimension(:), allocatable :: effective_radius ! m |
---|
| 98 | |
---|
| 99 | ! Matrix mapping optical properties in the file to values per |
---|
| 100 | ! g-point or band, such that in the thin-averaging case, |
---|
| 101 | ! this%mass_ext=matmul(mapping,file%mass_ext), so mapping is |
---|
| 102 | ! dimensioned (ngpoint,nwav) |
---|
| 103 | real(jprb), dimension(:,:), allocatable :: mapping |
---|
| 104 | |
---|
| 105 | ! The NetCDF file containing the coefficients |
---|
| 106 | type(netcdf_file) :: file |
---|
| 107 | |
---|
| 108 | real(jprb) :: diff_spread |
---|
| 109 | integer :: iverb |
---|
| 110 | integer :: nre ! Number of effective radii |
---|
| 111 | integer :: nwav ! Number of wavenumbers describing cloud |
---|
| 112 | |
---|
| 113 | logical :: use_bands_local, use_thick_averaging_local |
---|
| 114 | |
---|
| 115 | real(jphook) :: hook_handle |
---|
| 116 | |
---|
| 117 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:setup',0,hook_handle) |
---|
| 118 | |
---|
| 119 | ! Set local values of optional inputs |
---|
| 120 | if (present(iverbose)) then |
---|
| 121 | iverb = iverbose |
---|
| 122 | else |
---|
| 123 | iverb = 2 |
---|
| 124 | end if |
---|
| 125 | |
---|
| 126 | if (present(use_bands)) then |
---|
| 127 | use_bands_local = use_bands |
---|
| 128 | else |
---|
| 129 | use_bands_local = .false. |
---|
| 130 | end if |
---|
| 131 | |
---|
| 132 | if (present(use_thick_averaging)) then |
---|
| 133 | use_thick_averaging_local = use_thick_averaging |
---|
| 134 | else |
---|
| 135 | use_thick_averaging_local = .false. |
---|
| 136 | end if |
---|
| 137 | |
---|
| 138 | ! Open the scattering file and configure the way it is read |
---|
| 139 | call file%open(trim(file_name), iverbose=iverb) |
---|
| 140 | !call file%transpose_matrices() |
---|
| 141 | |
---|
| 142 | ! Read coordinate variables |
---|
| 143 | call file%get('wavenumber', wavenumber) |
---|
| 144 | call file%get('effective_radius', effective_radius) |
---|
| 145 | |
---|
| 146 | ! Read the band-specific coefficients |
---|
| 147 | call file%get('mass_extinction_coefficient', mass_ext) |
---|
| 148 | call file%get('single_scattering_albedo', ssa) |
---|
| 149 | call file%get('asymmetry_factor', asymmetry) |
---|
| 150 | |
---|
| 151 | ! Close scattering file |
---|
| 152 | call file%close() |
---|
| 153 | |
---|
| 154 | ! Check effective radius is evenly spaced |
---|
| 155 | nre = size(effective_radius) |
---|
| 156 | ! Fractional range of differences, should be near zero for evenly |
---|
| 157 | ! spaced data |
---|
| 158 | diff_spread = (maxval(effective_radius(2:nre)-effective_radius(1:nre-1)) & |
---|
| 159 | & -minval(effective_radius(2:nre)-effective_radius(1:nre-1))) & |
---|
| 160 | & / minval(abs(effective_radius(2:nre)-effective_radius(1:nre-1))) |
---|
| 161 | if (diff_spread > 0.01_jprb) then |
---|
| 162 | write(nulerr, '(a,a,a)') '*** Error: effective_radius in ', & |
---|
| 163 | & trim(file_name), ', is not evenly spaced to 1%' |
---|
| 164 | call radiation_abort('Radiation configuration error') |
---|
| 165 | end if |
---|
| 166 | |
---|
| 167 | ! Set up effective radius coordinate variable |
---|
| 168 | this%n_effective_radius = nre |
---|
| 169 | this%effective_radius_0 = effective_radius(1) |
---|
| 170 | this%d_effective_radius = effective_radius(2) - effective_radius(1) |
---|
| 171 | |
---|
| 172 | nwav = size(wavenumber) |
---|
| 173 | |
---|
| 174 | ! Define the mapping matrix |
---|
| 175 | call specdef%calc_mapping(wavenumber, mapping, & |
---|
| 176 | weighting_temperature=weighting_temperature, use_bands=use_bands) |
---|
| 177 | |
---|
| 178 | ! Thick averaging should be performed on delta-Eddington scaled |
---|
| 179 | ! quantities (it makes no difference to thin averaging) |
---|
| 180 | call delta_eddington(mass_ext, ssa, asymmetry) |
---|
| 181 | |
---|
| 182 | ! Thin averaging |
---|
| 183 | this%mass_ext = matmul(mapping, mass_ext) |
---|
| 184 | this%ssa = matmul(mapping, mass_ext*ssa) / this%mass_ext |
---|
| 185 | this%asymmetry = matmul(mapping, mass_ext*ssa*asymmetry) / (this%mass_ext*this%ssa) |
---|
| 186 | |
---|
| 187 | if (use_thick_averaging_local) then |
---|
| 188 | ! Thick averaging as described by Edwards and Slingo (1996), |
---|
| 189 | ! modifying only the single-scattering albedo |
---|
| 190 | allocate(ref_inf(nwav, nre)) |
---|
| 191 | |
---|
| 192 | ! Eqs. 18 and 17 of Edwards & Slingo (1996) |
---|
| 193 | ref_inf = sqrt((1.0_jprb - ssa) / (1.0_jprb - ssa*asymmetry)) |
---|
| 194 | ref_inf = (1.0_jprb - ref_inf) / (1.0_jprb + ref_inf) |
---|
| 195 | ! Here the left-hand side is actually the averaged ref_inf |
---|
| 196 | this%ssa = matmul(mapping, ref_inf) |
---|
| 197 | ! Eq. 19 of Edwards and Slingo (1996) |
---|
| 198 | this%ssa = 4.0_jprb * this%ssa / ((1.0_jprb + this%ssa)**2 & |
---|
| 199 | & - this%asymmetry * (1.0_jprb - this%ssa)**2) |
---|
| 200 | |
---|
| 201 | deallocate(ref_inf) |
---|
| 202 | end if |
---|
| 203 | |
---|
| 204 | deallocate(mapping) |
---|
| 205 | |
---|
| 206 | ! Revert back to unscaled quantities |
---|
| 207 | call revert_delta_eddington(this%mass_ext, this%ssa, this%asymmetry) |
---|
| 208 | |
---|
| 209 | if (iverb >= 2) then |
---|
| 210 | write(nulout,'(a,a)') ' File: ', trim(file_name) |
---|
| 211 | if (present(weighting_temperature)) then |
---|
| 212 | write(nulout,'(a,f7.1,a)') ' Weighting temperature: ', weighting_temperature, ' K' |
---|
| 213 | else |
---|
| 214 | write(nulout,'(a,f7.1,a)') ' Weighting temperature: ', specdef%reference_temperature, ' K' |
---|
| 215 | end if |
---|
| 216 | if (use_thick_averaging_local) then |
---|
| 217 | write(nulout,'(a)') ' SSA averaging: optically thick limit' |
---|
| 218 | else |
---|
| 219 | write(nulout,'(a)') ' SSA averaging: optically thin limit' |
---|
| 220 | end if |
---|
| 221 | if (use_bands_local) then |
---|
| 222 | write(nulout,'(a,i0,a)') ' Spectral discretization: ', specdef%nband, ' bands' |
---|
| 223 | else |
---|
| 224 | write(nulout,'(a,i0,a)') ' Spectral discretization: ', specdef%ng, ' g-points' |
---|
| 225 | end if |
---|
| 226 | write(nulout,'(a,i0,a,f6.1,a,f6.1,a)') ' Effective radius look-up: ', nre, ' points in range ', & |
---|
| 227 | & effective_radius(1)*1.0e6_jprb, '-', effective_radius(nre)*1.0e6_jprb, ' um' |
---|
| 228 | write(nulout,'(a,i0,a,i0,a)') ' Wavenumber range: ', int(specdef%min_wavenumber()), '-', & |
---|
| 229 | & int(specdef%max_wavenumber()), ' cm-1' |
---|
| 230 | end if |
---|
| 231 | |
---|
| 232 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:setup',1,hook_handle) |
---|
| 233 | |
---|
| 234 | end subroutine setup_general_cloud_optics |
---|
| 235 | |
---|
| 236 | |
---|
| 237 | !--------------------------------------------------------------------- |
---|
| 238 | ! Add the optical properties of a particular cloud type to the |
---|
| 239 | ! accumulated optical properties of all cloud types |
---|
| 240 | subroutine add_optical_properties(this, ng, nlev, ncol, & |
---|
| 241 | & cloud_fraction, & |
---|
| 242 | & water_path, effective_radius, & |
---|
| 243 | & od, scat_od, scat_asymmetry) |
---|
| 244 | |
---|
| 245 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 246 | |
---|
| 247 | class(general_cloud_optics_type), intent(in) :: this |
---|
| 248 | |
---|
| 249 | ! Number of g points, levels and columns |
---|
| 250 | integer, intent(in) :: ng, nlev, ncol |
---|
| 251 | |
---|
| 252 | ! Properties of present cloud type, dimensioned (ncol,nlev) |
---|
| 253 | real(jprb), intent(in) :: cloud_fraction(:,:) |
---|
| 254 | real(jprb), intent(in) :: water_path(:,:) ! kg m-2 |
---|
| 255 | real(jprb), intent(in) :: effective_radius(:,:) ! m |
---|
| 256 | |
---|
| 257 | ! Optical properties which are additive per cloud type, |
---|
| 258 | ! dimensioned (ng,nlev,ncol) |
---|
| 259 | real(jprb), intent(inout), dimension(ng,nlev,ncol) & |
---|
| 260 | & :: od ! Optical depth of layer |
---|
| 261 | real(jprb), intent(inout), dimension(ng,nlev,ncol), optional & |
---|
| 262 | & :: scat_od, & ! Scattering optical depth of layer |
---|
| 263 | & scat_asymmetry ! Scattering optical depth x asymmetry factor |
---|
| 264 | |
---|
| 265 | real(jprb) :: od_local(ng) |
---|
| 266 | |
---|
| 267 | real(jprb) :: re_index, weight1, weight2 |
---|
| 268 | integer :: ire |
---|
| 269 | |
---|
| 270 | integer :: jcol, jlev |
---|
| 271 | |
---|
| 272 | real(jphook) :: hook_handle |
---|
| 273 | |
---|
| 274 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:add_optical_properties',0,hook_handle) |
---|
| 275 | |
---|
| 276 | if (present(scat_od)) then |
---|
| 277 | do jcol = 1,ncol |
---|
| 278 | do jlev = 1,nlev |
---|
| 279 | if (cloud_fraction(jcol, jlev) > 0.0_jprb) then |
---|
| 280 | re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) & |
---|
| 281 | & / this%d_effective_radius, this%n_effective_radius-0.0001_jprb)) |
---|
| 282 | ire = int(re_index) |
---|
| 283 | weight2 = re_index - ire |
---|
| 284 | weight1 = 1.0_jprb - weight2 |
---|
| 285 | od_local = water_path(jcol, jlev) * (weight1*this%mass_ext(:,ire) & |
---|
| 286 | & +weight2*this%mass_ext(:,ire+1)) |
---|
| 287 | od(:,jlev,jcol) = od(:,jlev,jcol) + od_local |
---|
| 288 | od_local = od_local * (weight1*this%ssa(:,ire) & |
---|
| 289 | & +weight2*this%ssa(:,ire+1)) |
---|
| 290 | scat_od(:,jlev,jcol) = scat_od(:,jlev,jcol) + od_local |
---|
| 291 | scat_asymmetry(:,jlev,jcol) = scat_asymmetry(:,jlev,jcol) & |
---|
| 292 | & + od_local * (weight1*this%asymmetry(:,ire) & |
---|
| 293 | & +weight2*this%asymmetry(:,ire+1)) |
---|
| 294 | end if |
---|
| 295 | end do |
---|
| 296 | end do |
---|
| 297 | else |
---|
| 298 | ! No scattering: return the absorption optical depth |
---|
| 299 | do jcol = 1,ncol |
---|
| 300 | do jlev = 1,nlev |
---|
| 301 | if (water_path(jcol, jlev) > 0.0_jprb) then |
---|
| 302 | re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) & |
---|
| 303 | & / this%d_effective_radius, this%n_effective_radius-0.0001_jprb)) |
---|
| 304 | ire = int(re_index) |
---|
| 305 | weight2 = re_index - ire |
---|
| 306 | weight1 = 1.0_jprb - weight2 |
---|
| 307 | od(:,jlev,jcol) = od(:,jlev,jcol) & |
---|
| 308 | & + water_path(jcol, jlev) * (weight1*this%mass_ext(:,ire) & |
---|
| 309 | & +weight2*this%mass_ext(:,ire+1)) & |
---|
| 310 | & * (1.0_jprb - (weight1*this%ssa(:,ire)+weight2*this%ssa(:,ire+1))) |
---|
| 311 | end if |
---|
| 312 | end do |
---|
| 313 | end do |
---|
| 314 | end if |
---|
| 315 | |
---|
| 316 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:add_optical_properties',1,hook_handle) |
---|
| 317 | |
---|
| 318 | end subroutine add_optical_properties |
---|
| 319 | |
---|
| 320 | |
---|
| 321 | !--------------------------------------------------------------------- |
---|
| 322 | ! Return the Planck function (in W m-2 (cm-1)-1) for a given |
---|
| 323 | ! wavenumber (cm-1) and temperature (K), ensuring double precision |
---|
| 324 | ! for internal calculation |
---|
| 325 | elemental function calc_planck_function_wavenumber(wavenumber, temperature) |
---|
| 326 | |
---|
| 327 | use parkind1, only : jprb, jprd |
---|
| 328 | use radiation_constants, only : SpeedOfLight, BoltzmannConstant, PlanckConstant |
---|
| 329 | |
---|
| 330 | real(jprb), intent(in) :: wavenumber ! cm-1 |
---|
| 331 | real(jprb), intent(in) :: temperature ! K |
---|
| 332 | real(jprb) :: calc_planck_function_wavenumber |
---|
| 333 | |
---|
| 334 | real(jprd) :: freq ! Hz |
---|
| 335 | real(jprd) :: planck_fn_freq ! W m-2 Hz-1 |
---|
| 336 | |
---|
| 337 | freq = 100.0_jprd * real(SpeedOfLight,jprd) * real(wavenumber,jprd) |
---|
| 338 | planck_fn_freq = 2.0_jprd * real(PlanckConstant,jprd) * freq**3 & |
---|
| 339 | & / (real(SpeedOfLight,jprd)**2 * (exp(real(PlanckConstant,jprd)*freq & |
---|
| 340 | & /(real(BoltzmannConstant,jprd)*real(temperature,jprd))) - 1.0_jprd)) |
---|
| 341 | calc_planck_function_wavenumber = real(planck_fn_freq * 100.0_jprd * real(SpeedOfLight,jprd), jprb) |
---|
| 342 | |
---|
| 343 | end function calc_planck_function_wavenumber |
---|
| 344 | |
---|
| 345 | !--------------------------------------------------------------------- |
---|
| 346 | ! Save cloud optical properties in the named file |
---|
| 347 | subroutine save_general_cloud_optics_data(this, file_name, iverbose) |
---|
| 348 | |
---|
| 349 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 350 | use easy_netcdf, only : netcdf_file |
---|
| 351 | |
---|
| 352 | class(general_cloud_optics_type), intent(in) :: this |
---|
| 353 | character(len=*), intent(in) :: file_name |
---|
| 354 | integer, optional, intent(in) :: iverbose |
---|
| 355 | |
---|
| 356 | ! Object for output NetCDF file |
---|
| 357 | type(netcdf_file) :: out_file |
---|
| 358 | |
---|
| 359 | real(jprb) :: effective_radius(this%n_effective_radius) |
---|
| 360 | integer :: ire |
---|
| 361 | |
---|
| 362 | real(jphook) :: hook_handle |
---|
| 363 | |
---|
| 364 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:save',0,hook_handle) |
---|
| 365 | |
---|
| 366 | ! Create the file |
---|
| 367 | call out_file%create(trim(file_name), iverbose=iverbose) |
---|
| 368 | |
---|
| 369 | ! Define dimensions |
---|
| 370 | call out_file%define_dimension("band", size(this%mass_ext,1)) |
---|
| 371 | call out_file%define_dimension("effective_radius", this%n_effective_radius) |
---|
| 372 | |
---|
| 373 | ! Put global attributes |
---|
| 374 | call out_file%put_global_attributes( & |
---|
| 375 | & title_str="Optical properties of "//trim(this%type_name) & |
---|
| 376 | & //" hydrometeors using the spectral intervals of ecRad", & |
---|
| 377 | & source_str="ecRad offline radiation model") |
---|
| 378 | |
---|
| 379 | ! Define variables |
---|
| 380 | call out_file%define_variable("effective_radius", units_str="m", & |
---|
| 381 | & long_name="Effective radius", dim1_name="effective_radius") |
---|
| 382 | call out_file%define_variable("mass_extinction_coefficient", units_str="m2 kg-1", & |
---|
| 383 | & long_name="Mass-extinction coefficient", & |
---|
| 384 | & dim2_name="effective_radius", dim1_name="band") |
---|
| 385 | call out_file%define_variable("single_scattering_albedo", units_str="1", & |
---|
| 386 | & long_name="Single scattering albedo", & |
---|
| 387 | & dim2_name="effective_radius", dim1_name="band") |
---|
| 388 | call out_file%define_variable("asymmetry_factor", units_str="1", & |
---|
| 389 | & long_name="Asymmetry factor", & |
---|
| 390 | & dim2_name="effective_radius", dim1_name="band") |
---|
| 391 | |
---|
| 392 | ! Define effective radius |
---|
| 393 | do ire = 1,this%n_effective_radius |
---|
| 394 | effective_radius(ire) = this%effective_radius_0 + this%d_effective_radius*(ire-1) |
---|
| 395 | end do |
---|
| 396 | |
---|
| 397 | ! Write variables |
---|
| 398 | call out_file%put("effective_radius", effective_radius) |
---|
| 399 | call out_file%put("mass_extinction_coefficient", this%mass_ext) |
---|
| 400 | call out_file%put("single_scattering_albedo", this%ssa) |
---|
| 401 | call out_file%put("asymmetry_factor", this%asymmetry) |
---|
| 402 | |
---|
| 403 | call out_file%close() |
---|
| 404 | |
---|
| 405 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:save',1,hook_handle) |
---|
| 406 | |
---|
| 407 | end subroutine save_general_cloud_optics_data |
---|
| 408 | |
---|
| 409 | end module radiation_general_cloud_optics_data |
---|