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 |
---|
38 | use setup_aerosol_optics_lmdz_m, only: setup_aerosol_optics_lmdz |
---|
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 |
---|
49 | call setup_aerosol_optics_lmdz(config%aerosol_optics, & |
---|
50 | trim(config%aerosol_optics_file_name)) |
---|
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 |
---|
91 | USE phys_local_var_mod, ONLY: rhcl |
---|
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 | |
---|
190 | !AI juin 2022 |
---|
191 | !call gas%get(IH2O, IMassMixingRatio, h2o_mmr, istartcol=istartcol) |
---|
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 |
---|
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 |
---|
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 |
---|