Line data Source code
1 : module aer_rad_props
2 :
3 : !------------------------------------------------------------------------------------------------
4 : ! Converts aerosol masses to bulk optical properties for sw and lw radiation
5 : ! computations.
6 : !------------------------------------------------------------------------------------------------
7 :
8 : use shr_kind_mod, only: r8 => shr_kind_r8
9 : use ppgrid, only: pcols, pver
10 : use physconst, only: rga
11 : use physics_types, only: physics_state
12 :
13 : use physics_buffer, only: physics_buffer_desc
14 : use radconstants, only: nswbands, nlwbands, idx_sw_diag
15 : use phys_prop, only: nrh, ot_length
16 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, &
17 : rad_cnst_get_aer_props
18 : use wv_saturation, only: qsat
19 : use aerosol_optics_cam,only: aerosol_optics_cam_init, aerosol_optics_cam_sw, aerosol_optics_cam_lw
20 : use cam_history, only: fieldname_len, addfld, outfld, add_default, horiz_only
21 : use cam_history_support, only : fillvalue
22 : ! Placed here due to PGI bug.
23 : use ref_pres, only: clim_modal_aero_top_lev
24 :
25 : use cam_abortutils, only: endrun
26 :
27 : implicit none
28 : private
29 : save
30 :
31 : integer :: top_lev = 1
32 :
33 : public :: &
34 : aer_rad_props_init, &
35 : aer_rad_props_sw, & ! return SW optical props of aerosols
36 : aer_rad_props_lw ! return LW optical props of aerosols
37 :
38 : ! Private data
39 : character(len=fieldname_len), pointer :: odv_names(:) ! outfld names for visible OD
40 :
41 :
42 : !==============================================================================
43 : contains
44 : !==============================================================================
45 :
46 1536 : subroutine aer_rad_props_init()
47 : use phys_control, only: phys_getopts
48 :
49 :
50 : integer :: i
51 : integer :: numaerosols ! number of aerosols
52 1536 : character(len=64), pointer :: aernames(:) ! aerosol names
53 : logical :: history_amwg ! output the variables used by the AMWG diag package
54 : logical :: history_aero_optics ! Output aerosol optics diagnostics
55 : logical :: history_dust ! Output dust diagnostics
56 : logical :: prog_modal_aero ! Prognostic modal aerosols present
57 : integer :: nmodes ! number of aerosol modes
58 : integer :: nbins ! number of aerosol bins
59 :
60 : !----------------------------------------------------------------------------
61 :
62 : call phys_getopts( history_aero_optics_out = history_aero_optics, &
63 : history_amwg_out = history_amwg, &
64 : history_dust_out = history_dust, &
65 1536 : prog_modal_aero_out = prog_modal_aero )
66 :
67 : ! Limit modal aerosols with top_lev here.
68 1536 : if (prog_modal_aero) top_lev = clim_modal_aero_top_lev
69 :
70 : call addfld ('AEROD_v ', horiz_only, 'A', '1', &
71 1536 : 'Total Aerosol Optical Depth in visible band', flag_xyfill=.true.)
72 :
73 : call addfld ('AODvstrt', horiz_only, 'A', '1', &
74 1536 : 'Stratospheric Aerosol Optical Depth in visible band', flag_xyfill=.true.)
75 :
76 : ! Contributions to AEROD_v from individual aerosols (climate species).
77 :
78 : ! number of bulk aerosols in climate list
79 1536 : call rad_cnst_get_info(0, naero=numaerosols)
80 :
81 : ! get names of bulk aerosols
82 3072 : allocate(aernames(numaerosols))
83 1536 : call rad_cnst_get_info(0, aernames=aernames, nmodes=nmodes, nbins=nbins)
84 :
85 : ! diagnostic output for bulk aerosols
86 : ! create outfld names for visible OD
87 3072 : allocate(odv_names(numaerosols))
88 1536 : do i = 1, numaerosols
89 0 : odv_names(i) = 'ODV_'//trim(aernames(i))
90 0 : call addfld (odv_names(i), horiz_only, 'A', '1', &
91 1536 : trim(aernames(i))//' optical depth in visible band', flag_xyfill=.true.)
92 : end do
93 :
94 : ! Determine default fields
95 1536 : if (history_amwg .or. history_dust ) then
96 1536 : call add_default ('AEROD_v', 1, ' ')
97 : endif
98 :
99 1536 : if ( history_aero_optics ) then
100 0 : call add_default ('AEROD_v', 1, ' ')
101 0 : do i = 1, numaerosols
102 0 : odv_names(i) = 'ODV_'//trim(aernames(i))
103 0 : call add_default (odv_names(i), 1, ' ')
104 : end do
105 : endif
106 :
107 1536 : if (nmodes>0 .or. nbins>0) then
108 1536 : call aerosol_optics_cam_init()
109 : end if
110 :
111 1536 : deallocate(aernames)
112 :
113 1536 : end subroutine aer_rad_props_init
114 :
115 : !==============================================================================
116 :
117 38400 : subroutine aer_rad_props_sw(list_idx, state, pbuf, nnite, idxnite, &
118 : tau, tau_w, tau_w_g, tau_w_f)
119 :
120 : ! Return bulk layer tau, omega, g, f for all spectral intervals.
121 :
122 : use physics_buffer, only : physics_buffer_desc
123 : use tropopause, only : tropopause_find_cam
124 : ! Arguments
125 : integer, intent(in) :: list_idx ! index of the climate or a diagnostic list
126 : type(physics_state), intent(in), target :: state
127 :
128 : type(physics_buffer_desc), pointer :: pbuf(:)
129 : integer, intent(in) :: nnite ! number of night columns
130 : integer, intent(in) :: idxnite(:) ! local column indices of night columns
131 :
132 : real(r8), intent(out) :: tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth
133 : real(r8), intent(out) :: tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau
134 : real(r8), intent(out) :: tau_w_g(pcols,0:pver,nswbands) ! aerosol asymmetry parameter * tau * w
135 : real(r8), intent(out) :: tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * tau * w
136 :
137 : ! Local variables
138 :
139 : integer :: ncol
140 : integer :: lchnk
141 : integer :: k ! index
142 : integer :: troplev(pcols)
143 :
144 : ! optical props for each aerosol
145 : ! hygroscopic
146 38400 : real(r8), pointer :: h_ext(:,:)
147 38400 : real(r8), pointer :: h_ssa(:,:)
148 38400 : real(r8), pointer :: h_asm(:,:)
149 : ! non-hygroscopic
150 38400 : real(r8), pointer :: n_ext(:)
151 38400 : real(r8), pointer :: n_ssa(:)
152 38400 : real(r8), pointer :: n_asm(:)
153 38400 : real(r8), pointer :: n_scat(:)
154 38400 : real(r8), pointer :: n_ascat(:)
155 : ! radius-dependent
156 38400 : real(r8), pointer :: r_ext(:,:) ! radius-dependent mass-specific extinction
157 38400 : real(r8), pointer :: r_scat(:,:)
158 38400 : real(r8), pointer :: r_ascat(:,:)
159 38400 : real(r8), pointer :: r_mu(:) ! log(radius) domain variable for r_ext, r_scat, r_ascat
160 :
161 : ! radiative properties for each aerosol
162 : real(r8) :: ta (pcols,pver,nswbands)
163 : real(r8) :: tw (pcols,pver,nswbands)
164 : real(r8) :: twf(pcols,pver,nswbands)
165 : real(r8) :: twg(pcols,pver,nswbands)
166 :
167 : ! aerosol masses
168 38400 : real(r8), pointer :: aermmr(:,:) ! mass mixing ratio of aerosols
169 : real(r8) :: mmr_to_mass(pcols,pver) ! conversion factor for mmr to mass
170 : real(r8) :: aermass(pcols,pver) ! mass of aerosols
171 :
172 : ! for table lookup into rh grid
173 : real(r8) :: es(pcols,pver) ! saturation vapor pressure
174 : real(r8) :: qs(pcols,pver) ! saturation specific humidity
175 : real(r8) :: rh(pcols,pver)
176 : real(r8) :: rhtrunc(pcols,pver)
177 : real(r8) :: wrh(pcols,pver)
178 : integer :: krh(pcols,pver)
179 :
180 : integer :: numaerosols ! number of bulk aerosols in climate/diagnostic list
181 : integer :: nmodes ! number of aerosol modes in climate/diagnostic list
182 : integer :: nbins ! number of aerosol bins in climate/diagnostic list
183 : integer :: iaerosol ! index into bulk aerosol list
184 :
185 : character(len=ot_length) :: opticstype ! hygro or nonhygro
186 : character(len=16) :: pbuf_fld
187 : !-----------------------------------------------------------------------------
188 :
189 38400 : ncol = state%ncol
190 38400 : lchnk = state%lchnk
191 :
192 : ! compute mixing ratio to mass conversion
193 1267200 : do k = 1, pver
194 18961920 : mmr_to_mass(:ncol,k) = rga * state%pdeldry(:ncol,k)
195 : enddo
196 :
197 : ! initialize to conditions that would cause failure
198 302169600 : tau (:,:,:) = -100._r8
199 302169600 : tau_w (:,:,:) = -100._r8
200 302169600 : tau_w_g (:,:,:) = -100._r8
201 302169600 : tau_w_f (:,:,:) = -100._r8
202 :
203 : ! top layer (ilev = 0) has no aerosol (ie tau = 0)
204 : ! also initialize rest of layers to accumulate od's
205 273784320 : tau (1:ncol,:,:) = 0._r8
206 273784320 : tau_w (1:ncol,:,:) = 0._r8
207 273784320 : tau_w_g(1:ncol,:,:) = 0._r8
208 273784320 : tau_w_f(1:ncol,:,:) = 0._r8
209 :
210 : ! calculate relative humidity for table lookup into rh grid
211 1267200 : do k = 1, pver
212 1267200 : call qsat(state%t(1:ncol,k), state%pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol)
213 : end do
214 18961920 : rh(1:ncol,1:pver) = state%q(1:ncol,1:pver,1) / qs(1:ncol,1:pver)
215 :
216 18961920 : rhtrunc(1:ncol,1:pver) = min(rh(1:ncol,1:pver),1._r8)
217 18961920 : krh(1:ncol,1:pver) = min(floor( rhtrunc(1:ncol,1:pver) * nrh ) + 1, nrh - 1) ! index into rh mesh
218 18961920 : wrh(1:ncol,1:pver) = rhtrunc(1:ncol,1:pver) * nrh - krh(1:ncol,1:pver) ! (-) weighting on left side values
219 :
220 : ! get number of bulk aerosols and number of modes in current list
221 38400 : call rad_cnst_get_info(list_idx, naero=numaerosols, nmodes=nmodes, nbins=nbins)
222 :
223 : ! Contributions from modal and bin aerosols.
224 38400 : if (nmodes>0 .or. nbins>0) then
225 : call aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, &
226 38400 : tau, tau_w, tau_w_g, tau_w_f)
227 : else
228 0 : tau (1:ncol,:,:) = 0._r8
229 0 : tau_w (1:ncol,:,:) = 0._r8
230 0 : tau_w_g(1:ncol,:,:) = 0._r8
231 0 : tau_w_f(1:ncol,:,:) = 0._r8
232 : end if
233 :
234 : !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
235 38400 : troplev = 0
236 : !REMOVECAM_END
237 38400 : call tropopause_find_cam(state, troplev)
238 :
239 : ! Contributions from bulk aerosols.
240 38400 : do iaerosol = 1, numaerosols
241 :
242 : ! get bulk aerosol mass mixing ratio
243 0 : call rad_cnst_get_aer_mmr(list_idx, iaerosol, state, pbuf, aermmr)
244 0 : aermass(1:ncol,1:top_lev-1) = 0._r8
245 0 : aermass(1:ncol,top_lev:pver) = aermmr(1:ncol,top_lev:pver) * mmr_to_mass(1:ncol,top_lev:pver)
246 :
247 : ! get optics type
248 0 : call rad_cnst_get_aer_props(list_idx, iaerosol, opticstype=opticstype)
249 :
250 0 : select case (trim(opticstype))
251 : case('hygro','hygroscopic','hygroscopi')
252 : ! get optical properties for hygroscopic aerosols
253 0 : call rad_cnst_get_aer_props(list_idx, iaerosol, sw_hygro_ext=h_ext, sw_hygro_ssa=h_ssa, sw_hygro_asm=h_asm)
254 0 : call get_hygro_rad_props(ncol, krh, wrh, aermass, h_ext, h_ssa, h_asm, ta, tw, twg, twf)
255 0 : tau (1:ncol,1:pver,:) = tau (1:ncol,1:pver,:) + ta (1:ncol,:,:)
256 0 : tau_w (1:ncol,1:pver,:) = tau_w (1:ncol,1:pver,:) + tw (1:ncol,:,:)
257 0 : tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:)
258 0 : tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:)
259 :
260 : case('nonhygro','insoluble ')
261 : ! get optical properties for non-hygroscopic aerosols
262 : call rad_cnst_get_aer_props(list_idx, iaerosol, sw_nonhygro_ext=n_ext, sw_nonhygro_ssa=n_ssa, &
263 0 : sw_nonhygro_asm=n_asm)
264 :
265 0 : call get_nonhygro_rad_props(ncol, aermass, n_ext, n_ssa, n_asm, ta, tw, twg, twf)
266 0 : tau (1:ncol,1:pver,:) = tau (1:ncol,1:pver,:) + ta (1:ncol,:,:)
267 0 : tau_w (1:ncol,1:pver,:) = tau_w (1:ncol,1:pver,:) + tw (1:ncol,:,:)
268 0 : tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:)
269 0 : tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:)
270 :
271 : case('volcanic')
272 : ! get optical properties for volcanic aerosols
273 : call rad_cnst_get_aer_props(list_idx, iaerosol, sw_nonhygro_ext=n_ext, sw_nonhygro_scat=n_scat, &
274 0 : sw_nonhygro_ascat=n_ascat)
275 :
276 0 : call get_volcanic_rad_props(ncol, aermass, n_ext, n_scat, n_ascat, ta, tw, twg, twf)
277 0 : tau (1:ncol,1:pver,:) = tau (1:ncol,1:pver,:) + ta (1:ncol,:,:)
278 0 : tau_w (1:ncol,1:pver,:) = tau_w (1:ncol,1:pver,:) + tw (1:ncol,:,:)
279 0 : tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:)
280 0 : tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:)
281 :
282 : case('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3')
283 0 : pbuf_fld = 'VOLC_RAD_GEOM '
284 0 : if (len_trim(opticstype)>15) then
285 0 : pbuf_fld = trim(pbuf_fld)//opticstype(16:16)
286 : endif
287 : ! get optical properties for volcanic aerosols
288 0 : call rad_cnst_get_aer_props(list_idx, iaerosol, r_sw_ext=r_ext, r_sw_scat=r_scat, r_sw_ascat=r_ascat, mu=r_mu)
289 0 : call get_volcanic_radius_rad_props(ncol, aermass, pbuf_fld, pbuf, r_ext, r_scat, r_ascat, r_mu, ta, tw, twg, twf)
290 0 : tau (1:ncol,1:pver,:) = tau (1:ncol,1:pver,:) + ta (1:ncol,:,:)
291 0 : tau_w (1:ncol,1:pver,:) = tau_w (1:ncol,1:pver,:) + tw (1:ncol,:,:)
292 0 : tau_w_g(1:ncol,1:pver,:) = tau_w_g(1:ncol,1:pver,:) + twg(1:ncol,:,:)
293 0 : tau_w_f(1:ncol,1:pver,:) = tau_w_f(1:ncol,1:pver,:) + twf(1:ncol,:,:)
294 :
295 : case('zero')
296 : ! no effect of "zero" aerosols, so update nothing
297 : case default
298 0 : call endrun('aer_rad_props_sw: unsupported opticstype :'//trim(opticstype)//':')
299 : end select
300 :
301 : ! diagnostic output of individual aerosol optical properties
302 : ! currently implemented for climate list only
303 38400 : call aer_vis_diag_out(lchnk, ncol, nnite, idxnite, iaerosol, ta(:,:,idx_sw_diag), list_idx, troplev)
304 :
305 : enddo
306 :
307 : ! diagnostic output of total aerosol optical properties
308 : ! currently implemented for climate list only
309 38400 : call aer_vis_diag_out(lchnk, ncol, nnite, idxnite, 0, tau(:,:,idx_sw_diag), list_idx, troplev)
310 :
311 76800 : end subroutine aer_rad_props_sw
312 :
313 : !==============================================================================
314 :
315 38400 : subroutine aer_rad_props_lw(list_idx, state, pbuf, odap_aer)
316 :
317 : ! Purpose: Compute aerosol transmissions needed in absorptivity/
318 : ! emissivity calculations
319 :
320 : ! lw extinction is the same representation for all
321 : ! species. If this changes, this routine will need to do something
322 : ! similar to the sw with routines like get_hygro_lw_abs
323 :
324 38400 : use physics_buffer, only : pbuf_get_field, pbuf_get_index, physics_buffer_desc
325 :
326 : ! Arguments
327 : integer, intent(in) :: list_idx ! index of the climate or a diagnostic list
328 : type(physics_state), intent(in), target :: state
329 :
330 : type(physics_buffer_desc), pointer :: pbuf(:)
331 : real(r8), intent(out) :: odap_aer(pcols,pver,nlwbands) ! [fraction] absorption optical depth, per layer
332 :
333 : ! Local variables
334 :
335 : integer :: bnd_idx ! LW band index
336 : integer :: i ! column index
337 : integer :: k ! lev index
338 : integer :: ncol ! number of columns
339 : integer :: numaerosols ! number of bulk aerosols in climate/diagnostic list
340 : integer :: nmodes ! number of aerosol modes in climate/diagnostic list
341 : integer :: nbins ! number of aerosol bins in climate/diagnostic list
342 : integer :: iaerosol ! index into bulk aerosol list
343 : character(len=ot_length) :: opticstype ! hygro or nonhygro
344 :
345 : ! optical props for each aerosol
346 38400 : real(r8), pointer :: lw_abs(:)
347 38400 : real(r8), pointer :: lw_hygro_abs(:,:)
348 38400 : real(r8), pointer :: geometric_radius(:,:)
349 :
350 : ! volcanic lookup table
351 38400 : real(r8), pointer :: r_lw_abs(:,:) ! radius dependent mass-specific absorption coefficient
352 38400 : real(r8), pointer :: r_mu(:) ! log(geometric_mean_radius) domain samples of r_lw_abs(:,:)
353 : integer :: idx ! index to pbuf for geometric radius
354 : real(r8) :: mu(pcols,pver) ! log(geometric_radius)
355 : real(r8) :: r_mu_min, r_mu_max, wmu, mutrunc
356 : integer :: nmu, kmu
357 :
358 : ! for table lookup into rh grid
359 : real(r8) :: es(pcols,pver) ! saturation vapor pressure
360 : real(r8) :: qs(pcols,pver) ! saturation specific humidity
361 : real(r8) :: rh(pcols,pver)
362 : real(r8) :: rhtrunc(pcols,pver)
363 : real(r8) :: wrh(pcols,pver)
364 : integer :: krh(pcols,pver)
365 :
366 : ! aerosol (vertical) mass path and extinction
367 : ! aerosol masses
368 38400 : real(r8), pointer :: aermmr(:,:) ! mass mixing ratio of aerosols
369 : real(r8) :: mmr_to_mass(pcols,pver) ! conversion factor for mmr to mass
370 : real(r8) :: aermass(pcols,pver) ! mass of aerosols
371 :
372 : character(len=16) :: pbuf_fld
373 : !-----------------------------------------------------------------------------
374 :
375 38400 : ncol = state%ncol
376 :
377 : ! get number of bulk aerosols and number of modes in current list
378 38400 : call rad_cnst_get_info(list_idx, naero=numaerosols, nmodes=nmodes, nbins=nbins)
379 :
380 : ! Contributions from modal and sectional aerosols.
381 38400 : if (nmodes>0 .or. nbins>0) then
382 38400 : call aerosol_optics_cam_lw(list_idx, state, pbuf, odap_aer)
383 : else
384 0 : odap_aer = 0._r8
385 : end if
386 :
387 : ! Contributions from bulk aerosols.
388 38400 : if (numaerosols > 0) then
389 :
390 : ! compute mixing ratio to mass conversion
391 0 : do k = 1, pver
392 0 : mmr_to_mass(:ncol,k) = rga * state%pdeldry(:ncol,k)
393 : end do
394 :
395 : ! calculate relative humidity for table lookup into rh grid
396 0 : do k = 1, pver
397 0 : call qsat(state%t(1:ncol,k), state%pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol)
398 : end do
399 0 : rh(1:ncol,1:pver) = state%q(1:ncol,1:pver,1) / qs(1:ncol,1:pver)
400 :
401 0 : rhtrunc(1:ncol,1:pver) = min(rh(1:ncol,1:pver),1._r8)
402 0 : krh(1:ncol,1:pver) = min(floor( rhtrunc(1:ncol,1:pver) * nrh ) + 1, nrh - 1) ! index into rh mesh
403 0 : wrh(1:ncol,1:pver) = rhtrunc(1:ncol,1:pver) * nrh - krh(1:ncol,1:pver) ! (-) weighting on left side values
404 :
405 : end if
406 :
407 : ! Loop over bulk aerosols in list.
408 38400 : do iaerosol = 1, numaerosols
409 :
410 : ! get aerosol mass mixing ratio
411 0 : call rad_cnst_get_aer_mmr(list_idx, iaerosol, state, pbuf, aermmr)
412 0 : aermass(1:ncol,1:top_lev-1) = 0._r8
413 0 : aermass(1:ncol,top_lev:pver) = aermmr(1:ncol,top_lev:pver) * mmr_to_mass(1:ncol,top_lev:pver)
414 :
415 : ! get optics type
416 0 : call rad_cnst_get_aer_props(list_idx, iaerosol, opticstype=opticstype)
417 38400 : select case (trim(opticstype))
418 : case('hygroscopic')
419 : ! get optical properties for hygroscopic aerosols
420 0 : call rad_cnst_get_aer_props(list_idx, iaerosol, lw_hygro_ext=lw_hygro_abs)
421 0 : do bnd_idx = 1, nlwbands
422 0 : do k = 1, pver
423 0 : do i = 1, ncol
424 0 : odap_aer(i, k, bnd_idx) = odap_aer(i, k, bnd_idx) + &
425 : aermass(i, k) * &
426 0 : ((1 + wrh(i,k)) * lw_hygro_abs(krh(i,k)+1,bnd_idx) &
427 0 : - wrh(i,k) * lw_hygro_abs(krh(i,k), bnd_idx))
428 : end do
429 : end do
430 : end do
431 : case('insoluble','nonhygro','hygro','volcanic')
432 : ! get optical properties for hygroscopic aerosols
433 0 : call rad_cnst_get_aer_props(list_idx, iaerosol, lw_ext=lw_abs)
434 0 : do bnd_idx = 1, nlwbands
435 0 : do k = 1, pver
436 0 : do i = 1, ncol
437 0 : odap_aer(i,k,bnd_idx) = odap_aer(i,k,bnd_idx) + lw_abs(bnd_idx)*aermass(i,k)
438 : end do
439 : end do
440 : end do
441 :
442 : case('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3')
443 0 : pbuf_fld = 'VOLC_RAD_GEOM '
444 0 : if (len_trim(opticstype)>15) then
445 0 : pbuf_fld = trim(pbuf_fld)//opticstype(16:16)
446 : endif
447 :
448 : ! get optical properties for hygroscopic aerosols
449 0 : call rad_cnst_get_aer_props(list_idx, iaerosol, r_lw_abs=r_lw_abs, mu=r_mu)
450 : ! get microphysical properties for volcanic aerosols
451 0 : idx = pbuf_get_index(pbuf_fld)
452 0 : call pbuf_get_field(pbuf, idx, geometric_radius )
453 :
454 : ! interpolate in radius
455 : ! caution: clip the table with no warning when outside bounds
456 0 : nmu = size(r_mu)
457 0 : r_mu_max = r_mu(nmu)
458 0 : r_mu_min = r_mu(1)
459 0 : do i = 1, ncol
460 0 : do k = 1, pver
461 0 : if(geometric_radius(i,k) > 0._r8) then
462 0 : mu(i,k) = log(geometric_radius(i,k))
463 : else
464 0 : mu(i,k) = 0._r8
465 : endif
466 0 : mutrunc = max(min(mu(i,k),r_mu_max),r_mu_min)
467 0 : kmu = max(min(1 + (mutrunc-r_mu_min)/(r_mu_max-r_mu_min)*(nmu-1),nmu-1._r8),1._r8)
468 0 : wmu = max(min( (mutrunc -r_mu(kmu)) / (r_mu(kmu+1) - r_mu(kmu)) ,1._r8),0._r8)
469 0 : do bnd_idx = 1, nlwbands
470 0 : odap_aer(i,k,bnd_idx) = odap_aer(i,k,bnd_idx) + &
471 : aermass(i,k) * &
472 0 : ((1._r8 - wmu) * r_lw_abs(bnd_idx, kmu ) + &
473 0 : (wmu) * r_lw_abs(bnd_idx, kmu+1))
474 : end do
475 : end do
476 : end do
477 :
478 : case('zero')
479 : ! zero aerosols types have no optical effect, so do nothing.
480 : case default
481 0 : call endrun('aer_rad_props_lw: unsupported opticstype: '//trim(opticstype))
482 : end select
483 : end do
484 :
485 76800 : end subroutine aer_rad_props_lw
486 :
487 : !==============================================================================
488 : ! Private methods
489 : !==============================================================================
490 :
491 0 : subroutine get_hygro_rad_props(ncol, krh, wrh, mass, ext, ssa, asm, &
492 : tau, tau_w, tau_w_g, tau_w_f)
493 :
494 : ! Arguments
495 : integer, intent(in) :: ncol
496 : integer, intent(in) :: krh(pcols,pver) ! index for linear interpolation of optics on rh
497 : real(r8), intent(in) :: wrh(pcols,pver) ! weight for linear interpolation of optics on rh
498 : real(r8), intent(in) :: mass(pcols,pver)
499 : real(r8), intent(in) :: ext(:,:)
500 : real(r8), intent(in) :: ssa(:,:)
501 : real(r8), intent(in) :: asm(:,:)
502 :
503 : real(r8), intent(out) :: tau (pcols,pver,nswbands)
504 : real(r8), intent(out) :: tau_w (pcols,pver,nswbands)
505 : real(r8), intent(out) :: tau_w_g(pcols,pver,nswbands)
506 : real(r8), intent(out) :: tau_w_f(pcols,pver,nswbands)
507 :
508 : ! Local variables
509 : real(r8) :: ext1, ssa1, asm1
510 : integer :: icol, ilev, iswband
511 : !-----------------------------------------------------------------------------
512 :
513 0 : do iswband = 1, nswbands
514 0 : do icol = 1, ncol
515 0 : do ilev = 1, pver
516 0 : ext1 = (1 + wrh(icol,ilev)) * ext(krh(icol,ilev)+1,iswband) &
517 0 : - wrh(icol,ilev) * ext(krh(icol,ilev), iswband)
518 0 : ssa1 = (1 + wrh(icol,ilev)) * ssa(krh(icol,ilev)+1,iswband) &
519 0 : - wrh(icol,ilev) * ssa(krh(icol,ilev), iswband)
520 0 : asm1 = (1 + wrh(icol,ilev)) * asm(krh(icol,ilev)+1,iswband) &
521 0 : - wrh(icol,ilev) * asm(krh(icol,ilev), iswband)
522 :
523 0 : tau (icol, ilev, iswband) = mass(icol, ilev) * ext1
524 0 : tau_w (icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1
525 0 : tau_w_g(icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1 * asm1
526 0 : tau_w_f(icol, ilev, iswband) = mass(icol, ilev) * ext1 * ssa1 * asm1 * asm1
527 : enddo
528 : enddo
529 : enddo
530 :
531 38400 : end subroutine get_hygro_rad_props
532 :
533 : !==============================================================================
534 :
535 0 : subroutine get_nonhygro_rad_props(ncol, mass, ext, ssa, asm, &
536 : tau, tau_w, tau_w_g, tau_w_f)
537 :
538 : ! Arguments
539 : integer, intent(in) :: ncol
540 : real(r8), intent(in) :: mass(pcols, pver)
541 : real(r8), intent(in) :: ext(:)
542 : real(r8), intent(in) :: ssa(:)
543 : real(r8), intent(in) :: asm(:)
544 :
545 : real(r8), intent(out) :: tau (pcols, pver, nswbands)
546 : real(r8), intent(out) :: tau_w (pcols, pver, nswbands)
547 : real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands)
548 : real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands)
549 :
550 : ! Local variables
551 : integer :: iswband
552 : real(r8) :: ext1, ssa1, asm1
553 : !-----------------------------------------------------------------------------
554 :
555 0 : do iswband = 1, nswbands
556 0 : ext1 = ext(iswband)
557 0 : ssa1 = ssa(iswband)
558 0 : asm1 = asm(iswband)
559 0 : tau (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1
560 0 : tau_w (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1
561 0 : tau_w_g(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1 * asm1
562 0 : tau_w_f(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext1 * ssa1 * asm1 * asm1
563 : enddo
564 :
565 0 : end subroutine get_nonhygro_rad_props
566 :
567 : !==============================================================================
568 :
569 0 : subroutine get_volcanic_radius_rad_props(ncol, mass, pbuf_radius_name, pbuf, r_ext, r_scat, r_ascat, r_mu, &
570 : tau, tau_w, tau_w_g, tau_w_f)
571 :
572 :
573 : use physics_buffer, only : pbuf_get_field, pbuf_get_index
574 :
575 : ! Arguments
576 : integer, intent(in) :: ncol
577 : real(r8), intent(in) :: mass(pcols, pver)
578 : character(len=*) :: pbuf_radius_name
579 : type(physics_buffer_desc), pointer :: pbuf(:)
580 : real(r8), intent(in) :: r_ext(:,:)
581 : real(r8), intent(in) :: r_scat(:,:)
582 : real(r8), intent(in) :: r_ascat(:,:)
583 : real(r8), intent(in) :: r_mu(:) ! log(radius) domain of mass-specific optics
584 :
585 : real(r8), intent(out) :: tau (pcols, pver, nswbands)
586 : real(r8), intent(out) :: tau_w (pcols, pver, nswbands)
587 : real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands)
588 : real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands)
589 :
590 : ! Local variables
591 : integer :: iswband
592 : real(r8) :: g
593 :
594 : integer :: idx ! index to radius in physics buffer
595 0 : real(r8), pointer :: geometric_radius(:,:) ! geometric mean radius of volcanic aerosol
596 : real(r8) :: mu(pcols,pver) ! log(geometric mean radius of volcanic aerosol)
597 : integer :: kmu, nmu
598 : real(r8) :: wmu, mutrunc, r_mu_max, r_mu_min
599 :
600 : ! interpolated values from table
601 : real(r8) :: ext(nswbands)
602 : real(r8) :: scat(nswbands)
603 : real(r8) :: ascat(nswbands)
604 :
605 : integer :: i, k ! column level iterator
606 : !-----------------------------------------------------------------------------
607 :
608 0 : tau =0._r8
609 0 : tau_w =0._r8
610 0 : tau_w_g=0._r8
611 0 : tau_w_f=0._r8
612 :
613 : ! get microphysical properties for volcanic aerosols
614 0 : idx = pbuf_get_index(pbuf_radius_name)
615 0 : call pbuf_get_field(pbuf, idx, geometric_radius )
616 :
617 : ! interpolate in radius
618 : ! caution: clip the table with no warning when outside bounds
619 0 : nmu = size(r_mu)
620 0 : r_mu_max = r_mu(nmu)
621 0 : r_mu_min = r_mu(1)
622 0 : do i = 1, ncol
623 0 : do k = 1, pver
624 0 : if(geometric_radius(i,k) > 0._r8) then
625 0 : mu(i,k) = log(geometric_radius(i,k))
626 : else
627 0 : mu(i,k) = 0._r8
628 : endif
629 0 : mutrunc = max(min(mu(i,k),r_mu_max),r_mu_min)
630 0 : kmu = max(min(1 + (mutrunc-r_mu_min)/(r_mu_max-r_mu_min)*(nmu-1),nmu-1._r8),1._r8)
631 0 : wmu = max(min( (mutrunc -r_mu(kmu)) / (r_mu(kmu+1) - r_mu(kmu)) ,1._r8),0._r8)
632 0 : do iswband = 1, nswbands
633 0 : ext(iswband) = &
634 0 : ((1._r8 - wmu) * r_ext(iswband, kmu ) + &
635 0 : (wmu) * r_ext(iswband, kmu+1))
636 : scat(iswband) = &
637 0 : ((1._r8 - wmu) * r_scat(iswband, kmu ) + &
638 0 : (wmu) * r_scat(iswband, kmu+1))
639 : ascat(iswband) = &
640 0 : ((1._r8 - wmu) * r_ascat(iswband, kmu ) + &
641 0 : (wmu) * r_ascat(iswband, kmu+1))
642 0 : if (scat(iswband).gt.0._r8) then
643 0 : g = ascat(iswband)/scat(iswband)
644 : else
645 : g=0._r8
646 : endif
647 0 : tau (i,k,iswband) = mass(i,k) * ext(iswband)
648 0 : tau_w (i,k,iswband) = mass(i,k) * scat(iswband)
649 0 : tau_w_g(i,k,iswband) = mass(i,k) * ascat(iswband)
650 0 : tau_w_f(i,k,iswband) = mass(i,k) * g * ascat(iswband)
651 : end do
652 : enddo
653 : enddo
654 :
655 0 : end subroutine get_volcanic_radius_rad_props
656 :
657 : !==============================================================================
658 :
659 0 : subroutine get_volcanic_rad_props(ncol, mass, ext, scat, ascat, &
660 : tau, tau_w, tau_w_g, tau_w_f)
661 :
662 : ! Arguments
663 : integer, intent(in) :: ncol
664 : real(r8), intent(in) :: mass(pcols, pver)
665 : real(r8), intent(in) :: ext(:)
666 : real(r8), intent(in) :: scat(:)
667 : real(r8), intent(in) :: ascat(:)
668 :
669 : real(r8), intent(out) :: tau (pcols, pver, nswbands)
670 : real(r8), intent(out) :: tau_w (pcols, pver, nswbands)
671 : real(r8), intent(out) :: tau_w_g(pcols, pver, nswbands)
672 : real(r8), intent(out) :: tau_w_f(pcols, pver, nswbands)
673 :
674 : ! Local variables
675 : integer :: iswband
676 : real(r8) :: g
677 : !-----------------------------------------------------------------------------
678 :
679 0 : do iswband = 1, nswbands
680 0 : if (scat(iswband).gt.0._r8) then
681 0 : g = ascat(iswband)/scat(iswband)
682 : else
683 : g=0._r8
684 : endif
685 0 : tau (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ext(iswband)
686 0 : tau_w (1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * scat(iswband)
687 0 : tau_w_g(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * ascat(iswband)
688 0 : tau_w_f(1:ncol,1:pver,iswband) = mass(1:ncol,1:pver) * g * ascat(iswband)
689 : enddo
690 :
691 0 : end subroutine get_volcanic_rad_props
692 :
693 : !==============================================================================
694 :
695 38400 : subroutine aer_vis_diag_out(lchnk, ncol, nnite, idxnite, iaer, tau, diag_idx, troplev)
696 :
697 : ! output aerosol optical depth for the visible band
698 :
699 : integer, intent(in) :: lchnk
700 : integer, intent(in) :: ncol ! number of columns
701 : integer, intent(in) :: nnite ! number of night columns
702 : integer, intent(in) :: idxnite(:) ! local column indices of night columns
703 : integer, intent(in) :: iaer ! aerosol index -- if 0 then tau is a total for all aerosols
704 : real(r8), intent(in) :: tau(:,:) ! aerosol optical depth for the visible band
705 : integer, intent(in) :: diag_idx ! identifies whether the aerosol optics
706 : ! is for the climate calc or a diagnostic calc
707 : integer, intent(in) :: troplev(:) ! tropopause level
708 :
709 : ! Local variables
710 : integer :: i
711 : real(r8) :: tmp(pcols), tmp2(pcols)
712 : !-----------------------------------------------------------------------------
713 :
714 : ! currently only implemented for climate calc
715 38400 : if (diag_idx > 0) return
716 :
717 : ! compute total column aerosol optical depth
718 18839040 : tmp(1:ncol) = sum(tau(1:ncol,:), 2)
719 : ! use fillvalue to indicate night columns
720 314880 : do i = 1, nnite
721 314880 : tmp(idxnite(i)) = fillvalue
722 : end do
723 :
724 38400 : if (iaer > 0) then
725 0 : call outfld(odv_names(iaer), tmp, pcols, lchnk)
726 : else
727 38400 : call outfld('AEROD_v', tmp, pcols, lchnk)
728 591360 : do i = 1, ncol
729 8685989 : tmp2(i) = sum(tau(i,:troplev(i)))
730 : end do
731 38400 : call outfld('AODvstrt', tmp2, pcols, lchnk)
732 : end if
733 :
734 : end subroutine aer_vis_diag_out
735 :
736 : !==============================================================================
737 :
738 : end module aer_rad_props
|
