1 module module_mosaic_ext
5 ! scm_temp - subr dumpxx is temporary for scm testing
6 subroutine dumpxx( dumpch, dtchem, t_k, p_atm, ah2o, &
7 jaerosolstate, jphase, jhyst_leg, &
8 aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
11 use module_data_mosaic_kind, only: r8
13 use module_data_mosaic_aero, only: &
14 iso4_a, ino3_a, icl_a, inh4_a, &
15 ina_a, ioin_a, ioc_a, ibc_a, &
16 ipcg2_b_c_a, ipcg2_b_o_a, ipcg1_b_c_a, ipcg1_b_o_a, &
17 iopcg1_b_c_a, iopcg1_b_o_a, ipcg2_f_c_a, ipcg2_f_o_a, &
18 ipcg1_f_c_a, ipcg1_f_o_a, iopcg1_f_c_a, iopcg1_f_o_a, &
19 iant1_c_a, iant1_o_a, ibiog1_c_a, ibiog1_o_a, &
21 jtotal, naer, nbin_a, nbin_a_max, ngas_aerchtot, ngas_volatile, &
26 integer, intent(in), dimension(nbin_a_max) :: jaerosolstate, jphase, jhyst_leg
28 real(r8), intent(in) :: dtchem, t_k, p_atm, ah2o
29 real(r8), intent(in), dimension(naer,3,nbin_a_max) :: aer
30 real(r8), intent(in), dimension(ngas_volatile) :: gas
31 real(r8), intent(in), dimension(nbin_a_max) :: num_a, water_a, water_a_hyst, dp_dry_a
33 character(len=*), intent(in) :: dumpch
35 type (mosaic_vars_aa_type), intent(in) :: mosaic_vars_aa
37 integer, parameter :: lunaa = 131
38 integer, parameter :: naer_dump = 24
39 integer :: ii, iv, ix, jj, kk, ktau
40 integer :: iaer_dump(naer_dump) = 24
41 character(len=80) :: fmtaa
42 character(len=2) :: tmpch2
43 character(len=256) :: tmpchaa
44 character(len= 32) :: tmpchbb
47 ! return ! scm_temp - for non-scm might want to activate this line
49 if ( mosaic_vars_aa%idiagbb_host < 100 ) return
51 ii = mosaic_vars_aa%hostgridinfo(2)
52 jj = mosaic_vars_aa%hostgridinfo(3)
53 kk = mosaic_vars_aa%hostgridinfo(4)
54 if ( ii*jj*kk /= 1 ) return
56 ktau = mosaic_vars_aa%it_host
58 (/ iso4_a, ino3_a, icl_a, inh4_a, &
59 ina_a, ioin_a, ioc_a, ibc_a, &
60 ipcg2_b_c_a, ipcg2_b_o_a, ipcg1_b_c_a, ipcg1_b_o_a, &
61 iopcg1_b_c_a, iopcg1_b_o_a, ipcg2_f_c_a, ipcg2_f_o_a, &
62 ipcg1_f_c_a, ipcg1_f_o_a, iopcg1_f_c_a, iopcg1_f_o_a, &
63 iant1_c_a, iant1_o_a, ibiog1_c_a, ibiog1_o_a /)
68 write(lunaa,'(/2a,i5,3i3)') tmpch2, 'dump', &
70 write(lunaa,'(a,f11.2,f11.4,f11.4 )') 't p a ', t_k, p_atm, ah2o
71 ! if (tmpch2 == 'aa') &
72 ! write(lunaa,'(a,f11.2 )') 'cairm3', cair_mol_m3
74 ! write(lunaa,'(a,1p,8e11.3)') 'gas1:4 ', gas(1:4)
77 if (nbin_a > 8) fmtaa = '(a,8i11/(10x,8i11))'
78 write(lunaa,fmtaa) 'jstate ', jaerosolstate(1:nbin_a)
79 write(lunaa,fmtaa) 'jphase ', jphase(1:nbin_a)
80 write(lunaa,fmtaa) 'jhyst ', jhyst_leg(1:nbin_a)
82 fmtaa = '(a,1p,8e11.3)'
83 if (nbin_a > 8) fmtaa = '(a,1p,8e11.3/(10x,1p,8e11.3))'
84 write(lunaa,fmtaa) 'num ', num_a(1:nbin_a)
85 write(lunaa,fmtaa) 'dpdry ', dp_dry_a(1:nbin_a)
86 write(lunaa,fmtaa) 'water ', water_a(1:nbin_a)
87 write(lunaa,fmtaa) 'hyswtr ', water_a_hyst(1:nbin_a)
89 ! write(lunaa,'(a,1p,8e11.3)') 'so4 ', aer(iso4_a,jtotal,1:nbin_a)
90 ! write(lunaa,'(a,1p,8e11.3)') 'no3 ', aer(ino3_a,jtotal,1:nbin_a)
91 ! write(lunaa,'(a,1p,8e11.3)') 'cl ', aer(icl_a, jtotal,1:nbin_a)
92 ! write(lunaa,'(a,1p,8e11.3)') 'nh4 ', aer(inh4_a,jtotal,1:nbin_a)
93 !!! write(lunaa,'(a,1p,8e11.3)') 'msa ', aer(imsa_a,jtotal,1:nbin_a)
94 !!! write(lunaa,'(a,1p,8e11.3)') 'co3 ', aer(ico3_a,jtotal,1:nbin_a)
95 ! write(lunaa,'(a,1p,8e11.3)') 'na ', aer(ina_a, jtotal,1:nbin_a)
96 !!! write(lunaa,'(a,1p,8e11.3)') 'ca ', aer(ica_a, jtotal,1:nbin_a)
97 ! write(lunaa,'(a,1p,8e11.3)') 'oin ', aer(ioin_a,jtotal,1:nbin_a)
98 ! write(lunaa,'(a,1p,8e11.3)') 'oc ', aer(ioc_a, jtotal,1:nbin_a)
99 ! write(lunaa,'(a,1p,8e11.3)') 'bc ', aer(ibc_a, jtotal,1:nbin_a)
102 !need to dump soa species
105 if (iv < 1 .or. iv > naer) cycle
106 fmtaa = '(a,1p,8e11.3)'
107 write(tmpchaa,fmtaa) aer_name(iv)(1:10), aer(iv,jtotal,1:min(nbin_a,8))
108 if (iv <= ngas_aerchtot) then
109 write(tmpchbb,fmtaa) gas_name(iv)(1:10), gas(iv)
110 tmpchaa = trim(tmpchaa) // ' ' // trim(tmpchbb)
112 write(lunaa,'(a)') trim(tmpchaa)
113 fmtaa = '(10x,1p,8e11.3)'
114 if (nbin_a > 8) write(lunaa,fmtaa) aer(iv,jtotal,9:nbin_a)
119 end subroutine dumpxx
125 !***********************************************************************
126 ! determines phase state of an aerosol bin. includes kelvin effect.
128 ! author: Rahul A. Zaveri
130 !-----------------------------------------------------------------------
131 subroutine aerosol_phase_state( ibin, jaerosolstate, jphase, aer, &
132 jhyst_leg, electrolyte, epercent, kel, activity, mc, num_a, mass_wet_a, mass_dry_a, &
133 mass_soluble_a, vol_dry_a, vol_wet_a, water_a, water_a_hyst, water_a_up, aH2O_a, &
134 aH2O, ma, gam, log_gamZ, zc, za, gam_ratio, &
135 xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, partial_molar_vol, sigma_soln, T_K, & ! RAZ deleted a_zsr
136 RH_pc, mw_aer_mac, dens_aer_mac, sigma_water, Keq_ll, Keq_sl, MW_a, MW_c, &
137 growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, jsalt_present, jsalt_index, &
138 jsulf_poor, jsulf_rich, phi_salt_old, &
139 kappa_nonelectro, mosaic_vars_aa )
141 use module_data_mosaic_aero, only: r8, nbin_a_max, &
142 ngas_aerchtot, ngas_volatile, nelectrolyte, &!Parameters
143 Ncation, naer, jtotal, all_solid, jhyst_up, all_liquid, Nanion, nrxn_aer_ll, &
144 nrxn_aer_sl, nsalt, MDRH_T_NUM, jsulf_poor_NUM, jsulf_rich_NUM, &!Parameters
145 inh4_a, ina_a, ica_a, ico3_a, imsa_a, icl_a, ino3_a, iso4_a, & ! TBD
146 a_zsr, b_zsr, aw_min, &! RAZ added a_zsr, b_zsr, aw_min
153 integer, intent(in):: ibin
154 integer, intent(in), dimension(nsalt) :: jsalt_index
155 integer, intent(in), dimension(jsulf_poor_NUM) :: jsulf_poor
156 integer, intent(in), dimension(jsulf_rich_NUM) :: jsulf_rich
158 real(r8), intent(in) :: aH2O,T_K,RH_pc,rtol_mesa
159 real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
160 real(r8), intent(in), dimension(Ncation) :: zc,MW_c
161 real(r8), intent(in), dimension(Nanion) :: za,MW_a
162 real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
163 real(r8), intent(in), dimension(ngas_aerchtot) :: partial_molar_vol
164 real(r8), intent(in), dimension(naer) :: kappa_nonelectro
167 integer, intent(inout), dimension(nsalt) :: jsalt_present
168 integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg
170 real(r8), intent(inout) :: sigma_water
171 real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
172 real(r8), intent(inout), dimension(Nanion) :: xeq_a,na_Ma
173 real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_wet_a,mass_dry_a
174 real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,gam_ratio
175 real(r8), intent(inout), dimension(nbin_a_max) :: vol_dry_a,vol_wet_a,water_a
176 real(r8), intent(inout), dimension(nbin_a_max) :: water_a_hyst,water_a_up,aH2O_a
177 real(r8), intent(inout), dimension(nbin_a_max) :: sigma_soln,growth_factor,MDRH
178 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 ! RAZ 5/20/2014
179 real(r8), intent(inout), dimension (nrxn_aer_ll) :: Keq_ll
180 real(r8), intent(inout), dimension (nrxn_aer_sl) :: Keq_sl
181 real(r8), intent(inout), dimension(MDRH_T_NUM) :: MDRH_T
182 real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kel
183 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
184 real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
185 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
186 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
187 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
188 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
189 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
190 real(r8), intent(inout), dimension(nsalt) :: phi_salt_old
192 type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
195 character(len=256) :: errmsg
196 integer, parameter :: aer_pha_sta_diagaa = -1 !BALLI- changed from 100 to -1
197 integer, parameter :: iter_kelvin_method = 3
198 ! iter_kelvin_method = 1 - use rahuls original iteration method
199 ! iter_kelvin_method = 2 - use bisection
200 ! iter_kelvin_method = 3 - start with rahuls original iteration method, but if it fails, switch to bisection
201 integer, parameter :: iter_kelvin_meth1_max = 10
202 integer, parameter :: iter_kelvin_meth2_max = 100
203 integer :: iaer, iv, itmpa
204 integer :: iter_kelvin, iter_kelvin_meth1, iter_kelvin_state
208 real(r8):: aH2O_range_bisect_toler
209 real(r8) :: aH2O_a_new, aH2O_a_old, aH2O_a_oldn, aH2O_a_oldp, aH2O_a_del_state3
210 real(r8), dimension(nbin_a_max) :: DpmV
211 real(r8), dimension(nbin_a_max) :: kelvin
212 real(r8) :: kelvin_old, kelvin_oldn, kelvin_oldp
213 real(r8) :: kelvin_toler
214 real(r8) :: rel_err, rel_err_old, rel_err_old2, rel_err_oldn, rel_err_oldp
215 real(r8) :: term, tmpa
216 real(r8) :: water_a_old, water_a_oldn, water_a_oldp
219 if (aer_pha_sta_diagaa >= 3) &
220 write(*,'(/a,5i5,2f12.8,1p,2e11.3)') 'aer_pha_sta_a', ibin, jhyst_leg(ibin), jaerosolstate(ibin), -1, 0, aH2O, aH2O_a(ibin)
221 !aH2O = RH_pc*0.01 !**BALLI, this is already done in init subr
224 do iv = 1, ngas_aerchtot
228 ! if(RH_pc .le. 97.0)then ! RAZ
229 ! kelvin_toler = 1.e-4
231 ! kelvin_toler = 1.e-10 ! RAZ
233 ! define error tolerances become stricter as aH2O approaches 1.0
234 kelvin_toler = 1.e-6_r8 * max( 1.0_r8-aH2O, 1.0e-4_r8 )
235 aH2O_range_bisect_toler = 1.e-6_r8 * max( 1.0_r8-aH2O, 1.0e-4_r8 )
238 ! calculate dry mass and dry volume of a bin
239 mass_dry_a(ibin) = 0.0 ! initialize to 0.0
240 vol_dry_a(ibin) = 0.0 ! initialize to 0.0
242 aer_H = (2.*aer(iso4_a,jtotal,ibin) + &
243 aer(ino3_a,jtotal,ibin) + &
244 aer(icl_a,jtotal,ibin) + &
245 aer(imsa_a,jtotal,ibin) + &
246 2.*aer(ico3_a,jtotal,ibin))- &
247 (2.*aer(ica_a,jtotal,ibin) + &
248 aer(ina_a,jtotal,ibin) + &
249 aer(inh4_a,jtotal,ibin))
250 aer_H = max(aer_H, 0.0d0)
253 mass_dry_a(ibin) = mass_dry_a(ibin) + &
254 aer(iaer,jtotal,ibin)*mw_aer_mac(iaer) ! ng/m^3(air)
255 vol_dry_a(ibin) = vol_dry_a(ibin) + &
256 aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer) ! ncc/m^3(air)
258 mass_dry_a(ibin) = mass_dry_a(ibin) + aer_H
259 vol_dry_a(ibin) = vol_dry_a(ibin) + aer_H
261 mass_dry_a(ibin) = mass_dry_a(ibin)*1.e-15 ! g/cc(air)
262 vol_dry_a(ibin) = vol_dry_a(ibin)*1.e-15 ! cc(aer)/cc(air) or m^3/m^3(air)
264 ! wet mass and wet volume
265 mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3 ! g/cc(air)
266 vol_wet_a(ibin) = vol_dry_a(ibin) + water_a(ibin)*1.e-3 ! cc(aer)/cc(air) or m^3/m^3(air)
269 water_a_up(ibin) = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr) ! for hysteresis curve determination
272 iter_kelvin_meth1 = 0
274 iter_kelvin_state = 0
275 if (iter_kelvin_method == 2) iter_kelvin_state = 2
279 rel_err_old = 1.0e30_r8
280 rel_err_old2 = 1.0e30_r8
283 aH2O_a_del_state3 = 1.0e-3_r8
288 rel_err_oldn = 1.0e30_r8
289 rel_err_oldp = 1.0e30_r8
290 water_a_oldp = 0.0_r8
291 water_a_oldn = 0.0_r8
295 10 iter_kelvin = iter_kelvin + 1
296 aH2O_a(ibin) = aH2O_a_new
298 ! RAZ uncommented the next 3 lines
299 do je = 1, nelectrolyte
300 molality0(je,ibin) = bin_molality(je,ibin,aH2O_a,b_zsr,a_zsr,aw_min) ! compute aH2O dependent binary molalities EFFI
302 call MESA( ibin, jaerosolstate, jphase, aer, jhyst_leg, &
303 electrolyte, epercent, activity, mc, num_a, mass_wet_a, mass_dry_a, &
304 mass_soluble_a, vol_dry_a, vol_wet_a, water_a, water_a_hyst, water_a_up, aH2O_a, &
305 aH2O, ma, gam, log_gamZ, zc, za, gam_ratio, &
306 xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_ll, &
307 Keq_sl, MW_c, MW_a, growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, &
308 jsalt_present, jsalt_index, jsulf_poor, jsulf_rich, phi_salt_old, &
309 kappa_nonelectro, mosaic_vars_aa )
311 if(jaerosolstate(ibin) .eq. all_solid)then
312 if (aer_pha_sta_diagaa >= 2) &
313 write(*,'(a,5i5,2f12.8,1p,2e11.3)') 'aer_pha_sta_b', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
314 iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin)
317 ! new wet mass and wet volume
318 mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3 ! g/cc(air)
319 vol_wet_a(ibin) = vol_dry_a(ibin) + water_a(ibin)*1.e-3 ! cc(aer)/cc(air) or m^3/m^3(air)
321 call calculate_kelvin(ibin,num_a,vol_wet_a,aH2O_a,DpmV,kelvin,sigma_soln,T_K, &
324 kelvin(ibin) = max( kelvin(ibin), 1.0_r8 )
325 if (water_a(ibin) <= 0.0_r8) kelvin(ibin) = 1.0_r8
327 aH2O_a_new = aH2O/kelvin(ibin)
329 ! if(RH_pc .le. 97.0)then
330 ! rel_err = abs( (aH2O_a_new - aH2O_a(ibin))/aH2O_a(ibin))
332 ! if(water_a(ibin) .gt. 0.0)then
333 ! rel_err = abs( (water_a(ibin) - water_a_old)/water_a(ibin))
335 ! rel_err = 0.0 ! no soluble material is present
338 ! the above rel_err involve differences between current and previous
339 ! iteration values, and is not suitable for bisection
340 ! this rel_err below uses error from the exact solution, and is suitable for bisection
341 rel_err = (aH2O_a(ibin)*kelvin(ibin) - aH2O) / max( aH2O, 0.01_r8 )
343 if (aer_pha_sta_diagaa >= 10) &
344 write(*,'(a,2i5, 1p,e10.2, 0p,f14.10, 2x,2f14.10, 2x,1p,2e18.10)') &
345 'iter_kelvin', iter_kelvin_state, iter_kelvin, rel_err, kelvin(ibin), &
346 aH2O_a(ibin), aH2O_a_new, water_a_old, water_a(ibin)
348 if (abs(rel_err) <= kelvin_toler) then
349 iter_kelvin_state = iter_kelvin_state + 100
353 if (iter_kelvin_state <= 0) then
354 ! doing rahuls original iteration method
356 if (iter_kelvin >= iter_kelvin_meth1_max) then
358 else if (iter_kelvin >= iter_kelvin_meth1_max) then
359 tmpa = min( rel_err_old, rel_err_old2 )
360 if (tmpa < 0.0_r8 .and. rel_err <= tmpa) itmpa = 1
361 tmpa = max( rel_err_old, rel_err_old2 )
362 if (tmpa > 0.0_r8 .and. rel_err >= tmpa) itmpa = 1
366 if (iter_kelvin_method <= 1) then
367 ! quit if number of iterations is too large OR
368 ! rel_err is outside the range of the previous two rel_err values,
369 ! and one previous rel_err is positive, and one previous rel_err is negative
370 aH2O_a(ibin) = aH2O_a_new ! do this to get same output as prev version
371 if (aer_pha_sta_diagaa >= 1) &
372 write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_err', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
373 iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
374 iter_kelvin_state = 100
377 ! switch to method 2 but do not iterate yet
378 iter_kelvin_state = 1
379 iter_kelvin_meth1 = iter_kelvin
382 ! save current values to old then do next iteration
383 aH2O_a_old = aH2O_a(ibin)
384 kelvin_old = kelvin(ibin)
385 rel_err_old2 = rel_err_old
386 rel_err_old = rel_err
387 water_a_old = water_a(ibin)
389 ! call MTEM_compute_log_gamZ ! recompute activity coeffs (for surface tension and solid-liquid equilibria)
394 if (iter_kelvin_state == 1) then
395 ! rahuls original iteration method failed, so do some things before switching to bisection
396 iter_kelvin_state = 2
397 if (rel_err < 0.0_r8) then
398 ! current aH2O_a has negative rel_err so must start at the beginning
402 ! current aH2O_a has positive rel_err and can be used in bisection
408 if (iter_kelvin_state == 2) then
409 ! this is first "setup" step of bisection, and the algorithm is expecting that
410 ! the current aH2O_a has hel_err be > 0, and can be used as one of the 2 bisection points
411 if (rel_err < 0.0_r8) then
412 ! error should be positive, so this is a fatal error
413 if (aer_pha_sta_diagaa >= 1) &
414 write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er2', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
415 iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
416 iter_kelvin_state = 100
419 ! current aH2O_a will work as one of the two initial bisection points
420 ! (the one with a positive error)
421 aH2O_a_oldp = aH2O_a(ibin)
422 kelvin_oldp = kelvin(ibin)
423 rel_err_oldp = rel_err
424 water_a_oldp = water_a(ibin)
425 aH2O_a_new = min( aH2O/kelvin(ibin), 0.999999_r8 ) ! is this needed, or should it be 1.0, or ???
426 iter_kelvin_state = 3
430 if (iter_kelvin_state == 3) then
431 ! this is the second "setup" step of bisection, and the algorithm is looking for an aH2O_a
432 ! that has rel_err < 0, so that the "root" will be bracketed and bisection can begin
433 if (rel_err < 0.0_r8) then
434 ! current aH2O_a will work as one of the two initial bisection points
435 ! (the one with a negative error)
436 aH2O_a_oldn = aH2O_a(ibin)
437 kelvin_oldn = kelvin(ibin)
438 rel_err_oldn = rel_err
439 water_a_oldn = water_a(ibin)
440 aH2O_a_new = 0.5_r8*(aH2O_a_oldn + aH2O_a_oldp)
441 iter_kelvin_state = 4
444 ! need to find a point with a negative error
445 if ( (rel_err >= rel_err_oldp) .or. &
446 (aH2O_a_del_state3 >= 0.999_r8) ) then
447 ! cannot find such a point -- this is a fatal error
448 if (aer_pha_sta_diagaa >= 1) &
449 write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er3', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
450 iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
451 iter_kelvin_state = 200
454 ! save current aH2O_a as the initial bisection point with positive error
455 ! then calc aH2O_a_new = aH2O_a(ibin) - aH2O_a_del_state3
456 ! which will hopefully have a negative error
457 aH2O_a_oldp = aH2O_a(ibin)
458 kelvin_oldp = kelvin(ibin)
459 rel_err_oldp = rel_err
460 water_a_oldp = water_a(ibin)
461 aH2O_a_new = aH2O_a(ibin) - aH2O_a_del_state3
462 aH2O_a_del_state3 = aH2O_a_del_state3*1.5_r8
463 if (aH2O_a_new .le. 0.01_r8) then
465 aH2O_a_del_state3 = 1.0_r8
472 if (iter_kelvin_state == 4) then
473 ! at this point, the algorithm is doing bisection
474 if ( iter_kelvin >= iter_kelvin_meth2_max + iter_kelvin_meth1 ) then
475 ! maximum iterations is exceeded
476 if (aer_pha_sta_diagaa >= 1) &
477 write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er4', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
478 iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
479 iter_kelvin_state = 301
481 else if ( abs(aH2O_a_oldp - aH2O_a_oldn) <= aH2O_range_bisect_toler ) then
482 ! the aH2O_a_oldp to aH2O_a_oldn range is very small, which is treated as convergence
483 ! if (aer_pha_sta_diagaa >= 1) &
484 ! write(*,'(a,5i5,2f12.8,1p,4e11.3)') 'iter_kelv_er5', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
485 ! iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler, &
486 ! aH2O_range_bisect_toler
487 iter_kelvin_state = 302
490 ! decide if the current aH2O_a should replace the old negative-error point
491 ! or the old positive-error point
492 if (rel_err >= 0.0_r8) then
493 if (rel_err >= rel_err_oldp) then
494 ! current aH2O_a has positive error, but the error is not smaller
495 ! than the old positive-error point -- this is a fatal error
496 if (aer_pha_sta_diagaa >= 1) &
497 write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er6', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
498 iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
499 iter_kelvin_state = 303
502 ! current aH2O_a has positive error and replaces the the old positive-error point
503 aH2O_a_oldp = aH2O_a(ibin)
504 kelvin_oldp = kelvin(ibin)
505 rel_err_oldp = rel_err
506 water_a_oldp = water_a(ibin)
509 if (rel_err <= rel_err_oldn) then
510 ! current aH2O_a has negative error, but the error is not smaller
511 ! than the old negative-error point -- this is a fatal error
512 if (aer_pha_sta_diagaa >= 1) &
513 write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er7', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
514 iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
515 iter_kelvin_state = 304
518 ! current aH2O_a has negative error and replaces the the old negative-error point
519 aH2O_a_oldn = aH2O_a(ibin)
520 kelvin_oldn = kelvin(ibin)
521 rel_err_oldn = rel_err
522 water_a_oldn = water_a(ibin)
525 aH2O_a_new = 0.5_r8*(aH2O_a_oldn + aH2O_a_oldp)
529 write(errmsg,'(a,4i5)') 'iter_kelv fatal err 1', ibin, iter_kelvin, iter_kelvin_state
530 call wrf_error_fatal(trim(adjustl(errmsg)))
533 ! kelvin iterations completed
534 90 if (iter_kelvin_state == 200) then
535 ! select aH2O_a(ibin) or aH2O_a_oldp, whichever has lowest error
536 if (abs(rel_err_oldp) < abs(rel_err)) then
537 aH2O_a(ibin) = aH2O_a_oldp
538 rel_err = rel_err_oldp
540 else if (iter_kelvin_state >= 300 .and. iter_kelvin_state <= 304) then
541 ! select aH2O_a(ibin) or aH2O_a_oldp or aH2O_a_oldn, whichever has lowest error
542 tmpa = min( abs(rel_err_oldn), abs(rel_err_oldp), abs(rel_err) )
543 if (abs(rel_err_oldp) == tmpa) then
544 aH2O_a(ibin) = aH2O_a_oldp
545 rel_err = rel_err_oldp
546 else if (abs(rel_err_oldn) == tmpa) then
547 aH2O_a(ibin) = aH2O_a_oldn
548 rel_err = rel_err_oldn
552 if(jaerosolstate(ibin) .eq. all_liquid)jhyst_leg(ibin) = jhyst_up
554 ! now compute kelvin effect terms for condensing species (nh3, hno3, and hcl)
555 do iv = 1, ngas_aerchtot
556 term = 4.*sigma_soln(ibin)*partial_molar_vol(iv)/ &
557 (8.3144e7*T_K*DpmV(ibin))
558 kel(iv,ibin) = 1. + term*(1. + 0.5*term*(1. + term/3.))
561 if (aer_pha_sta_diagaa >= 2) &
562 write(*,'(a,5i5,2f12.8,1p,e11.3,e14.5)') 'aer_pha_sta_c', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
563 iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, water_a(ibin)
565 end subroutine aerosol_phase_state
569 !**********************************`*************************************
570 ! MESA: Multicomponent Equilibrium Solver for Aerosols.
571 ! Computes equilibrum solid and liquid phases by integrating
572 ! pseudo-transient dissolution and precipitation reactions
574 ! author: Rahul A. Zaveri
577 ! 9/3/2015 RAZ: Bugfix - fixed phase state calculations for aerosols that dont contain any salts,
578 ! but can still contain water due to presence of BC, OC, SOA, and OIN, which are now
579 ! allowed to absorb some water.
580 !-----------------------------------------------------------------------
581 subroutine MESA( ibin, jaerosolstate, jphase, aer, jhyst_leg, &
582 electrolyte, epercent, activity, mc, num_a, mass_wet_a, mass_dry_a, mass_soluble_a, &
583 vol_dry_a, vol_wet_a, water_a, water_a_hyst, water_a_up, aH2O_a, aH2O, &
584 ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, &
585 na_Ma, nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_ll, Keq_sl, MW_c, &
586 MW_a, growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, jsalt_present, jsalt_index, &
587 jsulf_poor, jsulf_rich, phi_salt_old, &
588 kappa_nonelectro, mosaic_vars_aa )
590 use module_data_mosaic_aero, only: r8, nbin_a_max, nelectrolyte, Ncation, naer, &!Parameters
591 jtotal, all_solid, jsolid, all_liquid, jliquid, jhyst_lo, mhyst_uporlo_jhyst, &!Parameters
592 jhyst_up, mhyst_uporlo_waterhyst, nsoluble, nsalt, Nanion, nrxn_aer_sl, &
593 nrxn_aer_ll, MDRH_T_NUM, jsulf_poor_NUM, jsulf_rich_NUM, &!Parameters
594 ptol_mol_astem, mhyst_force_lo, mhyst_force_up, &!Input
595 jcacl2, jcano3, mhyst_method, ioin_a, ibc_a, jcaco3, jcaso4, & !TBD
603 integer, intent(in) :: ibin
604 integer, intent(inout), dimension(nbin_a_max) :: jhyst_leg
605 integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase
606 integer, intent(in), dimension(nsalt) :: jsalt_index
607 integer, intent(inout), dimension(nsalt) :: jsalt_present
608 integer, intent(in), dimension(jsulf_poor_NUM) :: jsulf_poor
609 integer, intent(in), dimension(jsulf_rich_NUM) :: jsulf_rich
611 real(r8), intent(in) :: aH2O,rtol_mesa
612 real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
613 real(r8), intent(in), dimension(Ncation) :: zc,MW_c
614 real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
615 real(r8), intent(in), dimension(Nanion) :: za,MW_a
616 real(r8), intent(inout), dimension(Nanion) :: xeq_a,na_Ma
617 real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_wet_a,mass_dry_a
618 real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,vol_dry_a
619 real(r8), intent(inout), dimension(nbin_a_max) :: vol_wet_a,gam_ratio
620 real(r8), intent(inout), dimension(nbin_a_max) :: water_a,water_a_hyst,water_a_up
621 real(r8), intent(inout), dimension(nbin_a_max) :: aH2O_a,growth_factor,MDRH
622 real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
623 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
624 real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
625 real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
626 real(r8), intent(inout), dimension(MDRH_T_NUM) :: MDRH_T
627 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
628 real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
629 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
630 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
631 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
632 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
633 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
634 real(r8), intent(inout), dimension(nsalt) :: phi_salt_old
635 real(r8), intent(in), dimension(naer) :: kappa_nonelectro
637 type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
640 integer :: idissolved, j_index, jsalt_dum, jdum, js, je ! 9/3/2015 RAZ: added jsalt_dum
641 real(r8) :: CRH, solids, sum_soluble, sum_insoluble, XT !BALLI** XT, should it be subr arg?
642 !real(r8) :: aerosol_water ! mosaic func
643 !real(r8) :: drh_mutual ! mosaic func
644 real(r8) :: H_ion, sum_dum
648 !! calculate percent composition
650 do je = 1, nelectrolyte
651 sum_dum = sum_dum + electrolyte(je,jtotal,ibin)
654 if(sum_dum .eq. 0.)sum_dum = 1.0
656 do je = 1, nelectrolyte
657 epercent(je,jtotal,ibin) = 100.*electrolyte(je,jtotal,ibin)/sum_dum
661 call calculate_XT(ibin,jtotal,XT,aer)
665 !! begin new algorithm - 6/3/2015 RAZ
666 jsalt_dum = 0 ! 9/3/2015 RAZ
668 jsalt_present(js) = 0 ! default value - salt absent
670 if(epercent(js,jtotal,ibin) .gt. ptol_mol_astem)then
671 jsalt_present(js) = 1 ! salt present
672 jsalt_dum = jsalt_dum + jsalt_index(js) ! 9/3/2015 RAZ
677 if( (epercent(jcano3,jtotal,ibin) .gt. ptol_mol_astem) .or. &
678 (epercent(jcacl2,jtotal,ibin) .gt. ptol_mol_astem) )then
679 CRH = 0.0 ! no crystrallization or efflorescence point
681 CRH = 0.35 ! default value
685 if(jsalt_dum .eq. 0)then ! no salts or acids are present ! 9/3/2015 RAZ: updated algorithm for jsalt_dum = 0
689 jaerosolstate(ibin) = all_solid
690 jphase(ibin) = jsolid
691 jhyst_leg(ibin) = jhyst_lo
692 call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
693 water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg, & ! 9/3/2015 RAZ: water due to nonelectrolytes (OC, BC, SOA, OIN)
694 electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
697 elseif(XT .lt. 1. .and. XT .gt. 0.0)then ! excess sulfate, always liquid, MDRH=0.0
699 elseif(XT .ge. 2.0 .or. XT .lt. 0.0)then ! sulfate poor
700 j_index = jsulf_poor(jsalt_dum) ! 9/3/2015 RAZ
701 MDRH(ibin) = MDRH_T(j_index)
703 j_index = jsulf_rich(jsalt_dum) ! 9/3/2015 RAZ
704 MDRH(ibin) = MDRH_T(j_index)
707 CRH = min(CRH, MDRH(ibin)/100.0) ! 6/3/2015 RAZ
709 !! end new algorithm - 6/3/2015 RAZ
712 ! modified step 1: 9/3/2015 RAZ
713 ! step 1: check if aH2O is below CRH (crystallization or efflorescence point)
714 if( aH2O_a(ibin).lt.CRH )then
715 jaerosolstate(ibin) = all_solid
716 jphase(ibin) = jsolid
717 jhyst_leg(ibin) = jhyst_lo
718 call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
719 water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg, & ! 9/3/2015 RAZ: water due to nonelectrolytes (OC, BC, SOA, OIN)
720 electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
725 ! step 2: check mhyst_method for supersaturation/metastable state
727 if (mhyst_method == mhyst_uporlo_jhyst) then ! BOX method/logic
728 if (jhyst_leg(ibin) == jhyst_up) jdum = 1
729 elseif (mhyst_method == mhyst_uporlo_waterhyst) then ! 3-D method/logic
730 if (water_a_hyst(ibin) > 0.5*water_a_up(ibin)) jdum = 1
731 !BSINGH - 05/28/2013(RCE updates)
732 elseif (mhyst_method == mhyst_force_lo) then
734 elseif (mhyst_method == mhyst_force_up) then
736 !BSINGH - 05/28/2013(RCE updates ENDS)
738 call wrf_error_fatal('*** MESA - bad mhyst_method')
740 if (jdum == 1) then ! the aerosol is fully deliquesced in metastable or subsaturated state
741 call do_full_deliquescence(ibin,aer,electrolyte)
743 ! call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,nc_Mc,xeq_c) ! for Li and Lu surface tension
744 ! call compute_activities(ibin,jphase,aer,jhyst_leg,electrolyte, &
745 !activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,log_gamZ,gam_ratio) ! for Li and Lu surface tension
750 ! MODIFIED LOGIC IF SOA, POA, BC, OIN ARE ASSUMED TO BE SLIGHTLY HYGROSCOPIC RAZ 4/16/2014
752 ! do js = 1, nsoluble
753 ! sum_soluble = sum_soluble + electrolyte(js,jtotal,ibin)
756 ! solids = electrolyte(jcaso4,jtotal,ibin) + &
757 ! electrolyte(jcaco3,jtotal,ibin) + &
758 ! aer(ioin_a,jtotal,ibin) + &
759 ! aer(ibc_a,jtotal,ibin)
762 ! if(sum_soluble .le. 0.0 .and. solids .gt. 0.0)then ! RAZ modified logic
765 ! jaerosolstate(ibin) = all_solid ! no soluble material present, so go back to solid state
766 ! jphase(ibin) = jsolid
767 ! call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
769 ! ! new wet mass and wet volume
770 ! mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3 ! g/cc(air)
771 ! vol_wet_a(ibin) = vol_dry_a(ibin) + water_a(ibin)*1.e-3 ! cc(aer)/cc(air) or m^3/m^3(air)
772 ! growth_factor(ibin) = mass_wet_a(ibin)/mass_dry_a(ibin) ! mass growth factor
776 ! elseif(sum_soluble .gt. 0.0)then ! RAZ modified logic
778 jaerosolstate(ibin) = all_liquid
779 jhyst_leg(ibin) = jhyst_up
780 jphase(ibin) = jliquid
781 water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg, &
782 electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
783 if(water_a(ibin) .le. 0.0)then ! one last attempt to catch bad input
785 jaerosolstate(ibin) = all_solid ! no soluble material present
786 jphase(ibin) = jsolid
787 jhyst_leg(ibin) = jhyst_lo
788 call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
790 call adjust_liquid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent)
791 call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg, &
792 electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a, &
793 aH2O,ma,gam,log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)
796 ! new wet mass and wet volume
797 mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3 ! g/cc(air)
798 vol_wet_a(ibin) = vol_dry_a(ibin) + water_a(ibin)*1.e-3 ! cc(aer)/cc(air) or m^3/m^3(air)
799 growth_factor(ibin) = mass_wet_a(ibin)/mass_dry_a(ibin) ! mass growth factor
809 ! step 3: diagnose phase state based on MDRH
810 if(aH2O_a(ibin)*100. .lt. MDRH(ibin)) then
811 jaerosolstate(ibin) = all_solid
812 jphase(ibin) = jsolid
813 jhyst_leg(ibin) = jhyst_lo
814 call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
819 ! step 4: none of the above means it must be sub-saturated or mixed-phase
820 10 call do_full_deliquescence(ibin,aer,electrolyte)
821 call MESA_PTC( ibin, jaerosolstate, jphase, aer, jhyst_leg, &
822 electrolyte, epercent, activity, mc, num_a, mass_dry_a, mass_wet_a, &
823 mass_soluble_a, vol_dry_a, vol_wet_a, water_a, aH2O, &
824 ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, &
825 mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_sl, MW_c, MW_a, Keq_ll, &
826 growth_factor, molality0, rtol_mesa, jsalt_present, phi_salt_old, &
827 kappa_nonelectro, mosaic_vars_aa ) ! determines jaerosolstate(ibin)
833 !***********************************************************************
834 ! computes kelvin effect term (kelvin => 1.0)
836 ! author: Rahul A. Zaveri
838 !-----------------------------------------------------------------------
839 subroutine calculate_kelvin(ibin,num_a,vol_wet_a,aH2O_a,DpmV,kelvin,sigma_soln, &
841 use module_data_mosaic_constants, only: pi
842 use module_data_mosaic_aero, only: r8,nbin_a_max !Parameters
847 integer, intent(in) :: ibin
848 real(r8), intent(in) :: T_K,sigma_water
849 real(r8), intent(in), dimension(nbin_a_max) :: num_a
850 real(r8), intent(inout), dimension(nbin_a_max) :: sigma_soln
851 real(r8), intent(inout), dimension(nbin_a_max) ::vol_wet_a,aH2O_a,DpmV,kelvin
854 real(r8) :: term, sum_dum
855 real(r8), dimension(nbin_a_max) :: volume_a
857 volume_a(ibin) = vol_wet_a(ibin) ! [cc/cc(air)]
858 DpmV(ibin)=(6.*volume_a(ibin)/(num_a(ibin)*pi))**(1./3.) ! [cm]
861 ! Li and Lu (2001) surface tension model:
863 ! do je = 1, nelectrolyte
864 ! sum_dum = sum_dum + G_MX(je)*
865 ! & alog(1./(1.+K_MX(je)*activity(je,ibin)))
867 ! sigma_soln(ibin) = sigma_water + 8.3144e7*T_K*sum_dum
870 ! simpler correlation for solution surface tension:
871 sigma_soln(ibin) = sigma_water + 49.0*(1. - aH2O_a(ibin)) ! [dyn/cm]
875 term = 72.*sigma_soln(ibin)/(8.3144e7*T_K*DpmV(ibin)) ! [-]
876 ! kelvin(ibin) = exp(term)
877 kelvin(ibin) = 1. + term*(1. + 0.5*term*(1. + term/3.))
881 end subroutine calculate_kelvin
885 !***********************************************************************
886 ! computes sulfate ratio
888 ! author: Rahul A. Zaveri
890 !-----------------------------------------------------------------------
891 subroutine calculate_XT(ibin,jp,XT,aer)
892 use module_data_mosaic_aero, only: r8,naer,nbin_a_max, &
893 imsa_a,iso4_a,ica_a,ina_a,inh4_a
898 integer, intent(in) :: ibin, jp
899 real(r8), intent(inout) :: XT
900 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
903 if( (aer(iso4_a,jp,ibin)+aer(imsa_a,jp,ibin)) .gt.0.0)then
904 XT = ( aer(inh4_a,jp,ibin) + &
905 aer(ina_a,jp,ibin) + &
906 2.*aer(ica_a,jp,ibin) )/ &
907 (aer(iso4_a,jp,ibin)+0.5*aer(imsa_a,jp,ibin))
914 end subroutine calculate_XT
918 !***********************************************************************
919 ! called when aerosol bin is completely solid.
921 ! author: Rahul A. Zaveri
923 !-----------------------------------------------------------------------
924 subroutine adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent, &
927 use module_data_mosaic_aero, only: r8,nbin_a_max,naer,nelectrolyte,jsolid, &!Parameters
928 jhyst_lo,jtotal,jliquid, &!Parameters
929 inh4_a,ino3_a,icl_a !TBD
934 integer, intent(in) :: ibin
935 integer, intent(inout), dimension(nbin_a_max) :: jphase,jhyst_leg
936 real(r8), intent(inout), dimension(nbin_a_max) :: water_a
937 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
938 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte,epercent
943 jphase(ibin) = jsolid
944 jhyst_leg(ibin) = jhyst_lo ! lower curve
947 ! transfer aer(jtotal) to aer(jsolid)
949 aer(iaer, jsolid, ibin) = aer(iaer,jtotal,ibin)
950 aer(iaer, jliquid,ibin) = 0.0
953 ! transfer electrolyte(jtotal) to electrolyte(jsolid)
954 do je = 1, nelectrolyte
955 electrolyte(je,jliquid,ibin) = 0.0
956 epercent(je,jliquid,ibin) = 0.0
957 electrolyte(je,jsolid,ibin) = electrolyte(je,jtotal,ibin)
958 epercent(je,jsolid,ibin) = epercent(je,jtotal,ibin)
961 ! update aer(jtotal) that may have been affected above
962 aer(inh4_a,jtotal,ibin) = aer(inh4_a,jsolid,ibin)
963 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid,ibin)
964 aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid,ibin)
968 end subroutine adjust_solid_aerosol
972 !***********************************************************************
973 ! called when aerosol bin is completely liquid.
975 ! author: Rahul A. Zaveri
977 !-----------------------------------------------------------------------
978 subroutine adjust_liquid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent) ! TOUCH
980 use module_data_mosaic_aero, only: r8,nbin_a_max,naer,nelectrolyte,jliquid, &!Parameters
981 jhyst_up,jsolid,jtotal, &!Parameters
983 inh4_a,ina_a,ica_a,imsa_a,icl_a,ino3_a,iso4_a
988 integer, intent(in) :: ibin
989 integer, intent(inout), dimension(nbin_a_max) :: jphase,jhyst_leg
991 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
992 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
993 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
997 jphase(ibin) = jliquid
998 jhyst_leg(ibin) = jhyst_up ! upper curve
1000 ! partition all electrolytes into liquid phase
1001 do je = 1, nelectrolyte
1002 electrolyte(je,jsolid,ibin) = 0.0
1003 epercent(je,jsolid,ibin) = 0.0
1004 electrolyte(je,jliquid,ibin) = electrolyte(je,jtotal,ibin)
1005 epercent(je,jliquid,ibin) = epercent(je,jtotal,ibin)
1007 ! except these electrolytes, which always remain in the solid phase
1008 electrolyte(jcaco3,jsolid,ibin) = electrolyte(jcaco3,jtotal,ibin)
1009 electrolyte(jcaso4,jsolid,ibin) = electrolyte(jcaso4,jtotal,ibin)
1010 epercent(jcaco3,jsolid,ibin) = epercent(jcaco3,jtotal,ibin)
1011 epercent(jcaso4,jsolid,ibin) = epercent(jcaso4,jtotal,ibin)
1012 electrolyte(jcaco3,jliquid,ibin)= 0.0
1013 electrolyte(jcaso4,jliquid,ibin)= 0.0
1014 epercent(jcaco3,jliquid,ibin) = 0.0
1015 epercent(jcaso4,jliquid,ibin) = 0.0
1018 ! partition all the aer species into
1021 aer(iaer,jsolid,ibin) = aer(iaer,jtotal,ibin)
1023 aer(iso4_a,jsolid,ibin) = electrolyte(jcaso4,jsolid,ibin)
1024 aer(ino3_a,jsolid,ibin) = 0.0
1025 aer(icl_a,jsolid,ibin) = 0.0
1026 aer(inh4_a,jsolid,ibin) = 0.0
1027 aer(imsa_a,jsolid,ibin) = 0.0
1028 aer(ina_a,jsolid,ibin) = 0.0
1029 aer(ica_a,jsolid,ibin) = electrolyte(jcaco3,jsolid,ibin) + &
1030 electrolyte(jcaso4,jsolid,ibin)
1031 ! aer(ico3_a,jsolid,ibin) = aer(ico3_a,jtotal,ibin)
1032 ! aer(ioc_a,jsolid,ibin) = aer(ioc_a,jtotal,ibin)
1033 ! aer(ibc_a, jsolid,ibin)= aer(ibc_a, jtotal,ibin)
1034 ! aer(ioin_a, jsolid,ibin)= aer(ioin_a, jtotal,ibin)
1035 ! aer(iaro1_a,jsolid,ibin)= aer(iaro1_a,jtotal,ibin)
1036 ! aer(iaro2_a,jsolid,ibin)= aer(iaro2_a,jtotal,ibin)
1037 ! aer(ialk1_a,jsolid,ibin)= aer(ialk1_a,jtotal,ibin)
1038 ! aer(iole1_a,jsolid,ibin)= aer(iole1_a,jtotal,ibin)
1039 ! aer(iapi1_a,jsolid,ibin)= aer(iapi1_a,jtotal,ibin)
1040 ! aer(iapi2_a,jsolid,ibin)= aer(iapi2_a,jtotal,ibin)
1041 ! aer(ilim1_a,jsolid,ibin)= aer(ilim1_a,jtotal,ibin)
1042 ! aer(ilim2_a,jsolid,ibin)= aer(ilim2_a,jtotal,ibin)
1046 aer(iaer,jliquid,ibin) = 0.0
1048 aer(iso4_a,jliquid,ibin) = aer(iso4_a,jtotal,ibin) - &
1049 aer(iso4_a,jsolid,ibin)
1050 aer(iso4_a,jliquid,ibin) = max(0.d0, aer(iso4_a,jliquid,ibin)) ! RAZ 4/16/2014
1051 aer(ino3_a,jliquid,ibin) = aer(ino3_a,jtotal,ibin)
1052 aer(icl_a,jliquid,ibin) = aer(icl_a,jtotal,ibin)
1053 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jtotal,ibin)
1054 aer(imsa_a,jliquid,ibin) = aer(imsa_a,jtotal,ibin)
1055 aer(ina_a,jliquid,ibin) = aer(ina_a,jtotal,ibin)
1056 aer(ica_a,jliquid,ibin) = aer(ica_a,jtotal,ibin) - &
1057 aer(ica_a,jsolid,ibin)
1058 aer(ica_a,jliquid,ibin) = max(0.d0, aer(ica_a,jliquid,ibin)) ! RAZ 4/16/2014
1059 ! aer(ico3_a,jliquid,ibin) = 0.0
1060 ! aer(ioc_a,jliquid,ibin) = 0.0
1061 ! aer(ibc_a, jliquid,ibin)= 0.0
1062 ! aer(ioin_a, jliquid,ibin)= 0.0
1063 ! aer(iaro1_a,jliquid,ibin)= 0.0
1064 ! aer(iaro2_a,jliquid,ibin)= 0.0
1065 ! aer(ialk1_a,jliquid,ibin)= 0.0
1066 ! aer(iole1_a,jliquid,ibin)= 0.0
1067 ! aer(iapi1_a,jliquid,ibin)= 0.0
1068 ! aer(iapi2_a,jliquid,ibin)= 0.0
1069 ! aer(ilim1_a,jliquid,ibin)= 0.0
1070 ! aer(ilim2_a,jliquid,ibin)= 0.0
1073 end subroutine adjust_liquid_aerosol
1077 !***********************************************************************
1078 ! this subroutine completely deliquesces an aerosol and partitions
1079 ! all the soluble electrolytes into the liquid phase and insoluble
1080 ! ones into the solid phase. It also calculates the corresponding
1081 ! aer(js,jliquid,ibin) and aer(js,jsolid,ibin) generic species
1084 ! author: Rahul A. Zaveri
1086 !-----------------------------------------------------------------------
1087 subroutine do_full_deliquescence(ibin,aer,electrolyte) ! TOUCH
1088 use module_data_mosaic_aero, only: r8,naer,nbin_a_max,nelectrolyte,jtotal,jsolid, &!Parameters
1089 jliquid, &!Parameters
1090 jcacl2,jcano3,ioin_a,jcaco3,jcaso4, &
1091 inh4_a,ina_a,ica_a,imsa_a,icl_a,ino3_a,iso4_a
1098 integer, intent(in) :: ibin
1099 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1100 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
1104 ! partition all electrolytes into liquid phase
1105 do js = 1, nelectrolyte
1106 electrolyte(js,jsolid,ibin) = 0.0
1107 electrolyte(js,jliquid,ibin) = electrolyte(js,jtotal,ibin)
1110 ! except these electrolytes, which always remain in the solid phase
1111 electrolyte(jcaco3,jsolid,ibin) = electrolyte(jcaco3,jtotal,ibin)
1112 electrolyte(jcaso4,jsolid,ibin) = electrolyte(jcaso4,jtotal,ibin)
1113 electrolyte(jcaco3,jliquid,ibin)= 0.0
1114 electrolyte(jcaso4,jliquid,ibin)= 0.0
1117 ! partition all the generic aer species into solid and liquid phases
1120 aer(iaer,jsolid,ibin)= aer(iaer,jtotal,ibin)
1122 aer(iso4_a,jsolid,ibin) = electrolyte(jcaso4,jsolid,ibin)
1123 aer(ino3_a,jsolid,ibin) = 0.0
1124 aer(icl_a, jsolid,ibin) = 0.0
1125 aer(inh4_a,jsolid,ibin) = 0.0
1126 aer(imsa_a,jsolid,ibin) = 0.0
1127 aer(ina_a, jsolid,ibin) = 0.0
1128 aer(ica_a, jsolid,ibin) = electrolyte(jcaco3,jsolid,ibin) + &
1129 electrolyte(jcaso4,jsolid,ibin)
1130 ! aer(ico3_a,jsolid,ibin) = aer(ico3_a,jtotal,ibin)
1131 ! aer(ioc_a, jsolid,ibin) = aer(ioc_a,jtotal,ibin)
1132 ! aer(ibc_a, jsolid,ibin) = aer(ibc_a,jtotal,ibin)
1133 ! aer(ioin_a,jsolid,ibin) = aer(ioin_a,jtotal,ibin)
1134 ! aer(iaro1_a,jsolid,ibin)= aer(iaro1_a,jtotal,ibin)
1135 ! aer(iaro2_a,jsolid,ibin)= aer(iaro2_a,jtotal,ibin)
1136 ! aer(ialk1_a,jsolid,ibin)= aer(ialk1_a,jtotal,ibin)
1137 ! aer(iole1_a,jsolid,ibin)= aer(iole1_a,jtotal,ibin)
1138 ! aer(iapi1_a,jsolid,ibin)= aer(iapi1_a,jtotal,ibin)
1139 ! aer(iapi2_a,jsolid,ibin)= aer(iapi2_a,jtotal,ibin)
1140 ! aer(ilim1_a,jsolid,ibin)= aer(ilim1_a,jtotal,ibin)
1141 ! aer(ilim2_a,jsolid,ibin)= aer(ilim2_a,jtotal,ibin)
1145 aer(iaer,jliquid,ibin) = 0.0
1147 aer(iso4_a,jliquid,ibin) = max(0.0_r8, aer(iso4_a,jtotal,ibin) - &
1148 electrolyte(jcaso4,jsolid,ibin)) ! added max() RAZ 4/16/2014
1149 aer(ino3_a,jliquid,ibin) = aer(ino3_a,jtotal,ibin)
1150 aer(icl_a, jliquid,ibin) = aer(icl_a,jtotal,ibin)
1151 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jtotal,ibin)
1152 aer(imsa_a,jliquid,ibin) = aer(imsa_a,jtotal,ibin)
1153 aer(ina_a, jliquid,ibin) = aer(ina_a,jtotal,ibin)
1154 aer(ica_a, jliquid,ibin) = electrolyte(jcano3,jtotal,ibin) + &
1155 electrolyte(jcacl2,jtotal,ibin)
1156 ! aer(ioc_a, jliquid,ibin) = 0.0
1157 ! aer(ico3_a,jliquid,ibin) = 0.0
1158 ! aer(ibc_a, jliquid,ibin) = 0.0
1159 ! aer(ioin_a,jliquid,ibin) = 0.0
1160 ! aer(iaro1_a,jliquid,ibin)= 0.0
1161 ! aer(iaro2_a,jliquid,ibin)= 0.0
1162 ! aer(ialk1_a,jliquid,ibin)= 0.0
1163 ! aer(iole1_a,jliquid,ibin)= 0.0
1164 ! aer(iapi1_a,jliquid,ibin)= 0.0
1165 ! aer(iapi2_a,jliquid,ibin)= 0.0
1166 ! aer(ilim1_a,jliquid,ibin)= 0.0
1167 ! aer(ilim2_a,jliquid,ibin)= 0.0
1170 end subroutine do_full_deliquescence
1174 !***********************************************************************
1175 ! MESA: Multicomponent Equilibrium Solver for Aerosol-phase
1176 ! computes equilibrum solid and liquid phases by integrating
1177 ! pseudo-transient dissolution and precipitation reactions
1179 ! author: Rahul A. Zaveri
1181 ! Reference: Zaveri R.A., R.C. Easter, and L.K. Peters, JGR, 2005b
1182 !-----------------------------------------------------------------------
1183 subroutine MESA_PTC(ibin, jaerosolstate, jphase, aer, jhyst_leg, &
1184 electrolyte, epercent, activity, mc, num_a, mass_dry_a, mass_wet_a, mass_soluble_a, &
1185 vol_dry_a, vol_wet_a, water_a, aH2O, ma, gam, &
1186 log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac, &
1187 dens_aer_mac, Keq_sl, MW_c, MW_a, Keq_ll, growth_factor, molality0, rtol_mesa, &
1188 jsalt_present, phi_salt_old, &
1189 kappa_nonelectro, mosaic_vars_aa ) ! TOUCH
1191 use module_data_mosaic_aero, only: r8, nbin_a_max, nelectrolyte, Ncation, naer, nsalt, &!Parameters
1192 jhyst_lo, mixed, all_liquid, jsolid, jliquid, jtotal, mYES, &!Parameters
1193 all_solid, Nanion, nrxn_aer_sl, nrxn_aer_ll, &!Parameters
1194 ino3_a, iso4_a, ioc_a, ilim1_a, ilim2_a, inh4_a, ina_a, ica_a, ico3_a, imsa_a, icl_a, &
1200 integer, intent(in) :: ibin
1201 integer, intent(inout), dimension(nsalt) :: jsalt_present
1202 integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg
1204 real(r8), intent(in) :: aH2O,rtol_mesa
1205 real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
1206 real(r8), intent(in), dimension(Ncation) :: zc,MW_c
1207 real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
1208 real(r8), intent(in), dimension(Nanion) :: za,MW_a
1209 real(r8), intent(inout), dimension(Nanion) :: xeq_a,na_Ma
1210 real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_dry_a,mass_wet_a
1211 real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,vol_dry_a
1212 real(r8), intent(inout), dimension(nbin_a_max) :: growth_factor
1213 real(r8), intent(inout), dimension(nbin_a_max) :: vol_wet_a,water_a,gam_ratio
1214 real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
1215 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
1216 real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
1217 real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1218 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
1219 real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
1220 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1221 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1222 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1223 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
1224 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
1225 real(r8), intent(inout), dimension(nsalt) :: phi_salt_old
1226 real(r8), intent(in), dimension(naer) :: kappa_nonelectro
1228 type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
1231 integer iaer, iconverge, iconverge_flux, iconverge_mass, &
1232 idissolved, itdum, js, je, jp
1234 real(r8) :: tau_p(nsalt), tau_d(nsalt)
1235 real(r8) :: frac_solid, sumflux, hsalt_min, alpha, XT, dumdum, &
1237 real(r8) :: phi_prod, alpha_fac, sum_dum
1238 real(r8) :: aer_H,hsalt_max
1239 real(r8), dimension(nelectrolyte) :: eleliquid
1240 real(r8), dimension(nbin_a_max) :: mass_dry_salt
1241 real(r8), dimension(nsalt) :: phi_salt,flux_sl,phi_bar,alpha_salt
1242 real(r8), dimension(nsalt) :: sat_ratio,hsalt
1245 !real(r8) :: aerosol_water
1248 itdum = 0 ! initialize time
1262 !! EFFI calculate percent composition
1264 do je = 1, nelectrolyte
1265 sum_dum = sum_dum + electrolyte(je,jtotal,ibin)
1268 if(sum_dum .eq. 0.)sum_dum = 1.0
1270 do je = 1, nelectrolyte
1271 epercent(je,jtotal,ibin) = 100.*electrolyte(je,jtotal,ibin)/sum_dum
1278 jsalt_present(js) = 0 ! default value - salt absent
1279 if(epercent(js,jtotal,ibin) .gt. 1.0)then
1280 jsalt_present(js) = 1 ! salt present
1285 mass_dry_a(ibin) = 0.0
1287 aer_H = (2.*aer(iso4_a,jtotal,ibin) + &
1288 aer(ino3_a,jtotal,ibin) + &
1289 aer(icl_a,jtotal,ibin) + &
1290 aer(imsa_a,jtotal,ibin) + &
1291 2.*aer(ico3_a,jtotal,ibin))- &
1292 (2.*aer(ica_a,jtotal,ibin) + &
1293 aer(ina_a,jtotal,ibin) + &
1294 aer(inh4_a,jtotal,ibin))
1295 aer_H = max(aer_H, 0.0d0)
1298 mass_dry_a(ibin) = mass_dry_a(ibin) + &
1299 aer(iaer,jtotal,ibin)*mw_aer_mac(iaer) ! [ng/m^3(air)]
1300 vol_dry_a(ibin) = vol_dry_a(ibin) + &
1301 aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer) ! ncc/m^3(air)
1303 mass_dry_a(ibin) = mass_dry_a(ibin) + aer_H
1304 vol_dry_a(ibin) = vol_dry_a(ibin) + aer_H
1306 mass_dry_a(ibin) = mass_dry_a(ibin)*1.e-15 ! [g/cc(air)]
1307 vol_dry_a(ibin) = vol_dry_a(ibin)*1.e-15 ! [cc(aer)/cc(air)]
1309 mass_dry_salt(ibin) = 0.0 ! soluble salts only
1311 mass_dry_salt(ibin) = mass_dry_salt(ibin) + &
1312 electrolyte(je,jtotal,ibin)*mw_electrolyte(je)*1.e-15 ! g/cc(air)
1315 mosaic_vars_aa%jMESA_call = mosaic_vars_aa%jMESA_call + 1
1317 !----begin pseudo time continuation loop-------------------------------
1319 do 500 itdum = 1, mosaic_vars_aa%Nmax_MESA
1322 ! compute new salt fluxes
1323 call MESA_flux_salt(ibin,jaerosolstate,jphase, aer,jhyst_leg,electrolyte, &
1324 epercent,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,&
1325 gam,log_gamZ,zc,za,gam_ratio,xeq_a,na_Ma,nc_Mc,xeq_c,mw_electrolyte, &
1326 Keq_sl,MW_c,MW_a,Keq_ll,eleliquid,flux_sl,phi_salt,sat_ratio, &
1327 molality0,jsalt_present,kappa_nonelectro)
1331 call MESA_convergence_criterion(ibin,iconverge_mass,iconverge_flux,idissolved, &
1332 aer,electrolyte,mass_dry_a,mass_dry_salt,mw_electrolyte,mw_aer_mac, &
1333 flux_sl,phi_salt,rtol_mesa)
1335 if(iconverge_mass .eq. mYES)then
1336 mosaic_vars_aa%iter_MESA(ibin) = mosaic_vars_aa%iter_MESA(ibin) + itdum
1337 mosaic_vars_aa%niter_MESA = mosaic_vars_aa%niter_MESA + float(itdum)
1338 mosaic_vars_aa%niter_MESA_max = max( mosaic_vars_aa%niter_MESA_max, itdum)
1339 jaerosolstate(ibin) = all_solid
1340 call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent, &
1342 jhyst_leg(ibin) = jhyst_lo
1343 growth_factor(ibin) = 1.0
1345 elseif(iconverge_flux .eq. mYES)then
1346 mosaic_vars_aa%iter_MESA(ibin) = mosaic_vars_aa%iter_MESA(ibin) + itdum
1347 mosaic_vars_aa%niter_MESA = mosaic_vars_aa%niter_MESA + float(itdum)
1348 mosaic_vars_aa%niter_MESA_max = max( mosaic_vars_aa%niter_MESA_max, itdum)
1349 jaerosolstate(ibin) = mixed
1350 vol_wet_a(ibin) = vol_dry_a(ibin) + water_a(ibin)*1.e-3 ! cc(aer)/cc(air) or m^3/m^3(air)
1351 growth_factor(ibin) = mass_wet_a(ibin)/mass_dry_a(ibin) ! mass growth factor
1353 if(idissolved .eq. myes)then
1354 jaerosolstate(ibin) = all_liquid
1355 ! jhyst_leg(ibin) = jhyst_up ! ! do this later (to avoid tripping kelvin iterations)
1357 jaerosolstate(ibin) = mixed
1358 jhyst_leg(ibin) = jhyst_lo
1361 ! calculate epercent(jsolid) composition in mixed-phase aerosol EFFI
1364 !! do je = 1, nelectrolyte
1365 !! electrolyte(je,jp,ibin) = max(0.d0,electrolyte(je,jp,ibin)) ! remove -ve
1366 !! sum_dum = sum_dum + electrolyte(je,jp,ibin)
1368 !! electrolyte_sum(jp,ibin) = sum_dum
1369 !! if(sum_dum .eq. 0.)sum_dum = 1.0
1370 !! do je = 1, nelectrolyte
1371 !! epercent(je,jp,ibin) = 100.*electrolyte(je,jp,ibin)/sum_dum
1377 ! calculate hsalt(js) ! time step
1382 phi_prod = phi_salt(js) * phi_salt_old(js)
1384 if(itdum .gt. 1 .and. phi_prod .gt. 0.0)then
1385 phi_bar(js) = (abs(phi_salt(js))-abs(phi_salt_old(js)))/ &
1388 phi_bar(js) = 0.0 ! oscillating, or phi_salt and/or phi_salt_old may be zero
1391 if(phi_bar(js) .lt. 0.0)then ! good. phi getting lower. maybe able to take bigger alphas
1392 phi_bar(js) = max(phi_bar(js), -10.0d0)
1393 alpha_fac = 3.0*exp(phi_bar(js))
1394 alpha_salt(js) = min(alpha_fac*abs(phi_salt(js)), 0.9d0)
1395 elseif(phi_bar(js) .gt. 0.0)then ! bad - phi is getting bigger. so be conservative with alpha
1396 alpha_salt(js) = min(abs(phi_salt(js)), 0.5d0)
1397 else ! very bad - phi is oscillating. be very conservative
1398 alpha_salt(js) = min(abs(phi_salt(js))/3.0d0, 0.5d0)
1401 ! alpha_salt(js) = max(alpha_salt(js), 0.01)
1403 phi_salt_old(js) = phi_salt(js) ! update old array
1406 if(flux_sl(js) .gt. 0.)then
1408 tau_p(js) = eleliquid(js)/flux_sl(js) ! precipitation time scale
1409 if(tau_p(js) .eq. 0.0)then
1414 hsalt(js) = alpha_salt(js)*tau_p(js)
1417 elseif(flux_sl(js) .lt. 0.)then
1419 tau_p(js) = -eleliquid(js)/flux_sl(js) ! precipitation time scale
1420 tau_d(js) = -electrolyte(js,jsolid,ibin)/flux_sl(js) ! dissolution time scale
1421 if(tau_p(js) .eq. 0.0)then
1422 hsalt(js) = alpha_salt(js)*tau_d(js)
1424 hsalt(js) = alpha_salt(js)*min(tau_p(js),tau_d(js))
1433 hsalt_min = min(hsalt(js), hsalt_min)
1437 !---------------------------------
1439 ! integrate electrolyte(solid)
1441 electrolyte(js,jsolid,ibin) = ( &
1442 (electrolyte(js,jsolid,ibin)) + &
1443 (hsalt(js)) * (flux_sl(js)) )
1447 ! compute aer(solid) from electrolyte(solid)
1448 call electrolytes_to_ions(jsolid,ibin,aer,electrolyte)
1451 ! compute new electrolyte(liquid) from mass balance
1453 aer(iaer,jliquid,ibin) = ( (aer(iaer,jtotal,ibin)) - &
1454 (aer(iaer,jsolid,ibin)) )
1457 !---------------------------------
1461 500 continue ! end time continuation loop
1462 !--------------------------------------------------------------------
1463 mosaic_vars_aa%jMESA_fail = mosaic_vars_aa%jMESA_fail + 1
1464 mosaic_vars_aa%iter_MESA(ibin) = mosaic_vars_aa%iter_MESA(ibin) + itdum
1465 mosaic_vars_aa%niter_MESA = mosaic_vars_aa%niter_MESA + float(itdum)
1466 jaerosolstate(ibin) = mixed
1467 jhyst_leg(ibin) = jhyst_lo
1468 mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3 ! g/cc(air)
1469 vol_wet_a(ibin) = vol_dry_a(ibin) + water_a(ibin)*1.e-3 ! cc(aer)/cc(air) or m^3/m^3(air)
1470 growth_factor(ibin) = mass_wet_a(ibin)/mass_dry_a(ibin) ! mass growth factor
1473 end subroutine MESA_PTC
1477 !***********************************************************************
1478 ! part of MESA: calculates solid-liquid fluxes of soluble salts
1480 ! author: Rahul A. Zaveri
1482 !-----------------------------------------------------------------------
1483 subroutine MESA_flux_salt(ibin, jaerosolstate,jphase,aer,jhyst_leg,electrolyte, &
1484 epercent,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam, &
1485 log_gamZ,zc,za,gam_ratio,xeq_a,na_Ma,nc_Mc,xeq_c,mw_electrolyte,Keq_sl,MW_c,&
1486 MW_a,Keq_ll,eleliquid,flux_sl,phi_salt,sat_ratio,molality0,jsalt_present, &
1487 kappa_nonelectro ) ! TOUCH
1489 use module_data_mosaic_aero, only: r8,nbin_a_max,nelectrolyte,Ncation,naer, &!Parameters
1490 jliquid,nsalt,jsolid,Nanion,nrxn_aer_sl,nrxn_aer_ll,nrxn_aer_sl, &!Parameter
1491 jna3hso4,ica_a,jcano3,jcacl2 !TBD
1496 integer, intent(in) :: ibin
1497 integer, intent(inout), dimension(nsalt) :: jsalt_present
1498 integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg
1500 real(r8), intent(in) :: aH2O
1501 real(r8), intent(inout), dimension(nsalt) :: flux_sl,phi_salt,sat_ratio
1502 real(r8), intent(in), dimension(Ncation) :: zc,MW_c
1503 real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
1504 real(r8), intent(in), dimension(Nanion) :: za,MW_a
1505 real(r8), intent(inout), dimension(Nanion) :: xeq_a,na_Ma
1506 real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_dry_a,gam_ratio
1507 real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,water_a
1508 real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
1509 real(r8), intent(inout), dimension(nelectrolyte) :: eleliquid
1510 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
1511 real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1512 real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
1513 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
1514 real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
1515 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1516 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1517 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1518 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte,epercent
1519 real(r8), intent(in), dimension(naer) :: kappa_nonelectro
1523 real(r8) :: XT, calcium, sum_salt, sum_dum !**BALLI XT should it be subr arg??
1524 real(r8), dimension(nsalt) :: frac_salt_liq,frac_salt_solid
1527 ! compute activities and water content
1528 call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma, &
1529 nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)
1530 call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,electrolyte, &
1531 activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,log_gamZ, &
1532 gam_ratio,Keq_ll,molality0,kappa_nonelectro)
1533 activity(jna3hso4,ibin) = 0.0
1535 if(water_a(ibin) .le. 0.0)then
1543 call MESA_estimate_eleliquid(ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,nc_Mc, &
1544 xeq_c,mw_electrolyte,MW_c,MW_a,eleliquid)
1546 calcium = aer(ica_a,jliquid,ibin)
1550 !! EFFI calculate percent composition
1552 do je = 1, nelectrolyte
1553 sum_dum = sum_dum + electrolyte(je,jliquid,ibin)
1556 if(sum_dum .eq. 0.)sum_dum = 1.0
1558 do je = 1, nelectrolyte
1559 epercent(je,jliquid,ibin) = 100.*electrolyte(je,jliquid,ibin)/sum_dum
1565 ! calculate % electrolyte composition in the solid and liquid phases
1568 sum_salt = sum_salt + electrolyte(js,jsolid,ibin)
1571 if(sum_salt .eq. 0.0)sum_salt = 1.0
1573 frac_salt_solid(js) = electrolyte(js,jsolid,ibin)/sum_salt
1574 frac_salt_liq(js) = epercent(js,jliquid,ibin)/100.
1577 ! compute salt fluxes
1578 do js = 1, nsalt ! soluble solid salts
1580 ! compute new saturation ratio
1581 sat_ratio(js) = activity(js,ibin)/Keq_sl(js)
1582 ! compute relative driving force
1583 phi_salt(js) = (sat_ratio(js) - 1.0)/max(sat_ratio(js),1.0d0)
1585 ! check if too little solid-phase salt is trying to dissolve
1586 if(sat_ratio(js) .lt. 1.00 .and. &
1587 frac_salt_solid(js) .lt. 0.01 .and. &
1588 frac_salt_solid(js) .gt. 0.0)then
1589 call MESA_dissolve_small_salt(ibin,js,aer,electrolyte)
1590 call MESA_estimate_eleliquid(ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma, &
1591 nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a,eleliquid)
1592 sat_ratio(js) = activity(js,ibin)/Keq_sl(js)
1596 flux_sl(js) = sat_ratio(js) - 1.0
1598 ! apply Heaviside function
1599 if( (sat_ratio(js) .lt. 1.0 .and. &
1600 electrolyte(js,jsolid,ibin) .eq. 0.0) .or. &
1601 (calcium .gt. 0.0 .and. frac_salt_liq(js).lt.0.01).or. &
1602 (calcium .gt. 0.0 .and. jsalt_present(js).eq.0) )then
1610 ! force cacl2 and cano3 fluxes to zero
1611 sat_ratio(jcano3) = 1.0
1612 phi_salt(jcano3) = 0.0
1613 flux_sl(jcano3) = 0.0
1615 sat_ratio(jcacl2) = 1.0
1616 phi_salt(jcacl2) = 0.0
1617 flux_sl(jcacl2) = 0.0
1621 end subroutine MESA_flux_salt
1623 !***********************************************************************
1624 ! computes activities
1626 ! author: Rahul A. Zaveri
1628 !-----------------------------------------------------------------------
1629 subroutine compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg, &
1630 electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,&
1631 log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)
1633 use module_data_mosaic_aero, only: r8,nbin_a_max,nelectrolyte,Ncation,naer, &
1634 jliquid,Nanion,nrxn_aer_ll, &
1635 iso4_a,ja_so4,ja_hso4,ino3_a,ja_no3,icl_a,ja_cl,imsa_a,ja_msa,ica_a,jc_ca,&
1636 inh4_a,jc_nh4,ina_a,jc_na,jc_h,jhcl,jhno3,jcacl2,jcano3,jnacl,jnano3, &
1637 jna2so4,jnh4so4,jnh4cl,jnh4no3,jlvcite,jnh4hso4,jnh4msa,jna3hso4,jnahso4, &
1638 jnamsa,jcamsa2,jh2so4,jhhso4,jmsa !TBD
1643 integer, intent(in) :: ibin
1644 integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg
1646 real(r8), intent(in) :: aH2O
1647 real(r8), intent(in), dimension(nbin_a_max) :: num_a
1648 real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_a,mass_soluble_a
1649 real(r8), intent(inout), dimension(nbin_a_max) :: water_a,gam_ratio
1650 real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1651 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
1652 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
1653 real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
1654 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1655 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1656 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
1657 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1658 real(r8), intent(in), dimension(naer) :: kappa_nonelectro
1661 real(r8), dimension(nelectrolyte) :: log_gam
1663 real(r8) :: XT, xmol(Nelectrolyte), sum_elec, dumK, c_bal, a_c !BALLI** should xt be subr arg??
1664 real(r8) :: quad, aq, bq, cq, xq, dum, mSULF
1665 !real(r8) :: aerosol_water ! mosaic function
1668 water_a(ibin) = aerosol_water(jliquid,ibin,jaerosolstate,jphase, &
1669 jhyst_leg,electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O, &
1670 molality0) ! Kg/m^3(air)
1671 if(water_a(ibin) .eq. 0.0)return
1674 call calculate_XT(ibin,jliquid,XT,aer)
1677 if(XT.ge.2.0 .or. XT.lt.0.)then ! changed .gt. to .ge. RAZ 4/16/2014
1678 ! SULFATE POOR: fully dissociated electrolytes
1681 ! anion molalities (mol/kg water)
1682 ma(ja_so4,ibin) = 1.e-9*aer(iso4_a,jliquid,ibin)/water_a(ibin)
1683 ma(ja_hso4,ibin) = 0.0
1684 ma(ja_no3,ibin) = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
1685 ma(ja_cl,ibin) = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)
1686 ma(ja_msa,ibin) = 1.e-9*aer(imsa_a,jliquid,ibin)/water_a(ibin)
1688 ! cation molalities (mol/kg water)
1689 mc(jc_ca,ibin) = 1.e-9*aer(ica_a, jliquid,ibin)/water_a(ibin)
1690 mc(jc_nh4,ibin) = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
1691 mc(jc_na,ibin) = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)
1693 (2.*ma(ja_so4,ibin)+ &
1696 ma(ja_msa,ibin)) - &
1697 (2.*mc(jc_ca,ibin) + &
1701 mc(jc_h,ibin) = 0.5*( (a_c) + &
1702 (sqrt(a_c**2 + 4.*Keq_ll(3))) )
1704 if(mc(jc_h,ibin) .le. 0.0)then ! changed .eq. to .le. RAZ 4/16/2014
1705 mc(jc_h,ibin) = 1.e-10
1712 sum_elec = 2.*electrolyte(jnh4no3,jp,ibin) + &
1713 2.*electrolyte(jnh4cl,jp,ibin) + &
1714 3.*electrolyte(jnh4so4,jp,ibin) + &
1715 3.*electrolyte(jna2so4,jp,ibin) + &
1716 2.*electrolyte(jnano3,jp,ibin) + &
1717 2.*electrolyte(jnacl,jp,ibin) + &
1718 3.*electrolyte(jcano3,jp,ibin) + &
1719 3.*electrolyte(jcacl2,jp,ibin) + &
1720 2.*electrolyte(jhno3,jp,ibin) + &
1721 2.*electrolyte(jhcl,jp,ibin)
1723 if(sum_elec .eq. 0.0)then
1724 do jA = 1, nelectrolyte
1731 ! ionic mole fractions
1732 xmol(jnh4no3) = 2.*electrolyte(jnh4no3,jp,ibin)/sum_elec
1733 xmol(jnh4cl) = 2.*electrolyte(jnh4cl,jp,ibin) /sum_elec
1734 xmol(jnh4so4) = 3.*electrolyte(jnh4so4,jp,ibin)/sum_elec
1735 xmol(jna2so4) = 3.*electrolyte(jna2so4,jp,ibin)/sum_elec
1736 xmol(jnano3) = 2.*electrolyte(jnano3,jp,ibin) /sum_elec
1737 xmol(jnacl) = 2.*electrolyte(jnacl,jp,ibin) /sum_elec
1738 xmol(jcano3) = 3.*electrolyte(jcano3,jp,ibin) /sum_elec
1739 xmol(jcacl2) = 3.*electrolyte(jcacl2,jp,ibin) /sum_elec
1740 xmol(jhno3) = 2.*electrolyte(jhno3,jp,ibin) /sum_elec
1741 xmol(jhcl) = 2.*electrolyte(jhcl,jp,ibin) /sum_elec
1745 if(xmol(jA).gt.0.0)then
1746 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1747 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1748 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1749 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1750 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1751 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1752 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1753 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1754 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1755 xmol(jhcl) *log_gamZ(jA,jhcl)
1756 gam(jA,ibin) = 10.**log_gam(jA)
1757 activity(jnh4so4,ibin) = mc(jc_nh4,ibin)**2 * ma(ja_so4,ibin) * &
1758 gam(jnh4so4,ibin)**3
1764 ! always calculate gam(jnh4no3), even if xmol(jnh4no3) = 0. this to calculate gam_ratio
1766 ! if(xmol(jA).gt.0.0)then
1767 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1768 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1769 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1770 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1771 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1772 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1773 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1774 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1775 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1776 xmol(jhcl) *log_gamZ(jA,jhcl)
1777 gam(jA,ibin) = 10.**log_gam(jA)
1778 activity(jnh4no3,ibin) = mc(jc_nh4,ibin) * ma(ja_no3,ibin) * &
1779 gam(jnh4no3,ibin)**2
1784 if(xmol(jA).gt.0.0)then
1785 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1786 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1787 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1788 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1789 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1790 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1791 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1792 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1793 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1794 xmol(jhcl) *log_gamZ(jA,jhcl)
1795 gam(jA,ibin) = 10.**log_gam(jA)
1796 activity(jnh4cl,ibin) = mc(jc_nh4,ibin) * ma(ja_cl,ibin) * &
1802 if(xmol(jA).gt.0.0)then
1803 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1804 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1805 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1806 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1807 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1808 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1809 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1810 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1811 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1812 xmol(jhcl) *log_gamZ(jA,jhcl)
1813 gam(jA,ibin) = 10.**log_gam(jA)
1814 activity(jna2so4,ibin) = mc(jc_na,ibin)**2 * ma(ja_so4,ibin) * &
1815 gam(jna2so4,ibin)**3
1820 if(xmol(jA).gt.0.0)then
1821 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1822 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1823 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1824 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1825 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1826 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1827 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1828 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1829 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1830 xmol(jhcl) *log_gamZ(jA,jhcl)
1831 gam(jA,ibin) = 10.**log_gam(jA)
1832 activity(jnano3,ibin) = mc(jc_na,ibin) * ma(ja_no3,ibin) * &
1839 if(xmol(jA).gt.0.0)then
1840 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1841 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1842 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1843 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1844 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1845 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1846 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1847 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1848 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1849 xmol(jhcl) *log_gamZ(jA,jhcl)
1850 gam(jA,ibin) = 10.**log_gam(jA)
1851 activity(jnacl,ibin) = mc(jc_na,ibin) * ma(ja_cl,ibin) * &
1858 !c if(xmol(jA).gt.0.0)then
1859 !c gam(jA,ibin) = 1.0
1860 !c activity(jcano3,ibin) = 1.0
1866 !c if(xmol(jA).gt.0.0)then
1867 !c gam(jA,ibin) = 1.0
1868 !c activity(jcacl2,ibin) = 1.0
1872 if(xmol(jA).gt.0.0)then
1873 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1874 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1875 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1876 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1877 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1878 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1879 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1880 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1881 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1882 xmol(jhcl) *log_gamZ(jA,jhcl)
1883 gam(jA,ibin) = 10.**log_gam(jA)
1884 activity(jcano3,ibin) = mc(jc_ca,ibin) * ma(ja_no3,ibin)**2 * &
1891 if(xmol(jA).gt.0.0)then
1892 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1893 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1894 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1895 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1896 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1897 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1898 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1899 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1900 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1901 xmol(jhcl) *log_gamZ(jA,jhcl)
1902 gam(jA,ibin) = 10.**log_gam(jA)
1903 activity(jcacl2,ibin) = mc(jc_ca,ibin) * ma(ja_cl,ibin)**2 * &
1909 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1910 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1911 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1912 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1913 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1914 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1915 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1916 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1917 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1918 xmol(jhcl) *log_gamZ(jA,jhcl)
1919 gam(jA,ibin) = 10.**log_gam(jA)
1920 activity(jhno3,ibin) = mc(jc_h,ibin) * ma(ja_no3,ibin) * &
1925 log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) + &
1926 xmol(jnh4cl) *log_gamZ(jA,jnh4cl) + &
1927 xmol(jnh4so4)*log_gamZ(jA,jnh4so4) + &
1928 xmol(jna2so4)*log_gamZ(jA,jna2so4) + &
1929 xmol(jnano3) *log_gamZ(jA,jnano3) + &
1930 xmol(jnacl) *log_gamZ(jA,jnacl) + &
1931 xmol(jcano3) *log_gamZ(jA,jcano3) + &
1932 xmol(jcacl2) *log_gamZ(jA,jcacl2) + &
1933 xmol(jhno3) *log_gamZ(jA,jhno3) + &
1934 xmol(jhcl) *log_gamZ(jA,jhcl)
1935 gam(jA,ibin) = 10.**log_gam(jA)
1936 activity(jhcl,ibin) = mc(jc_h,ibin) * ma(ja_cl,ibin) * &
1940 10 gam(jlvcite,ibin) = 1.0
1942 gam(jnh4hso4,ibin)= 1.0
1944 gam(jnh4msa,ibin) = 1.0
1946 gam(jna3hso4,ibin) = 1.0
1948 gam(jnahso4,ibin) = 1.0
1950 gam(jnamsa,ibin) = 1.0
1952 gam(jcamsa2,ibin) = 1.0
1954 activity(jlvcite,ibin) = 0.0
1956 activity(jnh4hso4,ibin)= 0.0
1958 activity(jnh4msa,ibin) = mc(jc_nh4,ibin) * ma(ja_msa,ibin) * &
1959 gam(jnh4msa,ibin)**2
1961 activity(jna3hso4,ibin)= 0.0
1963 activity(jnahso4,ibin) = 0.0
1965 activity(jnamsa,ibin) = mc(jc_na,ibin) * ma(ja_msa,ibin) * &
1968 activity(jcamsa2,ibin) = mc(jc_ca,ibin) * ma(ja_msa,ibin)**2 * &
1969 gam(jcamsa2,ibin)**3
1971 gam_ratio(ibin) = gam(jnh4no3,ibin)**2/gam(jhno3,ibin)**2
1975 ! SULFATE-RICH: solve for SO4= and HSO4- ions
1979 sum_elec = 3.*electrolyte(jh2so4,jp,ibin) + &
1980 2.*electrolyte(jnh4hso4,jp,ibin) + &
1981 5.*electrolyte(jlvcite,jp,ibin) + &
1982 3.*electrolyte(jnh4so4,jp,ibin) + &
1983 2.*electrolyte(jnahso4,jp,ibin) + &
1984 5.*electrolyte(jna3hso4,jp,ibin) + &
1985 3.*electrolyte(jna2so4,jp,ibin) + &
1986 2.*electrolyte(jhno3,jp,ibin) + &
1987 2.*electrolyte(jhcl,jp,ibin)
1990 if(sum_elec .eq. 0.0)then
1991 do jA = 1, nelectrolyte
1998 xmol(jh2so4) = 3.*electrolyte(jh2so4,jp,ibin)/sum_elec
1999 xmol(jnh4hso4)= 2.*electrolyte(jnh4hso4,jp,ibin)/sum_elec
2000 xmol(jlvcite) = 5.*electrolyte(jlvcite,jp,ibin)/sum_elec
2001 xmol(jnh4so4) = 3.*electrolyte(jnh4so4,jp,ibin)/sum_elec
2002 xmol(jnahso4) = 2.*electrolyte(jnahso4,jp,ibin)/sum_elec
2003 xmol(jna3hso4)= 5.*electrolyte(jna3hso4,jp,ibin)/sum_elec
2004 xmol(jna2so4) = 3.*electrolyte(jna2so4,jp,ibin)/sum_elec
2005 xmol(jhno3) = 2.*electrolyte(jhno3,jp,ibin)/sum_elec
2006 xmol(jhcl) = 2.*electrolyte(jhcl,jp,ibin)/sum_elec
2011 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2012 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2013 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2014 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2015 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2016 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2017 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2018 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2019 xmol(jhcl) *log_gamZ(jA,jhcl)
2020 gam(jA,ibin) = 10.**log_gam(jA)
2025 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2026 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2027 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2028 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2029 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2030 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2031 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2032 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2033 xmol(jhcl) *log_gamZ(jA,jhcl)
2034 gam(jA,ibin) = 10.**log_gam(jA)
2039 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2040 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2041 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2042 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2043 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2044 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2045 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2046 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2047 xmol(jhcl) *log_gamZ(jA,jhcl)
2048 gam(jA,ibin) = 10.**log_gam(jA)
2053 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2054 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2055 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2056 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2057 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2058 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2059 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2060 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2061 xmol(jhcl) *log_gamZ(jA,jhcl)
2062 gam(jA,ibin) = 10.**log_gam(jA)
2067 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2068 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2069 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2070 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2071 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2072 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2073 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2074 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2075 xmol(jhcl) *log_gamZ(jA,jhcl)
2076 gam(jA,ibin) = 10.**log_gam(jA)
2081 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2082 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2083 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2084 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2085 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2086 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2087 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2088 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2089 xmol(jhcl) *log_gamZ(jA,jhcl)
2090 gam(jA,ibin) = 10.**log_gam(jA)
2095 ! log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) +
2096 ! & xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+
2097 ! & xmol(jlvcite) *log_gamZ(jA,jlvcite) +
2098 ! & xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +
2099 ! & xmol(jnahso4) *log_gamZ(jA,jnahso4) +
2100 ! & xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+
2101 ! & xmol(jna2so4) *log_gamZ(jA,jna2so4) +
2102 ! & xmol(jhno3) *log_gamZ(jA,jhno3) +
2103 ! & xmol(jhcl) *log_gamZ(jA,jhcl)
2104 ! gam(jA,ibin) = 10.**log_gam(jA)
2110 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2111 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2112 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2113 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2114 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2115 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2116 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2117 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2118 xmol(jhcl) *log_gamZ(jA,jhcl)
2119 gam(jA,ibin) = 10.**log_gam(jA)
2124 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2125 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2126 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2127 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2128 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2129 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2130 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2131 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2132 xmol(jhcl) *log_gamZ(jA,jhcl)
2133 gam(jA,ibin) = 10.**log_gam(jA)
2138 log_gam(jA) = xmol(jh2so4) *log_gamZ(jA,jh2so4) + &
2139 xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+ &
2140 xmol(jlvcite) *log_gamZ(jA,jlvcite) + &
2141 xmol(jnh4so4) *log_gamZ(jA,jnh4so4) + &
2142 xmol(jnahso4) *log_gamZ(jA,jnahso4) + &
2143 xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+ &
2144 xmol(jna2so4) *log_gamZ(jA,jna2so4) + &
2145 xmol(jhno3) *log_gamZ(jA,jhno3) + &
2146 xmol(jhcl) *log_gamZ(jA,jhcl)
2147 gam(jA,ibin) = 10.**log_gam(jA)
2150 20 gam(jnh4no3,ibin) = 1.0
2151 gam(jnh4cl,ibin) = 1.0
2152 gam(jnano3,ibin) = 1.0
2153 gam(jnacl,ibin) = 1.0
2154 gam(jcano3,ibin) = 1.0
2155 gam(jcacl2,ibin) = 1.0
2157 gam(jnh4msa,ibin) = 1.0
2158 gam(jnamsa,ibin) = 1.0
2159 gam(jcamsa2,ibin) = 1.0
2163 ! compute equilibrium pH
2164 ! cation molalities (mol/kg water)
2165 mc(jc_ca,ibin) = 1.e-9*aer(ica_a,jliquid,ibin)/water_a(ibin)
2166 mc(jc_nh4,ibin) = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
2167 mc(jc_na,ibin) = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)
2169 ! anion molalities (mol/kg water)
2170 mSULF = 1.e-9*aer(iso4_a,jliquid,ibin)/water_a(ibin)
2171 ma(ja_hso4,ibin) = 0.0
2172 ma(ja_so4,ibin) = 0.0
2173 ma(ja_no3,ibin) = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
2174 ma(ja_cl,ibin) = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)
2175 ma(ja_msa,ibin) = 1.e-9*aer(imsa_a,jliquid,ibin)/water_a(ibin)
2177 gam_ratio(ibin) = gam(jnh4hso4,ibin)**2/gam(jhhso4,ibin)**2
2178 dumK = Keq_ll(1)*gam(jhhso4,ibin)**2/gam(jh2so4,ibin)**3
2180 c_bal = mc(jc_nh4,ibin) + mc(jc_na,ibin) + 2.*mc(jc_ca,ibin) &
2181 - ma(ja_no3,ibin) - ma(ja_cl,ibin) - mSULF - ma(ja_msa,ibin)
2185 cq = dumK*(c_bal - mSULF)
2188 !--quadratic solution
2190 xq = 4.*(1./bq)*(cq/bq)
2195 if(abs(xq) .lt. 1.e-6)then
2196 dum = xq*(0.5 + xq*(0.125 + xq*0.0625))
2197 quad = (-0.5*bq/aq)*dum
2198 if(quad .lt. 0.)then
2199 quad = -bq/aq - quad
2202 quad = 0.5*(-bq+sqrt(bq*bq - 4.*cq))
2204 !--end of quadratic solution
2206 mc(jc_h,ibin) = max(quad, 1.d-7)
2207 ma(ja_so4,ibin) = mSULF*dumK/(mc(jc_h,ibin) + dumK)
2208 ma(ja_hso4,ibin)= mSULF - ma(ja_so4,ibin)
2210 activity(jcamsa2,ibin) = mc(jc_ca,ibin) * ma(ja_msa,ibin)**2 * &
2211 gam(jcamsa2,ibin)**3
2213 activity(jnh4so4,ibin) = mc(jc_nh4,ibin)**2 * ma(ja_so4,ibin) * &
2214 gam(jnh4so4,ibin)**3
2216 activity(jlvcite,ibin) = mc(jc_nh4,ibin)**3 * ma(ja_hso4,ibin) * &
2217 ma(ja_so4,ibin) * gam(jlvcite,ibin)**5
2219 activity(jnh4hso4,ibin)= mc(jc_nh4,ibin) * ma(ja_hso4,ibin) * &
2220 gam(jnh4hso4,ibin)**2
2222 activity(jnh4msa,ibin) = mc(jc_nh4,ibin) * ma(ja_msa,ibin) * &
2223 gam(jnh4msa,ibin)**2
2225 activity(jna2so4,ibin) = mc(jc_na,ibin)**2 * ma(ja_so4,ibin) * &
2226 gam(jna2so4,ibin)**3
2228 activity(jnahso4,ibin) = mc(jc_na,ibin) * ma(ja_hso4,ibin) * &
2229 gam(jnahso4,ibin)**2
2231 activity(jnamsa,ibin) = mc(jc_na,ibin) * ma(ja_msa,ibin) * &
2234 ! activity(jna3hso4,ibin)= mc(jc_na,ibin)**3 * ma(ja_hso4,ibin) *
2235 ! & ma(ja_so4,ibin) * gam(jna3hso4,ibin)**5
2237 activity(jna3hso4,ibin)= 0.0
2239 activity(jhno3,ibin) = mc(jc_h,ibin) * ma(ja_no3,ibin) * &
2242 activity(jhcl,ibin) = mc(jc_h,ibin) * ma(ja_cl,ibin) * &
2245 activity(jmsa,ibin) = mc(jc_h,ibin) * ma(ja_msa,ibin) * &
2249 ! sulfate-poor species
2250 activity(jnh4no3,ibin) = 0.0
2252 activity(jnh4cl,ibin) = 0.0
2254 activity(jnano3,ibin) = 0.0
2256 activity(jnacl,ibin) = 0.0
2258 activity(jcano3,ibin) = 0.0
2260 activity(jcacl2,ibin) = 0.0
2265 end subroutine compute_activities
2269 !***********************************************************************
2270 ! part of MESA: checks MESA convergence
2272 ! author: Rahul A. Zaveri
2274 !-----------------------------------------------------------------------
2275 subroutine MESA_convergence_criterion(ibin,iconverge_mass,iconverge_flux, &
2276 idissolved,aer,electrolyte,mass_dry_a,mass_dry_salt, &
2277 mw_electrolyte,mw_aer_mac,flux_sl,phi_salt,rtol_mesa) ! TOUCH
2279 use module_data_mosaic_aero, only: r8,nbin_a_max,naer,nelectrolyte,nsalt, &!Parameters
2280 jsolid, mYES, xhyst_up_crustal_thresh, &!Parameters
2281 mno,ioin_a,jcaso4,jcaco3 !TBD
2286 integer, intent(in) :: ibin
2287 integer, intent(inout) :: iconverge_mass, iconverge_flux, idissolved
2288 real(r8), intent(in) :: rtol_mesa
2289 real(r8), intent(inout), dimension(nsalt) :: flux_sl,phi_salt
2290 real(r8), intent(in), dimension(nbin_a_max) :: mass_dry_a
2291 real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_salt
2292 real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
2293 real(r8), intent(in), dimension(naer) :: mw_aer_mac
2294 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2295 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2297 integer je, js, iaer
2298 real(r8) :: mass_solid, mass_solid_salt,frac_solid, XT, H_ion, &
2299 crustal_solids, sumflux
2302 idissolved = mno ! default = not completely dissolved
2304 ! check mass convergence
2305 iconverge_mass = mNO ! default value = no convergence
2307 ! call electrolytes_to_ions(jsolid,ibin,aer,electrolyte)
2310 ! mass_solid = mass_solid +
2311 ! & aer(iaer,jsolid,ibin)*mw_aer_mac(iaer)*1.e-15 ! g/cc(air)
2314 mass_solid_salt = 0.0
2316 mass_solid_salt = mass_solid_salt + &
2317 electrolyte(je,jsolid,ibin)*mw_electrolyte(je)*1.e-15 ! g/cc(air)
2322 ! frac_solid = mass_solid/mass_dry_a(ibin)
2325 if(mass_dry_salt(ibin) .le. 0.0)then
2328 frac_solid = mass_solid_salt/mass_dry_salt(ibin)
2332 if(frac_solid .ge. 0.98)then
2333 iconverge_mass = mYES
2339 ! check relative driving force convergence
2340 iconverge_flux = mYES
2342 if(abs(phi_salt(js)).gt. rtol_mesa)then
2343 iconverge_flux = mNO
2350 ! check if all the fluxes are zero
2354 sumflux = sumflux + abs(flux_sl(js))
2357 ! crustal_solids = electrolyte(jcaco3,jsolid,ibin) + &
2358 ! electrolyte(jcaso4,jsolid,ibin) + &
2359 ! aer(ioin_a,jsolid,ibin)
2360 crustal_solids = electrolyte(jcaco3,jsolid,ibin)*mw_electrolyte(jcaco3) + &
2361 electrolyte(jcaso4,jsolid,ibin)*mw_electrolyte(jcaso4) + &
2362 aer(ioin_a,jsolid,ibin)*mw_aer_mac(ioin_a)
2364 ! if(sumflux .eq. 0.0 .and. crustal_solids .eq. 0.0)then
2365 if ( sumflux .eq. 0.0 .and. &
2366 crustal_solids .le. xhyst_up_crustal_thresh*(mass_dry_a(ibin)*1.0e15) ) then
2367 ! crustal_solids is ng/m^3, mass_dry_a is g/cm^3
2374 end subroutine MESA_convergence_criterion
2378 !***********************************************************************
2379 ! computes ions from electrolytes
2381 ! author: Rahul A. Zaveri
2383 !-----------------------------------------------------------------------
2384 subroutine electrolytes_to_ions(jp,ibin,aer,electrolyte)
2386 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &!Parameters
2387 jh2so4,jnh4hso4,jlvcite,jnh4so4,jnahso4,jna3hso4,jna2so4,jcaso4,iso4_a, &!TBD
2388 jhno3,jnh4no3,jcano3,jnano3,ino3_a,jhcl,jnh4cl,jcacl2,jnacl,icl_a,jmsa, &!TBD
2389 jcamsa2,jnamsa,jnh4msa,imsa_a,jcaco3,ico3_a,ica_a,ina_a,inh4_a !TBD
2393 integer, intent(in) :: jp, ibin
2394 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2395 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2400 aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
2401 electrolyte(jna2so4,jp,ibin) + &
2402 2.*electrolyte(jna3hso4,jp,ibin)+ &
2403 electrolyte(jnahso4,jp,ibin) + &
2404 electrolyte(jnh4so4,jp,ibin) + &
2405 2.*electrolyte(jlvcite,jp,ibin) + &
2406 electrolyte(jnh4hso4,jp,ibin)+ &
2407 electrolyte(jh2so4,jp,ibin)
2409 aer(ino3_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
2410 2.*electrolyte(jcano3,jp,ibin) + &
2411 electrolyte(jnh4no3,jp,ibin) + &
2412 electrolyte(jhno3,jp,ibin)
2414 aer(icl_a,jp,ibin) = electrolyte(jnacl,jp,ibin) + &
2415 2.*electrolyte(jcacl2,jp,ibin) + &
2416 electrolyte(jnh4cl,jp,ibin) + &
2417 electrolyte(jhcl,jp,ibin)
2419 aer(imsa_a,jp,ibin) = electrolyte(jnh4msa,jp,ibin) + &
2420 electrolyte(jnamsa,jp,ibin) + &
2421 2.*electrolyte(jcamsa2,jp,ibin) + &
2422 electrolyte(jmsa,jp,ibin)
2424 aer(ico3_a,jp,ibin) = electrolyte(jcaco3,jp,ibin)
2426 aer(ica_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
2427 electrolyte(jcano3,jp,ibin) + &
2428 electrolyte(jcacl2,jp,ibin) + &
2429 electrolyte(jcaco3,jp,ibin) + &
2430 electrolyte(jcamsa2,jp,ibin)
2432 aer(ina_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
2433 electrolyte(jnacl,jp,ibin) + &
2434 2.*electrolyte(jna2so4,jp,ibin) + &
2435 3.*electrolyte(jna3hso4,jp,ibin)+ &
2436 electrolyte(jnahso4,jp,ibin) + &
2437 electrolyte(jnamsa,jp,ibin)
2439 aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) + &
2440 electrolyte(jnh4cl,jp,ibin) + &
2441 2.*electrolyte(jnh4so4,jp,ibin) + &
2442 3.*electrolyte(jlvcite,jp,ibin) + &
2443 electrolyte(jnh4hso4,jp,ibin)+ &
2444 electrolyte(jnh4msa,jp,ibin)
2448 end subroutine electrolytes_to_ions
2452 !***********************************************************************
2453 ! combinatorial method for computing electrolytes from ions
2456 ! - to be used for liquid-phase or total-phase only
2457 ! - transfers caso4 and caco3 from liquid to solid phase
2459 ! author: Rahul A. Zaveri (based on code provided by A.S. Wexler)
2461 !-----------------------------------------------------------------------
2462 subroutine ions_to_electrolytes(jp,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma, &
2463 nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)
2465 use module_data_mosaic_aero, only: r8,naer,nbin_a_max,nelectrolyte,ncation, &!Parameters
2466 nanion,jliquid,jsolid, &!Parameters
2467 ica_a,iso4_a,jcaso4,imsa_a,ina_a,inh4_a,ja_hso4,ja_so4,ino3_a,ja_no3, &
2468 icl_a,ja_cl,ja_msa,jc_ca,jc_na,jc_nh4,jc_h,jna2so4,jnahso4,jnamsa,jnano3, &
2469 jnacl,jnh4so4,jnh4hso4,jnh4msa,jnh4no3,jnh4cl,jcano3,jcacl2,jcamsa2, &
2470 jh2so4,jhno3,jhcl,jmsa,jlvcite,jna3hso4 !TBD
2474 integer, intent(in) :: ibin, jp
2475 real(r8), intent(inout) :: XT
2476 real(r8), intent(in), dimension(Ncation) :: zc,MW_c
2477 real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
2478 real(r8), intent(in), dimension(Nanion) :: za,MW_a
2479 real(r8), intent(inout), dimension(Nanion) :: xeq_a,na_Ma
2480 real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
2481 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2482 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2484 character(len=256) :: errmsg
2485 integer iaer, je, jc, ja, icase
2486 real(r8) :: store(naer), sum_dum, sum_naza, sum_nczc, sum_na_nh4, &
2487 f_nh4, f_na, xh, xb, xl, xs, cat_net, rem_nh4, rem_na
2488 real(r8) :: nc(ncation), na(nanion)
2493 if(jp .ne. jliquid)then
2494 call wrf_message(' jp must be jliquid')
2495 call wrf_message(' in ions_to_electrolytes sub')
2496 write(errmsg,*)' wrong jp = ', jp
2497 call wrf_error_fatal(trim(adjustl(errmsg)))
2501 ! remove negative concentrations, if any
2503 ! aer(iaer,jp,ibin) = max(0.0d0, aer(iaer,jp,ibin)) ! EFFI
2507 ! first transfer caso4 from liquid to solid phase (caco3 should not be present here)
2508 store(ica_a) = aer(ica_a, jp,ibin)
2509 store(iso4_a) = aer(iso4_a,jp,ibin)
2511 call form_caso4(store,jp,ibin,electrolyte)
2513 if(jp .eq. jliquid)then ! transfer caso4 from liquid to solid phase
2514 aer(ica_a,jliquid,ibin) = aer(ica_a,jliquid,ibin) - &
2515 electrolyte(jcaso4,jliquid,ibin)
2517 aer(iso4_a,jliquid,ibin)= aer(iso4_a,jliquid,ibin)- &
2518 electrolyte(jcaso4,jliquid,ibin)
2520 aer(ica_a,jsolid,ibin) = aer(ica_a,jsolid,ibin) + &
2521 electrolyte(jcaso4,jliquid,ibin)
2523 aer(iso4_a,jsolid,ibin) = aer(iso4_a,jsolid,ibin) + &
2524 electrolyte(jcaso4,jliquid,ibin)
2526 electrolyte(jcaso4,jsolid,ibin)=electrolyte(jcaso4,jsolid,ibin) &
2527 +electrolyte(jcaso4,jliquid,ibin)
2528 electrolyte(jcaso4,jliquid,ibin)= 0.0
2532 ! calculate sulfate ratio
2533 ! call calculate_XT(ibin,jp,XT,aer) ! EFFI
2535 if( (aer(iso4_a,jp,ibin)+aer(imsa_a,jp,ibin)) .gt.0.0)then
2536 XT = ( aer(inh4_a,jp,ibin) + &
2537 aer(ina_a,jp,ibin) + &
2538 2.*aer(ica_a,jp,ibin) )/ &
2539 (aer(iso4_a,jp,ibin)+0.5*aer(imsa_a,jp,ibin))
2547 ! if(XT.ge.1.9999 .or. XT.lt.0.)then ! commented out by RAZ 4/16/2014
2548 if(XT.ge.2.0 .or. XT.lt.0.)then ! Slightly different logic, consistent with that in compute_activities subr. RAZ 4/16/2014
2549 icase = 1 ! sulfate poor: near neutral (acidity is caused by HCl and/or HNO3)
2551 icase = 2 ! sulfate rich: acidic (acidity is caused by excess SO4)
2555 ! initialize to zero
2556 do je = 1, nelectrolyte
2557 electrolyte(je,jp,ibin) = 0.0
2561 !---------------------------------------------------------
2562 ! initialize moles of ions depending on the sulfate domain
2564 if(icase.eq.1)then ! XT >= 2 or XT < 0: SULFATE POOR (OR NO SULFATE) DOMAIN. RAZ 4/16/2014
2567 na(ja_so4) = aer(iso4_a,jp,ibin)
2568 na(ja_no3) = aer(ino3_a,jp,ibin)
2569 na(ja_cl) = aer(icl_a, jp,ibin)
2570 na(ja_msa) = aer(imsa_a,jp,ibin)
2572 nc(jc_ca) = aer(ica_a, jp,ibin)
2573 nc(jc_na) = aer(ina_a, jp,ibin)
2574 nc(jc_nh4) = aer(inh4_a,jp,ibin)
2577 (2.*na(ja_so4)+na(ja_no3)+na(ja_cl)+na(ja_msa)) - &
2578 (2.*nc(jc_ca) +nc(jc_nh4)+nc(jc_na)) )
2580 if(cat_net .lt. 0.0)then
2584 else ! cat_net must be 0.0 or positive
2591 ! now compute equivalent fractions
2594 sum_naza = sum_naza + na(ja)*za(ja)
2599 sum_nczc = sum_nczc + nc(jc)*zc(jc)
2602 if(sum_naza .eq. 0. .or. sum_nczc .eq. 0.)then ! it's ok. this may happen if the aerosol is assumed to be composed of hygroscopic SOA, POA, BC, OIN, but does not contain any inorganic electrolytes
2603 ! write(6,*)'ionic concentrations are zero in ibin', ibin ! commented out by RAZ 4/16/2014
2604 ! write(6,*)'sum_naza = ', sum_naza ! commented out by RAZ 4/16/2014
2605 ! write(6,*)'sum_nczc = ', sum_nczc ! commented out by RAZ 4/16/2014
2610 xeq_a(ja) = na(ja)*za(ja)/sum_naza
2614 xeq_c(jc) = nc(jc)*zc(jc)/sum_nczc
2617 na_Ma(ja_so4) = na(ja_so4) *MW_a(ja_so4)
2618 na_Ma(ja_no3) = na(ja_no3) *MW_a(ja_no3)
2619 na_Ma(ja_cl) = na(ja_cl) *MW_a(ja_cl)
2620 na_Ma(ja_msa) = na(ja_msa) *MW_a(ja_msa)
2621 na_Ma(ja_hso4)= na(ja_hso4)*MW_a(ja_hso4)
2623 nc_Mc(jc_ca) = nc(jc_ca) *MW_c(jc_ca)
2624 nc_Mc(jc_na) = nc(jc_na) *MW_c(jc_na)
2625 nc_Mc(jc_nh4) = nc(jc_nh4)*MW_c(jc_nh4)
2626 nc_Mc(jc_h) = nc(jc_h) *MW_c(jc_h)
2629 ! now compute electrolyte moles
2630 if(xeq_c(jc_na) .gt. 0. .and. xeq_a(ja_so4) .gt. 0.)then
2631 electrolyte(jna2so4,jp,ibin) = (xeq_c(jc_na) *na_Ma(ja_so4) + &
2632 xeq_a(ja_so4)*nc_Mc(jc_na))/ &
2633 mw_electrolyte(jna2so4)
2636 electrolyte(jnahso4,jp,ibin) = 0.0
2638 if(xeq_c(jc_na) .gt. 0. .and. xeq_a(ja_msa) .gt. 0.)then
2639 electrolyte(jnamsa,jp,ibin) = (xeq_c(jc_na) *na_Ma(ja_msa) + &
2640 xeq_a(ja_msa)*nc_Mc(jc_na))/ &
2641 mw_electrolyte(jnamsa)
2644 if(xeq_c(jc_na) .gt. 0. .and. xeq_a(ja_no3) .gt. 0.)then
2645 electrolyte(jnano3,jp,ibin) = (xeq_c(jc_na) *na_Ma(ja_no3) + &
2646 xeq_a(ja_no3)*nc_Mc(jc_na))/ &
2647 mw_electrolyte(jnano3)
2650 if(xeq_c(jc_na) .gt. 0. .and. xeq_a(ja_cl) .gt. 0.)then
2651 electrolyte(jnacl,jp,ibin) = (xeq_c(jc_na) *na_Ma(ja_cl) + &
2652 xeq_a(ja_cl) *nc_Mc(jc_na))/ &
2653 mw_electrolyte(jnacl)
2656 if(xeq_c(jc_nh4) .gt. 0. .and. xeq_a(ja_so4) .gt. 0.)then
2657 electrolyte(jnh4so4,jp,ibin) = (xeq_c(jc_nh4)*na_Ma(ja_so4) + &
2658 xeq_a(ja_so4)*nc_Mc(jc_nh4))/ &
2659 mw_electrolyte(jnh4so4)
2662 electrolyte(jnh4hso4,jp,ibin)= 0.0
2664 if(xeq_c(jc_nh4) .gt. 0. .and. xeq_a(ja_msa) .gt. 0.)then
2665 electrolyte(jnh4msa,jp,ibin) = (xeq_c(jc_nh4)*na_Ma(ja_msa) + &
2666 xeq_a(ja_msa)*nc_Mc(jc_nh4))/ &
2667 mw_electrolyte(jnh4msa)
2670 if(xeq_c(jc_nh4) .gt. 0. .and. xeq_a(ja_no3) .gt. 0.)then
2671 electrolyte(jnh4no3,jp,ibin) = (xeq_c(jc_nh4)*na_Ma(ja_no3) + &
2672 xeq_a(ja_no3)*nc_Mc(jc_nh4))/ &
2673 mw_electrolyte(jnh4no3)
2676 if(xeq_c(jc_nh4) .gt. 0. .and. xeq_a(ja_cl) .gt. 0.)then
2677 electrolyte(jnh4cl,jp,ibin) = (xeq_c(jc_nh4)*na_Ma(ja_cl) + &
2678 xeq_a(ja_cl) *nc_Mc(jc_nh4))/ &
2679 mw_electrolyte(jnh4cl)
2682 if(xeq_c(jc_ca) .gt. 0. .and. xeq_a(ja_no3) .gt. 0.0)then
2683 electrolyte(jcano3, jp,ibin) = (xeq_c(jc_ca) *na_Ma(ja_no3) + &
2684 xeq_a(ja_no3)*nc_Mc(jc_ca))/ &
2685 mw_electrolyte(jcano3)
2688 if(xeq_c(jc_ca) .gt. 0. .and. xeq_a(ja_cl) .gt. 0.)then
2689 electrolyte(jcacl2,jp,ibin) = (xeq_c(jc_ca) *na_Ma(ja_cl) + &
2690 xeq_a(ja_cl) *nc_Mc(jc_ca))/ &
2691 mw_electrolyte(jcacl2)
2694 if(xeq_c(jc_ca) .gt. 0. .and. xeq_a(ja_msa) .gt. 0.)then
2695 electrolyte(jcamsa2,jp,ibin) = (xeq_c(jc_ca) *na_Ma(ja_msa) + &
2696 xeq_a(ja_msa) *nc_Mc(jc_ca))/ &
2697 mw_electrolyte(jcamsa2)
2700 electrolyte(jh2so4, jp,ibin) = 0.0
2702 if(xeq_c(jc_h) .gt. 0. .and. xeq_a(ja_no3) .gt. 0.)then
2703 electrolyte(jhno3,jp,ibin) = (xeq_c(jc_h) *na_Ma(ja_no3) + &
2704 xeq_a(ja_no3)*nc_Mc(jc_h))/ &
2705 mw_electrolyte(jhno3)
2708 if(xeq_c(jc_h) .gt. 0. .and. xeq_a(ja_cl) .gt. 0.)then
2709 electrolyte(jhcl,jp,ibin) = (xeq_c(jc_h) *na_Ma(ja_cl) + &
2710 xeq_a(ja_cl)*nc_Mc(jc_h))/ &
2711 mw_electrolyte(jhcl)
2714 if(xeq_c(jc_h) .gt. 0. .and. xeq_a(ja_msa) .gt. 0.)then
2715 electrolyte(jmsa,jp,ibin) = (xeq_c(jc_h) *na_Ma(ja_msa) + &
2716 xeq_a(ja_msa)*nc_Mc(jc_h))/ &
2717 mw_electrolyte(jmsa)
2720 !--------------------------------------------------------------------
2722 elseif(icase.eq.2)then ! XT < 2 : SULFATE RICH DOMAIN
2724 store(imsa_a) = aer(imsa_a,jp,ibin)
2725 store(ica_a) = aer(ica_a, jp,ibin)
2727 call form_camsa2(store,jp,ibin,electrolyte)
2729 sum_na_nh4 = aer(ina_a,jp,ibin) + aer(inh4_a,jp,ibin)
2731 if(sum_na_nh4 .gt. 0.0)then
2732 f_na = aer(ina_a,jp,ibin)/sum_na_nh4
2733 f_nh4 = aer(inh4_a,jp,ibin)/sum_na_nh4
2739 ! first form msa electrolytes
2740 if(sum_na_nh4 .gt. store(imsa_a))then
2741 electrolyte(jnamsa,jp,ibin) = f_na *store(imsa_a)
2742 electrolyte(jnh4msa,jp,ibin) = f_nh4*store(imsa_a)
2743 rem_na = max(0.0_r8, aer(ina_a,jp,ibin) - electrolyte(jnamsa,jp,ibin)) ! remaining na RAZ 4/16/2014
2744 rem_nh4= max(0.0_r8, aer(inh4_a,jp,ibin)- electrolyte(jnh4msa,jp,ibin)) ! remaining nh4 RAZ 4/16/2014
2746 electrolyte(jnamsa,jp,ibin) = aer(ina_a,jp,ibin)
2747 electrolyte(jnh4msa,jp,ibin) = aer(inh4_a,jp,ibin)
2748 electrolyte(jmsa,jp,ibin) = max(0.0_r8, store(imsa_a) - sum_na_nh4) ! RAZ 4/16/2014
2749 rem_nh4 = 0.0 ! remaining nh4
2750 rem_na = 0.0 ! remaining na
2755 if(aer(iso4_a,jp,ibin).gt.0.0)then
2756 XT = (rem_nh4 + rem_na)/aer(iso4_a,jp,ibin)
2761 if(XT .le. 1.0)then ! h2so4 + bisulfate
2762 xh = max(0.0_r8, (1.0_r8 - XT)) ! RAZ 4/16/2014
2764 electrolyte(jh2so4,jp,ibin) = xh*aer(iso4_a,jp,ibin)
2765 electrolyte(jnh4hso4,jp,ibin) = xb*f_nh4*aer(iso4_a,jp,ibin)
2766 electrolyte(jnahso4,jp,ibin) = xb*f_na *aer(iso4_a,jp,ibin)
2767 elseif(XT .le. 1.5)then ! bisulfate + letovicite
2768 xb = max(0.0_r8, 3.0_r8 - 2.0_r8*XT) ! RAZ 4/16/2014
2769 xl = max(0.0_r8, XT - 1.0_r8) ! RAZ 4/16/2014
2770 electrolyte(jnh4hso4,jp,ibin) = xb*f_nh4*aer(iso4_a,jp,ibin)
2771 electrolyte(jnahso4,jp,ibin) = xb*f_na *aer(iso4_a,jp,ibin)
2772 electrolyte(jlvcite,jp,ibin) = xl*f_nh4*aer(iso4_a,jp,ibin)
2773 electrolyte(jna3hso4,jp,ibin) = xl*f_na *aer(iso4_a,jp,ibin)
2774 else ! letovicite + sulfate
2775 xl = max(0.0_r8, 2.0_r8 - XT) ! RAZ 4/16/2014
2776 xs = max(0.0_r8, 2.0_r8*XT - 3.0_r8) ! RAZ 4/16/2014
2777 electrolyte(jlvcite,jp,ibin) = xl*f_nh4*aer(iso4_a,jp,ibin)
2778 electrolyte(jna3hso4,jp,ibin) = xl*f_na *aer(iso4_a,jp,ibin)
2779 electrolyte(jnh4so4,jp,ibin) = xs*f_nh4*aer(iso4_a,jp,ibin)
2780 electrolyte(jna2so4,jp,ibin) = xs*f_na *aer(iso4_a,jp,ibin)
2783 electrolyte(jhno3,jp,ibin) = aer(ino3_a,jp,ibin)
2784 electrolyte(jhcl,jp,ibin) = aer(icl_a,jp,ibin)
2787 !---------------------------------------------------------
2789 ! calculate % composition EFFI
2791 !! do je = 1, nelectrolyte
2792 !! sum_dum = sum_dum + electrolyte(je,jp,ibin)
2795 !! if(sum_dum .eq. 0.)sum_dum = 1.0
2796 !! electrolyte_sum(jp,ibin) = sum_dum
2798 !! do je = 1, nelectrolyte
2799 !! epercent(je,jp,ibin) = 100.*electrolyte(je,jp,ibin)/sum_dum
2804 end subroutine ions_to_electrolytes
2808 !***********************************************************************
2809 ! part of MESA: calculates liquid electrolytes from ions
2812 ! - this subroutine is to be used for liquid-phase or total-phase only
2813 ! - this sub transfers caso4 and caco3 from liquid to solid phase
2815 ! author: Rahul A. Zaveri
2817 !-----------------------------------------------------------------------
2818 subroutine MESA_estimate_eleliquid(ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma, &
2819 nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a,eleliquid) ! TOUCH
2820 use module_data_mosaic_aero, only: r8,naer,nbin_a_max,nelectrolyte,ncation, &!Parameters
2821 nanion,jliquid, &!Parameters
2822 jh2so4,jhno3,jhcl,jmsa,jlvcite,jnh4no3,jnh4cl,jcamsa2,jcano3,jcacl2, &
2823 jnano3,jnacl,jnh4so4,jnh4hso4,jnh4msa,jna2so4,jnahso4,jnamsa,iso4_a, &
2824 ja_so4,ja_no3,ja_cl,imsa_a,ja_msa,jc_ca,ina_a,jc_na,inh4_a,jc_nh4,jc_h, &
2825 ica_a,ino3_a,icl_a,ja_hso4
2830 integer, intent(in) :: ibin
2831 real(r8), intent(inout) :: XT
2832 real(r8), intent(in), dimension(Ncation) :: zc,MW_c
2833 real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
2834 real(r8), intent(in), dimension(Nanion) :: za,MW_a
2835 real(r8), intent(inout), dimension(Nanion) :: xeq_a,na_Ma
2836 real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
2837 real(r8), intent(inout), dimension(nelectrolyte) :: eleliquid
2838 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2839 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2841 integer iaer, je, jc, ja, icase, jp
2842 real(r8) :: store(naer), sum_dum, sum_naza, sum_nczc, sum_na_nh4, &
2843 f_nh4, f_na, xh, xb, xl, xs, XT_d, XNa_d, XNH4_d, &
2845 real(r8) :: nc(ncation), na(nanion)
2846 real(r8) :: dum_ca, dum_no3, dum_cl, cano3, cacl2
2848 !nc(:) = 0.0_r8!BSINGH - initialized to zero
2850 ! remove negative concentrations, if any
2852 aer(iaer,jliquid,ibin) = max(0.0d0, aer(iaer,jliquid,ibin))
2856 ! calculate sulfate ratio
2857 call calculate_XT(ibin,jliquid,XT,aer)
2859 if(XT .ge. 2.0 .or. XT.lt.0.)then
2860 icase = 1 ! near neutral (acidity is caused by HCl and/or HNO3)
2862 icase = 2 ! acidic (acidity is caused by excess SO4)
2866 ! initialize to zero
2867 do je = 1, nelectrolyte
2872 !---------------------------------------------------------
2873 ! initialize moles of ions depending on the sulfate domain
2877 if(icase.eq.1)then ! XT >= 2 : SULFATE POOR DOMAIN
2879 dum_ca = aer(ica_a,jp,ibin)
2880 dum_no3 = aer(ino3_a,jp,ibin)
2881 dum_cl = aer(icl_a,jp,ibin)
2883 cano3 = min(dum_ca, 0.5*dum_no3)
2884 dum_ca = max(0.d0, dum_ca - cano3)
2885 dum_no3 = max(0.d0, dum_no3 - 2.*cano3)
2887 cacl2 = min(dum_ca, 0.5*dum_cl)
2888 dum_ca = max(0.d0, dum_ca - cacl2)
2889 dum_cl = max(0.d0, dum_cl - 2.*cacl2)
2892 na(ja_so4) = aer(iso4_a,jp,ibin)
2893 na(ja_no3) = aer(ino3_a,jp,ibin)
2894 na(ja_cl) = aer(icl_a, jp,ibin)
2895 na(ja_msa) = aer(imsa_a,jp,ibin)
2897 nc(jc_ca) = aer(ica_a, jp,ibin)
2898 nc(jc_na) = aer(ina_a, jp,ibin)
2899 nc(jc_nh4) = aer(inh4_a,jp,ibin)
2902 (2.*na(ja_so4)+na(ja_no3)+na(ja_cl)+na(ja_msa)) - &
2903 (2.*nc(jc_ca) +nc(jc_nh4)+nc(jc_na)) ) ! RAZ 11/11/2014: bug fix. remove nc(jc_h)
2905 if(cat_net .lt. 0.0)then
2909 else ! cat_net must be 0.0 or positive
2916 ! now compute equivalent fractions
2919 sum_naza = sum_naza + na(ja)*za(ja)
2924 sum_nczc = sum_nczc + nc(jc)*zc(jc)
2927 if(sum_naza .eq. 0. .or. sum_nczc .eq. 0.)then
2928 write(6,*)'ionic concentrations are zero in ibin', ibin
2929 write(6,*)'sum_naza = ', sum_naza
2930 write(6,*)'sum_nczc = ', sum_nczc
2935 xeq_a(ja) = na(ja)*za(ja)/sum_naza
2939 xeq_c(jc) = nc(jc)*zc(jc)/sum_nczc
2942 na_Ma(ja_so4) = na(ja_so4) *MW_a(ja_so4)
2943 na_Ma(ja_no3) = na(ja_no3) *MW_a(ja_no3)
2944 na_Ma(ja_cl) = na(ja_cl) *MW_a(ja_cl)
2945 na_Ma(ja_hso4)= na(ja_hso4)*MW_a(ja_hso4)
2946 na_Ma(ja_msa) = na(ja_msa) *MW_a(ja_msa)
2948 nc_Mc(jc_ca) = nc(jc_ca) *MW_c(jc_ca)
2949 nc_Mc(jc_na) = nc(jc_na) *MW_c(jc_na)
2950 nc_Mc(jc_nh4) = nc(jc_nh4)*MW_c(jc_nh4)
2951 nc_Mc(jc_h) = nc(jc_h) *MW_c(jc_h)
2954 ! now compute electrolyte moles
2955 eleliquid(jna2so4) = (xeq_c(jc_na) *na_Ma(ja_so4) + &
2956 xeq_a(ja_so4)*nc_Mc(jc_na))/ &
2957 mw_electrolyte(jna2so4)
2959 eleliquid(jnahso4) = (xeq_c(jc_na) *na_Ma(ja_hso4) + &
2960 xeq_a(ja_hso4)*nc_Mc(jc_na))/ &
2961 mw_electrolyte(jnahso4)
2963 eleliquid(jnamsa) = (xeq_c(jc_na) *na_Ma(ja_msa) + &
2964 xeq_a(ja_msa)*nc_Mc(jc_na))/ &
2965 mw_electrolyte(jnamsa)
2967 eleliquid(jnano3) = (xeq_c(jc_na) *na_Ma(ja_no3) + &
2968 xeq_a(ja_no3)*nc_Mc(jc_na))/ &
2969 mw_electrolyte(jnano3)
2971 eleliquid(jnacl) = (xeq_c(jc_na) *na_Ma(ja_cl) + &
2972 xeq_a(ja_cl) *nc_Mc(jc_na))/ &
2973 mw_electrolyte(jnacl)
2975 eleliquid(jnh4so4) = (xeq_c(jc_nh4)*na_Ma(ja_so4) + &
2976 xeq_a(ja_so4)*nc_Mc(jc_nh4))/ &
2977 mw_electrolyte(jnh4so4)
2979 eleliquid(jnh4hso4)= (xeq_c(jc_nh4)*na_Ma(ja_hso4) + &
2980 xeq_a(ja_hso4)*nc_Mc(jc_nh4))/ &
2981 mw_electrolyte(jnh4hso4)
2983 eleliquid(jnh4msa) = (xeq_c(jc_nh4) *na_Ma(ja_msa) + &
2984 xeq_a(ja_msa)*nc_Mc(jc_nh4))/ &
2985 mw_electrolyte(jnh4msa)
2987 eleliquid(jnh4no3) = (xeq_c(jc_nh4)*na_Ma(ja_no3) + &
2988 xeq_a(ja_no3)*nc_Mc(jc_nh4))/ &
2989 mw_electrolyte(jnh4no3)
2991 eleliquid(jnh4cl) = (xeq_c(jc_nh4)*na_Ma(ja_cl) + &
2992 xeq_a(ja_cl) *nc_Mc(jc_nh4))/ &
2993 mw_electrolyte(jnh4cl)
2995 eleliquid(jcamsa2) = (xeq_c(jc_ca) *na_Ma(ja_msa) + &
2996 xeq_a(ja_msa)*nc_Mc(jc_ca))/ &
2997 mw_electrolyte(jcamsa2)
2999 eleliquid(jcano3) = (xeq_c(jc_ca) *na_Ma(ja_no3) + &
3000 xeq_a(ja_no3)*nc_Mc(jc_ca))/ &
3001 mw_electrolyte(jcano3)
3003 eleliquid(jcacl2) = (xeq_c(jc_ca) *na_Ma(ja_cl) + &
3004 xeq_a(ja_cl) *nc_Mc(jc_ca))/ &
3005 mw_electrolyte(jcacl2)
3007 eleliquid(jh2so4) = (xeq_c(jc_h) *na_Ma(ja_hso4) + &
3008 xeq_a(ja_hso4)*nc_Mc(jc_h))/ &
3009 mw_electrolyte(jh2so4)
3011 eleliquid(jhno3) = (xeq_c(jc_h) *na_Ma(ja_no3) + &
3012 xeq_a(ja_no3)*nc_Mc(jc_h))/ &
3013 mw_electrolyte(jhno3)
3015 eleliquid(jhcl) = (xeq_c(jc_h) *na_Ma(ja_cl) + &
3016 xeq_a(ja_cl)*nc_Mc(jc_h))/ &
3017 mw_electrolyte(jhcl)
3019 eleliquid(jmsa) = (xeq_c(jc_h) *na_Ma(ja_msa) + &
3020 xeq_a(ja_msa)*nc_Mc(jc_h))/ &
3021 mw_electrolyte(jmsa)
3023 !--------------------------------------------------------------------
3025 elseif(icase.eq.2)then ! XT < 2 : SULFATE RICH DOMAIN
3029 store(iso4_a) = aer(iso4_a,jp,ibin)
3030 store(imsa_a) = aer(imsa_a,jp,ibin)
3031 store(inh4_a) = aer(inh4_a,jp,ibin)
3032 store(ina_a) = aer(ina_a, jp,ibin)
3033 store(ica_a) = aer(ica_a, jp,ibin)
3035 call form_camsa2(store,jp,ibin,electrolyte)
3037 sum_na_nh4 = store(ina_a) + store(inh4_a)
3038 if(sum_na_nh4 .gt. 0.0)then
3039 f_nh4 = store(inh4_a)/sum_na_nh4
3040 f_na = store(ina_a)/sum_na_nh4
3046 ! first form msa electrolytes
3047 if(sum_na_nh4 .gt. store(imsa_a))then
3048 eleliquid(jnh4msa) = f_nh4*store(imsa_a)
3049 eleliquid(jnamsa) = f_na *store(imsa_a)
3050 store(inh4_a)= store(inh4_a)-eleliquid(jnh4msa) ! remaining nh4
3051 store(ina_a) = store(ina_a) -eleliquid(jnamsa) ! remaining na
3053 eleliquid(jnh4msa) = store(inh4_a)
3054 eleliquid(jnamsa) = store(ina_a)
3055 eleliquid(jmsa) = store(imsa_a) - sum_na_nh4
3056 store(inh4_a)= 0.0 ! remaining nh4
3057 store(ina_a) = 0.0 ! remaining na
3060 if(store(iso4_a).eq.0.0)goto 10
3063 XNa_d = 1. + 0.5*store(ina_a)/store(iso4_a)
3064 xdum = store(iso4_a) - store(inh4_a)
3066 dum = ( (2.*store(iso4_a)) - &
3068 if(store(inh4_a) .gt. 0.0 .and. dum .gt. 0.0)then
3069 XNH4_d = 2.*store(inh4_a)/ &
3070 (2.*store(iso4_a) - store(ina_a))
3076 IF(store(inh4_a) .gt. 0.0)THEN
3077 if(XT_d .ge. XNa_d)then
3078 eleliquid(jna2so4) = 0.5*store(ina_a)
3080 if(XNH4_d .ge. 5./3.)then
3081 eleliquid(jnh4so4) = 1.5*store(ina_a) &
3082 - 3.*xdum - store(inh4_a)
3083 eleliquid(jlvcite) = 2.*xdum + store(inh4_a) &
3085 elseif(XNH4_d .ge. 1.5)then
3086 eleliquid(jnh4so4) = store(inh4_a)/5.
3087 eleliquid(jlvcite) = store(inh4_a)/5.
3088 elseif(XNH4_d .ge. 1.0)then
3089 eleliquid(jnh4so4) = store(inh4_a)/6.
3090 eleliquid(jlvcite) = store(inh4_a)/6.
3091 eleliquid(jnh4hso4)= store(inh4_a)/6.
3094 elseif(XT_d .gt. 1.0)then
3095 eleliquid(jnh4so4) = store(inh4_a)/6.
3096 eleliquid(jlvcite) = store(inh4_a)/6.
3097 eleliquid(jnh4hso4) = store(inh4_a)/6.
3098 eleliquid(jna2so4) = store(ina_a)/3.
3099 eleliquid(jnahso4) = store(ina_a)/3.
3100 elseif(XT_d .le. 1.0)then
3101 eleliquid(jna2so4) = store(ina_a)/4.
3102 eleliquid(jnahso4) = store(ina_a)/2.
3103 eleliquid(jlvcite) = store(inh4_a)/6.
3104 eleliquid(jnh4hso4) = store(inh4_a)/2.
3109 if(XT_d .gt. 1.0)then
3110 eleliquid(jna2so4) = store(ina_a) - store(iso4_a)
3111 eleliquid(jnahso4) = 2.*store(iso4_a) - &
3114 eleliquid(jna2so4) = store(ina_a)/4.
3115 eleliquid(jnahso4) = store(ina_a)/2.
3124 !---------------------------------------------------------
3128 end subroutine MESA_estimate_eleliquid
3132 !***********************************************************************
3133 ! part of MESA: completely dissolves small amounts of soluble salts
3135 ! author: Rahul A. Zaveri
3137 !-----------------------------------------------------------------------
3138 subroutine MESA_dissolve_small_salt(ibin,js,aer,electrolyte)
3140 use module_data_mosaic_aero, only:r8,naer,nbin_a_max,nelectrolyte,jsolid, &!Parameters
3141 jliquid, &!Parameters
3142 jh2so4,jhno3,jhcl,jlvcite,jnh4no3,jnh4cl,jcamsa2,jcano3,jcacl2,jnano3, &!TBD
3143 jnacl,jnh4so4,jnh4hso4,jnh4msa,jna2so4,jnahso4,jnamsa,iso4_a,ina_a, &!TBD
3144 inh4_a,jna3hso4,jcaso4,jcaco3,ica_a,ino3_a,icl_a
3149 integer, intent(in) :: ibin, js
3150 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
3151 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
3158 if(js .eq. jnh4so4)then
3159 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + &
3160 2.*electrolyte(js,jsolid,ibin)
3161 aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) + &
3162 electrolyte(js,jsolid,ibin)
3164 electrolyte(js,jsolid,ibin) = 0.0
3166 aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) + &
3167 electrolyte(jnh4cl,jp,ibin) + &
3168 2.*electrolyte(jnh4so4,jp,ibin) + &
3169 3.*electrolyte(jlvcite,jp,ibin) + &
3170 electrolyte(jnh4hso4,jp,ibin)+ &
3171 electrolyte(jnh4msa,jp,ibin)
3173 aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
3174 electrolyte(jna2so4,jp,ibin) + &
3175 2.*electrolyte(jna3hso4,jp,ibin)+ &
3176 electrolyte(jnahso4,jp,ibin) + &
3177 electrolyte(jnh4so4,jp,ibin) + &
3178 2.*electrolyte(jlvcite,jp,ibin) + &
3179 electrolyte(jnh4hso4,jp,ibin)+ &
3180 electrolyte(jh2so4,jp,ibin)
3185 if(js .eq. jlvcite)then
3186 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + &
3187 3.*electrolyte(js,jsolid,ibin)
3188 aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) + &
3189 2.*electrolyte(js,jsolid,ibin)
3191 electrolyte(js,jsolid,ibin) = 0.0
3193 aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) + &
3194 electrolyte(jnh4cl,jp,ibin) + &
3195 2.*electrolyte(jnh4so4,jp,ibin) + &
3196 3.*electrolyte(jlvcite,jp,ibin) + &
3197 electrolyte(jnh4hso4,jp,ibin)+ &
3198 electrolyte(jnh4msa,jp,ibin)
3200 aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
3201 electrolyte(jna2so4,jp,ibin) + &
3202 2.*electrolyte(jna3hso4,jp,ibin)+ &
3203 electrolyte(jnahso4,jp,ibin) + &
3204 electrolyte(jnh4so4,jp,ibin) + &
3205 2.*electrolyte(jlvcite,jp,ibin) + &
3206 electrolyte(jnh4hso4,jp,ibin)+ &
3207 electrolyte(jh2so4,jp,ibin)
3212 if(js .eq. jnh4hso4)then
3213 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + &
3214 electrolyte(js,jsolid,ibin)
3215 aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) + &
3216 electrolyte(js,jsolid,ibin)
3218 electrolyte(js,jsolid,ibin) = 0.0
3220 aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) + &
3221 electrolyte(jnh4cl,jp,ibin) + &
3222 2.*electrolyte(jnh4so4,jp,ibin) + &
3223 3.*electrolyte(jlvcite,jp,ibin) + &
3224 electrolyte(jnh4hso4,jp,ibin)+ &
3225 electrolyte(jnh4msa,jp,ibin)
3227 aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
3228 electrolyte(jna2so4,jp,ibin) + &
3229 2.*electrolyte(jna3hso4,jp,ibin)+ &
3230 electrolyte(jnahso4,jp,ibin) + &
3231 electrolyte(jnh4so4,jp,ibin) + &
3232 2.*electrolyte(jlvcite,jp,ibin) + &
3233 electrolyte(jnh4hso4,jp,ibin)+ &
3234 electrolyte(jh2so4,jp,ibin)
3239 if(js .eq. jna2so4)then
3240 aer(ina_a,jliquid,ibin) = aer(ina_a,jliquid,ibin) + &
3241 2.*electrolyte(js,jsolid,ibin)
3242 aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) + &
3243 electrolyte(js,jsolid,ibin)
3245 electrolyte(js,jsolid,ibin) = 0.0
3247 aer(ina_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
3248 electrolyte(jnacl,jp,ibin) + &
3249 2.*electrolyte(jna2so4,jp,ibin) + &
3250 3.*electrolyte(jna3hso4,jp,ibin)+ &
3251 electrolyte(jnahso4,jp,ibin) + &
3252 electrolyte(jnamsa,jp,ibin)
3254 aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
3255 electrolyte(jna2so4,jp,ibin) + &
3256 2.*electrolyte(jna3hso4,jp,ibin)+ &
3257 electrolyte(jnahso4,jp,ibin) + &
3258 electrolyte(jnh4so4,jp,ibin) + &
3259 2.*electrolyte(jlvcite,jp,ibin) + &
3260 electrolyte(jnh4hso4,jp,ibin)+ &
3261 electrolyte(jh2so4,jp,ibin)
3266 if(js .eq. jna3hso4)then
3267 aer(ina_a,jliquid,ibin) = aer(ina_a,jliquid,ibin) + &
3268 3.*electrolyte(js,jsolid,ibin)
3269 aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) + &
3270 2.*electrolyte(js,jsolid,ibin)
3272 electrolyte(js,jsolid,ibin) = 0.0
3274 aer(ina_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
3275 electrolyte(jnacl,jp,ibin) + &
3276 2.*electrolyte(jna2so4,jp,ibin) + &
3277 3.*electrolyte(jna3hso4,jp,ibin)+ &
3278 electrolyte(jnahso4,jp,ibin) + &
3279 electrolyte(jnamsa,jp,ibin)
3281 aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
3282 electrolyte(jna2so4,jp,ibin) + &
3283 2.*electrolyte(jna3hso4,jp,ibin)+ &
3284 electrolyte(jnahso4,jp,ibin) + &
3285 electrolyte(jnh4so4,jp,ibin) + &
3286 2.*electrolyte(jlvcite,jp,ibin) + &
3287 electrolyte(jnh4hso4,jp,ibin)+ &
3288 electrolyte(jh2so4,jp,ibin)
3293 if(js .eq. jnahso4)then
3294 aer(ina_a,jliquid,ibin) = aer(ina_a,jliquid,ibin) + &
3295 electrolyte(js,jsolid,ibin)
3296 aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) + &
3297 electrolyte(js,jsolid,ibin)
3299 electrolyte(js,jsolid,ibin) = 0.0
3301 aer(ina_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
3302 electrolyte(jnacl,jp,ibin) + &
3303 2.*electrolyte(jna2so4,jp,ibin) + &
3304 3.*electrolyte(jna3hso4,jp,ibin)+ &
3305 electrolyte(jnahso4,jp,ibin) + &
3306 electrolyte(jnamsa,jp,ibin)
3308 aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
3309 electrolyte(jna2so4,jp,ibin) + &
3310 2.*electrolyte(jna3hso4,jp,ibin)+ &
3311 electrolyte(jnahso4,jp,ibin) + &
3312 electrolyte(jnh4so4,jp,ibin) + &
3313 2.*electrolyte(jlvcite,jp,ibin) + &
3314 electrolyte(jnh4hso4,jp,ibin)+ &
3315 electrolyte(jh2so4,jp,ibin)
3320 if(js .eq. jnh4no3)then
3321 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + &
3322 electrolyte(js,jsolid,ibin)
3323 aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) + &
3324 electrolyte(js,jsolid,ibin)
3326 electrolyte(js,jsolid,ibin) = 0.0
3328 aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) + &
3329 electrolyte(jnh4cl,jp,ibin) + &
3330 2.*electrolyte(jnh4so4,jp,ibin) + &
3331 3.*electrolyte(jlvcite,jp,ibin) + &
3332 electrolyte(jnh4hso4,jp,ibin)+ &
3333 electrolyte(jnh4msa,jp,ibin)
3335 aer(ino3_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
3336 2.*electrolyte(jcano3,jp,ibin) + &
3337 electrolyte(jnh4no3,jp,ibin) + &
3338 electrolyte(jhno3,jp,ibin)
3343 if(js .eq. jnh4cl)then
3344 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + &
3345 electrolyte(js,jsolid,ibin)
3346 aer(icl_a,jliquid,ibin) = aer(icl_a,jliquid,ibin) + &
3347 electrolyte(js,jsolid,ibin)
3349 electrolyte(js,jsolid,ibin) = 0.0
3351 aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) + &
3352 electrolyte(jnh4cl,jp,ibin) + &
3353 2.*electrolyte(jnh4so4,jp,ibin) + &
3354 3.*electrolyte(jlvcite,jp,ibin) + &
3355 electrolyte(jnh4hso4,jp,ibin)+ &
3356 electrolyte(jnh4msa,jp,ibin)
3358 aer(icl_a,jp,ibin) = electrolyte(jnacl,jp,ibin) + &
3359 2.*electrolyte(jcacl2,jp,ibin) + &
3360 electrolyte(jnh4cl,jp,ibin) + &
3361 electrolyte(jhcl,jp,ibin)
3366 if(js .eq. jnano3)then
3367 aer(ina_a,jliquid,ibin) = aer(ina_a,jliquid,ibin) + &
3368 electrolyte(js,jsolid,ibin)
3369 aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) + &
3370 electrolyte(js,jsolid,ibin)
3372 electrolyte(js,jsolid,ibin) = 0.0
3374 aer(ina_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
3375 electrolyte(jnacl,jp,ibin) + &
3376 2.*electrolyte(jna2so4,jp,ibin) + &
3377 3.*electrolyte(jna3hso4,jp,ibin)+ &
3378 electrolyte(jnahso4,jp,ibin) + &
3379 electrolyte(jnamsa,jp,ibin)
3381 aer(ino3_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
3382 2.*electrolyte(jcano3,jp,ibin) + &
3383 electrolyte(jnh4no3,jp,ibin) + &
3384 electrolyte(jhno3,jp,ibin)
3389 if(js .eq. jnacl)then
3390 aer(ina_a,jliquid,ibin) = aer(ina_a,jliquid,ibin) + &
3391 electrolyte(js,jsolid,ibin)
3392 aer(icl_a,jliquid,ibin) = aer(icl_a,jliquid,ibin) + &
3393 electrolyte(js,jsolid,ibin)
3395 electrolyte(js,jsolid,ibin) = 0.0
3397 aer(ina_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
3398 electrolyte(jnacl,jp,ibin) + &
3399 2.*electrolyte(jna2so4,jp,ibin) + &
3400 3.*electrolyte(jna3hso4,jp,ibin)+ &
3401 electrolyte(jnahso4,jp,ibin) + &
3402 electrolyte(jnamsa,jp,ibin)
3404 aer(icl_a,jp,ibin) = electrolyte(jnacl,jp,ibin) + &
3405 2.*electrolyte(jcacl2,jp,ibin) + &
3406 electrolyte(jnh4cl,jp,ibin) + &
3407 electrolyte(jhcl,jp,ibin)
3412 if(js .eq. jcano3)then
3413 aer(ica_a,jliquid,ibin) = aer(ica_a,jliquid,ibin) + &
3414 electrolyte(js,jsolid,ibin)
3415 aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) + &
3416 2.*electrolyte(js,jsolid,ibin)
3418 electrolyte(js,jsolid,ibin) = 0.0
3420 aer(ica_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
3421 electrolyte(jcano3,jp,ibin) + &
3422 electrolyte(jcacl2,jp,ibin) + &
3423 electrolyte(jcaco3,jp,ibin) + &
3424 electrolyte(jcamsa2,jp,ibin)
3426 aer(ino3_a,jp,ibin) = electrolyte(jnano3,jp,ibin) + &
3427 2.*electrolyte(jcano3,jp,ibin) + &
3428 electrolyte(jnh4no3,jp,ibin) + &
3429 electrolyte(jhno3,jp,ibin)
3434 if(js .eq. jcacl2)then
3435 aer(ica_a,jliquid,ibin) = aer(ica_a,jliquid,ibin) + &
3436 electrolyte(js,jsolid,ibin)
3437 aer(icl_a,jliquid,ibin) = aer(icl_a,jliquid,ibin) + &
3438 2.*electrolyte(js,jsolid,ibin)
3440 electrolyte(js,jsolid,ibin) = 0.0
3442 aer(ica_a,jp,ibin) = electrolyte(jcaso4,jp,ibin) + &
3443 electrolyte(jcano3,jp,ibin) + &
3444 electrolyte(jcacl2,jp,ibin) + &
3445 electrolyte(jcaco3,jp,ibin) + &
3446 electrolyte(jcamsa2,jp,ibin)
3448 aer(icl_a,jp,ibin) = electrolyte(jnacl,jp,ibin) + &
3449 2.*electrolyte(jcacl2,jp,ibin) + &
3450 electrolyte(jnh4cl,jp,ibin) + &
3451 electrolyte(jhcl,jp,ibin)
3456 end subroutine MESA_dissolve_small_salt
3460 !***********************************************************************
3461 ! electrolyte formation subroutines
3463 ! author: Rahul A. Zaveri
3465 !-----------------------------------------------------------------------
3466 subroutine form_caso4(store,jp,ibin,electrolyte)
3467 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
3473 integer, intent(in) :: jp, ibin
3474 real(r8), intent(inout), dimension(naer) :: store
3475 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
3477 electrolyte(jcaso4,jp,ibin) = min(store(ica_a),store(iso4_a))
3478 store(ica_a) = ( (store(ica_a)) - &
3479 (electrolyte(jcaso4,jp,ibin)) )
3480 store(iso4_a) = ( (store(iso4_a)) - &
3481 (electrolyte(jcaso4,jp,ibin)) )
3482 store(ica_a) = max(0.d0, store(ica_a))
3483 store(iso4_a) = max(0.d0, store(iso4_a))
3486 end subroutine form_caso4
3490 subroutine form_camsa2(store,jp,ibin,electrolyte)
3491 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
3492 imsa_a,ica_a,jcamsa2
3496 integer, intent(in) :: jp, ibin
3497 real(r8), intent(inout), dimension(naer) :: store
3498 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
3500 electrolyte(jcamsa2,jp,ibin) = min(store(ica_a),0.5*store(imsa_a))
3501 store(ica_a) = ( (store(ica_a)) - &
3502 (electrolyte(jcamsa2,jp,ibin)) )
3503 store(imsa_a) = ( (store(imsa_a)) - &
3504 (2.*electrolyte(jcamsa2,jp,ibin)) )
3505 store(ica_a) = max(0.d0, store(ica_a))
3506 store(imsa_a) = max(0.d0, store(imsa_a))
3509 end subroutine form_camsa2
3513 !***********************************************************************
3514 ! computes mass transfer coefficients for each condensing species for
3515 ! all the aerosol bins
3517 ! author: Rahul A. Zaveri
3519 !-----------------------------------------------------------------------
3520 subroutine aerosolmtc( jaerosolstate, num_a, Dp_wet_a, sigmag_a, P_atm, T_K, kg )
3522 use module_data_mosaic_aero, only: r8, nbin_a_max, nbin_a, naer, naercomp, &!Parameters
3523 ngas_aerchtot, ngas_volatile, nelectrolyte, ngas_ioa, &
3524 mMODAL, no_aerosol, mUNSTRUCTURED, mSECTIONAL, mSIZE_FRAMEWORK, &!Input
3525 isoa_first, mw_gas, v_molar_gas, &!TBD
3526 i_gas2bin_uptk_flag, m_gas2bin_uptk_flag, &
3527 use_cam5mam_accom_coefs, mosaic_vars_aa_type
3532 !Subroutine Arguments
3533 integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate
3535 real(r8), intent(in) :: P_atm,T_K
3536 real(r8), intent(in), dimension(nbin_a_max) :: num_a
3537 real(r8), intent(in), dimension(nbin_a_max) :: sigmag_a
3538 real(r8), intent(inout), dimension(nbin_a_max) :: Dp_wet_a
3539 real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
3543 parameter (nghq = 2) ! gauss-hermite quadrature order
3544 integer ibin, iq, iv
3545 real(r8) :: tworootpi, root2, beta
3546 parameter (tworootpi = 3.5449077, root2 = 1.4142135, beta = 2.0)
3547 real(r8) :: cdum, Dp, Dp_avg, Fkn, Kn, lnsg, lnDpgn, lnDp, speed, &
3549 real(r8) :: xghq(nghq), wghq(nghq) ! quadrature abscissae and weights
3550 real(r8) :: accom(ngas_aerchtot) ! keep local
3551 real(r8) :: freepath(ngas_aerchtot) ! keep local
3552 real(r8) :: Dg(ngas_aerchtot) ! keep local
3553 !real(r8) :: fuchs_sutugin ! mosaic func
3554 !real(r8) :: gas_diffusivity ! mosaic func
3555 !real(r8) :: mean_molecular_speed ! mosaic func
3557 ! mass accommodation coefficients
3559 if ( use_cam5mam_accom_coefs > 0 ) tmpa = 0.65
3560 accom(1:ngas_aerchtot) = tmpa ! default
3561 ! accom(ih2so4_g) = tmpa
3562 ! accom(ihno3_g) = tmpa
3563 ! accom(ihcl_g) = tmpa
3564 ! accom(inh3_g) = tmpa
3565 ! accom(imsa_g) = tmpa
3566 ! accom(iaro1_g) = tmpa
3567 ! accom(iaro2_g) = tmpa
3568 ! accom(ialk1_g) = tmpa
3569 ! accom(iole1_g) = tmpa
3570 ! accom(iapi1_g) = tmpa
3571 ! accom(iapi2_g) = tmpa
3572 ! accom(ilim1_g) = tmpa
3573 ! accom(ilim2_g) = tmpa
3575 ! quadrature weights
3576 xghq(1) = 0.70710678
3577 xghq(2) = -0.70710678
3578 wghq(1) = 0.88622693
3579 wghq(2) = 0.88622693
3582 ! calculate gas diffusivity and mean free path for condensing gases
3585 speed = mean_molecular_speed(T_K,mw_gas(iv)) ! cm/s
3586 Dg(iv) = gas_diffusivity(T_K,P_atm,mw_gas(iv),v_molar_gas(iv)) ! cm^2/s
3587 freepath(iv) = 3.*Dg(iv)/speed ! cm
3591 do iv = isoa_first, ngas_volatile
3592 speed = mean_molecular_speed(T_K,mw_gas(iv)) ! cm/s
3593 Dg(iv) = 0.1 ! cm^2/s
3594 freepath(iv) = 3.*Dg(iv)/speed
3598 do iv = (ngas_volatile+1), ngas_aerchtot
3599 speed = mean_molecular_speed(t_k,mw_gas(iv)) ! cm/s
3600 dg(iv) = gas_diffusivity(t_k,p_atm,mw_gas(iv),v_molar_gas(iv)) ! cm^2/s
3601 freepath(iv) = 3.*dg(iv)/speed ! cm
3605 ! calc mass transfer coefficients for gases over various aerosol bins
3607 if (mSIZE_FRAMEWORK .eq. mMODAL) then
3609 ! for modal approach
3610 do 10 ibin = 1, nbin_a
3612 if(jaerosolstate(ibin) .eq. no_aerosol)goto 10
3614 lnsg = log(sigmag_a(ibin))
3616 ! following 2 lines were incorrect as Dp_wet_a is wet "average" Dp
3617 ! Dpgn_a(ibin) = Dp_wet_a(ibin) ! cm
3618 ! lnDpgn = log(Dpgn_a(ibin))
3619 ! do this instead which gives
3620 ! lnDpgn = ln( wet geometric-mean Dp of number distribution )
3621 lnDpgn = log(Dp_wet_a(ibin)) - 1.5*lnsg*lnsg
3623 cdum = tworootpi*num_a(ibin)* &
3624 exp(beta*lnDpgn + 0.5*(beta*lnsg)**2)
3626 do 20 iv = 1, ngas_aerchtot
3629 do 30 iq = 1, nghq ! sum over gauss-hermite quadrature points
3630 lnDp = lnDpgn + beta*lnsg**2 + root2*lnsg*xghq(iq)
3632 Kn = 2.*freepath(iv)/Dp
3633 Fkn = fuchs_sutugin(Kn,accom(iv))
3634 sumghq = sumghq + wghq(iq)*Dp*Fkn/(Dp**beta)
3637 kg(iv,ibin) = cdum*Dg(iv)*sumghq ! 1/s
3642 elseif ((mSIZE_FRAMEWORK .eq. mSECTIONAL ) .or. &
3643 (mSIZE_FRAMEWORK .eq. mUNSTRUCTURED)) then
3645 ! for sectional approach
3646 do 11 ibin = 1, nbin_a
3648 if(jaerosolstate(ibin) .eq. no_aerosol)goto 11
3650 cdum = 6.283185*Dp_wet_a(ibin)*num_a(ibin)
3652 do 21 iv = 1, ngas_aerchtot
3653 Kn = 2.*freepath(iv)/Dp_wet_a(ibin)
3654 Fkn = fuchs_sutugin(Kn,accom(iv))
3655 kg(iv,ibin) = cdum*Dg(iv)*Fkn ! 1/s
3661 call wrf_message('Error in the choice of mSIZE_FRAMEWORK')
3662 call wrf_error_fatal('Stopping in subr. aerosolmtc')
3666 if (m_gas2bin_uptk_flag <= 0) then
3667 ! when m_gas2bin_uptk_flag <= 0, some gases cannot condense onto every bin
3669 do iv = 1, ngas_aerchtot
3670 if (i_gas2bin_uptk_flag(iv,ibin) <= 0) kg(iv,ibin) = 0.0
3676 end subroutine aerosolmtc
3680 !***********************************************************************
3681 ! calculates dry and wet aerosol properties: density, refractive indices
3683 ! author: Rahul A. Zaveri
3685 !-----------------------------------------------------------------------
3686 subroutine calc_dry_n_wet_aerosol_props( &
3687 ibin, jaerosolstate, aer, electrolyte, water_a, num_a, & ! input
3688 dens_comp_a, mw_comp_a, dens_aer_mac, mw_aer_mac, ref_index_a, & ! input
3689 Dp_dry_a, Dp_wet_a, dp_core_a, & ! output
3690 area_dry_a, area_wet_a, mass_dry_a, mass_wet_a, & ! output
3691 vol_dry_a, vol_wet_a, dens_dry_a, dens_wet_a, & ! output
3692 ri_shell_a, ri_core_a, ri_avg_a ) ! output
3693 ! include 'v33com9a'
3695 use module_data_mosaic_constants, only: piover4,piover6,third
3696 use module_data_mosaic_aero, only: r8,nbin_a_max,naer,nelectrolyte,naercomp, &!Parameters
3697 ngas_soa,no_aerosol,msectional, &!Parameters
3698 maeroptic_aero,msize_framework, &!Input
3699 inh4_a,ina_a,ica_a,ico3_a,imsa_a,icl_a,ino3_a,jtotal,iso4_a,ioc_a,joc, &!TBD
3700 ibc_a,jbc,ioin_a,join,jh2o, isoa_first, jsoa_first
3702 use module_data_mosaic_asecthp, only: dcen_sect,isize_of_ibin,itype_of_ibin
3707 integer, intent(in) :: ibin
3708 integer, intent(in), dimension(nbin_a_max) :: jaerosolstate
3710 real(r8), intent(in), dimension(nbin_a_max) :: num_a
3711 real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
3712 real(r8), intent(in), dimension(naercomp) :: dens_comp_a,mw_comp_a
3713 real(r8), intent(inout), dimension(nbin_a_max) :: Dp_dry_a,Dp_wet_a
3714 real(r8), intent(inout), dimension(nbin_a_max) :: dp_core_a,vol_dry_a
3715 real(r8), intent(inout), dimension(nbin_a_max) :: vol_wet_a,dens_wet_a,water_a
3716 real(r8), intent(inout), dimension(nbin_a_max) :: area_dry_a,area_wet_a
3717 real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_a,mass_wet_a
3718 real(r8), intent(inout), dimension(nbin_a_max) :: dens_dry_a
3719 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
3720 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
3722 complex, intent(in), dimension(naercomp) :: ref_index_a
3723 complex, intent(inout), dimension(nbin_a_max) :: ri_avg_a,ri_core_a,ri_shell_a
3725 integer i, iaer, isize, itype, j, jc, je, k
3726 real(r8) :: aer_H, duma, vol_core, vol_shell, vol_dum
3727 real(r8),dimension(naercomp) :: comp_a
3728 complex rixvol_tot, rixvol_core, rixvol_shell
3731 ! calculate dry mass and dry volume of a bin
3732 mass_dry_a(ibin) = 0.0 ! initialize to 0.0
3733 vol_dry_a(ibin) = 0.0 ! initialize to 0.0
3734 area_dry_a(ibin) = 0.0 ! initialize to 0.0
3736 if(jaerosolstate(ibin) .ne. no_aerosol)then
3738 aer_H = (2.*aer(iso4_a,jtotal,ibin) + &
3739 aer(ino3_a,jtotal,ibin) + &
3740 aer(icl_a,jtotal,ibin) + &
3741 aer(imsa_a,jtotal,ibin) + &
3742 2.*aer(ico3_a,jtotal,ibin))- &
3743 (2.*aer(ica_a,jtotal,ibin) + &
3744 aer(ina_a,jtotal,ibin) + &
3745 aer(inh4_a,jtotal,ibin))
3746 aer_H = max(aer_H, 0.0d0)
3749 mass_dry_a(ibin) = mass_dry_a(ibin) + &
3750 aer(iaer,jtotal,ibin)*mw_aer_mac(iaer) ! ng/m^3(air)
3751 vol_dry_a(ibin) = vol_dry_a(ibin) + &
3752 aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer) ! ncc/m^3(air)
3754 mass_dry_a(ibin) = mass_dry_a(ibin) + aer_H
3755 vol_dry_a(ibin) = vol_dry_a(ibin) + aer_H
3757 mass_dry_a(ibin) = mass_dry_a(ibin)*1.e-15 ! g/cc(air)
3758 vol_dry_a(ibin) = vol_dry_a(ibin)*1.e-15 ! cc(aer)/cc(air)
3760 ! wet mass and wet volume
3761 mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3 ! g/cc(air)
3762 vol_wet_a(ibin) = vol_dry_a(ibin) + water_a(ibin)*1.e-3 ! cc(aer)/cc(air)
3764 ! calculate mean dry and wet particle densities
3765 dens_dry_a(ibin) = mass_dry_a(ibin)/vol_dry_a(ibin) ! g/cc(aerosol)
3766 dens_wet_a(ibin) = mass_wet_a(ibin)/vol_wet_a(ibin) ! g/cc(aerosol)
3768 ! calculate mean dry and wet particle diameters
3769 Dp_dry_a(ibin)=(vol_dry_a(ibin)/(piover6*num_a(ibin)))**third ! cm
3770 Dp_wet_a(ibin)=(vol_wet_a(ibin)/(piover6*num_a(ibin)))**third ! cm
3772 ! calculate mean dry and wet particle surface areas
3773 area_dry_a(ibin)= piover4*num_a(ibin)*Dp_dry_a(ibin)**2 ! cm^2/cc(air)
3774 area_wet_a(ibin)= piover4*num_a(ibin)*Dp_wet_a(ibin)**2 ! cm^2/cc(air)
3776 ! calculate volume average refractive index
3777 ! load comp_a array with component mass concentrations
3779 ! rahul had turned this off, but it is needed
3780 ! if(1 == 1)go to 100 ! TEMP
3781 if (maeroptic_aero <= 0) goto 100
3783 do je = 1, nelectrolyte
3784 comp_a(je)=electrolyte(je,jtotal,ibin)*mw_comp_a(je)*1.e-15 ! g/cc(air)
3786 comp_a(joc) = aer(ioc_a, jtotal,ibin)*mw_comp_a(joc )*1.e-15 ! g/cc(air)
3787 comp_a(jbc) = aer(ibc_a, jtotal,ibin)*mw_comp_a(jbc )*1.e-15 ! g/cc(air)
3788 comp_a(join) = aer(ioin_a, jtotal,ibin)*mw_comp_a(join )*1.e-15 ! g/cc(air)
3790 ! comp_a(jaro1)= aer(iaro1_a,jtotal,ibin)*mw_comp_a(jaro1)*1.e-15 ! g/cc(air)
3791 ! comp_a(jaro2)= aer(iaro2_a,jtotal,ibin)*mw_comp_a(jaro2)*1.e-15 ! g/cc(air)
3792 ! comp_a(jalk1)= aer(ialk1_a,jtotal,ibin)*mw_comp_a(jalk1)*1.e-15 ! g/cc(air)
3793 ! comp_a(jole1)= aer(iole1_a,jtotal,ibin)*mw_comp_a(jole1)*1.e-15 ! g/cc(air)
3794 ! comp_a(japi1)= aer(iapi1_a,jtotal,ibin)*mw_comp_a(japi1)*1.e-15 ! g/cc(air)
3795 ! comp_a(japi2)= aer(iapi2_a,jtotal,ibin)*mw_comp_a(japi2)*1.e-15 ! g/cc(air)
3796 ! comp_a(jlim1)= aer(ilim1_a,jtotal,ibin)*mw_comp_a(jlim1)*1.e-15 ! g/cc(air)
3797 ! comp_a(jlim2)= aer(ilim2_a,jtotal,ibin)*mw_comp_a(jlim2)*1.e-15 ! g/cc(air)
3799 j = jsoa_first + k - 1
3800 i = isoa_first + k - 1
3801 comp_a(j) = aer(i,jtotal,ibin)*mw_comp_a(j)*1.e-15 ! g/cc(air)
3804 comp_a(jh2o) = water_a(ibin)*1.e-3 ! g/cc(air)
3806 rixvol_tot = (0.0,0.0)
3808 comp_a(jc) = max( 0.0d0, comp_a(jc) )
3809 rixvol_tot = rixvol_tot &
3810 + ref_index_a(jc)*comp_a(jc)/dens_comp_a(jc)
3812 ri_avg_a(ibin) = rixvol_tot/vol_wet_a(ibin)
3815 ! shell/core calcs - first set values to default (corresponding to zero core)
3817 ri_shell_a(ibin) = ri_avg_a(ibin)
3818 ri_core_a(ibin) = (0.0,0.0)
3819 Dp_core_a(ibin) = 0.0
3821 ! sum ri*vol and vol for core species (bc and optionally oin=dust)
3822 ! currently just bc in core, but what about insoluble oin and dust species ???
3824 rixvol_core = ref_index_a(jc)*comp_a(jc)/dens_comp_a(jc)
3825 vol_core = comp_a(jc)/dens_comp_a(jc)
3826 vol_core = max( 0.0d0, min( vol_core, vol_wet_a(ibin) ) )
3828 ! neglect core if (core volume) < 1.0d-9*(total volume)
3829 ! or (core volume) < 1.0d-22 cm3 = (0.58 nm)**3
3830 ! neglect shell using similar criteria
3831 vol_dum = max( 1.0d-22, 1.0d-9*vol_wet_a(ibin) )
3832 vol_shell = vol_wet_a(ibin) - vol_core
3833 if (vol_core >= vol_dum) then
3834 if (vol_shell < vol_dum) then
3835 ri_shell_a(ibin) = (0.0,0.0)
3836 ri_core_a(ibin) = ri_avg_a(ibin)
3837 Dp_core_a(ibin) = Dp_wet_a(ibin)
3839 ri_core_a(ibin) = rixvol_core/vol_core
3840 Dp_core_a(ibin) = Dp_wet_a(ibin) &
3841 * (vol_core/vol_wet_a(ibin))**third
3843 if (vol_shell >= vol_dum) then
3844 rixvol_shell = rixvol_tot - rixvol_core
3845 ri_shell_a(ibin) = rixvol_shell/vol_shell
3847 ri_shell_a(ibin) = (0.0,0.0)
3853 ! use defaults when (jaerosolstate(ibin) .eq. no_aerosol)
3855 dens_dry_a(ibin) = 1.0 ! g/cc(aerosol)
3856 dens_wet_a(ibin) = 1.0 ! g/cc(aerosol)
3857 ! Dp_dry_a(ibin) = dcen_sect(ibin) ! cm
3858 ! Dp_wet_a(ibin) = dcen_sect(ibin) ! cm
3859 if (msize_framework == msectional) then
3860 isize = isize_of_ibin(ibin)
3861 itype = itype_of_ibin(ibin)
3862 Dp_dry_a(ibin) = dcen_sect(isize,itype)
3863 Dp_wet_a(ibin) = Dp_dry_a(ibin)
3866 ri_avg_a(ibin) = (1.5,0.0)
3867 ri_shell_a(ibin) = (1.5,0.0)
3868 ri_core_a(ibin) = (0.0,0.0)
3869 Dp_core_a(ibin) = 0.0
3871 endif ! if(jaerosolstate(ibin) .ne. no_aerosol)then
3877 end subroutine calc_dry_n_wet_aerosol_props
3881 !***********************************************************************
3882 ! forms electrolytes from ions
3884 ! author: Rahul A. Zaveri
3886 !-----------------------------------------------------------------------
3887 subroutine form_electrolytes(jp,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)
3888 use module_data_mosaic_aero, only: r8,naer,nbin_a_max, &
3889 ngas_aerchtot,ngas_volatile,nelectrolyte,jsolid, &
3890 imsa_a,iso4_a,ica_a,ina_a,inh4_a,ino3_a,icl_a,ico3_a
3895 integer, intent(in) :: ibin, jp
3896 real(r8), intent(inout) :: XT
3897 real(r8), intent(inout) :: tot_cl_in
3898 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
3899 real(r8), intent(inout), dimension(ngas_volatile) :: total_species
3900 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
3901 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
3903 integer i, iXT_case, j, je
3904 real(r8) :: sum_dum, XNa_prime, XNH4_prime, XT_prime
3905 real(r8) :: store(naer)
3907 ! remove negative concentrations, if any
3909 ! aer(i,jp,ibin) = max(0.0d0, aer(i,jp,ibin)) ! EFFI
3913 ! call calculate_XT(ibin,jp,XT,aer) ! EFFI
3915 if( (aer(iso4_a,jp,ibin)+aer(imsa_a,jp,ibin)) .gt.0.0)then
3916 XT = ( aer(inh4_a,jp,ibin) + &
3917 aer(ina_a,jp,ibin) + &
3918 2.*aer(ica_a,jp,ibin) )/ &
3919 (aer(iso4_a,jp,ibin)+0.5*aer(imsa_a,jp,ibin))
3927 ! if(XT .ge. 1.9999 .or. XT.lt.0.)then
3928 if(XT .ge. 2.0 .or. XT.lt.0.)then ! RAZ 11/10/2014
3929 iXT_case = 1 ! near neutral (acidity is caused by HCl and/or HNO3)
3931 iXT_case = 2 ! acidic (acidity is caused by excess SO4)
3936 ! put total aer(*) into store(*)
3937 store(iso4_a) = aer(iso4_a,jp,ibin)
3938 store(ino3_a) = aer(ino3_a,jp,ibin)
3939 store(icl_a) = aer(icl_a, jp,ibin)
3940 store(imsa_a) = aer(imsa_a,jp,ibin)
3941 store(ico3_a) = aer(ico3_a,jp,ibin)
3942 store(inh4_a) = aer(inh4_a,jp,ibin)
3943 store(ina_a) = aer(ina_a, jp,ibin)
3944 store(ica_a) = aer(ica_a, jp,ibin)
3947 electrolyte(j,jp,ibin) = 0.0
3951 !---------------------------------------------------------
3953 if(iXT_case.eq.1)then
3955 ! XT >= 2 : sulfate deficient
3956 call form_caso4(store,jp,ibin,electrolyte)
3957 call form_camsa2(store,jp,ibin,electrolyte)
3958 call form_na2so4(store,jp,ibin,electrolyte)
3959 call form_namsa(store,jp,ibin,electrolyte)
3960 call form_cano3(store,jp,ibin,electrolyte)
3961 call form_nano3(store,jp,ibin,electrolyte)
3962 call form_nacl(store,jp,ibin,aer,gas,electrolyte,total_species,tot_cl_in)
3963 call form_cacl2(store,jp,ibin,electrolyte)
3964 call form_caco3(store,jp,ibin,aer,electrolyte)
3965 call form_nh4so4(store,jp,ibin,electrolyte)
3966 call form_nh4msa(store,jp,ibin,electrolyte)
3967 call form_nh4no3(store,jp,ibin,electrolyte)
3968 call form_nh4cl(store,jp,ibin,electrolyte)
3969 call form_msa(store,jp,ibin,electrolyte)
3971 if(jp .eq. jsolid)then
3972 call degas_hno3(store,jp,ibin,aer,gas,electrolyte)
3973 call degas_hcl(store,jp,ibin,aer,gas,electrolyte)
3974 call degas_nh3(store,jp,ibin,aer,gas)
3976 call form_hno3(store,jp,ibin,electrolyte)
3977 call form_hcl(store,jp,ibin,electrolyte)
3978 call degas_nh3(store,jp,ibin,aer,gas)
3983 elseif(iXT_case.eq.2)then
3985 ! XT < 2 : sulfate enough or sulfate excess
3987 call form_caso4(store,jp,ibin,electrolyte)
3988 call form_camsa2(store,jp,ibin,electrolyte)
3989 call form_namsa(store,jp,ibin,electrolyte)
3990 call form_nh4msa(store,jp,ibin,electrolyte)
3991 call form_msa(store,jp,ibin,electrolyte)
3993 if(store(iso4_a).eq.0.0)goto 10
3996 XT_prime =(store(ina_a)+store(inh4_a))/ &
3998 XNa_prime=0.5*store(ina_a)/store(iso4_a) + 1.
4000 if(XT_prime.ge.XNa_prime)then
4001 call form_na2so4(store,jp,ibin,electrolyte)
4003 if(store(iso4_a).gt.1.e-15)then
4004 XNH4_prime = store(inh4_a)/store(iso4_a)
4007 if(XNH4_prime .ge. 1.5)then
4008 call form_nh4so4_lvcite(store,jp,ibin,electrolyte)
4010 call form_lvcite_nh4hso4(store,jp,ibin,electrolyte)
4013 elseif(XT_prime.ge.1.)then
4014 call form_nh4hso4(store,jp,ibin,electrolyte)
4015 call form_na2so4_nahso4(store,jp,ibin,electrolyte)
4016 elseif(XT_prime.lt.1.)then
4017 call form_nahso4(store,jp,ibin,electrolyte)
4018 call form_nh4hso4(store,jp,ibin,electrolyte)
4019 call form_h2so4(store,jp,ibin,electrolyte)
4022 10 if(jp .eq. jsolid)then
4023 call degas_hno3(store,jp,ibin,aer,gas,electrolyte)
4024 call degas_hcl(store,jp,ibin,aer,gas,electrolyte)
4025 call degas_nh3(store,jp,ibin,aer,gas)
4027 call form_hno3(store,jp,ibin,electrolyte)
4028 call form_hcl(store,jp,ibin,electrolyte)
4029 call degas_nh3(store,jp,ibin,aer,gas)
4035 ! re-calculate ions to eliminate round-off errors
4036 call electrolytes_to_ions(jp, ibin,aer,electrolyte)
4037 !---------------------------------------------------------
4039 ! calculate % composition EFFI
4041 !! do je = 1, nelectrolyte
4042 !! electrolyte(je,jp,ibin) = max(0.d0,electrolyte(je,jp,ibin)) ! remove -ve EFFI
4043 !! sum_dum = sum_dum + electrolyte(je,jp,ibin)
4046 !! if(sum_dum .eq. 0.)sum_dum = 1.0
4047 !! electrolyte_sum(jp,ibin) = sum_dum
4049 !! do je = 1, nelectrolyte
4050 !! epercent(je,jp,ibin) = 100.*electrolyte(je,jp,ibin)/sum_dum
4055 end subroutine form_electrolytes
4059 subroutine form_na2so4(store,jp,ibin,electrolyte)
4060 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4061 iso4_a,ina_a,jna2so4
4065 integer, intent(in) :: jp, ibin
4066 real(r8), intent(inout), dimension(naer) :: store(naer)
4067 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4068 electrolyte(jna2so4,jp,ibin) = min(.5*store(ina_a), &
4070 store(ina_a) =( (store(ina_a)) - &
4071 (2.*electrolyte(jna2so4,jp,ibin)) )
4072 store(iso4_a)=( (store(iso4_a)) - &
4073 (electrolyte(jna2so4,jp,ibin)) )
4074 store(ina_a) =max(0.d0, store(ina_a))
4075 store(iso4_a)=max(0.d0, store(iso4_a))
4078 end subroutine form_na2so4
4082 subroutine form_nahso4(store,jp,ibin,electrolyte)
4083 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4084 iso4_a,ina_a,jnahso4
4089 integer, intent(in) :: jp, ibin
4090 real(r8), intent(inout), dimension(naer) :: store
4091 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4093 electrolyte(jnahso4,jp,ibin) = min(store(ina_a), &
4095 store(ina_a) = ( (store(ina_a)) - &
4096 (electrolyte(jnahso4,jp,ibin)) )
4097 store(iso4_a) = ( (store(iso4_a)) - &
4098 (electrolyte(jnahso4,jp,ibin)) )
4099 store(ina_a) = max(0.d0, store(ina_a))
4100 store(iso4_a) = max(0.d0, store(iso4_a))
4103 end subroutine form_nahso4
4107 subroutine form_namsa(store,jp,ibin,electrolyte)
4108 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4113 integer, intent(in) :: jp, ibin
4114 real(r8), intent(inout), dimension(naer) :: store
4115 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4117 electrolyte(jnamsa,jp,ibin) = min(store(ina_a), &
4119 store(ina_a) = ( (store(ina_a)) - &
4120 (electrolyte(jnamsa,jp,ibin)) )
4121 store(imsa_a) = ( (store(imsa_a)) - &
4122 (electrolyte(jnamsa,jp,ibin)) )
4123 store(ina_a) = max(0.d0, store(ina_a))
4124 store(imsa_a) = max(0.d0, store(imsa_a))
4127 end subroutine form_namsa
4131 subroutine form_nano3(store,jp,ibin,electrolyte)
4132 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4137 integer, intent(in) :: jp, ibin
4138 real(r8), intent(inout), dimension(naer) :: store
4139 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4141 electrolyte(jnano3,jp,ibin)=min(store(ina_a),store(ino3_a))
4142 store(ina_a) = ( (store(ina_a)) - &
4143 (electrolyte(jnano3,jp,ibin)) )
4144 store(ino3_a) = ( (store(ino3_a)) - &
4145 (electrolyte(jnano3,jp,ibin)) )
4146 store(ina_a) = max(0.d0, store(ina_a))
4147 store(ino3_a) = max(0.d0, store(ino3_a))
4150 end subroutine form_nano3
4154 subroutine form_cano3(store,jp,ibin,electrolyte) ! Ca(NO3)2
4155 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4160 integer, intent(in) :: jp, ibin
4161 real(r8), intent(inout), dimension(naer) :: store
4162 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4164 electrolyte(jcano3,jp,ibin) = min(store(ica_a),0.5*store(ino3_a))
4166 store(ica_a) = ( (store(ica_a)) - &
4167 (electrolyte(jcano3,jp,ibin)) )
4168 store(ino3_a) = ( (store(ino3_a)) - &
4169 (2.*electrolyte(jcano3,jp,ibin)) )
4170 store(ica_a) = max(0.d0, store(ica_a))
4171 store(ino3_a) = max(0.d0, store(ino3_a))
4174 end subroutine form_cano3
4178 subroutine form_cacl2(store,jp,ibin,electrolyte)
4179 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4185 integer, intent(in) :: jp, ibin
4186 real(r8), intent(inout), dimension(naer) :: store
4187 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4189 electrolyte(jcacl2,jp,ibin) = min(store(ica_a),0.5*store(icl_a))
4191 store(ica_a) = ( (store(ica_a)) - &
4192 (electrolyte(jcacl2,jp,ibin)) )
4193 store(icl_a) = ( (store(icl_a)) - &
4194 (2.*electrolyte(jcacl2,jp,ibin)) )
4195 store(ica_a) = max(0.d0, store(ica_a))
4196 store(icl_a) = max(0.d0, store(icl_a))
4199 end subroutine form_cacl2
4202 subroutine form_caco3(store,jp,ibin,aer,electrolyte)
4203 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,jsolid, &
4210 integer, intent(in) :: jp, ibin
4211 real(r8), intent(inout), dimension(naer) :: store
4212 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4213 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4215 if(jp.eq.jtotal .or. jp.eq.jsolid)then
4216 electrolyte(jcaco3,jp,ibin) = store(ica_a)
4218 aer(ico3_a,jp,ibin)= electrolyte(jcaco3,jp,ibin) ! force co3 = caco3
4225 end subroutine form_caco3
4229 subroutine form_nacl(store,jp,ibin,aer,gas,electrolyte,total_species,tot_cl_in)
4230 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4231 ngas_aerchtot,ngas_volatile,jtotal,jsolid,jliquid, &
4232 ina_a,jnacl,icl_a,ihcl_g
4237 integer, intent(in) :: jp, ibin
4238 real(r8), intent(inout), dimension(naer) :: store
4239 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4240 real(r8), intent(inout), dimension(ngas_volatile) :: total_species
4241 real(r8), intent(inout) :: tot_cl_in
4242 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4243 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4246 electrolyte(jnacl,jp,ibin) = store(ina_a)
4249 store(icl_a) = ( (store(icl_a)) - &
4250 (electrolyte(jnacl,jp,ibin)) )
4252 if(store(icl_a) .lt. 0.)then ! cl deficit in aerosol. take some from gas
4253 aer(icl_a,jp,ibin)= aer(icl_a,jp,ibin)- store(icl_a) ! update aer(icl_a)
4255 if(jp .ne. jtotal)then
4256 aer(icl_a,jtotal,ibin)= aer(icl_a,jliquid,ibin)+ & ! update for jtotal
4257 aer(icl_a,jsolid,ibin)
4260 gas(ihcl_g) = gas(ihcl_g) + store(icl_a) ! update gas(ihcl_g)
4262 if(gas(ihcl_g) .lt. 0.0)then
4263 total_species(ihcl_g) = total_species(ihcl_g) - gas(ihcl_g) ! update total_species
4264 tot_cl_in = tot_cl_in - gas(ihcl_g) ! update tot_cl_in
4267 gas(ihcl_g) = max(0.d0, gas(ihcl_g)) ! restrict gas(ihcl_g) to >= 0.
4268 store(icl_a) = 0. ! force store(icl_a) to 0.
4272 store(icl_a) = max(0.d0, store(icl_a))
4275 end subroutine form_nacl
4279 subroutine form_nh4so4(store,jp,ibin,electrolyte) ! (nh4)2so4
4280 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4281 iso4_a,inh4_a,jnh4so4
4285 integer, intent(in) :: jp, ibin
4286 real(r8), intent(inout), dimension(naer) :: store
4287 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4289 electrolyte(jnh4so4,jp,ibin)= min(.5*store(inh4_a), &
4291 store(inh4_a)= ( (store(inh4_a)) - &
4292 (2.*electrolyte(jnh4so4,jp,ibin)) )
4293 store(iso4_a)= ( (store(iso4_a)) - &
4294 (electrolyte(jnh4so4,jp,ibin)) )
4295 store(inh4_a) = max(0.d0, store(inh4_a))
4296 store(iso4_a) = max(0.d0, store(iso4_a))
4299 end subroutine form_nh4so4
4303 subroutine form_nh4hso4(store,jp,ibin,electrolyte) ! nh4hso4
4304 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4305 iso4_a,inh4_a,jnh4hso4
4309 integer, intent(in) :: jp, ibin
4310 real(r8), intent(inout), dimension(naer) :: store
4311 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4313 electrolyte(jnh4hso4,jp,ibin) = min(store(inh4_a), &
4315 store(inh4_a)= ( (store(inh4_a)) - &
4316 (electrolyte(jnh4hso4,jp,ibin)) )
4317 store(iso4_a)= ( (store(iso4_a)) - &
4318 (electrolyte(jnh4hso4,jp,ibin)) )
4319 store(inh4_a) = max(0.d0, store(inh4_a))
4320 store(iso4_a) = max(0.d0, store(iso4_a))
4323 end subroutine form_nh4hso4
4327 subroutine form_nh4msa(store,jp,ibin,electrolyte)
4328 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4329 imsa_a,inh4_a,jnh4msa
4333 integer, intent(in) :: jp, ibin
4334 real(r8), intent(inout), dimension(naer) :: store
4335 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4337 electrolyte(jnh4msa,jp,ibin) = min(store(inh4_a), &
4339 store(inh4_a) = ( (store(inh4_a)) - &
4340 (electrolyte(jnh4msa,jp,ibin)) )
4341 store(imsa_a) = ( (store(imsa_a)) - &
4342 (electrolyte(jnh4msa,jp,ibin)) )
4343 store(inh4_a) = max(0.d0, store(inh4_a))
4344 store(imsa_a) = max(0.d0, store(imsa_a))
4347 end subroutine form_nh4msa
4351 subroutine form_nh4cl(store,jp,ibin,electrolyte)
4352 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4357 integer, intent(in) :: jp, ibin
4358 real(r8), intent(inout), dimension(naer) :: store
4359 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4361 electrolyte(jnh4cl,jp,ibin) = min(store(inh4_a), &
4363 store(inh4_a) = ( (store(inh4_a)) - &
4364 (electrolyte(jnh4cl,jp,ibin)) )
4365 store(icl_a) = ( (store(icl_a)) - &
4366 (electrolyte(jnh4cl,jp,ibin)) )
4367 store(inh4_a) = max(0.d0, store(inh4_a))
4368 store(icl_a) = max(0.d0, store(icl_a))
4371 end subroutine form_nh4cl
4375 subroutine form_nh4no3(store,jp,ibin,electrolyte)
4376 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4377 ino3_a,inh4_a,jnh4no3
4381 integer, intent(in) :: jp, ibin
4382 real(r8), intent(inout), dimension(naer) :: store
4383 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4385 electrolyte(jnh4no3,jp,ibin) = min(store(inh4_a), &
4387 store(inh4_a) = ( (store(inh4_a)) - &
4388 (electrolyte(jnh4no3,jp,ibin)) )
4389 store(ino3_a) = ( (store(ino3_a)) - &
4390 (electrolyte(jnh4no3,jp,ibin)) )
4391 store(inh4_a) = max(0.d0, store(inh4_a))
4392 store(ino3_a) = max(0.d0, store(ino3_a))
4395 end subroutine form_nh4no3
4399 subroutine form_nh4so4_lvcite(store,jp,ibin,electrolyte) ! (nh4)2so4 + (nh4)3h(so4)2
4400 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4401 iso4_a,inh4_a,jnh4so4,jlvcite
4405 integer, intent(in) :: jp, ibin
4406 real(r8), intent(inout), dimension(naer) :: store
4407 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4409 electrolyte(jnh4so4,jp,ibin)= ( (2.*store(inh4_a)) - &
4410 (3.*store(iso4_a)) )
4411 electrolyte(jlvcite,jp,ibin)= ( (2.*store(iso4_a)) - &
4413 electrolyte(jnh4so4,jp,ibin)= max(0.d0, &
4414 electrolyte(jnh4so4,jp,ibin))
4415 electrolyte(jlvcite,jp,ibin)= max(0.d0, &
4416 electrolyte(jlvcite,jp,ibin))
4421 end subroutine form_nh4so4_lvcite
4425 subroutine form_lvcite_nh4hso4(store,jp,ibin,electrolyte) ! (nh4)3h(so4)2 + nh4hso4
4426 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4427 iso4_a,inh4_a,jlvcite,jnh4hso4
4431 integer, intent(in) :: jp, ibin
4432 real(r8), intent(inout), dimension(naer) :: store
4433 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4435 electrolyte(jlvcite,jp,ibin) = ( (store(inh4_a)) - &
4437 electrolyte(jnh4hso4,jp,ibin)= ( (3.*store(iso4_a)) - &
4438 (2.*store(inh4_a)) )
4439 electrolyte(jlvcite,jp,ibin) = max(0.d0, &
4440 electrolyte(jlvcite,jp,ibin))
4441 electrolyte(jnh4hso4,jp,ibin)= max(0.d0, &
4442 electrolyte(jnh4hso4,jp,ibin))
4447 end subroutine form_lvcite_nh4hso4
4451 subroutine form_na2so4_nahso4(store,jp,ibin,electrolyte) ! na2so4 + nahso4
4452 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4453 iso4_a,ina_a,jna2so4,jnahso4
4457 integer, intent(in) :: jp, ibin
4458 real(r8), intent(inout), dimension(naer) :: store
4459 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4461 electrolyte(jna2so4,jp,ibin)= ( (store(ina_a)) - &
4463 electrolyte(jnahso4,jp,ibin)= ( (2.*store(iso4_a))- &
4465 electrolyte(jna2so4,jp,ibin)= max(0.d0, &
4466 electrolyte(jna2so4,jp,ibin))
4467 electrolyte(jnahso4,jp,ibin)= max(0.d0, &
4468 electrolyte(jnahso4,jp,ibin))
4472 ! write(6,*)'na2so4 + nahso4'
4475 end subroutine form_na2so4_nahso4
4479 subroutine form_h2so4(store,jp,ibin,electrolyte)
4480 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4485 integer, intent(in) :: jp, ibin
4486 real(r8), intent(inout), dimension(naer) :: store
4487 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4489 electrolyte(jh2so4,jp,ibin) = max(0.0d0, store(iso4_a))
4493 end subroutine form_h2so4
4497 subroutine form_msa(store,jp,ibin,electrolyte)
4498 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4503 integer, intent(in) :: jp, ibin
4504 real(r8), intent(inout), dimension(naer) :: store
4505 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4507 electrolyte(jmsa,jp,ibin) = max(0.0d0, store(imsa_a))
4511 end subroutine form_msa
4515 subroutine form_hno3(store,jp,ibin,electrolyte)
4516 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4521 integer, intent(in) :: jp, ibin
4522 real(r8), intent(inout), dimension(naer) :: store
4523 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4525 electrolyte(jhno3,jp,ibin) = max(0.0d0, store(ino3_a))
4529 end subroutine form_hno3
4533 subroutine form_hcl(store,jp,ibin,electrolyte)
4534 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4539 integer, intent(in) :: jp, ibin
4540 real(r8), intent(inout), dimension(naer) :: store
4541 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4543 electrolyte(jhcl,jp,ibin) = max(0.0d0, store(icl_a))
4547 end subroutine form_hcl
4551 subroutine degas_hno3(store,jp,ibin,aer,gas,electrolyte)
4552 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4553 ngas_aerchtot,jtotal,jliquid,jsolid, &
4554 ino3_a,ihno3_g,jhno3
4559 integer, intent(in) :: jp, ibin
4560 real(r8), intent(inout), dimension(naer) :: store
4561 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4562 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4563 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4565 store(ino3_a) = max(0.0d0, store(ino3_a))
4566 gas(ihno3_g) = gas(ihno3_g) + store(ino3_a)
4567 aer(ino3_a,jp,ibin) = ( (aer(ino3_a,jp,ibin)) - &
4569 aer(ino3_a,jp,ibin) = max(0.0d0,aer(ino3_a,jp,ibin))
4571 ! also do it for jtotal
4572 if(jp .ne. jtotal)then
4573 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid, ibin) + &
4574 aer(ino3_a,jliquid,ibin)
4577 electrolyte(jhno3,jp,ibin) = 0.0
4581 end subroutine degas_hno3
4585 subroutine degas_hcl(store,jp,ibin,aer,gas,electrolyte)
4586 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4587 ngas_aerchtot,jtotal,jliquid,jsolid, &
4592 integer, intent(in) :: jp, ibin
4593 real(r8), intent(inout) :: store(naer)
4594 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4595 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4596 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4598 store(icl_a) = max(0.0d0, store(icl_a))
4599 gas(ihcl_g) = gas(ihcl_g) + store(icl_a)
4600 aer(icl_a,jp,ibin) = ( (aer(icl_a,jp,ibin)) - &
4602 aer(icl_a,jp,ibin) = max(0.0d0,aer(icl_a,jp,ibin))
4604 ! also do it for jtotal
4605 if(jp .ne. jtotal)then
4606 aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid, ibin) + &
4607 aer(icl_a,jliquid,ibin)
4610 electrolyte(jhcl,jp,ibin) = 0.0
4614 end subroutine degas_hcl
4618 subroutine degas_nh3(store,jp,ibin,aer,gas)
4619 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4620 ngas_aerchtot,jtotal,jliquid,jsolid, &
4626 integer, intent(in) :: jp, ibin
4627 real(r8), intent(inout) :: store(naer)
4628 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4629 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4631 store(inh4_a) = max(0.0d0, store(inh4_a))
4632 gas(inh3_g) = gas(inh3_g) + store(inh4_a)
4633 aer(inh4_a,jp,ibin) = ( (aer(inh4_a,jp,ibin)) - &
4635 aer(inh4_a,jp,ibin) = max(0.0d0,aer(inh4_a,jp,ibin))
4637 ! also do it for jtotal
4638 if(jp .ne. jtotal)then
4639 aer(inh4_a,jtotal,ibin)= aer(inh4_a,jsolid, ibin) + &
4640 aer(inh4_a,jliquid,ibin)
4646 end subroutine degas_nh3
4650 !***********************************************************************
4651 ! subroutines to absorb and degas small amounts of volatile species
4653 ! author: Rahul A. Zaveri
4655 !-----------------------------------------------------------------------
4658 subroutine absorb_tiny_nh4no3(ibin,aer,gas,electrolyte,delta_nh3_max, &
4659 delta_hno3_max,electrolyte_sum)
4660 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4661 ngas_aerchtot,jtotal,jliquid,jsolid, &
4662 inh4_a,ino3_a,inh3_g,ihno3_g
4666 integer, intent(in) :: ibin
4667 real(r8), intent(inout), dimension(nbin_a_max) :: delta_nh3_max,delta_hno3_max
4668 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4669 real(r8), intent(inout), dimension(3,nbin_a_max) :: electrolyte_sum
4670 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4671 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4674 real(r8) :: small_aer, small_gas, small_amt
4679 electrolyte_sum(jtotal,ibin) = 0.0
4680 do je = 1, nelectrolyte
4681 electrolyte_sum(jtotal,ibin) = electrolyte_sum(jtotal,ibin) + &
4682 electrolyte(je,jtotal,ibin)
4687 small_gas = 0.01 * min(delta_nh3_max(ibin),delta_hno3_max(ibin))
4688 small_aer = 0.01 * electrolyte_sum(jtotal,ibin)
4689 if(small_aer .eq. 0.0)small_aer = small_gas
4691 small_amt = min(small_gas, small_aer)
4693 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + small_amt
4694 aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) + small_amt
4697 aer(inh4_a,jtotal,ibin) = aer(inh4_a,jsolid,ibin) + &
4698 aer(inh4_a,jliquid,ibin)
4699 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid,ibin) + &
4700 aer(ino3_a,jliquid,ibin)
4703 gas(inh3_g) = ((gas(inh3_g)) - (small_amt))
4704 gas(ihno3_g) = ((gas(ihno3_g)) - (small_amt))
4707 end subroutine absorb_tiny_nh4no3
4711 !--------------------------------------------------------------------
4713 subroutine absorb_tiny_nh4cl(ibin,aer,gas,electrolyte,delta_nh3_max, &
4714 delta_hcl_max,electrolyte_sum)
4715 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4716 ngas_aerchtot,jtotal,jliquid,jsolid, &
4717 inh4_a,icl_a,inh3_g,ihcl_g
4721 integer, intent(in) :: ibin
4722 real(r8), intent(inout), dimension(nbin_a_max) :: delta_nh3_max,delta_hcl_max
4723 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4724 real(r8), intent(inout), dimension(3,nbin_a_max) :: electrolyte_sum
4725 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4726 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4729 real(r8) :: small_aer, small_gas, small_amt
4733 electrolyte_sum(jtotal,ibin) = 0.0
4734 do je = 1, nelectrolyte
4735 electrolyte_sum(jtotal,ibin) = electrolyte_sum(jtotal,ibin) + &
4736 electrolyte(je,jtotal,ibin)
4742 small_gas = 0.01 * min(delta_nh3_max(ibin), delta_hcl_max(ibin))
4743 small_aer = 0.01 * electrolyte_sum(jtotal,ibin)
4744 if(small_aer .eq. 0.0)small_aer = small_gas
4746 small_amt = min(small_gas, small_aer)
4748 aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + small_amt
4749 aer(icl_a,jliquid,ibin) = aer(icl_a,jliquid,ibin) + small_amt
4752 aer(inh4_a,jtotal,ibin) = aer(inh4_a,jsolid,ibin) + &
4753 aer(inh4_a,jliquid,ibin)
4754 aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid,ibin) + &
4755 aer(icl_a,jliquid,ibin)
4758 gas(inh3_g) = ((gas(inh3_g)) - (small_amt))
4759 gas(ihcl_g) = ((gas(ihcl_g)) - (small_amt))
4762 end subroutine absorb_tiny_nh4cl
4766 !--------------------------------------------------------------------
4768 subroutine absorb_tiny_hno3(ibin,aer,gas,delta_hno3_max) ! and degas tiny hcl
4769 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4770 ngas_aerchtot,jliquid,jsolid,jtotal, &
4771 icl_a,ino3_a,ihno3_g,ihcl_g
4775 integer, intent(in) :: ibin
4776 real(r8), intent(inout), dimension(nbin_a_max) :: delta_hno3_max
4777 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4778 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4780 real(r8) :: small_aer, small_amt, small_gas
4782 small_gas = 0.01 * delta_hno3_max(ibin)
4783 small_aer = 0.01 * aer(icl_a,jliquid,ibin)
4785 small_amt = min(small_gas, small_aer)
4788 aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) + small_amt
4789 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid,ibin) + &
4790 aer(ino3_a,jliquid,ibin)
4791 gas(ihno3_g) = ((gas(ihno3_g))-(small_amt))
4794 aer(icl_a,jliquid,ibin) = ((aer(icl_a,jliquid,ibin))- &
4796 aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid,ibin) + &
4797 aer(icl_a,jliquid,ibin)
4800 gas(ihcl_g) = gas(ihcl_g) + small_amt
4803 end subroutine absorb_tiny_hno3
4807 !--------------------------------------------------------------
4809 subroutine absorb_tiny_hcl(ibin,aer,gas,delta_hcl_max) ! and degas tiny hno3
4810 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4811 ngas_aerchtot,jliquid,jtotal,jsolid, &
4812 ino3_a,icl_a,ihcl_g,ihno3_g
4816 integer, intent(in) :: ibin
4817 real(r8), intent(inout), dimension(nbin_a_max) :: delta_hcl_max
4818 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4819 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4821 real(r8) :: small_aer, small_amt, small_gas
4823 small_gas = 0.01 * delta_hcl_max(ibin)
4824 small_aer = 0.01 * aer(ino3_a,jliquid,ibin)
4826 small_amt = min(small_gas, small_aer)
4829 aer(icl_a,jliquid,ibin)= aer(icl_a,jliquid,ibin) + small_amt
4830 aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid,ibin) + &
4831 aer(icl_a,jliquid,ibin)
4832 gas(ihcl_g) = ((gas(ihcl_g))-(small_amt))
4835 aer(ino3_a,jliquid,ibin) = ((aer(ino3_a,jliquid,ibin))- &
4837 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid,ibin) + &
4838 aer(ino3_a,jliquid,ibin)
4841 gas(ihno3_g) = gas(ihno3_g) + small_amt
4844 end subroutine absorb_tiny_hcl
4848 !--------------------------------------------------------------
4850 subroutine degas_tiny_nh4no3(ibin,aer,gas,electrolyte)
4851 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4852 ngas_aerchtot,jliquid,jsolid,jtotal, &
4853 jnh4no3,inh4_a,ino3_a,inh3_g,ihno3_g
4857 integer, intent(in) :: ibin
4858 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4859 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4860 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4862 real(r8) :: small_amt
4864 small_amt = 0.01 * electrolyte(jnh4no3,jliquid,ibin)
4866 aer(inh4_a,jliquid,ibin) = ((aer(inh4_a,jliquid,ibin))- &
4868 aer(ino3_a,jliquid,ibin) = ((aer(ino3_a,jliquid,ibin))- &
4872 aer(inh4_a,jtotal,ibin) = aer(inh4_a,jsolid,ibin) + &
4873 aer(inh4_a,jliquid,ibin)
4874 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid,ibin) + &
4875 aer(ino3_a,jliquid,ibin)
4878 gas(inh3_g) = gas(inh3_g) + small_amt
4879 gas(ihno3_g) = gas(ihno3_g) + small_amt
4882 end subroutine degas_tiny_nh4no3
4887 !--------------------------------------------------------------------
4889 subroutine degas_tiny_nh4cl(ibin,aer,gas,electrolyte)
4890 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4891 ngas_aerchtot,jliquid,jsolid,jtotal, &
4892 jnh4cl,inh4_a,icl_a,inh3_g,ihcl_g
4896 integer, intent(in) :: ibin
4897 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4898 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4899 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4901 real(r8) :: small_amt
4904 small_amt = 0.01 * electrolyte(jnh4cl,jliquid,ibin)
4906 aer(inh4_a,jliquid,ibin) = ((aer(inh4_a,jliquid,ibin))- &
4908 aer(icl_a,jliquid,ibin) = ((aer(icl_a,jliquid,ibin))- &
4912 aer(inh4_a,jtotal,ibin) = aer(inh4_a,jsolid,ibin) + &
4913 aer(inh4_a,jliquid,ibin)
4914 aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid,ibin) + &
4915 aer(icl_a,jliquid,ibin)
4918 gas(inh3_g) = gas(inh3_g) + small_amt
4919 gas(ihcl_g) = gas(ihcl_g) + small_amt
4922 end subroutine degas_tiny_nh4cl
4926 !***********************************************************************
4927 ! subroutines to equilibrate volatile acids
4929 ! author: Rahul A. Zaveri
4931 !-----------------------------------------------------------------------
4932 subroutine equilibrate_acids(ibin,aer,gas,electrolyte,activity,mc,water_a, &
4933 total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
4934 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4935 ngas_aerchtot,ngas_volatile,Ncation,Nanion,nrxn_aer_gl,nrxn_aer_ll, &
4940 integer, intent(in) :: ibin
4941 real(r8), intent(inout), dimension(nbin_a_max) :: water_a
4942 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4943 real(r8), intent(inout), dimension(ngas_volatile) :: total_species
4944 real(r8), intent(inout) :: tot_cl_in
4945 real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
4946 real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
4947 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
4948 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
4949 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
4950 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4951 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4956 if(gas(ihcl_g)*gas(ihno3_g) .gt. 0.)then
4957 call equilibrate_hcl_and_hno3(ibin,aer,gas,electrolyte,activity,mc,water_a, &
4958 total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
4959 elseif(gas(ihcl_g) .gt. 0.)then
4960 call equilibrate_hcl(ibin,aer,gas,electrolyte,activity,mc,water_a, &
4961 total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
4962 elseif(gas(ihno3_g) .gt. 0.)then
4963 call equilibrate_hno3(ibin,aer,gas,electrolyte,activity,mc,water_a, &
4964 total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
4969 end subroutine equilibrate_acids
4974 subroutine equilibrate_hcl(ibin,aer,gas,electrolyte,activity,mc,water_a, &
4975 total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
4976 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
4977 ngas_aerchtot,ngas_volatile,Ncation,jliquid,jsolid,jtotal,Nanion, &
4978 nrxn_aer_gl,nrxn_aer_ll, &
4979 ja_so4,ja_hso4,ihcl_g,icl_a,jhcl,ino3_a,ica_a,inh4_a,ina_a,jc_h,jc_ca, &
4980 jc_nh4,jc_na,ja_cl,ja_no3,jhno3,jnh4cl
4985 integer, intent(in) :: ibin
4986 real(r8), intent(inout), dimension(nbin_a_max) :: water_a
4987 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
4988 real(r8), intent(inout), dimension(ngas_volatile) :: total_species
4989 real(r8), intent(inout) :: tot_cl_in
4990 real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
4991 real(r8), intent(inout), dimension(nrxn_aer_gl) ::Keq_gl
4992 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
4993 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
4994 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
4995 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
4996 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
4998 real(r8) :: a, aerH, aerHSO4, aerSO4, b, c, dum, Kdash_hcl, mH, Tcl, &
5000 !real(r8) :: quadratic ! mosaic func
5002 aerSO4 = ma(ja_so4,ibin)*water_a(ibin)*1.e+9
5003 aerHSO4= ma(ja_hso4,ibin)*water_a(ibin)*1.e+9
5005 Tcl = aer(icl_a,jliquid,ibin) + gas(ihcl_g) ! nmol/m^3(air)
5006 Kdash_hcl = Keq_gl(4)*1.e+18/gam(jhcl,ibin)**2 ! (nmol^2/kg^2)/(nmol/m^3(air))
5007 Z = ( aer(ina_a, jliquid,ibin) + & ! nmol/m^3(air)
5008 aer(inh4_a,jliquid,ibin) + &
5009 2.*aer(ica_a, jliquid,ibin) ) - &
5012 aer(ino3_a,jliquid,ibin) )
5015 W = water_a(ibin) ! kg/m^3(air)
5017 Kdash_hcl = Keq_gl(4)*1.e+18/gam(jhcl,ibin)**2 ! (nmol^2/kg^2)/(nmol/m^3(air))
5019 b = ((Kdash_hcl*W) + (Z/W))*1.e-9
5020 c = Kdash_hcl*(Z - Tcl)*1.e-18
5023 dum = ((b*b)-(4.*a*c))
5024 if (dum .lt. 0.) return ! no real root
5028 mH = quadratic(a,b,c) ! mol/kg(water)
5030 aer(icl_a,jliquid,ibin) = ((aerH) + (Z))
5032 mH = sqrt(Keq_ll(3))
5035 call form_electrolytes(jliquid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)
5037 ! update gas phase concentration
5038 gas(ihcl_g) = ( (Tcl) - (aer(icl_a,jliquid,ibin)) )
5041 ! update the following molalities
5042 ma(ja_so4,ibin) = 1.e-9*aerSO4/water_a(ibin)
5043 ma(ja_hso4,ibin) = 1.e-9*aerHSO4/water_a(ibin)
5044 ma(ja_no3,ibin) = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
5045 ma(ja_cl,ibin) = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)
5048 mc(jc_ca,ibin) = 1.e-9*aer(ica_a, jliquid,ibin)/water_a(ibin)
5049 mc(jc_nh4,ibin) = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
5050 mc(jc_na,ibin) = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)
5053 ! update the following activities
5054 activity(jhcl,ibin) = mc(jc_h,ibin) *ma(ja_cl,ibin) * &
5057 activity(jhno3,ibin) = mc(jc_h,ibin) *ma(ja_no3,ibin) * &
5060 activity(jnh4cl,ibin) = mc(jc_nh4,ibin)*ma(ja_cl,ibin) * &
5064 ! also update xyz(jtotal)
5065 aer(icl_a,jtotal,ibin) = aer(icl_a,jliquid,ibin) + &
5066 aer(icl_a,jsolid,ibin)
5068 electrolyte(jhcl,jtotal,ibin) = electrolyte(jhcl,jliquid,ibin)
5071 end subroutine equilibrate_hcl
5076 subroutine equilibrate_hno3(ibin,aer,gas,electrolyte,activity,mc,water_a, &
5077 total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
5078 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
5079 ngas_aerchtot,ngas_volatile,Ncation,jliquid,jsolid,jtotal,Nanion, &
5080 nrxn_aer_gl,nrxn_aer_ll, &
5081 ja_so4,ja_hso4,ihno3_g,ino3_a,jhno3,icl_a,ica_a,inh4_a,ina_a,jc_h,jc_ca, &
5082 jc_nh4,jc_na,ja_cl,jhcl,ja_no3,jnh4no3
5087 integer, intent(in) :: ibin
5088 real(r8), intent(inout), dimension(nbin_a_max) :: water_a
5089 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
5090 real(r8), intent(inout), dimension(ngas_volatile) :: total_species
5091 real(r8), intent(inout) :: tot_cl_in
5092 real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
5093 real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
5094 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
5095 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
5096 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
5097 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
5098 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
5100 real(r8) :: a, aerH, aerHSO4, aerSO4, b, c, dum, Kdash_hno3, mH, &
5102 !real(r8) :: quadratic ! mosaic func
5104 aerSO4 = ma(ja_so4,ibin)*water_a(ibin)*1.e+9
5105 aerHSO4= ma(ja_hso4,ibin)*water_a(ibin)*1.e+9
5107 Tno3 = aer(ino3_a,jliquid,ibin) + gas(ihno3_g) ! nmol/m^3(air)
5108 Kdash_hno3 = Keq_gl(3)*1.e+18/gam(jhno3,ibin)**2 ! (nmol^2/kg^2)/(nmol/m^3(air))
5109 Z = ( aer(ina_a, jliquid,ibin) + & ! nmol/m^3(air)
5110 aer(inh4_a,jliquid,ibin) + &
5111 2.*aer(ica_a, jliquid,ibin) ) - &
5114 aer(icl_a,jliquid,ibin) )
5117 W = water_a(ibin) ! kg/m^3(air)
5119 Kdash_hno3 = Keq_gl(3)*1.e+18/gam(jhno3,ibin)**2 ! (nmol^2/kg^2)/(nmol/m^3(air))
5121 b = ((Kdash_hno3*W) + (Z/W))*1.e-9
5122 c = Kdash_hno3*(Z - Tno3)*1.e-18
5124 dum = ((b*b)-(4.*a*c))
5125 if (dum .lt. 0.) return ! no real root
5130 mH = quadratic(a,b,c) ! mol/kg(water)
5132 aer(ino3_a,jliquid,ibin) = ((aerH) + (Z))
5134 mH = sqrt(Keq_ll(3))
5137 call form_electrolytes(jliquid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)
5139 ! update gas phase concentration
5140 gas(ihno3_g)= ( (Tno3) - (aer(ino3_a,jliquid,ibin)) )
5143 ! update the following molalities
5144 ma(ja_so4,ibin) = 1.e-9*aerSO4/water_a(ibin)
5145 ma(ja_hso4,ibin) = 1.e-9*aerHSO4/water_a(ibin)
5146 ma(ja_no3,ibin) = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
5147 ma(ja_cl,ibin) = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)
5150 mc(jc_ca,ibin) = 1.e-9*aer(ica_a, jliquid,ibin)/water_a(ibin)
5151 mc(jc_nh4,ibin) = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
5152 mc(jc_na,ibin) = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)
5155 ! update the following activities
5156 activity(jhcl,ibin) = mc(jc_h,ibin) *ma(ja_cl,ibin) * &
5159 activity(jhno3,ibin) = mc(jc_h,ibin) *ma(ja_no3,ibin) * &
5162 activity(jnh4no3,ibin) = mc(jc_nh4,ibin)*ma(ja_no3,ibin) * &
5163 gam(jnh4no3,ibin)**2
5166 ! also update xyz(jtotal)
5167 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jliquid,ibin) + &
5168 aer(ino3_a,jsolid,ibin)
5170 electrolyte(jhno3,jtotal,ibin) = electrolyte(jhno3,jliquid,ibin)
5173 end subroutine equilibrate_hno3
5178 subroutine equilibrate_hcl_and_hno3(ibin,aer,gas,electrolyte,activity,mc, &
5179 water_a,total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
5180 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
5181 ngas_aerchtot,ngas_volatile,Ncation,jliquid,jsolid,jtotal,Nanion, &
5182 nrxn_aer_gl,nrxn_aer_ll, &
5183 ja_so4,ja_hso4,ihcl_g,icl_a,ihno3_g,ino3_a,jhcl,jhno3, &
5184 ica_a,inh4_a,ina_a,jc_h,jc_ca,jc_nh4,jc_na,ja_cl,ja_no3,jnh4no3, &
5190 integer, intent(in) :: ibin
5191 real(r8), intent(inout), dimension(nbin_a_max) :: water_a
5192 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
5193 real(r8), intent(inout), dimension(ngas_volatile) :: total_species
5194 real(r8), intent(inout) :: tot_cl_in
5195 real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
5196 real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
5197 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
5198 real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
5199 real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
5200 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
5201 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
5203 real(r8) :: aerH, aerHSO4, aerSO4, Kdash_hcl, Kdash_hno3, &
5204 mH, p, q, r, Tcl, Tno3, W, XT, Z
5205 !real(r8) :: cubic ! mosaic func
5208 aerSO4 = ma(ja_so4,ibin)*water_a(ibin)*1.e+9
5209 aerHSO4= ma(ja_hso4,ibin)*water_a(ibin)*1.e+9
5211 Tcl = aer(icl_a,jliquid,ibin) + gas(ihcl_g) ! nmol/m^3(air)
5212 Tno3 = aer(ino3_a,jliquid,ibin) + gas(ihno3_g) ! nmol/m^3(air)
5214 Kdash_hcl = Keq_gl(4)*1.e+18/gam(jhcl,ibin)**2 ! (nmol^2/kg^2)/(nmol/m^3(air))
5215 Kdash_hno3 = Keq_gl(3)*1.e+18/gam(jhno3,ibin)**2 ! (nmol^2/kg^2)/(nmol/m^3(air))
5217 Z = ( aer(ina_a, jliquid,ibin) + & ! nmol/m^3(air)
5218 aer(inh4_a,jliquid,ibin) + &
5219 2.*aer(ica_a, jliquid,ibin) ) - &
5220 (2.*aerSO4 + aerHSO4 )
5225 Kdash_hcl = Keq_gl(4)*1.e+18/gam(jhcl,ibin)**2 ! (nmol^2/kg^2)/(nmol/m^3(air))
5226 Kdash_hno3 = Keq_gl(3)*1.e+18/gam(jhno3,ibin)**2 ! (nmol^2/kg^2)/(nmol/m^3(air))
5228 p = (Z/W + W*(Kdash_hcl + Kdash_hno3))*1.e-9
5230 q = 1.e-18*Kdash_hcl*Kdash_hno3*W**2 + &
5231 1.e-18*Z*(Kdash_hcl + Kdash_hno3) - &
5232 1.e-18*Kdash_hcl*Tcl - &
5233 1.e-18*Kdash_hno3*Tno3
5235 r = 1.e-18*Kdash_hcl*Kdash_hno3*W*(Z - Tcl - Tno3)*1.e-9
5241 aer(ino3_a,jliquid,ibin) = Kdash_hno3*W*W*Tno3/ &
5242 (aerH + Kdash_hno3*W*W)
5243 aer(icl_a, jliquid,ibin) = Kdash_hcl*W*W*Tcl/ &
5244 (aerH + Kdash_hcl*W*W)
5246 mH = sqrt(Keq_ll(3))
5249 call form_electrolytes(jliquid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)
5251 ! update gas phase concentration
5252 gas(ihno3_g)= ( (Tno3) - (aer(ino3_a,jliquid,ibin)) )
5253 gas(ihcl_g) = ( (Tcl) - (aer(icl_a,jliquid,ibin)) )
5256 ! update the following molalities
5257 ma(ja_so4,ibin) = 1.e-9*aerSO4/water_a(ibin)
5258 ma(ja_hso4,ibin) = 1.e-9*aerHSO4/water_a(ibin)
5259 ma(ja_no3,ibin) = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
5260 ma(ja_cl,ibin) = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)
5263 mc(jc_ca,ibin) = 1.e-9*aer(ica_a, jliquid,ibin)/water_a(ibin)
5264 mc(jc_nh4,ibin) = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
5265 mc(jc_na,ibin) = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)
5268 ! update the following activities
5269 activity(jhcl,ibin) = mc(jc_h,ibin)*ma(ja_cl,ibin) * &
5272 activity(jhno3,ibin) = mc(jc_h,ibin)*ma(ja_no3,ibin) * &
5275 activity(jnh4no3,ibin) = mc(jc_nh4,ibin)*ma(ja_no3,ibin)* &
5276 gam(jnh4no3,ibin)**2
5278 activity(jnh4cl,ibin) = mc(jc_nh4,ibin)*ma(ja_cl,ibin) * &
5282 ! also update xyz(jtotal)
5283 aer(icl_a,jtotal,ibin) = aer(icl_a,jliquid,ibin) + &
5284 aer(icl_a,jsolid,ibin)
5286 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jliquid,ibin) + &
5287 aer(ino3_a,jsolid,ibin)
5289 electrolyte(jhno3,jtotal,ibin) = electrolyte(jhno3,jliquid,ibin)
5290 electrolyte(jhcl, jtotal,ibin) = electrolyte(jhcl, jliquid,ibin)
5293 end subroutine equilibrate_hcl_and_hno3
5297 !***********************************************************************
5298 ! subroutines to evaporate solid volatile species
5300 ! author: Rahul A. Zaveri
5302 !-----------------------------------------------------------------------
5305 subroutine degas_solid_nh4no3(ibin,aer,gas,electrolyte,Keq_sg)
5306 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
5307 ngas_aerchtot,jsolid,jliquid,jtotal,nrxn_aer_sg, &
5308 ihno3_g,inh3_g,jnh4no3,inh4_a,ino3_a
5313 integer, intent(in) :: ibin
5314 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
5315 real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
5316 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
5317 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
5320 real(r8) :: a, b, c, xgas, XT
5321 !real(r8) :: quadratic ! mosaic func
5327 b = gas(inh3_g) + gas(ihno3_g)
5328 c = gas(inh3_g)*gas(ihno3_g) - Keq_sg(1)
5329 xgas = quadratic(a,b,c)
5331 if(xgas .ge. electrolyte(jnh4no3,jp,ibin))then ! degas all nh4no3
5333 gas(inh3_g) = gas(inh3_g) + electrolyte(jnh4no3,jp,ibin)
5334 gas(ihno3_g)= gas(ihno3_g) + electrolyte(jnh4no3,jp,ibin)
5335 aer(inh4_a,jp,ibin) = aer(inh4_a,jp,ibin) - &
5336 electrolyte(jnh4no3,jp,ibin)
5337 aer(ino3_a,jp,ibin) = aer(ino3_a,jp,ibin) - &
5338 electrolyte(jnh4no3,jp,ibin)
5340 else ! degas only xgas amount of nh4no3
5342 gas(inh3_g) = gas(inh3_g) + xgas
5343 gas(ihno3_g)= gas(ihno3_g) + xgas
5344 aer(inh4_a,jp,ibin) = aer(inh4_a,jp,ibin) - xgas
5345 aer(ino3_a,jp,ibin) = aer(ino3_a,jp,ibin) - xgas
5350 aer(inh4_a,jtotal,ibin) = aer(inh4_a,jsolid,ibin) + &
5351 aer(inh4_a,jliquid,ibin)
5352 aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid,ibin) + &
5353 aer(ino3_a,jliquid,ibin)
5356 end subroutine degas_solid_nh4no3
5361 subroutine degas_solid_nh4cl(ibin,aer,gas,electrolyte,Keq_sg)
5362 use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max, &
5363 ngas_aerchtot,jsolid,jliquid,jtotal,nrxn_aer_sg, &
5364 ihcl_g,inh3_g,jnh4cl,inh4_a,icl_a
5368 integer, intent(in) :: ibin
5369 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
5370 real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
5371 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
5372 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
5375 real(r8) :: a, b, c, xgas, XT
5376 !real(r8) :: quadratic ! mosaic func
5382 b = gas(inh3_g) + gas(ihcl_g)
5383 c = gas(inh3_g)*gas(ihcl_g) - Keq_sg(2)
5384 xgas = quadratic(a,b,c)
5386 if(xgas .ge. electrolyte(jnh4cl,jp,ibin))then ! degas all nh4cl
5388 gas(inh3_g) = gas(inh3_g) + electrolyte(jnh4cl,jp,ibin)
5389 gas(ihcl_g) = gas(ihcl_g) + electrolyte(jnh4cl,jp,ibin)
5390 aer(inh4_a,jp,ibin) = aer(inh4_a,jp,ibin) - &
5391 electrolyte(jnh4cl,jp,ibin)
5392 aer(icl_a,jp,ibin) = aer(icl_a,jp,ibin) - &
5393 electrolyte(jnh4cl,jp,ibin)
5395 else ! degas only xgas amount of nh4cl
5397 gas(inh3_g) = gas(inh3_g) + xgas
5398 gas(ihcl_g) = gas(ihcl_g) + xgas
5399 aer(inh4_a,jp,ibin) = aer(inh4_a,jp,ibin) - xgas
5400 aer(icl_a,jp,ibin) = aer(icl_a,jp,ibin) - xgas
5406 aer(inh4_a,jtotal,ibin) = aer(inh4_a,jsolid,ibin) + &
5407 aer(inh4_a,jliquid,ibin)
5408 aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid,ibin) + &
5409 aer(icl_a,jliquid,ibin)
5412 end subroutine degas_solid_nh4cl
5416 !***********************************************************************
5417 ! conforms aerosol generic species to a valid electrolyte composition
5419 ! author: Rahul A. Zaveri
5421 !-----------------------------------------------------------------------
5422 subroutine conform_electrolytes(jp,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)
5424 use module_data_mosaic_aero, only: r8,naer,nbin_a_max, &
5425 ngas_aerchtot, ngas_volatile, nelectrolyte, &
5426 imsa_a,iso4_a,ica_a,ina_a,inh4_a,ino3_a,icl_a,ico3_a
5431 integer, intent(in) :: ibin, jp
5432 real(r8), intent(inout) :: XT
5433 real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
5434 real(r8), intent(inout), dimension(ngas_volatile) :: total_species
5435 real(r8), intent(inout) :: tot_cl_in
5436 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
5437 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
5439 integer i, iXT_case, je
5440 real(r8) :: sum_dum, XNa_prime, XNH4_prime, XT_prime
5441 real(r8) :: store(naer)
5443 ! remove negative concentrations, if any
5445 ! aer(i,jp,ibin) = max(0.0d0, aer(i,jp,ibin)) ! EFFI
5449 ! call calculate_XT(ibin,jp,XT,aer) ! EFFI
5451 if( (aer(iso4_a,jp,ibin)+aer(imsa_a,jp,ibin)) .gt.0.0)then
5452 XT = ( aer(inh4_a,jp,ibin) + &
5453 aer(ina_a,jp,ibin) + &
5454 2.*aer(ica_a,jp,ibin) )/ &
5455 (aer(iso4_a,jp,ibin)+0.5*aer(imsa_a,jp,ibin))
5461 ! if(XT .ge. 1.9999 .or. XT.lt.0.)then ! RAZ 11/10/2014
5462 if(XT .ge. 2.0 .or. XT.lt.0.)then ! RAZ 11/10/2014
5463 iXT_case = 1 ! near neutral (acidity is caused by HCl and/or HNO3)
5465 iXT_case = 2 ! acidic (acidity is caused by excess SO4)
5470 ! put total aer(*) into store(*)
5471 store(iso4_a) = aer(iso4_a,jp,ibin)
5472 store(ino3_a) = aer(ino3_a,jp,ibin)
5473 store(icl_a) = aer(icl_a, jp,ibin)
5474 store(imsa_a) = aer(imsa_a,jp,ibin)
5475 store(ico3_a) = aer(ico3_a,jp,ibin)
5476 store(inh4_a) = aer(inh4_a,jp,ibin)
5477 store(ina_a) = aer(ina_a, jp,ibin)
5478 store(ica_a) = aer(ica_a, jp,ibin)
5480 do je=1,nelectrolyte
5481 electrolyte(je,jp,ibin) = 0.0
5485 !---------------------------------------------------------
5487 if(iXT_case.eq.1)then
5489 ! XT >= 2 : sulfate deficient
5491 call form_caso4(store,jp,ibin,electrolyte)
5492 call form_camsa2(store,jp,ibin,electrolyte)
5493 call form_na2so4(store,jp,ibin,electrolyte)
5494 call form_namsa(store,jp,ibin,electrolyte)
5495 call form_cano3(store,jp,ibin,electrolyte)
5496 call form_nano3(store,jp,ibin,electrolyte)
5497 call form_nacl(store,jp,ibin,aer,gas,electrolyte,total_species,tot_cl_in)
5498 call form_cacl2(store,jp,ibin,electrolyte)
5499 call form_caco3(store,jp,ibin,aer,electrolyte)
5500 call form_nh4so4(store,jp,ibin,electrolyte)
5501 call form_nh4msa(store,jp,ibin,electrolyte)
5502 call form_nh4no3(store,jp,ibin,electrolyte)
5503 call form_nh4cl(store,jp,ibin,electrolyte)
5504 call form_msa(store,jp,ibin,electrolyte)
5505 call degas_hno3(store,jp,ibin,aer,gas,electrolyte)
5506 call degas_hcl(store,jp,ibin,aer,gas,electrolyte)
5507 call degas_nh3(store,jp,ibin,aer,gas)
5509 elseif(iXT_case.eq.2)then
5511 ! XT < 2 : sulfate enough or sulfate excess
5513 call form_caso4(store,jp,ibin,electrolyte)
5514 call form_camsa2(store,jp,ibin,electrolyte)
5515 call form_namsa(store,jp,ibin,electrolyte)
5516 call form_nh4msa(store,jp,ibin,electrolyte)
5517 call form_msa(store,jp,ibin,electrolyte)
5519 if(store(iso4_a).eq.0.0)goto 10
5522 XT_prime =(store(ina_a)+store(inh4_a))/ &
5524 XNa_prime=0.5*store(ina_a)/store(iso4_a) + 1.
5526 if(XT_prime.ge.XNa_prime)then
5527 call form_na2so4(store,jp,ibin,electrolyte)
5529 if(store(iso4_a).gt.1.e-15)then
5530 XNH4_prime = store(inh4_a)/store(iso4_a)
5533 if(XNH4_prime .ge. 1.5)then
5534 call form_nh4so4_lvcite(store,jp,ibin,electrolyte)
5536 call form_lvcite_nh4hso4(store,jp,ibin,electrolyte)
5539 elseif(XT_prime.ge.1.)then
5540 call form_nh4hso4(store,jp,ibin,electrolyte)
5541 call form_na2so4_nahso4(store,jp,ibin,electrolyte)
5542 elseif(XT_prime.lt.1.)then
5543 call form_nahso4(store,jp,ibin,electrolyte)
5544 call form_nh4hso4(store,jp,ibin,electrolyte)
5545 call form_h2so4(store,jp,ibin,electrolyte)
5548 10 call degas_hno3(store,jp,ibin,aer,gas,electrolyte)
5549 call degas_hcl(store,jp,ibin,aer,gas,electrolyte)
5550 call degas_nh3(store,jp,ibin,aer,gas)
5555 ! re-calculate ions to eliminate round-off errors
5556 call electrolytes_to_ions(jp, ibin,aer,electrolyte)
5557 !---------------------------------------------------------
5559 ! calculate % composition EFFI
5561 !! do je = 1, nelectrolyte
5562 !! electrolyte(je,jp,ibin) = max(0.d0,electrolyte(je,jp,ibin)) ! remove -ve
5563 !! sum_dum = sum_dum + electrolyte(je,jp,ibin)
5566 !! if(sum_dum .eq. 0.)sum_dum = 1.0
5567 !! electrolyte_sum(jp,ibin) = sum_dum
5569 !! do je = 1, nelectrolyte
5570 !! epercent(je,jp,ibin) = 100.*electrolyte(je,jp,ibin)/sum_dum
5575 end subroutine conform_electrolytes
5579 !----------------------------------------------------------
5580 ! solution to x^3 + px^2 + qx + r = 0
5582 function cubic( psngl, qsngl, rsngl )
5583 use module_data_mosaic_kind, only: r8
5587 real(r8) :: psngl, qsngl, rsngl
5589 real(r8) :: p, q, r, A, B, D, M, N, third, y
5590 real(r8) :: k, phi, thesign, x(3), duma
5599 A = (1.d0/3.d0)*((3.d0*q) - (p*p))
5600 B = (1.d0/27.d0)*((2.d0*p*p*p) - (9.d0*p*q) + (27.d0*r))
5602 D = ( ((A*A*A)/27.d0) + ((B*B)/4.d0) )
5604 if(D .gt. 0.)then ! => 1 real and 2 complex roots
5606 elseif(D .eq. 0.)then ! => 3 real roots, atleast 2 identical
5608 else ! D < 0 => 3 distinct real roots
5622 M = thesign*((-B/2.d0) + (sqrt(D)))**(third)
5623 N = thesign*((-B/2.d0) - (sqrt(D)))**(third)
5625 cubic = ( (M) + (N) - (p/3.d0) )
5635 M = thesign*(-B/2.d0)**third
5638 x(1) = ( (M) + (N) - (p/3.d0) )
5639 x(2) = ( (-M/2.d0) + (-N/2.d0) - (p/3.d0) )
5640 x(2) = ( (-M/2.d0) + (-N/2.d0) - (p/3.d0) )
5644 if(x(kk).gt.cubic) cubic = x(kk)
5655 ! rce 18-nov-2004 -- make sure that acos argument is between +/-1.0
5656 ! phi = acos(thesign*sqrt( (B*B/4.d0)/(-A*A*A/27.d0) )) ! radians
5657 duma = thesign*sqrt( (B*B/4.d0)/(-A*A*A/27.d0) )
5658 duma = min( duma, +1.0d0 )
5659 duma = max( duma, -1.0d0 )
5660 phi = acos( duma ) ! radians
5666 y = 2.*Sqrt(-A/3.)*cos(phi + 120.*k*0.017453293)
5667 x(kk) = ((y) - (p/3.d0))
5668 if(x(kk).gt.cubic) cubic = x(kk)
5673 !----------------------------------------------------------
5676 !----------------------------------------------------------
5677 function quadratic(a,b,c)
5678 use module_data_mosaic_kind, only: r8
5680 real(r8) :: quadratic
5684 real(r8) :: x, dum, quad1, quad2
5693 if(abs(x) .lt. 1.e-6)then
5698 quadratic = (-0.5*b/a)*dum
5700 if(quadratic .lt. 0.)then
5701 quadratic = -b/a - quadratic
5705 quad1 = ((-b)+sqrt((b*b)-(4.*a*c)))/ &
5707 quad2 = ((-b)-sqrt((b*b)-(4.*a*c)))/ &
5710 quadratic = max(quad1, quad2)
5714 end function quadratic
5715 !----------------------------------------------------------
5718 !----------------------------------------------------------
5719 function mean_molecular_speed(T, MW) ! in cm/s
5720 use module_data_mosaic_kind, only: r8
5722 real(r8) :: mean_molecular_speed
5724 real(r8) :: T, MW ! T(K)
5726 mean_molecular_speed = 1.455e4 * sqrt(T/MW)
5729 end function mean_molecular_speed
5730 !----------------------------------------------------------
5732 !----------------------------------------------------------
5733 function gas_diffusivity(T, P, MW, Vm) ! in cm^2/s
5734 use module_data_mosaic_kind, only: r8
5735 use module_data_mosaic_constants, only: third
5737 real(r8) :: gas_diffusivity
5739 real(r8) :: MW, Vm, T, P ! T(K), P(atm)
5742 gas_diffusivity = (1.0e-3 * T**1.75 * sqrt(1./MW + 0.035))/ &
5743 (P * (Vm**third + 2.7189)**2)
5747 end function gas_diffusivity
5748 !----------------------------------------------------------
5751 !----------------------------------------------------------
5752 function fuchs_sutugin(rkn,a)
5753 use module_data_mosaic_kind, only: r8
5755 real(r8) :: fuchs_sutugin
5759 real(r8) :: rnum, denom
5762 rnum = 0.75*a*(1. + rkn)
5763 denom = rkn**2 + rkn + 0.283*rkn*a + 0.75*a
5764 fuchs_sutugin = rnum/denom
5767 end function fuchs_sutugin
5768 !----------------------------------------------------------
5771 !----------------------------------------------------------
5772 ! ZSR method at 60% RH
5774 function aerosol_water_up( ibin, electrolyte, aer, kappa_nonelectro, a_zsr ) ! kg (water)/m^3 (air)
5776 use module_data_mosaic_aero, only: r8, nelectrolyte, naer, nbin_a_max, jtotal, &
5777 nsalt, ioc_a, ibc_a, ilim2_a, ioin_a, dens_aer_mac, mw_aer_mac ! RAZ 4/16/2014
5779 use module_data_mosaic_asecthp, only: dens_water_aer
5783 real(r8) :: aerosol_water_up
5785 integer, intent(in) :: ibin
5786 real(r8), intent(in), dimension (6,nelectrolyte) :: a_zsr
5787 real(r8), intent(in), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
5788 real(r8), intent(in), dimension(naer,3,nbin_a_max) :: aer
5789 real(r8), intent(in), dimension(naer) :: kappa_nonelectro
5792 integer :: iaer, jp, je
5793 real(r8) :: tmpa, tmpb, aH2O_60 ! RAZ 4/16/2014
5795 !real(r8) :: bin_molality_60
5803 do je = 1, (nsalt+4) ! include hno3 and hcl in water calculation
5804 tmpa = tmpa + electrolyte(je,jp,ibin)/bin_molality_60(je,a_zsr)
5808 ! ( (aer(ilim2_a,jp,ibin)/dens_aer_mac(ilim2_a))*kappa_nonelectro(ilim2_a) + & ! RCE 5/20/2015
5809 ! (aer(ioin_a, jp,ibin)/dens_aer_mac(ioin_a ))*kappa_nonelectro(ioin_a ) + & ! " "
5810 ! (aer(ioc_a, jp,ibin)/dens_aer_mac(ioc_a ))*kappa_nonelectro(ioc_a ) + & ! " "
5811 ! (aer(ibc_a, jp,ibin)/dens_aer_mac(ibc_a ))*kappa_nonelectro(ibc_a ) )*aH2O_60/(1.0-aH2O_60) ! RCE 5/20/2015
5815 if (kappa_nonelectro(iaer) > 0.0_r8) then
5816 tmpb = tmpb + (aer(iaer,jp,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer))*kappa_nonelectro(iaer)
5819 tmpa = tmpa + tmpb * dens_water_aer * 1.0e-3 * aH2O_60/(1.0-aH2O_60) ! RCE 1/6/2016 added 1.0e-3 factor
5821 aerosol_water_up = tmpa*1.e-9 ! kg(water)/m^3(air)
5824 end function aerosol_water_up
5825 !----------------------------------------------------------
5828 !----------------------------------------------------------
5829 function bin_molality_60(je,a_zsr) ! TOUCH
5830 use module_data_mosaic_aero, only: r8,nelectrolyte
5834 real(r8) :: bin_molality_60
5836 integer, intent(in) :: je
5837 real(r8), intent(in), dimension (6,nelectrolyte) :: a_zsr
5844 xm = a_zsr(1,je) + &
5849 aw* a_zsr(6,je) ))))
5851 bin_molality_60 = 55.509_r8*xm/(1. - xm)
5854 end function bin_molality_60
5855 !----------------------------------------------------------
5858 !----------------------------------------------------------
5860 function aerosol_water( jp, ibin, jaerosolstate, jphase, jhyst_leg, electrolyte, aer, &
5861 kappa_nonelectro, num_a, mass_dry_a, mass_soluble_a, aH2O, molality0 ) ! kg (water)/m^3 (air). RAZ added aer
5863 use module_data_mosaic_aero, only: r8, nbin_a_max, nelectrolyte, nsoluble, naer, &
5864 all_solid, jsolid, jhyst_lo, ioc_a, ibc_a, ilim2_a, ioin_a, & ! RAZ 4/16/2014
5865 jtotal, ah2o_max, dens_aer_mac, mw_aer_mac, ename
5867 use module_data_mosaic_asecthp, only: dens_water_aer
5871 real(r8) :: aerosol_water
5873 integer, intent(in) :: jp, ibin
5874 integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg
5876 real(r8), intent(in) :: aH2O
5877 real(r8), intent(in), dimension(nbin_a_max) :: num_a,mass_dry_a,mass_soluble_a
5878 real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
5879 real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
5880 real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
5881 real(r8), intent(in), dimension(naer) :: kappa_nonelectro
5884 integer :: iaer, iclm_aer, jclm_aer, je
5885 real(r8) :: tmpa, tmpb
5887 real(r8) :: bin_molality
5892 if (jaerosolstate(ibin) .ne. all_solid) then ! RAZ 5/24/2016 - added this "if" check
5893 ! should this be "do je = 1, nsalt+4" as in aerosol_water_up ?
5894 ! or "do je = 1, jxxxx" where jxxxx is the index of the last electrolyte used here ?
5895 do je = 1, 19 ! include hno3 and hcl in water calculation
5896 tmpa = tmpa + electrolyte(je,jp,ibin)/molality0(je,ibin) ! RAZ 5/20/2014
5900 ! note that this only considers the ilim2_a soa species
5903 ! ( (aer(ioc_a, jtotal,ibin)/dens_aer_mac(ioc_a ))*kappa_nonelectro(ioc_a ) + & ! RCE 5/20/2015
5904 ! (aer(ilim2_a,jtotal,ibin)/dens_aer_mac(ilim2_a))*kappa_nonelectro(ilim2_a) + & ! RCE 5/20/2015
5905 ! (aer(ioin_a, jtotal,ibin)/dens_aer_mac(ioin_a ))*kappa_nonelectro(ioin_a ) + & ! RCE 5/20/2015
5906 ! (aer(ibc_a, jtotal,ibin)/dens_aer_mac(ibc_a ))*kappa_nonelectro(ibc_a ) ) & ! RCE 5/20/2015
5907 ! * 1.0e-3 * aH2O/(1.0-min(ah2o,ah2o_max)) ! RCE 5/20/2015 - need 1.0e-3 factor
5911 if (kappa_nonelectro(iaer) > 0.0_r8) then
5912 tmpb = tmpb + (aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer))*kappa_nonelectro(iaer)
5915 tmpa = tmpa + tmpb * dens_water_aer * 1.0e-3 * ah2o/(1.0-min(ah2o,ah2o_max)) ! ug(water)/m^3(air)
5917 aerosol_water = tmpa*1.e-9 ! kg(water)/m^3(air)
5920 if (aerosol_water .le. 0.0) then !BALLI- Commented out to avoid slow runtime.
5921 iclm_aer = 0 !BSINGH- THIS IS WRONG!!!
5922 jclm_aer = 0 !BSINGH- THIS IS WRONG!!!
5923 !write(6,*)'iclm jclm ibin jp = ', &
5924 ! iclm_aer, jclm_aer, ibin, jp !BSINGH- iclm_aer and jclm_aer are never set but they are used here.***
5926 !write(6,*)'aH2O, water = ', aH2O, aerosol_water
5927 !write(6,*)'dry mass = ', mass_dry_a(ibin)
5928 !write(6,*)'soluble mass = ', mass_soluble_a(ibin)
5929 !write(6,*)'number = ', num_a(ibin)
5930 !do je = 1, nsoluble
5931 ! write(6,44)ename(je), electrolyte(je,jp,ibin)
5933 !write(6,*)'Error in water calculation'
5934 !write(6,*)'ibin = ', ibin
5935 !write(6,*)'water content cannot be negative or zero'
5936 !write(6,*)'setting jaerosolstate to all_solid'
5940 jaerosolstate(ibin) = all_solid
5941 jphase(ibin) = jsolid
5942 jhyst_leg(ibin) = jhyst_lo
5946 44 format(a7, 2x, e11.3)
5950 end function aerosol_water
5956 !----------------------------------------------------------
5957 function bin_molality(je,ibin,aH2O_a,b_zsr,a_zsr,aw_min)
5958 use module_data_mosaic_aero, only:r8, nbin_a_max, nelectrolyte
5962 real(r8) :: bin_molality
5964 integer, intent(in) :: je, ibin
5965 real(r8), intent(in), dimension(nbin_a_max) :: aH2O_a
5966 real(r8), intent(in), dimension(nelectrolyte) :: b_zsr,aw_min
5967 real(r8), intent(in), dimension (6,nelectrolyte) :: a_zsr
5972 aw = max(aH2O_a(ibin), aw_min(je))
5973 aw = min(aw, 0.999999_r8)
5976 if(aw .lt. 0.97_r8)then
5978 xm = a_zsr(1,je) + &
5983 aw* a_zsr(6,je) ))))
5985 bin_molality = 55.509_r8*xm/(1. - xm)
5989 bin_molality = -b_zsr(je)*log(aw)
5995 end function bin_molality
5996 !----------------------------------------------------------
6004 end module module_mosaic_ext