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