Line data Source code
1 : module ndrop_bam
2 :
3 : !---------------------------------------------------------------------------------
4 : !
5 : ! CAM Interface for droplet activation by bulk aerosols.
6 : !
7 : !---------------------------------------------------------------------------------
8 :
9 : use shr_kind_mod, only: r8=>shr_kind_r8
10 : use spmd_utils, only: masterproc
11 : use ppgrid, only: pcols, pver, pverp
12 : use physconst, only: gravit, rair, tmelt, cpair, rh2o, &
13 : r_universal, mwh2o, rhoh2o, latvap
14 :
15 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_props
16 :
17 : use shr_spfn_mod, only: erf => shr_spfn_erf, &
18 : erfc => shr_spfn_erfc
19 : use wv_saturation, only: qsat
20 : use cam_history, only: addfld, add_default, outfld
21 : use cam_logfile, only: iulog
22 : use cam_abortutils, only: endrun
23 : use ref_pres, only: top_lev=>trop_cloud_top_lev
24 :
25 : implicit none
26 : private
27 : save
28 :
29 : public :: ndrop_bam_init, ndrop_bam_run, ndrop_bam_ccn
30 :
31 : ! activate parameters
32 :
33 : real(r8) :: pi ! pi
34 : integer, parameter :: psat = 6 ! number of supersaturations to calc ccn concentration
35 : real(r8), parameter :: supersat(psat) = & ! supersaturation (%) to determine ccn concentration
36 : (/0.02_r8,0.05_r8,0.1_r8,0.2_r8,0.5_r8,1.0_r8/)
37 : real(r8) :: super(psat)
38 : real(r8), allocatable :: ccnfact(:,:)
39 :
40 : real(r8), allocatable :: alogsig(:) ! natl log of geometric standard dev of aerosol
41 : real(r8), allocatable :: exp45logsig(:)
42 : real(r8), allocatable :: argfactor(:)
43 : real(r8), allocatable :: amcube(:) ! cube of dry mode radius (m)
44 : real(r8), allocatable :: smcrit(:) ! critical supersatuation for activation
45 : real(r8), allocatable :: lnsm(:) ! ln(smcrit)
46 : real(r8), allocatable :: amcubefactor(:) ! factors for calculating mode radius
47 : real(r8), allocatable :: smcritfactor(:) ! factors for calculating critical supersaturation
48 : real(r8), allocatable :: f1(:), f2(:) ! abdul-razzak functions of width
49 :
50 : real(r8) :: aten
51 : real(r8) :: third, sixth
52 : real(r8) :: sq2, sqpi
53 : real(r8) :: alogten, alog2, alog3, alogaten
54 :
55 : real(r8) :: pref = 1013.25e2_r8 ! reference pressure (Pa)
56 :
57 : ! aerosol properties
58 : character(len=20), allocatable :: aername(:)
59 : real(r8), allocatable :: dryrad_aer(:)
60 : real(r8), allocatable :: density_aer(:)
61 : real(r8), allocatable :: hygro_aer(:)
62 : real(r8), allocatable :: dispersion_aer(:)
63 : real(r8), allocatable :: num_to_mass_aer(:)
64 :
65 : integer :: naer_all ! number of aerosols affecting climate
66 : integer :: idxsul = -1 ! index in aerosol list for sulfate
67 :
68 : !===============================================================================
69 : contains
70 : !===============================================================================
71 :
72 0 : subroutine ndrop_bam_init
73 :
74 : use phys_control, only: phys_getopts
75 :
76 : !-----------------------------------------------------------------------
77 : !
78 : ! Initialize constants for droplet activation by bulk aerosols
79 : !
80 : !-----------------------------------------------------------------------
81 :
82 : integer :: l, m, iaer
83 : real(r8) :: surften ! surface tension of water w/respect to air (N/m)
84 : real(r8) :: arg
85 : logical :: history_amwg
86 : !-------------------------------------------------------------------------------
87 :
88 0 : call phys_getopts(history_amwg_out=history_amwg)
89 :
90 : ! Access the physical properties of the bulk aerosols that are affecting the climate
91 : ! by using routines from the rad_constituents module.
92 :
93 0 : call rad_cnst_get_info(0, naero=naer_all)
94 : allocate( &
95 0 : aername(naer_all), &
96 0 : dryrad_aer(naer_all), &
97 0 : density_aer(naer_all), &
98 0 : hygro_aer(naer_all), &
99 0 : dispersion_aer(naer_all), &
100 0 : num_to_mass_aer(naer_all) )
101 :
102 0 : do iaer = 1, naer_all
103 : call rad_cnst_get_aer_props(0, iaer, &
104 0 : aername = aername(iaer), &
105 0 : dryrad_aer = dryrad_aer(iaer), &
106 0 : density_aer = density_aer(iaer), &
107 0 : hygro_aer = hygro_aer(iaer), &
108 0 : dispersion_aer = dispersion_aer(iaer), &
109 0 : num_to_mass_aer = num_to_mass_aer(iaer) )
110 :
111 : ! Look for sulfate aerosol in this list (Bulk aerosol only)
112 0 : if (trim(aername(iaer)) == 'SULFATE') idxsul = iaer
113 :
114 : ! aerosol number concentration
115 0 : call addfld(trim(aername(iaer))//'_m3', (/ 'lev' /), 'A', 'm-3', 'aerosol number concentration')
116 :
117 : end do
118 :
119 0 : if (masterproc) then
120 0 : write(iulog,*) 'ndrop_bam_init: iaer, name, dryrad, density, hygro, dispersion, num_to_mass'
121 0 : do iaer = 1, naer_all
122 0 : write(iulog,*) iaer, aername(iaer), dryrad_aer(iaer), density_aer(iaer), hygro_aer(iaer), &
123 0 : dispersion_aer(iaer), num_to_mass_aer(iaer)
124 : end do
125 0 : if (idxsul < 1) then
126 0 : write(iulog,*) 'ndrop_bam_init: SULFATE aerosol properties NOT FOUND'
127 : else
128 0 : write(iulog,*) 'ndrop_bam_init: SULFATE aerosol properties FOUND at index ',idxsul
129 : end if
130 : end if
131 :
132 0 : call addfld ('CCN1',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.02%')
133 0 : call addfld ('CCN2',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.05%')
134 0 : call addfld ('CCN3',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.1%')
135 0 : call addfld ('CCN4',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.2%')
136 0 : call addfld ('CCN5',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.5%')
137 0 : call addfld ('CCN6',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=1.0%')
138 :
139 0 : if (history_amwg) then
140 0 : call add_default('CCN3', 1, ' ')
141 : end if
142 :
143 : ! set parameters for droplet activation, following abdul-razzak and ghan 2000, JGR
144 :
145 0 : third = 1._r8/3._r8
146 0 : sixth = 1._r8/6._r8
147 0 : sq2 = sqrt(2._r8)
148 0 : pi = 4._r8*atan(1.0_r8)
149 0 : sqpi = sqrt(pi)
150 0 : surften = 0.076_r8
151 0 : aten = 2._r8*mwh2o*surften/(r_universal*tmelt*rhoh2o)
152 0 : alogaten= log(aten)
153 0 : alog2 = log(2._r8)
154 0 : alog3 = log(3._r8)
155 0 : super(:)= 0.01_r8*supersat(:)
156 :
157 : allocate( &
158 0 : alogsig(naer_all), &
159 0 : exp45logsig(naer_all), &
160 0 : argfactor(naer_all), &
161 0 : f1(naer_all), &
162 0 : f2(naer_all), &
163 0 : amcubefactor(naer_all), &
164 0 : smcritfactor(naer_all), &
165 0 : amcube(naer_all), &
166 0 : smcrit(naer_all), &
167 0 : lnsm(naer_all), &
168 0 : ccnfact(psat,naer_all) )
169 :
170 :
171 0 : do m = 1, naer_all
172 :
173 : ! Skip aerosols that don't have a dispersion defined.
174 0 : if (dispersion_aer(m) == 0._r8) cycle
175 :
176 0 : alogsig(m) = log(dispersion_aer(m))
177 0 : exp45logsig(m) = exp(4.5_r8*alogsig(m)*alogsig(m))
178 0 : argfactor(m) = 2._r8/(3._r8*sqrt(2._r8)*alogsig(m))
179 0 : f1(m) = 0.5_r8*exp(2.5_r8*alogsig(m)*alogsig(m))
180 0 : f2(m) = 1._r8 + 0.25_r8*alogsig(m)
181 0 : amcubefactor(m)= 3._r8/(4._r8*pi*exp45logsig(m)*density_aer(m))
182 0 : smcritfactor(m)= 2._r8*aten*sqrt(aten/(27._r8*max(1.e-10_r8,hygro_aer(m))))
183 0 : amcube(m) = amcubefactor(m)/num_to_mass_aer(m)
184 :
185 0 : if (hygro_aer(m) .gt. 1.e-10_r8) then
186 0 : smcrit(m) = smcritfactor(m)/sqrt(amcube(m))
187 : else
188 0 : smcrit(m) = 100._r8
189 : endif
190 0 : lnsm(m) = log(smcrit(m))
191 :
192 0 : do l = 1, psat
193 0 : arg = argfactor(m)*log(smcrit(m)/super(l))
194 0 : if (arg < 2) then
195 0 : if (arg < -2) then
196 0 : ccnfact(l,m) = 1.e-6_r8
197 : else
198 0 : ccnfact(l,m) = 1.e-6_r8*0.5_r8*erfc(arg)
199 : endif
200 : else
201 0 : ccnfact(l,m) = 0._r8
202 : endif
203 : enddo
204 :
205 : end do
206 :
207 0 : end subroutine ndrop_bam_init
208 :
209 : !===============================================================================
210 :
211 0 : subroutine ndrop_bam_run( &
212 0 : wbar, tair, rhoair, na, pmode, &
213 0 : nmode, ma, nact)
214 :
215 : ! calculates number fraction of aerosols activated as CCN
216 : ! assumes an internal mixture within each of up to pmode multiple aerosol modes
217 : ! a gaussian spectrum of updrafts can be treated.
218 :
219 : ! mks units
220 :
221 : ! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
222 : ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
223 :
224 : ! input
225 : integer, intent(in) :: pmode ! dimension of modes
226 : integer, intent(in) :: nmode ! number of aerosol modes
227 : real(r8), intent(in) :: wbar ! grid cell mean vertical velocity (m/s)
228 : real(r8), intent(in) :: tair ! air temperature (K)
229 : real(r8), intent(in) :: rhoair ! air density (kg/m3)
230 : real(r8), intent(in) :: na(pmode) ! aerosol number concentration (1/m3)
231 : real(r8), intent(in) :: ma(pmode) ! aerosol mass concentration (kg/m3)
232 :
233 : ! output
234 : real(r8), intent(out) :: nact ! number fraction of aerosols activated
235 :
236 : ! local variables
237 : integer :: maxmodes
238 :
239 0 : real(r8), allocatable :: volc(:) ! total aerosol volume concentration (m3/m3)
240 0 : real(r8), allocatable :: eta(:)
241 0 : real(r8), allocatable :: smc(:)
242 0 : real(r8), allocatable :: etafactor2(:)
243 0 : real(r8), allocatable :: amcubeloc(:)
244 0 : real(r8), allocatable :: lnsmloc(:)
245 :
246 : real(r8) :: pres ! pressure (Pa)
247 : real(r8) :: diff0
248 : real(r8) :: conduct0 ! thermal conductivity (Joule/m/sec/deg)
249 : real(r8) :: qs ! water vapor saturation mixing ratio
250 : real(r8) :: dqsdt ! change in qs with temperature
251 : real(r8) :: gloc ! thermodynamic function (m2/s)
252 : real(r8) :: zeta
253 : real(r8) :: lnsmax ! ln(smax)
254 : real(r8) :: alpha
255 : real(r8) :: gammaloc
256 : real(r8) :: beta
257 : real(r8) :: sqrtg
258 : real(r8) :: wnuc
259 : real(r8) :: alw
260 : real(r8) :: sqrtalw
261 : real(r8) :: smax
262 : real(r8) :: x
263 : real(r8) :: etafactor1
264 : real(r8) :: etafactor2max
265 : real(r8) :: es
266 : integer :: m
267 : !-------------------------------------------------------------------------------
268 :
269 0 : maxmodes = naer_all
270 : allocate( &
271 : volc(maxmodes), &
272 : eta(maxmodes), &
273 : smc(maxmodes), &
274 : etafactor2(maxmodes), &
275 : amcubeloc(maxmodes), &
276 0 : lnsmloc(maxmodes) )
277 :
278 0 : if (maxmodes < pmode) then
279 0 : write(iulog,*)'ndrop_bam_run: maxmodes,pmode=',maxmodes,pmode
280 0 : call endrun('ndrop_bam_run')
281 : endif
282 :
283 0 : nact = 0._r8
284 :
285 0 : if (nmode .eq. 1 .and. na(1) .lt. 1.e-20_r8) return
286 :
287 0 : if (wbar .le. 0._r8) return
288 :
289 0 : pres = rair*rhoair*tair
290 0 : diff0 = 0.211e-4_r8*(pref/pres)*(tair/tmelt)**1.94_r8
291 0 : conduct0 = (5.69_r8 + 0.017_r8*(tair-tmelt))*4.186e2_r8*1.e-5_r8 ! convert to J/m/s/deg
292 0 : call qsat(tair, pres, es, qs)
293 0 : dqsdt = latvap/(rh2o*tair*tair)*qs
294 0 : alpha = gravit*(latvap/(cpair*rh2o*tair*tair) - 1._r8/(rair*tair))
295 0 : gammaloc = (1 + latvap/cpair*dqsdt)/(rhoair*qs)
296 : ! growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
297 : ! should depend on mean radius of mode to account for gas kinetic effects
298 : gloc = 1._r8/(rhoh2o/(diff0*rhoair*qs) &
299 0 : + latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair) - 1._r8))
300 0 : sqrtg = sqrt(gloc)
301 0 : beta = 4._r8*pi*rhoh2o*gloc*gammaloc
302 0 : etafactor2max = 1.e10_r8/(alpha*wbar)**1.5_r8 ! this should make eta big if na is very small.
303 :
304 0 : do m = 1, nmode
305 : ! skip aerosols with no dispersion, since they aren't meant to be CCN
306 0 : if (dispersion_aer(m) == 0._r8) then
307 0 : smc(m)=100._r8
308 0 : cycle
309 : end if
310 : ! internal mixture of aerosols
311 0 : volc(m) = ma(m)/(density_aer(m)) ! only if variable size dist
312 0 : if (volc(m) > 1.e-39_r8 .and. na(m) > 1.e-39_r8) then
313 0 : etafactor2(m) = 1._r8/(na(m)*beta*sqrtg) !fixed or variable size dist
314 : ! number mode radius (m)
315 0 : amcubeloc(m) = (3._r8*volc(m)/(4._r8*pi*exp45logsig(m)*na(m))) ! only if variable size dist
316 0 : smc(m) = smcrit(m) ! only for prescribed size dist
317 :
318 0 : if (hygro_aer(m) > 1.e-10_r8) then ! loop only if variable size dist
319 0 : smc(m) = 2._r8*aten*sqrt(aten/(27._r8*hygro_aer(m)*amcubeloc(m)))
320 : else
321 0 : smc(m) = 100._r8
322 : endif
323 : else
324 0 : smc(m) = 1._r8
325 0 : etafactor2(m) = etafactor2max ! this should make eta big if na is very small.
326 : end if
327 0 : lnsmloc(m) = log(smc(m)) ! only if variable size dist
328 : end do
329 :
330 : ! single updraft
331 0 : wnuc = wbar
332 0 : alw = alpha*wnuc
333 0 : sqrtalw = sqrt(alw)
334 0 : zeta = 2._r8*sqrtalw*aten/(3._r8*sqrtg)
335 0 : etafactor1 = 2._r8*alw*sqrtalw
336 :
337 0 : do m = 1, nmode
338 : ! skip aerosols with no dispersion, since they aren't meant to be CCN
339 0 : if (dispersion_aer(m) /= 0._r8) eta(m) = etafactor1*etafactor2(m)
340 : end do
341 :
342 0 : call maxsat(zeta, eta, nmode, smc, smax)
343 0 : lnsmax = log(smax)
344 :
345 0 : nact = 0._r8
346 0 : do m = 1, nmode
347 : ! skip aerosols with no dispersion, since they aren't meant to be CCN
348 0 : if (dispersion_aer(m) == 0._r8) cycle
349 0 : x = 2*(lnsmloc(m) - lnsmax)/(3*sq2*alogsig(m))
350 0 : nact = nact + 0.5_r8*(1._r8 - erf(x))*na(m)
351 : end do
352 0 : nact = nact/rhoair ! convert from #/m3 to #/kg
353 :
354 0 : deallocate( &
355 : volc, &
356 0 : eta, &
357 0 : smc, &
358 0 : etafactor2, &
359 0 : amcubeloc, &
360 0 : lnsmloc )
361 :
362 0 : end subroutine ndrop_bam_run
363 :
364 : !===============================================================================
365 :
366 0 : subroutine ndrop_bam_ccn(lchnk, ncol, maerosol, naer2)
367 :
368 : !-------------------------------------------------------------------------------
369 : !
370 : ! Write diagnostic bulk aerosol ccn concentration
371 : !
372 : !-------------------------------------------------------------------------------
373 :
374 : ! Input arguments
375 : integer, intent(in) :: lchnk ! chunk identifier
376 : integer, intent(in) :: ncol ! number of columns
377 : real(r8), intent(in) :: naer2(:,:,:) ! aerosol number concentration (1/m3)
378 : real(r8), intent(in) :: maerosol(:,:,:) ! aerosol mass conc (kg/m3)
379 :
380 : ! Local variables
381 : integer :: i, k, l, m
382 : real(r8) :: arg
383 :
384 : character*8, parameter :: ccn_name(psat)=(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
385 : real(r8) :: amcubesulfate(pcols) ! cube of dry mode radius (m) of sulfate
386 : real(r8) :: smcritsulfate(pcols) ! critical supersatuation for activation of sulfate
387 : real(r8) :: ccnfactsulfate
388 : real(r8) :: ccn(pcols,pver,psat) ! number conc of aerosols activated at supersat
389 : !-------------------------------------------------------------------------------
390 :
391 0 : ccn(:ncol,:,:) = 0._r8
392 :
393 0 : do k = top_lev, pver
394 :
395 0 : do m = 1, naer_all
396 :
397 0 : if (m == idxsul) then
398 : ! Lohmann treatment for sulfate has variable size distribution
399 0 : do i = 1, ncol
400 0 : if (naer2(i,k,m) > 0._r8) then
401 0 : amcubesulfate(i) = amcubefactor(m)*maerosol(i,k,m)/(naer2(i,k,m))
402 0 : smcritsulfate(i) = smcritfactor(m)/sqrt(amcubesulfate(i))
403 : else
404 0 : smcritsulfate(i) = 1._r8
405 : end if
406 : end do
407 : end if
408 :
409 0 : do l = 1, psat
410 :
411 0 : if (m == idxsul) then
412 : ! This code is modifying ccnfact for sulfate only.
413 0 : do i = 1, ncol
414 0 : arg = argfactor(m)*log(smcritsulfate(i)/super(l))
415 0 : if (arg < 2) then
416 0 : if (arg < -2) then
417 : ccnfactsulfate = 1.0e-6_r8
418 : else
419 0 : ccnfactsulfate = 0.5e-6_r8*erfc(arg)
420 : end if
421 : else
422 : ccnfactsulfate = 0.0_r8
423 : end if
424 0 : ccn(i,k,l) = ccn(i,k,l) + naer2(i,k,m)*ccnfactsulfate
425 : end do
426 : else
427 : ! Non-sulfate species use ccnfact computed by the init routine
428 0 : ccn(:ncol,k,l) = ccn(:ncol,k,l) + naer2(:ncol,k,m)*ccnfact(l,m)
429 : end if
430 :
431 : end do ! supersaturation
432 : end do ! bulk aerosol
433 : end do ! level
434 :
435 0 : do l = 1, psat
436 0 : call outfld(ccn_name(l), ccn(1,1,l), pcols, lchnk)
437 : end do
438 :
439 0 : do l = 1, naer_all
440 0 : call outfld(trim(aername(l))//'_m3', naer2(:,:,l), pcols, lchnk)
441 : end do
442 :
443 0 : end subroutine ndrop_bam_ccn
444 :
445 : !===============================================================================
446 :
447 0 : subroutine maxsat(zeta, eta, nmode, smc, smax)
448 :
449 : ! calculates maximum supersaturation for multiple
450 : ! competing aerosol modes.
451 :
452 : ! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
453 : ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
454 :
455 : real(r8), intent(in) :: zeta
456 : integer, intent(in) :: nmode ! number of modes
457 : real(r8), intent(in) :: smc(:) ! critical supersaturation for number mode radius
458 : real(r8), intent(in) :: eta(:)
459 :
460 : real(r8), intent(out) :: smax ! maximum supersaturation
461 :
462 : integer :: m ! mode index
463 : real(r8) :: sum, g1, g2
464 :
465 0 : do m=1,nmode
466 0 : if(zeta.gt.1.e5_r8*eta(m).or.smc(m)*smc(m).gt.1.e5_r8*eta(m))then
467 : ! weak forcing. essentially none activated
468 0 : smax=1.e-20_r8
469 : else
470 : ! significant activation of this mode. calc activation all modes.
471 : go to 1
472 : endif
473 : enddo
474 :
475 0 : return
476 :
477 : 1 continue
478 :
479 : sum=0
480 0 : do m=1,nmode
481 0 : if(eta(m).gt.1.e-20_r8)then
482 0 : g1=sqrt(zeta/eta(m))
483 0 : g1=g1*g1*g1
484 0 : g2=smc(m)/sqrt(eta(m)+3*zeta)
485 0 : g2=sqrt(g2)
486 0 : g2=g2*g2*g2
487 0 : sum=sum+(f1(m)*g1+f2(m)*g2)/(smc(m)*smc(m))
488 : else
489 : sum=1.e20_r8
490 : endif
491 : enddo
492 :
493 0 : smax=1._r8/sqrt(sum)
494 :
495 : end subroutine maxsat
496 :
497 : !===============================================================================
498 :
499 : end module ndrop_bam
|
