Line data Source code
1 : module mo_aerosols
2 : !-----------------------------------------------------------------
3 : !
4 : ! this module computes the production of ammonium nitrate
5 : ! using the formulation by Seinfeld and Pandis (p531, 1998)
6 : ! with the simplification of activity coefficients and
7 : ! aerosol molality using the parameterizations
8 : ! from Metzger et al. (JGR, ACH-16, 107(D16), 2002)
9 : !
10 : ! written by Jean-Francois Lamarque (April 2004)
11 : ! adapted for CAM (May 2004)
12 : !
13 : !-----------------------------------------------------------------
14 :
15 : use shr_kind_mod, only : r8 => shr_kind_r8
16 : use ppgrid, only : pver
17 : use cam_logfile, only: iulog
18 :
19 : private
20 : public :: aerosols_inti,aerosols_formation
21 : public :: has_aerosols
22 :
23 : save
24 :
25 : integer, target :: spc_ndx(5)
26 : integer, pointer :: nh3_ndx, nh4no3_ndx, nh4_ndx
27 : integer, pointer :: so4_ndx, hno3_ndx
28 : integer :: xhno3_ndx, xnh4no3_ndx
29 : integer :: nu_i(2)
30 : real(r8) :: zeta_inv
31 : real(r8) :: z_i(2)
32 : logical :: has_aerosols = .true.
33 :
34 : contains
35 :
36 1536 : subroutine aerosols_inti()
37 :
38 : use mo_chem_utls, only : get_spc_ndx
39 : use cam_history, only : addfld
40 : use spmd_utils, only : masterproc
41 :
42 : implicit none
43 :
44 : !-----------------------------------------------------------------
45 : ! ... local variables
46 : !-----------------------------------------------------------------
47 : integer :: m
48 :
49 1536 : nh3_ndx => spc_ndx(1)
50 1536 : nh4no3_ndx => spc_ndx(2)
51 1536 : so4_ndx => spc_ndx(3)
52 1536 : hno3_ndx => spc_ndx(4)
53 1536 : nh4_ndx => spc_ndx(5)
54 :
55 : !-----------------------------------------------------------------
56 : ! ... set species index
57 : !-----------------------------------------------------------------
58 1536 : nh3_ndx = get_spc_ndx( 'NH3' )
59 1536 : nh4no3_ndx = get_spc_ndx( 'NH4NO3' )
60 1536 : so4_ndx = get_spc_ndx( 'SO4' )
61 1536 : hno3_ndx = get_spc_ndx( 'HNO3' )
62 1536 : nh4_ndx = get_spc_ndx( 'NH4' )
63 1536 : xnh4no3_ndx = get_spc_ndx( 'XNH4NO3' )
64 1536 : xhno3_ndx = get_spc_ndx( 'XHNO3' )
65 :
66 1536 : has_aerosols = all( spc_ndx(:) > 0 )
67 1536 : if( .not. has_aerosols ) then
68 1536 : if (masterproc) then
69 2 : write(iulog,*) '-----------------------------------------'
70 2 : write(iulog,*) 'mozart will NOT do nh4no3'
71 2 : write(iulog,*) 'following species are missing'
72 12 : do m = 1,size(spc_ndx)
73 12 : if( spc_ndx(m) < 1 ) then
74 10 : write(iulog,*) m
75 : end if
76 : end do
77 2 : write(iulog,*) '-----------------------------------------'
78 : endif
79 1536 : return
80 : else
81 0 : if (masterproc) then
82 0 : write(iulog,*) '-----------------------------------------'
83 0 : write(iulog,*) 'mozart will do nh4no3'
84 0 : write(iulog,*) '-----------------------------------------'
85 : end if
86 : end if
87 :
88 : !
89 : ! define parameters
90 : !
91 : ! ammonium nitrate (NH4NO3)
92 : !
93 0 : zeta_inv = 1._r8/4._r8
94 0 : nu_i(1) = 4
95 0 : z_i (1) = 1._r8
96 : !
97 : ! ammonium sulfate
98 : !
99 0 : nu_i(2) = 4
100 0 : z_i (2) = 0.5_r8
101 :
102 0 : call addfld ('TSO4_VMR', (/ 'lev' /), 'A', 'mol/mol','total sulfate in mo_aerosols')
103 0 : call addfld ('THNO3_VMR',(/ 'lev' /), 'A','mol/mol','total nitric acid in mo_aerosols')
104 :
105 0 : return
106 1536 : end subroutine aerosols_inti
107 :
108 0 : subroutine aerosols_formation( ncol, lchnk, tfld, rh, qin)
109 :
110 1536 : use ppgrid, only : pcols, pver
111 : use chem_mods, only : gas_pcnst, adv_mass
112 : use cam_history, only : outfld
113 :
114 : implicit none
115 : !
116 : ! input arguments
117 : !
118 : !
119 : ! input arguments
120 : !
121 : integer, intent(in) :: ncol ! number columns in chunk
122 : integer, intent(in) :: lchnk ! chunk index
123 : real(r8), intent(in) :: tfld(:,:) ! temperature
124 : real(r8), intent(in) :: rh(:,:) ! relative humidity
125 : real(r8), intent(inout) :: qin(:,:,:) ! xported species ( vmr )
126 :
127 : !
128 : ! local variables
129 : !
130 : integer :: i,j,k,n
131 : integer :: domain_number ! concentration domain
132 : real(r8) :: sulfate_state ! fraction of sulfate neutralized by ammonia
133 0 : real(r8), dimension(ncol,pver) :: tso4, & ! total sulfate
134 0 : thno3,& ! total nitric acid
135 0 : txhno3 ! total nitric acid ( XHNO3 )
136 : real(r8) :: tnh3, & ! total ammonia
137 : fnh3, & ! free ammonia
138 : rhd, & ! relative humidity of deliquescence
139 : gamma, & ! activity coefficient
140 : ssm_nh4no3, & ! single solute molality for NH4NO3
141 : ta, & ! total ammonia
142 : tn, & ! total nitrate
143 : kp, & ! equilibrium constant
144 : nh4no3 ! ammonium nitrate produced
145 : real(r8) :: log_t
146 : real(r8) :: ti
147 : real(r8) :: xnh4no3
148 :
149 0 : do k=1,pver
150 0 : do i=1,ncol
151 :
152 : !
153 : ! compute total concentrations
154 : !
155 0 : tnh3 = (qin(i,k, nh3_ndx)+qin(i,k,nh4_ndx))
156 0 : tso4(i,k) = qin(i,k,so4_ndx)
157 : !
158 : ! define concentration domain
159 : !
160 0 : if ( tnh3 < tso4(i,k) ) then
161 : domain_number = 4
162 : sulfate_state = 1.0_r8
163 0 : elseif ( tnh3 < 2._r8*tso4(i,k) ) then
164 : domain_number = 3
165 : sulfate_state = 1.5_r8
166 : else
167 0 : domain_number = 2
168 0 : sulfate_state = 2.0_r8
169 : endif
170 : !
171 : ! define free ammonia (ammonia available for ammonium nitrate production)
172 : !
173 0 : fnh3 = tnh3 - sulfate_state * tso4(i,k)
174 0 : fnh3 = max(0._r8,fnh3)
175 : !
176 : ! convert initial concentrations to ppbv
177 : !
178 0 : tso4(i,k) = tso4(i,k) * 1.e9_r8
179 0 : tnh3 = tnh3 * 1.e9_r8
180 0 : fnh3 = fnh3 * 1.e9_r8
181 0 : thno3(i,k) = (qin(i,k,hno3_ndx)+qin(i,k,nh4no3_ndx)) * 1.e9_r8
182 0 : if ( xhno3_ndx > 0 .and. xnh4no3_ndx > 0 ) then
183 0 : txhno3(i,k) = (qin(i,k,xhno3_ndx)+qin(i,k,xnh4no3_ndx)) * 1.e9_r8
184 : endif
185 : !
186 : ! compute relative humidity of deliquescence (%) for NH4NO3
187 : ! (Seinfeld and Pandis, p532)
188 : !
189 0 : ti = 1._r8/tfld(i,k)
190 0 : rhd = 0.01_r8 * exp( 1.6954_r8 + 723.7_r8*ti )
191 0 : log_t = log( tfld(i,k)/298._r8 )
192 0 : if ( rh(i,k) < rhd ) then
193 : !
194 : ! crystalline ammonium nitrate
195 : !
196 : ! compute equilibrium constant
197 : !
198 0 : kp = exp( 84.6_r8 - 24220._r8*ti - 6.1_r8*log_t )
199 : !
200 : else
201 : !
202 : ! aqueous phase ammonium nitrate
203 : !
204 : ! compute activity coefficients (from Menzger et al.)
205 : !
206 0 : n = domain_number
207 0 : gamma = (rh(i,k)**n/(1000._r8/n*(1._r8-rh(i,k))+n))**zeta_inv
208 : !
209 : ! compute single solute molality for NH4NO3
210 : !
211 0 : ssm_nh4no3 = (1000._r8 * 0.81_r8 * nu_i(1) * (1._r8/rh(i,k)-1._r8)/80._r8)**z_i(1)
212 : !
213 : ! compute equilibrium constant
214 : !
215 0 : kp = (gamma*ssm_nh4no3)**2 * exp( 53.19_r8 - 15850.62_r8*ti + 11.51_r8*log_t )
216 :
217 : endif
218 : !
219 : ! calculate production of NH4NO3 (in ppbv) using Seinfeld and Pandis (p534, 1998)
220 : !
221 0 : ta = fnh3
222 0 : tn = thno3(i,k)
223 0 : nh4no3 = 0.5_r8 * (ta + tn - sqrt(max(0._r8,(ta+tn)**2 - 4._r8*(ta*tn-kp))))
224 0 : nh4no3 = max(0._r8,nh4no3)
225 0 : if ( xhno3_ndx > 0 .and. xnh4no3_ndx > 0 ) then
226 0 : tn = txhno3(i,k)
227 0 : xnh4no3 = 0.5_r8 * (ta + tn - sqrt(max(0._r8,(ta+tn)**2 - 4._r8*(ta*tn-kp))))
228 0 : xnh4no3 = max(0._r8,xnh4no3)
229 : endif
230 : !
231 : ! reset concentrations according to equilibrium calculation
232 : !
233 0 : qin(i,k,nh4no3_ndx) = nh4no3
234 0 : if ( xhno3_ndx > 0 ) then
235 0 : qin(i,k,xnh4no3_ndx) = xnh4no3
236 : endif
237 0 : qin(i,k,nh3_ndx ) = max(0._r8,(fnh3-nh4no3))
238 0 : qin(i,k,nh4_ndx ) = max(0._r8,(tnh3-(fnh3-nh4no3)))
239 0 : qin(i,k,hno3_ndx ) = max(0._r8,(thno3(i,k)-nh4no3))
240 0 : if ( xhno3_ndx > 0 ) then
241 0 : qin(i,k,xhno3_ndx ) = max(0._r8,(txhno3(i,k)-xnh4no3))
242 : endif
243 0 : qin(i,k,so4_ndx ) = tso4(i,k)
244 : !
245 : ! convert from ppbv to vmr
246 : !
247 0 : qin(i,k,nh4no3_ndx) = qin(i,k,nh4no3_ndx) * 1.e-9_r8
248 0 : qin(i,k,nh3_ndx ) = qin(i,k,nh3_ndx ) * 1.e-9_r8
249 0 : qin(i,k,nh4_ndx ) = qin(i,k,nh4_ndx ) * 1.e-9_r8
250 0 : qin(i,k,hno3_ndx ) = qin(i,k,hno3_ndx ) * 1.e-9_r8
251 0 : qin(i,k,so4_ndx ) = qin(i,k,so4_ndx ) * 1.e-9_r8
252 0 : if ( xhno3_ndx > 0 ) then
253 0 : qin(i,k,xnh4no3_ndx) = qin(i,k,xnh4no3_ndx) * 1.e-9_r8
254 : endif
255 0 : if ( xhno3_ndx > 0 ) then
256 0 : qin(i,k,xhno3_ndx ) = qin(i,k,xhno3_ndx ) * 1.e-9_r8
257 : endif
258 :
259 : end do
260 : end do
261 : !
262 : ! outputs
263 : !
264 0 : call outfld ('TSO4_VMR' ,tso4 (:ncol,:), ncol, lchnk )
265 0 : call outfld ('THNO3_VMR',thno3(:ncol,:), ncol, lchnk )
266 :
267 0 : return
268 0 : end subroutine aerosols_formation
269 :
270 : end module mo_aerosols
|
