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