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 5969088 : 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 11938176 : real(r8) :: vol(ncol) ! volume concentration of aerosol species (m3/kg)
91 11938176 : real(r8) :: dryvol(ncol) ! volume concentration of aerosol mode (m3/kg)
92 : real(r8) :: specdens ! species density (kg/m3)
93 5969088 : real(r8), pointer :: specmmr(:,:) ! species mass mixing ratio
94 : real(r8) :: logsigma ! geometric standard deviation of number distribution
95 :
96 11938176 : real(r8) :: dgnumwet(ncol,nlev) ! aerosol wet number mode diameter (m)
97 11938176 : real(r8) :: qaerwat(ncol,nlev) ! aerosol water (g/g)
98 :
99 : real(r8), parameter :: rh2odens = 1._r8/rhoh2o
100 :
101 5969088 : allocate(newobj, stat=ierr)
102 5969088 : if (ierr/=0) then
103 5969088 : 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 5969088 : absplw=newobj%absplw, ncoef=newobj%ncoef, prefr=newobj%prefr, prefi=newobj%prefi)
113 :
114 23876352 : allocate(newobj%watervol(ncol,nlev),stat=ierr)
115 5969088 : if (ierr/=0) then
116 0 : nullify(newobj)
117 0 : return
118 : end if
119 17907264 : allocate(newobj%wetvol(ncol,nlev),stat=ierr)
120 5969088 : if (ierr/=0) then
121 0 : nullify(newobj)
122 0 : return
123 : end if
124 29845440 : allocate(newobj%cheb(newobj%ncoef,ncol,nlev),stat=ierr)
125 5969088 : if (ierr/=0) then
126 0 : nullify(newobj)
127 0 : return
128 : end if
129 17907264 : allocate(newobj%radsurf(ncol,nlev),stat=ierr)
130 5969088 : if (ierr/=0) then
131 0 : nullify(newobj)
132 0 : return
133 : end if
134 17907264 : allocate(newobj%logradsurf(ncol,nlev),stat=ierr)
135 5969088 : if (ierr/=0) then
136 0 : nullify(newobj)
137 0 : return
138 : end if
139 :
140 17907264 : allocate(newobj%crefwlw(nlw),stat=ierr)
141 5969088 : if (ierr/=0) then
142 0 : nullify(newobj)
143 0 : return
144 : end if
145 101474496 : newobj%crefwlw(:) = crefwlw(:)
146 :
147 17907264 : allocate(newobj%crefwsw(nsw),stat=ierr)
148 5969088 : if (ierr/=0) then
149 0 : nullify(newobj)
150 0 : return
151 : end if
152 89536320 : newobj%crefwsw(:) = crefwsw(:)
153 :
154 5969088 : call aero_state%water_uptake(aero_props, ilist, ibin, ncol, nlev, dgnumwet, qaerwat)
155 :
156 5969088 : nspec = aero_props%nspecies(ilist,ibin)
157 :
158 5969088 : 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 5969088 : newobj%radsurf, newobj%logradsurf, newobj%cheb)
163 :
164 561094272 : do ilev = 1, nlev
165 9269299584 : dryvol(:ncol) = 0._r8
166 2642813712 : do ispec = 1, nspec
167 2081719440 : call aero_state%get_ambient_mmr(ilist,ispec,ibin,specmmr)
168 2081719440 : call aero_props%get(ibin, ispec, list_ndx=ilist, density=specdens)
169 :
170 35314998624 : do icol = 1, ncol
171 32678154000 : vol(icol) = specmmr(icol,ilev)/specdens
172 32678154000 : dryvol(icol) = dryvol(icol) + vol(icol)
173 :
174 32678154000 : newobj%watervol(icol,ilev) = qaerwat(icol,ilev)*rh2odens
175 32678154000 : newobj%wetvol(icol,ilev) = newobj%watervol(icol,ilev) + dryvol(icol)
176 34759873440 : 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 5969088 : newobj%aero_state => aero_state
185 5969088 : newobj%aero_props => aero_props
186 5969088 : newobj%ilist = ilist
187 5969088 : newobj%ibin = ibin
188 :
189 5969088 : end function constructor
190 :
191 : !------------------------------------------------------------------------------
192 : ! returns short wave aerosol optics properties
193 : !------------------------------------------------------------------------------
194 3885876288 : 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 7771752576 : real(r8) :: refr(ncol) ! real part of refractive index
206 7771752576 : real(r8) :: refi(ncol) ! imaginary part of refractive index
207 7771752576 : real(r8) :: cext(self%ncoef,ncol), cabs(self%ncoef,ncol), casm(self%ncoef,ncol)
208 :
209 7771752576 : complex(r8) :: crefin(ncol) ! complex refractive index
210 : integer :: icol,icoef
211 :
212 7771752576 : type(table_interp_wghts) :: wghtsr(ncol)
213 7771752576 : type(table_interp_wghts) :: wghtsi(ncol)
214 :
215 64885097088 : crefin(:ncol) = self%aero_state%refractive_index_sw(ncol, ilev, self%ilist, self%ibin, iwav, self%aero_props)
216 :
217 64885097088 : do icol = 1, ncol
218 60999220800 : crefin(icol) = crefin(icol) + self%watervol(icol,ilev)*self%crefwsw(iwav)
219 60999220800 : crefin(icol) = crefin(icol)/max(self%wetvol(icol,ilev),1.e-60_r8)
220 60999220800 : refr(icol) = real(crefin(icol))
221 64885097088 : refi(icol) = abs(aimag(crefin(icol)))
222 : end do
223 :
224 : ! interpolate coefficients linear in refractive index
225 :
226 3885876288 : wghtsr = table_interp_calcwghts( self%prefr, self%refrtabsw(:,iwav), ncol, refr(:ncol) )
227 3885876288 : wghtsi = table_interp_calcwghts( self%prefi, self%refitabsw(:,iwav), ncol, refi(:ncol) )
228 :
229 3885876288 : cext(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%extpsw(:,:,:,iwav))
230 3885876288 : cabs(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%abspsw(:,:,:,iwav))
231 3885876288 : casm(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%asmpsw(:,:,:,iwav))
232 :
233 64885097088 : do icol = 1,ncol
234 :
235 60999220800 : if (self%logradsurf(icol,ilev) <= xrmax) then
236 60999220800 : pext(icol) = 0.5_r8*cext(1,icol)
237 >30499*10^7 : do icoef = 2, self%ncoef
238 >30499*10^7 : pext(icol) = pext(icol) + self%cheb(icoef,icol,ilev)*cext(icoef,icol)
239 : enddo
240 60999220800 : 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 60999220800 : pext(icol) = pext(icol)*self%wetvol(icol,ilev)*rhoh2o
247 60999220800 : pabs(icol) = 0.5_r8*cabs(1,icol)
248 60999220800 : pasm(icol) = 0.5_r8*casm(1,icol)
249 >30499*10^7 : do icoef = 2, self%ncoef
250 >24399*10^7 : pabs(icol) = pabs(icol) + self%cheb(icoef,icol,ilev)*cabs(icoef,icol)
251 >30499*10^7 : pasm(icol) = pasm(icol) + self%cheb(icoef,icol,ilev)*casm(icoef,icol)
252 : enddo
253 60999220800 : pabs(icol) = pabs(icol)*self%wetvol(icol,ilev)*rhoh2o
254 60999220800 : pabs(icol) = max(0._r8,pabs(icol))
255 60999220800 : pabs(icol) = min(pext(icol),pabs(icol))
256 :
257 64885097088 : palb(icol) = 1._r8-pabs(icol)/max(pext(icol),1.e-40_r8)
258 :
259 : end do
260 :
261 3885876288 : end subroutine sw_props
262 :
263 : !------------------------------------------------------------------------------
264 : ! returns long wave aerosol optics properties
265 : !------------------------------------------------------------------------------
266 4441001472 : 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 8882002944 : real(r8) :: refr(ncol) ! real part of refractive index
275 8882002944 : real(r8) :: refi(ncol) ! imaginary part of refractive index
276 8882002944 : real(r8) :: cabs(self%ncoef,ncol)
277 :
278 8882002944 : complex(r8) :: crefin(ncol) ! complex refractive index
279 : integer :: icol, icoef
280 :
281 8882002944 : type(table_interp_wghts) :: wghtsr(ncol)
282 8882002944 : type(table_interp_wghts) :: wghtsi(ncol)
283 :
284 74154396672 : crefin(:ncol) = self%aero_state%refractive_index_lw(ncol, ilev, self%ilist, self%ibin, iwav, self%aero_props)
285 :
286 74154396672 : do icol = 1, ncol
287 69713395200 : crefin(icol) = crefin(icol) + self%watervol(icol,ilev)*self%crefwlw(iwav)
288 69713395200 : crefin(icol) = crefin(icol)/max(self%wetvol(icol,ilev), 1.e-40_r8)
289 :
290 69713395200 : refr(icol) = real(crefin(icol))
291 74154396672 : refi(icol) = aimag(crefin(icol))
292 :
293 : end do
294 :
295 : ! interpolate coefficients linear in refractive index
296 :
297 4441001472 : wghtsr = table_interp_calcwghts( self%prefr, self%refrtablw(:,iwav), ncol, refr(:ncol) )
298 4441001472 : wghtsi = table_interp_calcwghts( self%prefi, self%refitablw(:,iwav), ncol, refi(:ncol) )
299 :
300 4441001472 : cabs(:,:ncol)= table_interp( self%ncoef,ncol, self%prefr,self%prefi, wghtsr,wghtsi, self%absplw(:,:,:,iwav))
301 :
302 74154396672 : do icol = 1,ncol
303 69713395200 : pabs(icol) = 0.5_r8*cabs(1,icol)
304 >34856*10^7 : do icoef = 2, self%ncoef
305 >34856*10^7 : pabs(icol) = pabs(icol) + self%cheb(icoef,icol,ilev)*cabs(icoef,icol)
306 : end do
307 69713395200 : pabs(icol) = pabs(icol)*self%wetvol(icol,ilev)*rhoh2o
308 74154396672 : pabs(icol) = max(0._r8,pabs(icol))
309 : end do
310 :
311 4441001472 : end subroutine lw_props
312 :
313 :
314 : !------------------------------------------------------------------------------
315 : !------------------------------------------------------------------------------
316 5969088 : subroutine destructor(self)
317 :
318 : type(refractive_aerosol_optics), intent(inout) :: self
319 :
320 5969088 : deallocate(self%watervol)
321 5969088 : deallocate(self%wetvol)
322 5969088 : deallocate(self%cheb)
323 5969088 : deallocate(self%radsurf)
324 5969088 : deallocate(self%logradsurf)
325 5969088 : deallocate(self%crefwsw)
326 5969088 : deallocate(self%crefwlw)
327 :
328 5969088 : nullify(self%aero_state)
329 5969088 : nullify(self%aero_props)
330 5969088 : nullify(self%extpsw)
331 5969088 : nullify(self%abspsw)
332 5969088 : nullify(self%asmpsw)
333 5969088 : nullify(self%absplw)
334 5969088 : nullify(self%refrtabsw)
335 5969088 : nullify(self%refitabsw)
336 5969088 : nullify(self%refrtablw)
337 5969088 : nullify(self%refitablw)
338 :
339 5969088 : end subroutine destructor
340 :
341 :
342 : ! Private routines
343 : !===============================================================================
344 :
345 : !===============================================================================
346 :
347 5969088 : 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 11938176 : real(r8) :: xrad(ncol) ! normalized aerosol radius
359 :
360 : !-------------------------------------------------------------------------------
361 :
362 5969088 : explnsigma = exp(2.0_r8*alnsg_amode*alnsg_amode)
363 :
364 561094272 : do k = 1, nlev
365 9275268672 : do i = 1, ncol
366 : ! convert from number mode diameter to surface area
367 8714174400 : radsurf(i,k) = max(0.5_r8*dgnumwet(i,k)*explnsigma,radmin)
368 8714174400 : logradsurf(i,k) = log(radsurf(i,k))
369 : ! normalize size parameter
370 8714174400 : xrad(i) = max(logradsurf(i,k),xrmin)
371 8714174400 : xrad(i) = min(xrad(i),xrmax)
372 8714174400 : xrad(i) = (2._r8*xrad(i)-xrmax-xrmin)/(xrmax-xrmin)
373 : ! chebyshev polynomials
374 8714174400 : cheb(1,i,k) = 1._r8
375 8714174400 : cheb(2,i,k) = xrad(i)
376 35411822784 : do nc = 3, ncoef
377 34856697600 : 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 5969088 : end subroutine modal_size_parameters
383 :
384 23876352 : end module refractive_aerosol_optics_mod
|
