Line data Source code
1 : module refractive_aerosol_optics_mod
2 : use shr_kind_mod, only: r8 => shr_kind_r8
3 : use aerosol_optics_mod, only: aerosol_optics
4 : use physconst, only: rhoh2o
5 : use aerosol_state_mod, only: aerosol_state
6 : use aerosol_properties_mod, only: aerosol_properties
7 :
8 : use table_interp_mod, only: table_interp, table_interp_wghts, table_interp_calcwghts
9 :
10 : implicit none
11 :
12 : private
13 : public :: refractive_aerosol_optics
14 :
15 : !> refractive_aerosol_optics
16 : !! Table look up implementation of aerosol_optics to parameterize aerosol radiative properties in terms of
17 : !! surface mode wet radius and wet refractive index using chebychev polynomials
18 : type, extends(aerosol_optics) :: refractive_aerosol_optics
19 :
20 : integer :: ibin, ilist
21 : class(aerosol_state), pointer :: aero_state ! aerosol_state object
22 : class(aerosol_properties), pointer :: aero_props ! aerosol_properties object
23 :
24 : real(r8), allocatable :: watervol(:,:) ! volume concentration of water in each mode (m3/kg)
25 : real(r8), allocatable :: wetvol(:,:) ! volume concentration of wet mode (m3/kg)
26 : real(r8), allocatable :: cheb(:,:,:) ! chebychev polynomials
27 : real(r8), allocatable :: radsurf(:,:) ! aerosol surface mode radius
28 : real(r8), allocatable :: logradsurf(:,:) ! log(aerosol surface mode radius)
29 :
30 : ! refractive index for water read in read_water_refindex
31 : complex(r8), allocatable :: crefwsw(:) ! complex refractive index for water visible
32 : complex(r8), allocatable :: crefwlw(:) ! complex refractive index for water infrared
33 :
34 : real(r8), pointer :: extpsw(:,:,:,:) => null() ! specific extinction
35 : real(r8), pointer :: abspsw(:,:,:,:) => null() ! specific absorption
36 : real(r8), pointer :: asmpsw(:,:,:,:) => null() ! asymmetry factor
37 : real(r8), pointer :: absplw(:,:,:,:) => null() ! specific absorption
38 :
39 : real(r8), pointer :: refrtabsw(:,:) => null() ! table of real refractive indices for aerosols
40 : real(r8), pointer :: refitabsw(:,:) => null() ! table of imag refractive indices for aerosols
41 : real(r8), pointer :: refrtablw(:,:) => null() ! table of real refractive indices for aerosols
42 : real(r8), pointer :: refitablw(:,:) => null() ! table of imag refractive indices for aerosols
43 :
44 : ! Dimension sizes in coefficient arrays used to parameterize aerosol radiative properties
45 : ! in terms of refractive index and wet radius
46 : integer :: ncoef = -1 ! number of chebychev coeficients
47 : integer :: prefr = -1 ! number of real refractive indices
48 : integer :: prefi = -1 ! number of imaginary refractive indices
49 :
50 : contains
51 :
52 : procedure :: sw_props
53 : procedure :: lw_props
54 :
55 : final :: destructor
56 :
57 : end type refractive_aerosol_optics
58 :
59 : interface refractive_aerosol_optics
60 : procedure :: constructor
61 : end interface refractive_aerosol_optics
62 :
63 : ! radius limits (m)
64 : real(r8), parameter :: radmin = 0.01e-6_r8 ! min aerosol surface mode radius (m)
65 : real(r8), parameter :: radmax = 25.e-6_r8 ! max aerosol surface mode radius (m)
66 : real(r8), parameter :: xrmin=log(radmin) ! min log(aerosol surface mode radius)
67 : real(r8), parameter :: xrmax=log(radmax) ! max log(aerosol surface mode radius)
68 :
69 : contains
70 :
71 : !------------------------------------------------------------------------------
72 : !------------------------------------------------------------------------------
73 247680 : function constructor(aero_props, aero_state, ilist, ibin, ncol, nlev, nsw, nlw, crefwsw, crefwlw) &
74 : result(newobj)
75 :
76 : class(aerosol_properties),intent(in), target :: aero_props ! aerosol_properties object
77 : class(aerosol_state),intent(in), target :: aero_state ! aerosol_state object
78 : integer, intent(in) :: ilist ! climate or a diagnostic list number
79 : integer, intent(in) :: ibin ! bin number
80 : integer, intent(in) :: ncol ! number of columns
81 : integer, intent(in) :: nlev ! number of levels
82 : integer, intent(in) :: nsw ! number of short wave lengths
83 : integer, intent(in) :: nlw ! number of long wave lengths
84 : complex(r8), intent(in) :: crefwsw(nsw) ! complex refractive index for water visible
85 : complex(r8), intent(in) :: crefwlw(nlw) ! complex refractive index for water infrared
86 :
87 : type(refractive_aerosol_optics), pointer :: newobj
88 :
89 : integer :: ierr, icol, ilev, ispec, nspec
90 495360 : real(r8) :: vol(ncol) ! volume concentration of aerosol species (m3/kg)
91 495360 : real(r8) :: dryvol(ncol) ! volume concentration of aerosol mode (m3/kg)
92 : real(r8) :: specdens ! species density (kg/m3)
93 247680 : real(r8), pointer :: specmmr(:,:) ! species mass mixing ratio
94 : real(r8) :: logsigma ! geometric standard deviation of number distribution
95 :
96 495360 : real(r8) :: dgnumwet(ncol,nlev) ! aerosol wet number mode diameter (m)
97 495360 : real(r8) :: qaerwat(ncol,nlev) ! aerosol water (g/g)
98 :
99 : real(r8), parameter :: rh2odens = 1._r8/rhoh2o
100 :
101 247680 : allocate(newobj, stat=ierr)
102 247680 : if (ierr/=0) then
103 247680 : nullify(newobj)
104 : return
105 : end if
106 :
107 : ! get mode properties
108 : call aero_props%optics_params(ilist, ibin, &
109 : refrtabsw=newobj%refrtabsw, refitabsw=newobj%refitabsw, &
110 : refrtablw=newobj%refrtablw, refitablw=newobj%refitablw,&
111 : extpsw=newobj%extpsw, abspsw=newobj%abspsw, asmpsw=newobj%asmpsw, &
112 247680 : absplw=newobj%absplw, ncoef=newobj%ncoef, prefr=newobj%prefr, prefi=newobj%prefi)
113 :
114 990720 : allocate(newobj%watervol(ncol,nlev),stat=ierr)
115 247680 : if (ierr/=0) then
116 0 : nullify(newobj)
117 0 : return
118 : end if
119 743040 : allocate(newobj%wetvol(ncol,nlev),stat=ierr)
120 247680 : if (ierr/=0) then
121 0 : nullify(newobj)
122 0 : return
123 : end if
124 1238400 : allocate(newobj%cheb(newobj%ncoef,ncol,nlev),stat=ierr)
125 247680 : if (ierr/=0) then
126 0 : nullify(newobj)
127 0 : return
128 : end if
129 743040 : allocate(newobj%radsurf(ncol,nlev),stat=ierr)
130 247680 : if (ierr/=0) then
131 0 : nullify(newobj)
132 0 : return
133 : end if
134 743040 : allocate(newobj%logradsurf(ncol,nlev),stat=ierr)
135 247680 : if (ierr/=0) then
136 0 : nullify(newobj)
137 0 : return
138 : end if
139 :
140 743040 : allocate(newobj%crefwlw(nlw),stat=ierr)
141 247680 : if (ierr/=0) then
142 0 : nullify(newobj)
143 0 : return
144 : end if
145 4210560 : newobj%crefwlw(:) = crefwlw(:)
146 :
147 743040 : allocate(newobj%crefwsw(nsw),stat=ierr)
148 247680 : if (ierr/=0) then
149 0 : nullify(newobj)
150 0 : return
151 : end if
152 3715200 : newobj%crefwsw(:) = crefwsw(:)
153 :
154 247680 : call aero_state%water_uptake(aero_props, ilist, ibin, ncol, nlev, dgnumwet, qaerwat)
155 :
156 247680 : nspec = aero_props%nspecies(ilist,ibin)
157 :
158 247680 : logsigma=aero_props%alogsig(ilist,ibin)
159 :
160 : ! calc size parameter for all columns
161 : call modal_size_parameters(newobj%ncoef, ncol, nlev, logsigma, dgnumwet, &
162 247680 : newobj%radsurf, newobj%logradsurf, newobj%cheb)
163 :
164 23281920 : do ilev = 1, nlev
165 384618240 : dryvol(:ncol) = 0._r8
166 109660320 : do ispec = 1, nspec
167 86378400 : call aero_state%get_ambient_mmr(ilist,ispec,ibin,specmmr)
168 86378400 : call aero_props%get(ibin, ispec, list_ndx=ilist, density=specdens)
169 :
170 1465352640 : do icol = 1, ncol
171 1355940000 : vol(icol) = specmmr(icol,ilev)/specdens
172 1355940000 : dryvol(icol) = dryvol(icol) + vol(icol)
173 :
174 1355940000 : newobj%watervol(icol,ilev) = qaerwat(icol,ilev)*rh2odens
175 1355940000 : newobj%wetvol(icol,ilev) = newobj%watervol(icol,ilev) + dryvol(icol)
176 1442318400 : if (newobj%watervol(icol,ilev) < 0._r8) then
177 0 : newobj%watervol(icol,ilev) = 0._r8
178 0 : newobj%wetvol(icol,ilev) = dryvol(icol)
179 : end if
180 : end do
181 : end do
182 : end do
183 :
184 247680 : newobj%aero_state => aero_state
185 247680 : newobj%aero_props => aero_props
186 247680 : newobj%ilist = ilist
187 247680 : newobj%ibin = ibin
188 :
189 247680 : end function constructor
190 :
191 : !------------------------------------------------------------------------------
192 : ! returns short wave aerosol optics properties
193 : !------------------------------------------------------------------------------
194 161239680 : subroutine sw_props(self, ncol, ilev, iwav, pext, pabs, palb, pasm)
195 :
196 : class(refractive_aerosol_optics), intent(in) :: self
197 : integer, intent(in) :: ncol ! number of columns
198 : integer, intent(in) :: ilev ! vertical level index
199 : integer, intent(in) :: iwav ! wave length index
200 : real(r8),intent(out) :: pext(ncol) ! parameterized specific extinction (m2/kg)
201 : real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg)
202 : real(r8),intent(out) :: palb(ncol) ! parameterized asymmetry factor
203 : real(r8),intent(out) :: pasm(ncol) ! parameterized single scattering albedo
204 :
205 322479360 : real(r8) :: refr(ncol) ! real part of refractive index
206 322479360 : real(r8) :: refi(ncol) ! imaginary part of refractive index
207 322479360 : real(r8) :: cext(self%ncoef,ncol), cabs(self%ncoef,ncol), casm(self%ncoef,ncol)
208 :
209 322479360 : complex(r8) :: crefin(ncol) ! complex refractive index
210 : integer :: icol,icoef
211 :
212 322479360 : type(table_interp_wghts) :: wghtsr(ncol)
213 322479360 : type(table_interp_wghts) :: wghtsi(ncol)
214 :
215 2692327680 : crefin(:ncol) = self%aero_state%refractive_index_sw(ncol, ilev, self%ilist, self%ibin, iwav, self%aero_props)
216 :
217 2692327680 : do icol = 1, ncol
218 2531088000 : crefin(icol) = crefin(icol) + self%watervol(icol,ilev)*self%crefwsw(iwav)
219 2531088000 : crefin(icol) = crefin(icol)/max(self%wetvol(icol,ilev),1.e-60_r8)
220 2531088000 : refr(icol) = real(crefin(icol))
221 2692327680 : refi(icol) = abs(aimag(crefin(icol)))
222 : end do
223 :
224 : ! interpolate coefficients linear in refractive index
225 :
226 161239680 : wghtsr = table_interp_calcwghts( self%prefr, self%refrtabsw(:,iwav), ncol, refr(:ncol) )
227 161239680 : wghtsi = table_interp_calcwghts( self%prefi, self%refitabsw(:,iwav), ncol, refi(:ncol) )
228 :
229 161239680 : cext(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%extpsw(:,:,:,iwav))
230 161239680 : cabs(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%abspsw(:,:,:,iwav))
231 161239680 : casm(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%asmpsw(:,:,:,iwav))
232 :
233 2692327680 : do icol = 1,ncol
234 :
235 2531088000 : if (self%logradsurf(icol,ilev) <= xrmax) then
236 2531088000 : pext(icol) = 0.5_r8*cext(1,icol)
237 12655440000 : do icoef = 2, self%ncoef
238 12655440000 : pext(icol) = pext(icol) + self%cheb(icoef,icol,ilev)*cext(icoef,icol)
239 : enddo
240 2531088000 : pext(icol) = exp(pext(icol))
241 : else
242 0 : pext(icol) = 1.5_r8/(self%radsurf(icol,ilev)*rhoh2o) ! geometric optics
243 : endif
244 :
245 : ! convert from m2/kg water to m2/kg aerosol
246 2531088000 : pext(icol) = pext(icol)*self%wetvol(icol,ilev)*rhoh2o
247 2531088000 : pabs(icol) = 0.5_r8*cabs(1,icol)
248 2531088000 : pasm(icol) = 0.5_r8*casm(1,icol)
249 12655440000 : do icoef = 2, self%ncoef
250 10124352000 : pabs(icol) = pabs(icol) + self%cheb(icoef,icol,ilev)*cabs(icoef,icol)
251 12655440000 : pasm(icol) = pasm(icol) + self%cheb(icoef,icol,ilev)*casm(icoef,icol)
252 : enddo
253 2531088000 : pabs(icol) = pabs(icol)*self%wetvol(icol,ilev)*rhoh2o
254 2531088000 : pabs(icol) = max(0._r8,pabs(icol))
255 2531088000 : pabs(icol) = min(pext(icol),pabs(icol))
256 :
257 2692327680 : palb(icol) = 1._r8-pabs(icol)/max(pext(icol),1.e-40_r8)
258 :
259 : end do
260 :
261 161239680 : end subroutine sw_props
262 :
263 : !------------------------------------------------------------------------------
264 : ! returns long wave aerosol optics properties
265 : !------------------------------------------------------------------------------
266 184273920 : subroutine lw_props(self, ncol, ilev, iwav, pabs)
267 :
268 : class(refractive_aerosol_optics), intent(in) :: self
269 : integer, intent(in) :: ncol ! number of columns
270 : integer, intent(in) :: ilev ! vertical level index
271 : integer, intent(in) :: iwav ! wave length index
272 : real(r8),intent(out) :: pabs(ncol) ! parameterized specific absorption (m2/kg)
273 :
274 368547840 : real(r8) :: refr(ncol) ! real part of refractive index
275 368547840 : real(r8) :: refi(ncol) ! imaginary part of refractive index
276 368547840 : real(r8) :: cabs(self%ncoef,ncol)
277 :
278 368547840 : complex(r8) :: crefin(ncol) ! complex refractive index
279 : integer :: icol, icoef
280 :
281 368547840 : type(table_interp_wghts) :: wghtsr(ncol)
282 368547840 : type(table_interp_wghts) :: wghtsi(ncol)
283 :
284 3076945920 : crefin(:ncol) = self%aero_state%refractive_index_lw(ncol, ilev, self%ilist, self%ibin, iwav, self%aero_props)
285 :
286 3076945920 : do icol = 1, ncol
287 2892672000 : crefin(icol) = crefin(icol) + self%watervol(icol,ilev)*self%crefwlw(iwav)
288 2892672000 : crefin(icol) = crefin(icol)/max(self%wetvol(icol,ilev), 1.e-40_r8)
289 :
290 2892672000 : refr(icol) = real(crefin(icol))
291 3076945920 : refi(icol) = aimag(crefin(icol))
292 :
293 : end do
294 :
295 : ! interpolate coefficients linear in refractive index
296 :
297 184273920 : wghtsr = table_interp_calcwghts( self%prefr, self%refrtablw(:,iwav), ncol, refr(:ncol) )
298 184273920 : wghtsi = table_interp_calcwghts( self%prefi, self%refitablw(:,iwav), ncol, refi(:ncol) )
299 :
300 184273920 : cabs(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%absplw(:,:,:,iwav))
301 :
302 3076945920 : do icol = 1,ncol
303 2892672000 : pabs(icol) = 0.5_r8*cabs(1,icol)
304 14463360000 : do icoef = 2, self%ncoef
305 14463360000 : pabs(icol) = pabs(icol) + self%cheb(icoef,icol,ilev)*cabs(icoef,icol)
306 : end do
307 2892672000 : pabs(icol) = pabs(icol)*self%wetvol(icol,ilev)*rhoh2o
308 3076945920 : pabs(icol) = max(0._r8,pabs(icol))
309 : end do
310 :
311 184273920 : end subroutine lw_props
312 :
313 :
314 : !------------------------------------------------------------------------------
315 : !------------------------------------------------------------------------------
316 247680 : subroutine destructor(self)
317 :
318 : type(refractive_aerosol_optics), intent(inout) :: self
319 :
320 247680 : deallocate(self%watervol)
321 247680 : deallocate(self%wetvol)
322 247680 : deallocate(self%cheb)
323 247680 : deallocate(self%radsurf)
324 247680 : deallocate(self%logradsurf)
325 247680 : deallocate(self%crefwsw)
326 247680 : deallocate(self%crefwlw)
327 :
328 247680 : nullify(self%aero_state)
329 247680 : nullify(self%aero_props)
330 247680 : nullify(self%extpsw)
331 247680 : nullify(self%abspsw)
332 247680 : nullify(self%asmpsw)
333 247680 : nullify(self%absplw)
334 247680 : nullify(self%refrtabsw)
335 247680 : nullify(self%refitabsw)
336 247680 : nullify(self%refrtablw)
337 247680 : nullify(self%refitablw)
338 :
339 247680 : end subroutine destructor
340 :
341 :
342 : ! Private routines
343 : !===============================================================================
344 :
345 : !===============================================================================
346 :
347 247680 : subroutine modal_size_parameters(ncoef,ncol,nlev, alnsg_amode, dgnumwet, radsurf, logradsurf, cheb)
348 :
349 : integer, intent(in) :: ncoef,ncol,nlev
350 : real(r8), intent(in) :: alnsg_amode ! geometric standard deviation of number distribution
351 : real(r8), intent(in) :: dgnumwet(:,:) ! aerosol wet number mode diameter (m)
352 : real(r8), intent(out) :: radsurf(:,:) ! aerosol surface mode radius
353 : real(r8), intent(out) :: logradsurf(:,:) ! log(aerosol surface mode radius)
354 : real(r8), intent(out) :: cheb(:,:,:)
355 :
356 : integer :: i, k, nc
357 : real(r8) :: explnsigma
358 495360 : real(r8) :: xrad(ncol) ! normalized aerosol radius
359 :
360 : !-------------------------------------------------------------------------------
361 :
362 247680 : explnsigma = exp(2.0_r8*alnsg_amode*alnsg_amode)
363 :
364 23281920 : do k = 1, nlev
365 384865920 : do i = 1, ncol
366 : ! convert from number mode diameter to surface area
367 361584000 : radsurf(i,k) = max(0.5_r8*dgnumwet(i,k)*explnsigma,radmin)
368 361584000 : logradsurf(i,k) = log(radsurf(i,k))
369 : ! normalize size parameter
370 361584000 : xrad(i) = max(logradsurf(i,k),xrmin)
371 361584000 : xrad(i) = min(xrad(i),xrmax)
372 361584000 : xrad(i) = (2._r8*xrad(i)-xrmax-xrmin)/(xrmax-xrmin)
373 : ! chebyshev polynomials
374 361584000 : cheb(1,i,k) = 1._r8
375 361584000 : cheb(2,i,k) = xrad(i)
376 1469370240 : do nc = 3, ncoef
377 1446336000 : cheb(nc,i,k) = 2._r8*xrad(i)*cheb(nc-1,i,k)-cheb(nc-2,i,k)
378 : end do
379 : end do
380 : end do
381 :
382 247680 : end subroutine modal_size_parameters
383 :
384 990720 : end module refractive_aerosol_optics_mod
|
