[3908] | 1 | ! radiation_aerosol_optics.F90 - Computing aerosol optical properties |
---|
| 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 | ! Modifications |
---|
| 16 | ! 2018-04-15 R Hogan Add "direct" option |
---|
| 17 | |
---|
| 18 | module radiation_aerosol_optics |
---|
| 19 | |
---|
| 20 | implicit none |
---|
| 21 | public |
---|
| 22 | |
---|
| 23 | contains |
---|
| 24 | |
---|
| 25 | ! Provides the elemental function "delta_eddington_extensive" |
---|
| 26 | #include "radiation_delta_eddington.h" |
---|
| 27 | |
---|
| 28 | !--------------------------------------------------------------------- |
---|
| 29 | ! Load aerosol scattering data; this subroutine delegates to one |
---|
| 30 | ! in radiation_aerosol_optics_data.F90 |
---|
| 31 | subroutine setup_aerosol_optics(config) |
---|
| 32 | |
---|
| 33 | use parkind1, only : jprb |
---|
| 34 | use yomhook, only : lhook, dr_hook |
---|
| 35 | use radiation_config, only : config_type |
---|
| 36 | use radiation_aerosol_optics_data, only : aerosol_optics_type |
---|
| 37 | use radiation_io, only : nulerr, radiation_abort |
---|
[4182] | 38 | use setup_aerosol_optics_lmdz_m, only: setup_aerosol_optics_lmdz |
---|
[3908] | 39 | |
---|
| 40 | type(config_type), intent(inout) :: config |
---|
| 41 | |
---|
| 42 | real(jprb) :: hook_handle |
---|
| 43 | |
---|
| 44 | if (lhook) call dr_hook('radiation_aerosol_optics:setup_aerosol_optics',0,hook_handle) |
---|
| 45 | |
---|
| 46 | if (config%n_aerosol_types > 0) then |
---|
| 47 | ! Load data from file and prepare to map config%n_aerosol_types |
---|
| 48 | ! aerosol types |
---|
[4182] | 49 | call setup_aerosol_optics_lmdz(config%aerosol_optics, & |
---|
| 50 | trim(config%aerosol_optics_file_name)) |
---|
[3908] | 51 | |
---|
| 52 | ! Check agreement in number of bands |
---|
| 53 | if (config%n_bands_lw /= config%aerosol_optics%n_bands_lw) then |
---|
| 54 | write(nulerr,'(a)') '*** Error: number of longwave bands does not match aerosol optics look-up table' |
---|
| 55 | call radiation_abort() |
---|
| 56 | end if |
---|
| 57 | if (config%n_bands_sw /= config%aerosol_optics%n_bands_sw) then |
---|
| 58 | write(nulerr,'(a)') '*** Error: number of shortwave bands does not match aerosol optics look-up table' |
---|
| 59 | call radiation_abort() |
---|
| 60 | end if |
---|
| 61 | |
---|
| 62 | ! Map aerosol types to those loaded from the data file |
---|
| 63 | call config%aerosol_optics%set_types(config%i_aerosol_type_map(1:config%n_aerosol_types)) |
---|
| 64 | end if |
---|
| 65 | |
---|
| 66 | call config%aerosol_optics%print_description(config%i_aerosol_type_map(1:config%n_aerosol_types)) |
---|
| 67 | |
---|
| 68 | if (lhook) call dr_hook('radiation_aerosol_optics:setup_aerosol_optics',1,hook_handle) |
---|
| 69 | |
---|
| 70 | end subroutine setup_aerosol_optics |
---|
| 71 | |
---|
| 72 | |
---|
| 73 | !--------------------------------------------------------------------- |
---|
| 74 | ! Compute aerosol optical properties and add to existing gas optical |
---|
| 75 | ! depth and scattering properties |
---|
| 76 | subroutine add_aerosol_optics(nlev,istartcol,iendcol, & |
---|
| 77 | & config, thermodynamics, gas, aerosol, & |
---|
| 78 | & od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw) |
---|
| 79 | |
---|
| 80 | use parkind1, only : jprb |
---|
| 81 | use radiation_io, only : nulout, nulerr, radiation_abort |
---|
| 82 | use yomhook, only : lhook, dr_hook |
---|
| 83 | use radiation_config, only : config_type |
---|
| 84 | use radiation_thermodynamics, only : thermodynamics_type |
---|
| 85 | use radiation_gas, only : gas_type, IH2O, IMassMixingRatio |
---|
| 86 | use radiation_aerosol, only : aerosol_type |
---|
| 87 | use radiation_constants, only : AccelDueToGravity |
---|
| 88 | use radiation_aerosol_optics_data, only : aerosol_optics_type, & |
---|
| 89 | & IAerosolClassUndefined, IAerosolClassIgnored, & |
---|
| 90 | & IAerosolClassHydrophobic, IAerosolClassHydrophilic |
---|
[4188] | 91 | USE phys_local_var_mod, ONLY: rhcl |
---|
[3908] | 92 | |
---|
| 93 | integer, intent(in) :: nlev ! number of model levels |
---|
| 94 | integer, intent(in) :: istartcol, iendcol ! range of columns to process |
---|
| 95 | type(config_type), intent(in), target :: config |
---|
| 96 | type(thermodynamics_type),intent(in) :: thermodynamics |
---|
| 97 | type(gas_type), intent(in) :: gas |
---|
| 98 | type(aerosol_type), intent(in) :: aerosol |
---|
| 99 | ! Optical depth, single scattering albedo and asymmetry factor of |
---|
| 100 | ! the atmosphere (gases on input, gases and aerosols on output) |
---|
| 101 | ! for each g point. Note that longwave ssa and asymmetry and |
---|
| 102 | ! shortwave asymmetry are all zero for gases, so are not yet |
---|
| 103 | ! defined on input and are therefore intent(out). |
---|
| 104 | real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), & |
---|
| 105 | & intent(inout) :: od_lw |
---|
| 106 | real(jprb), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol), & |
---|
| 107 | & intent(out) :: ssa_lw, g_lw |
---|
| 108 | real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), & |
---|
| 109 | & intent(inout) :: od_sw, ssa_sw |
---|
| 110 | real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), & |
---|
| 111 | & intent(out) :: g_sw |
---|
| 112 | |
---|
| 113 | ! Extinction optical depth, scattering optical depth and |
---|
| 114 | ! asymmetry-times-scattering-optical-depth for all the aerosols at |
---|
| 115 | ! a point in space for each spectral band of the shortwave and |
---|
| 116 | ! longwave spectrum |
---|
| 117 | real(jprb), dimension(config%n_bands_sw) & |
---|
| 118 | & :: od_sw_aerosol, scat_sw_aerosol, scat_g_sw_aerosol, local_od_sw |
---|
| 119 | real(jprb), dimension(config%n_bands_lw) :: od_lw_aerosol, local_od_lw |
---|
| 120 | real(jprb), dimension(config%n_bands_lw_if_scattering) & |
---|
| 121 | & :: scat_lw_aerosol, scat_g_lw_aerosol |
---|
| 122 | |
---|
| 123 | real(jprb) :: h2o_mmr(istartcol:iendcol,nlev) |
---|
| 124 | |
---|
| 125 | real(jprb) :: rh ! Relative humidity with respect to liquid water |
---|
| 126 | |
---|
| 127 | ! Factor (kg m-2) to convert mixing ratio (kg kg-1) to mass in |
---|
| 128 | ! path (kg m-2) |
---|
| 129 | real(jprb) :: factor |
---|
| 130 | |
---|
| 131 | ! Temporary extinction and scattering optical depths of aerosol |
---|
| 132 | ! plus gas |
---|
| 133 | real(jprb) :: local_od, local_scat |
---|
| 134 | |
---|
| 135 | ! Loop indices for column, level, g point, band and aerosol type |
---|
| 136 | integer :: jcol, jlev, jg, jtype |
---|
| 137 | |
---|
| 138 | ! Range of levels over which aerosols are present |
---|
| 139 | integer :: istartlev, iendlev |
---|
| 140 | |
---|
| 141 | ! Indices to spectral band and relative humidity look-up table |
---|
| 142 | integer :: iband, irh |
---|
| 143 | |
---|
| 144 | ! Pointer to the aerosol optics coefficients for brevity of access |
---|
| 145 | type(aerosol_optics_type), pointer :: ao |
---|
| 146 | |
---|
| 147 | real(jprb) :: hook_handle |
---|
| 148 | |
---|
| 149 | if (lhook) call dr_hook('radiation_aerosol_optics:add_aerosol_optics',0,hook_handle) |
---|
| 150 | |
---|
| 151 | if (aerosol%is_direct) then |
---|
| 152 | ! Aerosol optical properties have been provided in each band |
---|
| 153 | ! directly by the user |
---|
| 154 | call add_aerosol_optics_direct(nlev,istartcol,iendcol, & |
---|
| 155 | & config, aerosol, & |
---|
| 156 | & od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw) |
---|
| 157 | else |
---|
| 158 | ! Aerosol mixing ratios have been provided |
---|
| 159 | |
---|
| 160 | do jtype = 1,config%n_aerosol_types |
---|
| 161 | if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then |
---|
| 162 | write(nulerr,'(a)') '*** Error: not all aerosol types are defined' |
---|
| 163 | call radiation_abort() |
---|
| 164 | end if |
---|
| 165 | end do |
---|
| 166 | |
---|
| 167 | if (config%iverbose >= 2) then |
---|
| 168 | write(nulout,'(a)') 'Computing aerosol absorption/scattering properties' |
---|
| 169 | end if |
---|
| 170 | |
---|
| 171 | ao => config%aerosol_optics |
---|
| 172 | |
---|
| 173 | istartlev = lbound(aerosol%mixing_ratio,2) |
---|
| 174 | iendlev = ubound(aerosol%mixing_ratio,2) |
---|
| 175 | |
---|
| 176 | if (ubound(aerosol%mixing_ratio,3) /= config%n_aerosol_types) then |
---|
| 177 | write(nulerr,'(a,i0,a,i0)') '*** Error: aerosol%mixing_ratio contains ', & |
---|
| 178 | & ubound(aerosol%mixing_ratio,3), ' aerosol types, expected ', & |
---|
| 179 | & config%n_aerosol_types |
---|
| 180 | call radiation_abort() |
---|
| 181 | end if |
---|
| 182 | |
---|
| 183 | ! Set variables to zero that may not have been previously |
---|
| 184 | g_sw = 0.0_jprb |
---|
| 185 | if (config%do_lw_aerosol_scattering) then |
---|
| 186 | ssa_lw = 0.0_jprb |
---|
| 187 | g_lw = 0.0_jprb |
---|
| 188 | end if |
---|
| 189 | |
---|
[4188] | 190 | !AI juin 2022 |
---|
| 191 | !call gas%get(IH2O, IMassMixingRatio, h2o_mmr, istartcol=istartcol) |
---|
[3908] | 192 | |
---|
| 193 | ! Loop over position |
---|
| 194 | do jlev = istartlev,iendlev |
---|
| 195 | do jcol = istartcol,iendcol |
---|
| 196 | ! Compute relative humidity with respect to liquid |
---|
| 197 | ! saturation and the index to the relative-humidity index of |
---|
| 198 | ! hydrophilic-aerosol data |
---|
[4188] | 199 | ! AI juin 2022 |
---|
| 200 | ! rh = h2o_mmr(jcol,jlev) / thermodynamics%h2o_sat_liq(jcol,jlev) |
---|
| 201 | ! irh = ao%calc_rh_index(rh) |
---|
| 202 | irh = ao%calc_rh_index(rhcl(jcol,jlev)) |
---|
| 203 | ! print*,'irh=',irh |
---|
[3908] | 204 | |
---|
| 205 | factor = ( thermodynamics%pressure_hl(jcol,jlev+1) & |
---|
| 206 | & -thermodynamics%pressure_hl(jcol,jlev ) ) & |
---|
| 207 | & / AccelDueToGravity |
---|
| 208 | |
---|
| 209 | ! Reset temporary arrays |
---|
| 210 | od_sw_aerosol = 0.0_jprb |
---|
| 211 | scat_sw_aerosol = 0.0_jprb |
---|
| 212 | scat_g_sw_aerosol = 0.0_jprb |
---|
| 213 | od_lw_aerosol = 0.0_jprb |
---|
| 214 | scat_lw_aerosol = 0.0_jprb |
---|
| 215 | scat_g_lw_aerosol = 0.0_jprb |
---|
| 216 | |
---|
| 217 | do jtype = 1,config%n_aerosol_types |
---|
| 218 | ! Add the optical depth, scattering optical depth and |
---|
| 219 | ! scattering optical depth-weighted asymmetry factor for |
---|
| 220 | ! this aerosol type to the total for all aerosols. Note |
---|
| 221 | ! that the following expressions are array-wise, the |
---|
| 222 | ! dimension being spectral band. |
---|
| 223 | if (ao%iclass(jtype) == IAerosolClassHydrophobic) then |
---|
| 224 | local_od_sw = factor * aerosol%mixing_ratio(jcol,jlev,jtype) & |
---|
| 225 | & * ao%mass_ext_sw_phobic(:,ao%itype(jtype)) |
---|
| 226 | od_sw_aerosol = od_sw_aerosol + local_od_sw |
---|
| 227 | scat_sw_aerosol = scat_sw_aerosol & |
---|
| 228 | & + local_od_sw * ao%ssa_sw_phobic(:,ao%itype(jtype)) |
---|
| 229 | scat_g_sw_aerosol = scat_g_sw_aerosol & |
---|
| 230 | & + local_od_sw * ao%ssa_sw_phobic(:,ao%itype(jtype)) & |
---|
| 231 | & * ao%g_sw_phobic(:,ao%itype(jtype)) |
---|
| 232 | if (config%do_lw_aerosol_scattering) then |
---|
| 233 | local_od_lw = factor * aerosol%mixing_ratio(jcol,jlev,jtype) & |
---|
| 234 | & * ao%mass_ext_lw_phobic(:,ao%itype(jtype)) |
---|
| 235 | od_lw_aerosol = od_lw_aerosol + local_od_lw |
---|
| 236 | scat_lw_aerosol = scat_lw_aerosol & |
---|
| 237 | & + local_od_lw * ao%ssa_lw_phobic(:,ao%itype(jtype)) |
---|
| 238 | scat_g_lw_aerosol = scat_g_lw_aerosol & |
---|
| 239 | & + local_od_lw * ao%ssa_lw_phobic(:,ao%itype(jtype)) & |
---|
| 240 | & * ao%g_lw_phobic(:,ao%itype(jtype)) |
---|
| 241 | else |
---|
| 242 | ! If aerosol longwave scattering is not included then we |
---|
| 243 | ! weight the optical depth by the single scattering |
---|
| 244 | ! co-albedo |
---|
| 245 | od_lw_aerosol = od_lw_aerosol & |
---|
| 246 | & + factor * aerosol%mixing_ratio(jcol,jlev,jtype) & |
---|
| 247 | & * ao%mass_ext_lw_phobic(:,ao%itype(jtype)) & |
---|
| 248 | & * (1.0_jprb - ao%ssa_lw_phobic(:,ao%itype(jtype))) |
---|
| 249 | end if |
---|
| 250 | else if (ao%iclass(jtype) == IAerosolClassHydrophilic) then |
---|
| 251 | ! Hydrophilic aerosols require the look-up tables to |
---|
| 252 | ! be indexed with irh |
---|
| 253 | local_od_sw = factor * aerosol%mixing_ratio(jcol,jlev,jtype) & |
---|
| 254 | & * ao%mass_ext_sw_philic(:,irh,ao%itype(jtype)) |
---|
| 255 | od_sw_aerosol = od_sw_aerosol + local_od_sw |
---|
| 256 | scat_sw_aerosol = scat_sw_aerosol & |
---|
| 257 | & + local_od_sw * ao%ssa_sw_philic(:,irh,ao%itype(jtype)) |
---|
| 258 | scat_g_sw_aerosol = scat_g_sw_aerosol & |
---|
| 259 | & + local_od_sw * ao%ssa_sw_philic(:,irh,ao%itype(jtype)) & |
---|
| 260 | & * ao%g_sw_philic(:,irh,ao%itype(jtype)) |
---|
| 261 | if (config%do_lw_aerosol_scattering) then |
---|
| 262 | local_od_lw = factor * aerosol%mixing_ratio(jcol,jlev,jtype) & |
---|
| 263 | & * ao%mass_ext_lw_philic(:,irh,ao%itype(jtype)) |
---|
| 264 | od_lw_aerosol = od_lw_aerosol + local_od_lw |
---|
| 265 | scat_lw_aerosol = scat_lw_aerosol & |
---|
| 266 | & + local_od_lw * ao%ssa_lw_philic(:,irh,ao%itype(jtype)) |
---|
| 267 | scat_g_lw_aerosol = scat_g_lw_aerosol & |
---|
| 268 | & + local_od_lw * ao%ssa_lw_philic(:,irh,ao%itype(jtype)) & |
---|
| 269 | & * ao%g_lw_philic(:,irh,ao%itype(jtype)) |
---|
| 270 | else |
---|
| 271 | ! If aerosol longwave scattering is not included then we |
---|
| 272 | ! weight the optical depth by the single scattering |
---|
| 273 | ! co-albedo |
---|
| 274 | od_lw_aerosol = od_lw_aerosol & |
---|
| 275 | & + factor * aerosol%mixing_ratio(jcol,jlev,jtype) & |
---|
| 276 | & * ao%mass_ext_lw_philic(:,irh,ao%itype(jtype)) & |
---|
| 277 | & * (1.0_jprb - ao%ssa_lw_philic(:,irh,ao%itype(jtype))) |
---|
| 278 | end if |
---|
| 279 | end if |
---|
| 280 | ! Implicitly, if ao%iclass(jtype) == IAerosolClassNone, then |
---|
| 281 | ! no aerosol scattering properties are added |
---|
| 282 | |
---|
| 283 | end do ! Loop over aerosol type |
---|
| 284 | |
---|
| 285 | if (.not. config%do_sw_delta_scaling_with_gases) then |
---|
| 286 | ! Delta-Eddington scaling on aerosol only. Note that if |
---|
| 287 | ! do_sw_delta_scaling_with_gases==.true. then the delta |
---|
| 288 | ! scaling is done to the cloud-aerosol-gas mixture inside |
---|
| 289 | ! the solver |
---|
| 290 | call delta_eddington_extensive(od_sw_aerosol, scat_sw_aerosol, & |
---|
| 291 | & scat_g_sw_aerosol) |
---|
| 292 | end if |
---|
| 293 | |
---|
| 294 | ! Combine aerosol shortwave scattering properties with gas |
---|
| 295 | ! properties (noting that any gas scattering will have an |
---|
| 296 | ! asymmetry factor of zero) |
---|
| 297 | if (od_sw_aerosol(1) > 0.0_jprb) then |
---|
| 298 | do jg = 1,config%n_g_sw |
---|
| 299 | iband = config%i_band_from_reordered_g_sw(jg) |
---|
| 300 | local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband) |
---|
| 301 | local_scat = ssa_sw(jg,jlev,jcol) * od_sw(jg,jlev,jcol) & |
---|
| 302 | & + scat_sw_aerosol(iband) |
---|
| 303 | ! Note that asymmetry_sw of gases is zero so the following |
---|
| 304 | ! simply weights the aerosol asymmetry by the scattering |
---|
| 305 | ! optical depth |
---|
| 306 | g_sw(jg,jlev,jcol) = scat_g_sw_aerosol(iband) / local_scat |
---|
| 307 | ssa_sw(jg,jlev,jcol) = local_scat / local_od |
---|
| 308 | od_sw (jg,jlev,jcol) = local_od |
---|
| 309 | end do |
---|
| 310 | end if |
---|
| 311 | |
---|
| 312 | ! Combine aerosol longwave scattering properties with gas |
---|
| 313 | ! properties, noting that in the longwave, gases do not |
---|
| 314 | ! scatter at all |
---|
| 315 | if (config%do_lw_aerosol_scattering) then |
---|
| 316 | |
---|
| 317 | call delta_eddington_extensive(od_lw_aerosol, scat_lw_aerosol, & |
---|
| 318 | & scat_g_lw_aerosol) |
---|
| 319 | |
---|
| 320 | do jg = 1,config%n_g_lw |
---|
| 321 | iband = config%i_band_from_reordered_g_lw(jg) |
---|
| 322 | if (od_lw_aerosol(iband) > 0.0_jprb) then |
---|
| 323 | ! All scattering is due to aerosols, therefore the |
---|
| 324 | ! asymmetry factor is equal to the value for aerosols |
---|
| 325 | if (scat_lw_aerosol(iband) > 0.0_jprb) then |
---|
| 326 | g_lw(jg,jlev,jcol) = scat_g_lw_aerosol(iband) & |
---|
| 327 | & / scat_lw_aerosol(iband) |
---|
| 328 | else |
---|
| 329 | g_lw(jg,jlev,jcol) = 0.0_jprb |
---|
| 330 | end if |
---|
| 331 | local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband) |
---|
| 332 | ssa_lw(jg,jlev,jcol) = scat_lw_aerosol(iband) / local_od |
---|
| 333 | od_lw (jg,jlev,jcol) = local_od |
---|
| 334 | end if |
---|
| 335 | end do |
---|
| 336 | else |
---|
| 337 | do jg = 1,config%n_g_lw |
---|
| 338 | od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) & |
---|
| 339 | & + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg)) |
---|
| 340 | end do |
---|
| 341 | end if |
---|
| 342 | |
---|
| 343 | end do ! Loop over column |
---|
| 344 | end do ! Loop over level |
---|
| 345 | |
---|
| 346 | end if |
---|
| 347 | |
---|
| 348 | if (lhook) call dr_hook('radiation_aerosol_optics:add_aerosol_optics',1,hook_handle) |
---|
| 349 | |
---|
| 350 | end subroutine add_aerosol_optics |
---|
| 351 | |
---|
| 352 | |
---|
| 353 | !--------------------------------------------------------------------- |
---|
| 354 | ! Add precomputed optical properties to gas optical depth and |
---|
| 355 | ! scattering properties |
---|
| 356 | subroutine add_aerosol_optics_direct(nlev,istartcol,iendcol, & |
---|
| 357 | & config, aerosol, & |
---|
| 358 | & od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw) |
---|
| 359 | |
---|
| 360 | use parkind1, only : jprb |
---|
| 361 | use radiation_io, only : nulerr, radiation_abort |
---|
| 362 | use yomhook, only : lhook, dr_hook |
---|
| 363 | use radiation_config, only : config_type |
---|
| 364 | use radiation_aerosol, only : aerosol_type |
---|
| 365 | |
---|
| 366 | integer, intent(in) :: nlev ! number of model levels |
---|
| 367 | integer, intent(in) :: istartcol, iendcol ! range of columns to process |
---|
| 368 | type(config_type), intent(in), target :: config |
---|
| 369 | type(aerosol_type), intent(in) :: aerosol |
---|
| 370 | ! Optical depth, single scattering albedo and asymmetry factor of |
---|
| 371 | ! the atmosphere (gases on input, gases and aerosols on output) |
---|
| 372 | ! for each g point. Note that longwave ssa and asymmetry and |
---|
| 373 | ! shortwave asymmetry are all zero for gases, so are not yet |
---|
| 374 | ! defined on input and are therefore intent(out). |
---|
| 375 | real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), & |
---|
| 376 | & intent(inout) :: od_lw |
---|
| 377 | real(jprb), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol), & |
---|
| 378 | & intent(out) :: ssa_lw, g_lw |
---|
| 379 | real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), & |
---|
| 380 | & intent(inout) :: od_sw, ssa_sw |
---|
| 381 | real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), & |
---|
| 382 | & intent(out) :: g_sw |
---|
| 383 | |
---|
| 384 | ! Temporary extinction and scattering optical depths of aerosol |
---|
| 385 | ! plus gas |
---|
| 386 | real(jprb) :: local_od, local_scat |
---|
| 387 | |
---|
| 388 | ! Extinction optical depth, scattering optical depth and |
---|
| 389 | ! asymmetry-times-scattering-optical-depth for all the aerosols at |
---|
| 390 | ! a point in space for each spectral band of the shortwave and |
---|
| 391 | ! longwave spectrum |
---|
| 392 | real(jprb), dimension(config%n_bands_sw) & |
---|
| 393 | & :: od_sw_aerosol, scat_sw_aerosol, scat_g_sw_aerosol |
---|
| 394 | real(jprb), dimension(config%n_bands_lw) :: od_lw_aerosol |
---|
| 395 | real(jprb), dimension(config%n_bands_lw_if_scattering) & |
---|
| 396 | & :: scat_lw_aerosol, scat_g_lw_aerosol |
---|
| 397 | |
---|
| 398 | ! Loop indices for column, level, g point and band |
---|
| 399 | integer :: jcol, jlev, jg |
---|
| 400 | |
---|
| 401 | ! Range of levels over which aerosols are present |
---|
| 402 | integer :: istartlev, iendlev |
---|
| 403 | |
---|
| 404 | ! Indices to spectral band |
---|
| 405 | integer :: iband |
---|
| 406 | |
---|
| 407 | real(jprb) :: hook_handle |
---|
| 408 | |
---|
| 409 | if (lhook) call dr_hook('radiation_aerosol_optics:add_aerosol_optics_direct',0,hook_handle) |
---|
| 410 | |
---|
| 411 | if (config%do_sw) then |
---|
| 412 | ! Check array dimensions |
---|
| 413 | if (ubound(aerosol%od_sw,1) /= config%n_bands_sw) then |
---|
| 414 | write(nulerr,'(a,i0,a,i0)') '*** Error: aerosol%od_sw contains ', & |
---|
| 415 | & ubound(aerosol%od_sw,1), ' band, expected ', & |
---|
| 416 | & config%n_bands_sw |
---|
| 417 | call radiation_abort() |
---|
| 418 | end if |
---|
| 419 | |
---|
| 420 | istartlev = lbound(aerosol%od_sw,2) |
---|
| 421 | iendlev = ubound(aerosol%od_sw,2) |
---|
| 422 | |
---|
| 423 | ! Set variables to zero that may not have been previously |
---|
| 424 | g_sw = 0.0_jprb |
---|
| 425 | |
---|
| 426 | ! Loop over position |
---|
| 427 | do jcol = istartcol,iendcol |
---|
| 428 | do jlev = istartlev,iendlev |
---|
| 429 | od_sw_aerosol = aerosol%od_sw(:,jlev,jcol) |
---|
| 430 | scat_sw_aerosol = aerosol%ssa_sw(:,jlev,jcol) * od_sw_aerosol |
---|
| 431 | scat_g_sw_aerosol = aerosol%g_sw(:,jlev,jcol) * scat_sw_aerosol |
---|
| 432 | |
---|
| 433 | if (.not. config%do_sw_delta_scaling_with_gases) then |
---|
| 434 | ! Delta-Eddington scaling on aerosol only. Note that if |
---|
| 435 | ! do_sw_delta_scaling_with_gases==.true. then the delta |
---|
| 436 | ! scaling is done to the cloud-aerosol-gas mixture inside |
---|
| 437 | ! the solver |
---|
| 438 | call delta_eddington_extensive(od_sw_aerosol, scat_sw_aerosol, & |
---|
| 439 | & scat_g_sw_aerosol) |
---|
| 440 | end if |
---|
| 441 | |
---|
| 442 | ! Combine aerosol shortwave scattering properties with gas |
---|
| 443 | ! properties (noting that any gas scattering will have an |
---|
| 444 | ! asymmetry factor of zero) |
---|
| 445 | if (od_sw_aerosol(1) > 0.0_jprb) then |
---|
| 446 | do jg = 1,config%n_g_sw |
---|
| 447 | iband = config%i_band_from_reordered_g_sw(jg) |
---|
| 448 | local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband) |
---|
| 449 | local_scat = ssa_sw(jg,jlev,jcol) * od_sw(jg,jlev,jcol) & |
---|
| 450 | & + scat_sw_aerosol(iband) |
---|
| 451 | ! Note that asymmetry_sw of gases is zero so the following |
---|
| 452 | ! simply weights the aerosol asymmetry by the scattering |
---|
| 453 | ! optical depth |
---|
| 454 | g_sw(jg,jlev,jcol) = scat_g_sw_aerosol(iband) / local_scat |
---|
| 455 | ssa_sw(jg,jlev,jcol) = local_scat / local_od |
---|
| 456 | od_sw (jg,jlev,jcol) = local_od |
---|
| 457 | end do |
---|
| 458 | end if |
---|
| 459 | end do |
---|
| 460 | end do |
---|
| 461 | |
---|
| 462 | end if |
---|
| 463 | |
---|
| 464 | if (config%do_lw) then |
---|
| 465 | |
---|
| 466 | if (ubound(aerosol%od_lw,1) /= config%n_bands_lw) then |
---|
| 467 | write(nulerr,'(a,i0,a,i0)') '*** Error: aerosol%od_lw contains ', & |
---|
| 468 | & ubound(aerosol%od_lw,1), ' band, expected ', & |
---|
| 469 | & config%n_bands_lw |
---|
| 470 | call radiation_abort() |
---|
| 471 | end if |
---|
| 472 | |
---|
| 473 | istartlev = lbound(aerosol%od_lw,2) |
---|
| 474 | iendlev = ubound(aerosol%od_lw,2) |
---|
| 475 | |
---|
| 476 | if (config%do_lw_aerosol_scattering) then |
---|
| 477 | ssa_lw = 0.0_jprb |
---|
| 478 | g_lw = 0.0_jprb |
---|
| 479 | |
---|
| 480 | ! Loop over position |
---|
| 481 | do jcol = istartcol,iendcol |
---|
| 482 | do jlev = istartlev,iendlev |
---|
| 483 | od_lw_aerosol = aerosol%od_lw(:,jlev,jcol) |
---|
| 484 | scat_lw_aerosol = aerosol%ssa_lw(:,jlev,jcol) * od_lw_aerosol |
---|
| 485 | scat_g_lw_aerosol = aerosol%g_lw(:,jlev,jcol) * scat_lw_aerosol |
---|
| 486 | |
---|
| 487 | call delta_eddington_extensive(od_lw_aerosol, scat_lw_aerosol, & |
---|
| 488 | & scat_g_lw_aerosol) |
---|
| 489 | |
---|
| 490 | do jg = 1,config%n_g_lw |
---|
| 491 | iband = config%i_band_from_reordered_g_lw(jg) |
---|
| 492 | if (od_lw_aerosol(iband) > 0.0_jprb) then |
---|
| 493 | ! All scattering is due to aerosols, therefore the |
---|
| 494 | ! asymmetry factor is equal to the value for aerosols |
---|
| 495 | if (scat_lw_aerosol(iband) > 0.0_jprb) then |
---|
| 496 | g_lw(jg,jlev,jcol) = scat_g_lw_aerosol(iband) & |
---|
| 497 | & / scat_lw_aerosol(iband) |
---|
| 498 | else |
---|
| 499 | g_lw(jg,jlev,jcol) = 0.0_jprb |
---|
| 500 | end if |
---|
| 501 | local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband) |
---|
| 502 | ssa_lw(jg,jlev,jcol) = scat_lw_aerosol(iband) / local_od |
---|
| 503 | od_lw (jg,jlev,jcol) = local_od |
---|
| 504 | end if |
---|
| 505 | end do |
---|
| 506 | end do |
---|
| 507 | end do |
---|
| 508 | |
---|
| 509 | else ! No longwave scattering |
---|
| 510 | |
---|
| 511 | ! Loop over position |
---|
| 512 | do jcol = istartcol,iendcol |
---|
| 513 | do jlev = istartlev,iendlev |
---|
| 514 | ! If aerosol longwave scattering is not included then we |
---|
| 515 | ! weight the optical depth by the single scattering |
---|
| 516 | ! co-albedo |
---|
| 517 | od_lw_aerosol = aerosol%od_lw(:,jlev,jcol) & |
---|
| 518 | & * (1.0_jprb - aerosol%ssa_lw(:,jlev,jcol)) |
---|
| 519 | do jg = 1,config%n_g_lw |
---|
| 520 | od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) & |
---|
| 521 | & + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg)) |
---|
| 522 | end do |
---|
| 523 | end do |
---|
| 524 | end do |
---|
| 525 | |
---|
| 526 | end if |
---|
| 527 | end if |
---|
| 528 | |
---|
| 529 | |
---|
| 530 | if (lhook) call dr_hook('radiation_aerosol_optics:add_aerosol_optics_direct',1,hook_handle) |
---|
| 531 | |
---|
| 532 | end subroutine add_aerosol_optics_direct |
---|
| 533 | |
---|
| 534 | |
---|
| 535 | !--------------------------------------------------------------------- |
---|
| 536 | ! Sometimes it is useful to specify aerosol in terms of its optical |
---|
| 537 | ! depth at a particular wavelength. This function returns the dry |
---|
| 538 | ! shortwave mass-extinction coefficient, i.e. the extinction cross |
---|
| 539 | ! section per unit mass, for aerosol of type "itype" at shortwave |
---|
| 540 | ! band "iband". For hydrophilic types, the value at the first |
---|
| 541 | ! relative humidity bin is taken. |
---|
| 542 | function dry_aerosol_sw_mass_extinction(config, itype, iband) |
---|
| 543 | |
---|
| 544 | use parkind1, only : jprb |
---|
| 545 | use radiation_config, only : config_type |
---|
| 546 | use radiation_aerosol_optics_data, only : aerosol_optics_type, & |
---|
| 547 | & IAerosolClassUndefined, IAerosolClassIgnored, & |
---|
| 548 | & IAerosolClassHydrophobic, IAerosolClassHydrophilic |
---|
| 549 | |
---|
| 550 | type(config_type), intent(in), target :: config |
---|
| 551 | |
---|
| 552 | ! Aerosol type and shortwave band as indices to the array |
---|
| 553 | integer, intent(in) :: itype, iband |
---|
| 554 | |
---|
| 555 | real(jprb) dry_aerosol_sw_mass_extinction |
---|
| 556 | |
---|
| 557 | ! Pointer to the aerosol optics coefficients for brevity of access |
---|
| 558 | type(aerosol_optics_type), pointer :: ao |
---|
| 559 | |
---|
| 560 | ao => config%aerosol_optics |
---|
| 561 | |
---|
| 562 | if (ao%iclass(itype) == IAerosolClassHydrophobic) then |
---|
| 563 | dry_aerosol_sw_mass_extinction = ao%mass_ext_sw_phobic(iband,ao%itype(itype)) |
---|
| 564 | else if (ao%iclass(itype) == IAerosolClassHydrophilic) then |
---|
| 565 | ! Take the value at the first relative-humidity bin for the |
---|
| 566 | ! "dry" aerosol value |
---|
| 567 | dry_aerosol_sw_mass_extinction = ao%mass_ext_sw_philic(iband,1,ao%itype(itype)) |
---|
| 568 | else |
---|
| 569 | dry_aerosol_sw_mass_extinction = 0.0_jprb |
---|
| 570 | end if |
---|
| 571 | |
---|
| 572 | end function dry_aerosol_sw_mass_extinction |
---|
| 573 | |
---|
| 574 | |
---|
| 575 | !--------------------------------------------------------------------- |
---|
| 576 | ! Compute aerosol extinction coefficient at a particular shortwave |
---|
| 577 | ! band and a single height - this is useful for visibility |
---|
| 578 | ! diagnostics |
---|
| 579 | subroutine aerosol_sw_extinction(ncol,istartcol,iendcol, & |
---|
| 580 | & config, iband, mixing_ratio, relative_humidity, extinction) |
---|
| 581 | |
---|
| 582 | use parkind1, only : jprb |
---|
| 583 | use yomhook, only : lhook, dr_hook |
---|
| 584 | use radiation_io, only : nulerr, radiation_abort |
---|
| 585 | use radiation_config, only : config_type |
---|
| 586 | use radiation_aerosol_optics_data, only : aerosol_optics_type, & |
---|
| 587 | & IAerosolClassUndefined, IAerosolClassIgnored, & |
---|
| 588 | & IAerosolClassHydrophobic, IAerosolClassHydrophilic |
---|
| 589 | |
---|
| 590 | integer, intent(in) :: ncol ! number of columns |
---|
| 591 | integer, intent(in) :: istartcol, iendcol ! range of columns to process |
---|
| 592 | type(config_type), intent(in), target :: config |
---|
| 593 | integer, intent(in) :: iband ! Index of required spectral band |
---|
| 594 | real(jprb), intent(in) :: mixing_ratio(ncol,config%n_aerosol_types) |
---|
| 595 | real(jprb), intent(in) :: relative_humidity(ncol) |
---|
| 596 | real(jprb), intent(out) :: extinction(ncol) |
---|
| 597 | |
---|
| 598 | ! Local aerosol extinction |
---|
| 599 | real(jprb) :: ext |
---|
| 600 | |
---|
| 601 | ! Pointer to the aerosol optics coefficients for brevity of access |
---|
| 602 | type(aerosol_optics_type), pointer :: ao |
---|
| 603 | |
---|
| 604 | ! Loop indices for column and aerosol type |
---|
| 605 | integer :: jcol, jtype |
---|
| 606 | |
---|
| 607 | ! Relative humidity index |
---|
| 608 | integer :: irh |
---|
| 609 | |
---|
| 610 | real(jprb) :: hook_handle |
---|
| 611 | |
---|
| 612 | if (lhook) call dr_hook('radiation_aerosol_optics:aerosol_sw_extinction',0,hook_handle) |
---|
| 613 | |
---|
| 614 | do jtype = 1,config%n_aerosol_types |
---|
| 615 | if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then |
---|
| 616 | write(nulerr,'(a)') '*** Error: not all aerosol types are defined' |
---|
| 617 | call radiation_abort() |
---|
| 618 | end if |
---|
| 619 | end do |
---|
| 620 | |
---|
| 621 | ao => config%aerosol_optics |
---|
| 622 | |
---|
| 623 | ! Loop over position |
---|
| 624 | do jcol = istartcol,iendcol |
---|
| 625 | ext = 0.0_jprb |
---|
| 626 | ! Get relative-humidity index |
---|
| 627 | irh = ao%calc_rh_index(relative_humidity(jcol)) |
---|
| 628 | ! Add extinction coefficients from each aerosol type |
---|
| 629 | do jtype = 1,config%n_aerosol_types |
---|
| 630 | if (ao%iclass(jtype) == IAerosolClassHydrophobic) then |
---|
| 631 | ext = ext + mixing_ratio(jcol,jtype) & |
---|
| 632 | & * ao%mass_ext_sw_phobic(iband,ao%itype(jtype)) |
---|
| 633 | else if (ao%iclass(jtype) == IAerosolClassHydrophilic) then |
---|
| 634 | ext = ext + mixing_ratio(jcol,jtype) & |
---|
| 635 | & * ao%mass_ext_sw_philic(iband,irh,ao%itype(jtype)) |
---|
| 636 | end if |
---|
| 637 | end do |
---|
| 638 | |
---|
| 639 | extinction(jcol) = ext |
---|
| 640 | end do |
---|
| 641 | |
---|
| 642 | if (lhook) call dr_hook('radiation_aerosol_optics:aerosol_sw_extinction',1,hook_handle) |
---|
| 643 | |
---|
| 644 | end subroutine aerosol_sw_extinction |
---|
| 645 | |
---|
| 646 | end module radiation_aerosol_optics |
---|