Line data Source code
1 : module mo_sad
2 :
3 : use shr_kind_mod, only : r8 => shr_kind_r8
4 : use physconst, only : pi
5 : use ppgrid, only : pcols, pver
6 : use m_sad_data, only : a, b
7 : use cam_logfile, only : iulog
8 : use spmd_utils, only : masterproc
9 :
10 :
11 : implicit none
12 :
13 : private
14 : public :: sad_inti
15 : public :: sad_strat_calc
16 : public :: sad_top
17 :
18 : save
19 :
20 : real(r8), parameter :: four_thrd = 4._r8/3._r8
21 : real(r8), parameter :: one_thrd = 1._r8/3._r8
22 : real(r8), parameter :: two_thrd = 2._r8/3._r8
23 : real(r8), parameter :: four_pi = 4._r8*pi
24 :
25 : integer :: sad_top
26 : integer :: sad_topp
27 :
28 : contains
29 :
30 :
31 1536 : subroutine sad_inti(pbuf2d)
32 : !----------------------------------------------------------------------
33 : ! ... initialize the sad module
34 : !----------------------------------------------------------------------
35 :
36 : use ref_pres, only : pref_mid_norm
37 : use cam_history, only : addfld
38 : use physics_buffer, only : physics_buffer_desc, pbuf_set_field
39 :
40 : implicit none
41 :
42 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
43 :
44 : !----------------------------------------------------------------------
45 : ! ... Local variables
46 : !----------------------------------------------------------------------
47 : integer :: k
48 :
49 : !----------------------------------------------------------------------
50 : ! ... find level where etamids are all > 1 hPa
51 : !----------------------------------------------------------------------
52 1536 : sad_top = 0
53 130560 : do k = pver,1,-1
54 130560 : if( (pref_mid_norm(k)) < .001_r8 ) then
55 1536 : sad_top = k
56 1536 : exit
57 : end if
58 : end do
59 1536 : sad_topp = sad_top + 1
60 1536 : if (masterproc) then
61 2 : write(iulog,*) 'sad_inti: sad capped at level ',sad_top
62 2 : write(iulog,*) ' whose midpoint is ',pref_mid_norm(sad_topp)*1.e3_r8,' hPa'
63 : endif
64 :
65 3072 : call addfld( 'H2SO4M_C', (/ 'lev' /), 'I', 'ug/m3', 'chemical sulfate aerosol mass' )
66 :
67 1536 : end subroutine sad_inti
68 : !===============================================================================
69 : ! ROUTINE
70 : ! sad_strat_calc
71 : !
72 : ! Date...
73 : ! 14 October 1999
74 : !
75 : ! Programmed by...
76 : ! Douglas E. Kinnison
77 : ! Modified by
78 : ! Stacy Walters
79 : ! 1 November 1999
80 : ! Modified by
81 : ! Doug Kinnison
82 : ! 1 September 2004; Condensed phase H2O passed in from CAM
83 : ! 2 November 2004; New treatment of denitrificatoin (NAT)
84 : ! 14 November 2004; STS mode of operation.
85 : ! 27 March 2008; Using original NAT approach.
86 : ! 08 November 2010; STS Approach (HNO3 => STS; then HNO3 => NAT)
87 : ! 24 March 2011; updated mask logic and removed sm NAT logic
88 : ! 14 April 2011; updated EQUIL logic
89 : ! 19 December 2012; updated using Wegner et al., JGR, 2013a,b.
90 : ! 25 April 2013; Removed volcanic heating logic.
91 : ! 6 April 2022; update nat_part_dens from 1e-02 to 5.0e-4
92 : ! (Wilka et al., ACP, 2021, doi:10.5194/acp-21-15771-2021)
93 : !
94 : ! DESCRIPTION
95 : !
96 : ! This routine has the logic to derive the surface area density for
97 : ! three types of aerosols: Sulfate (LBS, STS); Nitric Acid Trihydrate (NAT);
98 : ! and Water-ICE. The surface area density is stored in sad_strat(3). The
99 : ! first, second, and third dimensions are SULFATE, NAT, and ICE SAD
100 : ! respectively. The effective radius of each particle is also stored
101 : ! in radius_strat(3).
102 : !
103 : ! NOTE1: For Sulfate and H2O ICE calculations
104 : ! The Surface Area and Volume Densities are calculated from the
105 : ! second and third moment of the LOG Normal Distribution. For an example
106 : ! see Binkowski et al., JGR, 100, 26191-26209, 1995. The Volume Density
107 : ! is substituted into the SAD equation so that the SAD is dependent on
108 : ! the # of particles cm-3, the width of the distribution (sigma), and
109 : ! the volume density of the aerosol. This approach is discussed in
110 : ! Considine et al., 2000 and Kinnison et al., 2007.
111 : !
112 : ! NOTE2: For the ternary solution calculation
113 : ! the total sulfate mass is derived from the SAGEII SAD data. This approach
114 : ! has been previously used in Considine et al., JGR, 1999. The thermodynamic
115 : ! models used in this routine are from A. Tabazedeh et al, 1994.
116 : !
117 : ! NOTE3: Updates to the PSC scheme is discussed in Wegner et al., 2013a.
118 : ! 80% of the total HNO3 is allowed to see STS, 20% NAT. The number density of
119 : ! the NAT and ICE particles are is set to 0.01 and 0.1 particle cm-3 respectively.
120 : !
121 : ! NOTE4: The HCl solubility (in STS) has been added and evalutede in Wegner et al., 2013b.
122 : ! This solubility is based on Carslaw et al., 1995.
123 : !
124 : ! REFERENCES for this PSC module:
125 : ! Considine, D. B., A. R. Douglass, P. S. Connell, D. E. Kinnison, and D. A., Rotman,
126 : ! A polar stratospheric cloud parameterization for the three dimensional model of
127 : ! the global modeling initiative and its response to stratospheric aircraft,
128 : ! J. Geophys. Res., 105, 3955-3975, 2000.
129 : !
130 : ! Kinnison, D. E.,et al., Sensitivity of chemical tracers to meteorological
131 : ! parameters in the MOZART-3 chemical transport model, J. Geophys. Res.,
132 : ! 112, D20302, doi:10.1029/2006JD007879, 2007.
133 : !
134 : ! Wegner, T, D. E. Kinnison, R. R. Garcia, S. Madronich, S. Solomon, and M. von Hobe,
135 : ! On the depletion of HCl in the Antarctic polar vortex,
136 : ! in review J. Geophys. Res., 2013.
137 : !
138 : ! Wegner, T, D. E. Kinnison, R. R. Garcia, S. Madronich, and S. Solomon,
139 : ! Polar Stratospheric Clouds in SD-WACCM4, in review J. Geophys. Res., 2013.
140 : !
141 : ! Tabazedeh, A., R. P. Turco, K. Drdla, M. Z. Jacobson, and O. B, Toon, A study
142 : ! of the type I polar stratosphere cloud formation,
143 : ! Geophys. Res. Lett., 21, 1619-1622, 1994.
144 : !
145 : ! Carslaw, K. S., S. L. Clegg, and P. Brimblecombe, A thermodynamic model of the
146 : ! system HCl-HNO3-H2SO4-H2O, including solubilities of HBr, from <200 to 328K,
147 : ! J. Phys. Chem., 99, 11,557-11,574, doi:1021/100029a039, 1995.
148 : !
149 : !
150 : ! ARGUMENTS
151 : ! INPUT:
152 : ! hno3_gas Nitric Acid gas phase abundance (mole fraction)
153 : ! hno3_cond(2) Nitric Acid cond. phase abundance (mole fraction)
154 : ! (1) in STS; (2) in NAT
155 : ! h2o_cond Water condensed phase (mole fraction)
156 : ! h2o_gas Water gas-phase abundance (mole fraction)
157 : !
158 : ! hcl_gas HCL gas-phase abundance (mole fraction)
159 : ! hcl_cond HCl condensed phase (STS) (mole fraction)
160 : !
161 : ! sage_sad SAGEII surface area density (cm2-aer cm-3 atm)
162 : ! m Airdensity (molecules cm-3)
163 : ! press Pressure (hPa)
164 : !
165 : !
166 : ! OUTPUT:
167 : !
168 : ! hno3_gas = Gas-phase HNO3 Used in chemical solver.
169 : ! hno3_cond(1) = Condensed HNO3 from STS Not used in mo_aero_settling.F90
170 : ! hno3_cond(2) = Condensed HNO3 from NAT Used in mo_aero_settling.F90
171 : !
172 : ! hcl_gas = Gas-phase HCL Used in chemical solver.
173 : ! hcl_cond = Condensed HCl from STS
174 : !
175 : ! SAD_strat(1) = Sulfate Aerosol... Used in mo_strato_rates.F90
176 : ! SAD_strat(2) = NAT Aerosol.... Used in mo_strato_rates.F90
177 : ! SAD_strat(3) = Water-Ice......... Used in mo_strato_rates.F90
178 : !
179 : ! RAD_strat(1) = Sulfate Aerosol... Used in mo_strato_rates.F90
180 : ! RAD_strat(2) = NAT large mode.... Used in mo_aero_settling.F90
181 : ! RAD_strat(3) = Water-Ice......... Not used in mo_aero_settling.F90
182 : !
183 : ! NOTE1: The sum of hno3_cond(1-2) will be added to hno3_gas for advection of HNO3 in
184 : ! mo_gas_phase_chemdr.F90.
185 : !
186 : ! NOTE2: The sum of hcl_cond will be added to hcl_gas for advection of HCl in
187 : ! mo_gas_phase_chemdr.F90.
188 : !
189 : ! NOTE3: This routine does not partition H2O.
190 : !
191 : !
192 : ! ROUTINES Called (in and below this routine):
193 : !
194 : ! sad_strat_calc
195 : ! nat_sat_temp ...... derives the NAT saturation temp
196 : !
197 : ! sulfate_sad_calc .. Calculates the sulfate SAD; HNO3, H2O cond. phase
198 : ! calc_radius_lbs ... Calculates the radius for a H2SO4 Binary soln. (T>200K)
199 : ! sad_to_h2so4 ...... Derives H2SO4 abundance (micrograms m-3)
200 : ! from SAGEII SAD.
201 : ! density............ A. Tabazedeh binary thermo model
202 : ! equil.............. A. Tabazedeh ternary thermo. model
203 : !
204 : ! nat_sad_calc....... Calculates the NAT SAD; HNO3, H2O cond. phase
205 : ! nat_cond........... Derives the NAT HNO3 gas/cond partitioning
206 : !
207 : ! ice_sad_calc....... derives the ICE SAD and H2O gas/cond partitioning
208 : !===============================================================================
209 :
210 0 : subroutine sad_strat_calc( lchnk, m, press, temper, hno3_gas, &
211 0 : hno3_cond, h2o_gas, h2o_cond, hcl_gas, hcl_cond, &
212 0 : sad_sage, radius_strat, sad_strat, ncol, pbuf )
213 :
214 1536 : use cam_history, only : outfld
215 : use physics_buffer, only : physics_buffer_desc
216 :
217 : implicit none
218 :
219 : !-------------------------------------------------------------------------------
220 : ! ... dummy arguments
221 : !-------------------------------------------------------------------------------
222 : integer, intent(in) :: lchnk ! chnk id
223 : integer, intent(in) :: ncol ! columns in chunk
224 : real(r8), intent(in) :: m (ncol,pver) ! Air density (molecules cm-3)
225 : real(r8), intent(in) :: sad_sage (ncol,pver) ! SAGEII surface area density (cm2 aer. cm-3 air)
226 : real(r8), intent(in) :: press (ncol,pver) ! Pressure, hPa
227 : real(r8), intent(in) :: temper (pcols,pver) ! Temperature (K)
228 : real(r8), intent(inout) :: h2o_gas (ncol,pver) ! H2O gas-phase (mole fraction)
229 : real(r8), intent(inout) :: h2o_cond (ncol,pver) ! H2O condensed phase (mole fraction)
230 : real(r8), intent(inout) :: hno3_gas (ncol,pver) ! HNO3 condensed phase (mole fraction)
231 : real(r8), intent(inout) :: hno3_cond (ncol,pver,2) ! HNO3 condensed phase (mole fraction)
232 : real(r8), intent(inout) :: hcl_gas (ncol,pver) ! HCL gas-phase (mole fraction)
233 : real(r8), intent(inout) :: hcl_cond (ncol,pver) ! HCL condensed phase (mole fraction)
234 : real(r8), intent(out) :: radius_strat(ncol,pver,3) ! Radius of Sulfate, NAT, and ICE (cm)
235 : real(r8), intent(out) :: sad_strat (ncol,pver,3) ! Surface area density of Sulfate, NAT, ICE (cm2 cm-3)
236 :
237 : !-------------------------------------------------------------------------------
238 : ! ... local variables
239 : !-------------------------------------------------------------------------------
240 : real(r8), parameter :: temp_floor = 0._r8
241 :
242 : integer :: i, k, n
243 : integer :: dims(1)
244 0 : real(r8) :: hno3_total (ncol,pver) ! HNO3 total = gas-phase + condensed
245 0 : real(r8) :: h2o_total (ncol,pver) ! H2O total = gas-phase + condensed
246 0 : real(r8) :: radius_lbs (ncol,pver) ! Radius of Liquid Binary Sulfate (cm)
247 0 : real(r8) :: radius_sulfate(ncol,pver) ! Radius of Sulfate aerosol (cm)
248 0 : real(r8) :: radius_nat (ncol,pver) ! Radius of NAT aerosol (cm)
249 0 : real(r8) :: radius_ice (ncol,pver) ! Radius of ICE aerosol (cm)
250 0 : real(r8) :: sad_nat (ncol,pver) ! SAD of NAT aerosol (cm2 cm-3)
251 0 : real(r8) :: sad_sulfate (ncol,pver) ! SAD of Sulfate aerosol (cm2 cm-3)
252 0 : real(r8) :: sad_ice (ncol,pver) ! SAD of ICE aerosol (cm2 cm-3)
253 0 : real(r8) :: tsat_nat (ncol,pver) ! Temperature for NAT saturation
254 0 : real(r8) :: h2o_avail (ncol,pver) ! H2O temporary arrays
255 0 : real(r8) :: hno3_avail (ncol,pver) ! HNO3 temporary array
256 0 : real(r8) :: hno3_gas_nat (ncol,pver) ! HNO3 after call to NAT routines
257 0 : real(r8) :: hno3_gas_sulf (ncol,pver) ! HNO3 after call to STS routines
258 0 : real(r8) :: hno3_cond_nat (ncol,pver) ! HNO3 condensed after call to NAT
259 0 : real(r8) :: hno3_cond_sulf(ncol,pver) ! HNO3 condensed after call to STS routines
260 0 : real(r8) :: hcl_total (ncol,pver) ! HCl total = gas-phase + condensed
261 0 : real(r8) :: hcl_avail (ncol,pver) ! HCL temporary arrays
262 0 : real(r8) :: hcl_gas_sulf (ncol,pver) ! HCL after call to STS routines
263 0 : real(r8) :: hcl_cond_sulf (ncol,pver) ! HCL condensed after call to STS routines
264 : real(r8) :: temp (pcols,pver) ! wrk temperature array
265 0 : real(r8) :: h2so4m (ncol,pver) ! wrk array
266 :
267 0 : logical :: z_val(ncol)
268 :
269 0 : logical :: mask_lbs(ncol,pver) ! LBS mask T: P>300hPa; P<2hPa2hPa; SAGE<1e-15; T>200K
270 0 : logical :: mask_ice(ncol,pver) ! ICE mask T: .not. mask_lbs; h2o_cond>0
271 0 : logical :: mask_sts(ncol,pver) ! STS mask T: .not. mask_sts
272 0 : logical :: mask_nat(ncol,pver) ! NAT mask T: mask_sts=T; T<Tsat_nat
273 : type(physics_buffer_desc), pointer :: pbuf(:)
274 :
275 : real(r8), parameter :: eighty_percent = 0.8_r8
276 : real(r8), parameter :: twenty_percent = 0.2_r8
277 :
278 : !----------------------------------------------------------------------
279 : ! ... initialize to zero
280 : !----------------------------------------------------------------------
281 0 : do n = 1,3
282 0 : do k = 1,pver
283 0 : radius_strat(:,k,n) = 0._r8
284 0 : sad_strat(:,k,n) = 0._r8
285 : end do
286 : end do
287 : !
288 0 : do n = 1,2
289 0 : do k = 1,pver
290 0 : hno3_cond(:,k,n) = 0._r8
291 : end do
292 : end do
293 : !
294 0 : do k = 1,pver
295 0 : h2o_total (:,k) = 0._r8
296 : !
297 0 : h2o_avail (:,k) = 0._r8
298 0 : hno3_avail (:,k) = 0._r8
299 0 : hcl_avail (:,k) = 0._r8
300 : !
301 0 : hno3_total (:,k) = 0._r8
302 0 : hno3_gas_nat (:,k) = 0._r8
303 0 : hno3_gas_sulf (:,k) = 0._r8
304 0 : hno3_cond_nat (:,k) = 0._r8
305 0 : hno3_cond_sulf(:,k) = 0._r8
306 : !
307 0 : hcl_total (:,k) = 0._r8
308 0 : hcl_cond_sulf (:,k) = 0._r8
309 0 : hcl_gas_sulf (:,k) = 0._r8
310 : end do
311 : !----------------------------------------------------------------------
312 : ! ... limit temperature
313 : !----------------------------------------------------------------------
314 0 : do k = 1,pver
315 0 : h2so4m(:ncol,k) = 0._r8
316 0 : temp(:ncol,k) = max( temp_floor,temper(:ncol,k) )
317 : end do
318 :
319 : !----------------------------------------------------------------------
320 : ! ... total HNO3 and H2O gas and condensed phases
321 : !----------------------------------------------------------------------
322 0 : do k = sad_topp,pver
323 0 : hno3_total(:,k) = hno3_gas(:,k) + hno3_cond(:,k,1) + hno3_cond(:,k,2)
324 0 : h2o_total(:,k) = h2o_gas(:,k) + h2o_cond(:,k)
325 0 : hcl_total (:,k) = hcl_gas(:,k) + hcl_cond(:,k)
326 : end do
327 :
328 : !======================================================================
329 : !======================================================================
330 : ! ... Set SAD to SAGEII if Temperature is GT 200K and Return
331 : !
332 : ! ... mask_lbs = true .... T > 200K or SAD_SULF <= 1e-15 or
333 : ! P < 2hPa or P > 300hPa
334 : ! ... mask_sts = false .... .not. mask_lbs
335 : ! ... mask_nat = false .... T <= TSAT_NAT
336 : ! ... mask_ice = false .... H2O_COND > 0.0
337 : !======================================================================
338 : !======================================================================
339 :
340 0 : do k = sad_topp,pver
341 0 : mask_lbs(:,k) = temp(:ncol,k) > 200._r8 .or. sad_sage(:,k) <= 1.e-15_r8 &
342 0 : .or. press(:ncol,k) < 2._r8 .or. press(:ncol,k) > 300._r8
343 : end do
344 :
345 : sage_sad : &
346 0 : if( any( mask_lbs(:,sad_topp:pver) ) ) then
347 0 : do k = sad_topp,pver
348 0 : where( mask_lbs(:,k) )
349 0 : sad_strat(:,k,1) = sad_sage(:,k)
350 : sad_strat(:,k,2) = 0._r8
351 : sad_strat(:,k,3) = 0._r8
352 : endwhere
353 : end do
354 : !----------------------------------------------------------------------
355 : ! ... Calculate Liquid Binary Sulfate Radius
356 : !----------------------------------------------------------------------
357 0 : call calc_radius_lbs( ncol, mask_lbs, sad_sage, radius_lbs )
358 0 : do k = sad_topp,pver
359 0 : where( mask_lbs(:,k) )
360 0 : radius_strat(:,k,1) = radius_lbs(:,k)
361 : radius_strat(:,k,2) = 0._r8
362 : radius_strat(:,k,3) = 0._r8
363 : hno3_gas (:,k) = hno3_total(:,k)
364 : hno3_cond (:,k,1) = 0._r8
365 : hno3_cond (:,k,2) = 0._r8
366 : hcl_gas (:,k) = hcl_total(:,k)
367 : hcl_cond (:,k) = 0._r8
368 : endwhere
369 : end do
370 0 : if( all( mask_lbs(:,sad_topp:pver) ) ) then
371 0 : call outfld( 'H2SO4M_C', h2so4m(:ncol,:), ncol, lchnk )
372 0 : return
373 : end if
374 : end if sage_sad
375 :
376 : !======================================================================
377 : !======================================================================
378 : ! ... Logic for deriving ICE
379 : ! Ice formation occurs here if condensed phase H2O exists.
380 : !
381 : ! mask_lbs = false.... T > 200K or SAD_SULF < 1e-15 or
382 : ! P >2hPa or P <300hPa
383 : ! mask_ice = true .... H2O_COND > 0.0
384 : !======================================================================
385 : !======================================================================
386 0 : do k = sad_topp,pver
387 0 : do i = 1,ncol
388 0 : if( .not. mask_lbs(i,k) ) then
389 0 : mask_ice(i,k) = h2o_cond(i,k) > 0._r8
390 : else
391 0 : mask_ice(i,k) = .false.
392 : end if
393 : end do
394 : end do
395 :
396 : all_ice : &
397 0 : if( any( mask_ice(:,sad_topp:pver) ) ) then
398 0 : do k = sad_topp,pver
399 0 : where( mask_ice(:,k) )
400 0 : h2o_avail(:,k) = h2o_cond(:,k)
401 : endwhere
402 : end do
403 : !----------------------------------------------------------------------
404 : ! ... ICE
405 : !----------------------------------------------------------------------
406 : call ice_sad_calc( ncol, press, temp, m, h2o_avail, &
407 0 : sad_ice, radius_ice, mask_ice )
408 :
409 0 : do k = sad_topp,pver
410 0 : where( mask_ice(:,k) )
411 0 : sad_strat (:,k,3) = sad_ice (:,k)
412 : radius_strat(:,k,3) = radius_ice (:,k)
413 : endwhere
414 : end do
415 : end if all_ice
416 :
417 : !======================================================================
418 : !======================================================================
419 : ! ... LOGIC for STS and NAT
420 : !
421 : ! mask_lbs = false .... T > 200K or SAD_SULF <= 1e-15 or
422 : ! P < 2hPa or P >300hPa
423 : ! mask_sts = true .... not mask_lbs
424 : ! mask_nat = true .... T <= TSAT_NAT and mask_sts = true
425 : !======================================================================
426 : !======================================================================
427 :
428 0 : do k = sad_topp,pver
429 0 : do i = 1,ncol
430 0 : if( .not. mask_lbs(i,k) ) then
431 0 : mask_sts(i,k) = .true.
432 : else
433 0 : mask_sts(i,k) = .false.
434 : end if
435 : end do
436 : end do
437 : !----------------------------------------------------------------------
438 : ! ... STS (80% of total HNO3 logic)
439 : !----------------------------------------------------------------------
440 : ! NOTE: STS only sees 80% of the total HNO3 (Wegner et al., JGR, 2013a)
441 : sts_nat_sad : &
442 0 : if( any( mask_sts(:,sad_topp:pver) ) ) then
443 0 : do k = sad_topp,pver
444 0 : where( mask_sts(:,k) )
445 0 : h2o_avail (:,k) = h2o_gas (:,k)
446 : hno3_avail(:,k) = hno3_total(:,k)*eighty_percent
447 : hcl_avail (:,k) = hcl_total (:,k)
448 : endwhere
449 0 : if( any(mask_sts(:,k)) ) then
450 0 : where( mask_sts(:,k) )
451 : z_val(:) = hno3_avail(:,k) == 0._r8
452 : elsewhere
453 : z_val(:) = .false.
454 : endwhere
455 0 : if( any( z_val(:) ) ) then
456 0 : write(iulog,*) 'sad_strat_calc: Before CHEM Sulfate_SAD_CALC_1 has zero hno3_avail at lchnk,k = ',lchnk,k
457 : end if
458 : end if
459 : end do
460 :
461 : call sulfate_sad_calc( ncol, press, temp, h2o_avail, hno3_avail, hcl_avail, &
462 : sad_sage, m, hno3_gas_sulf, hno3_cond_sulf, &
463 : hcl_gas_sulf, hcl_cond_sulf, sad_sulfate, &
464 0 : radius_sulfate, mask_sts, lchnk, 1, h2so4m, .true.)
465 0 : do k = sad_topp,pver
466 0 : where( mask_sts(:,k) )
467 0 : sad_strat (:,k,1) = sad_sulfate (:,k)
468 : radius_strat(:,k,1) = radius_sulfate(:,k)
469 : hno3_gas (:,k) = hno3_gas_sulf (:,k)
470 : hno3_cond (:,k,1) = hno3_cond_sulf(:,k)
471 : hcl_gas (:,k) = hcl_gas_sulf (:,k)
472 : hcl_cond (:,k) = hcl_cond_sulf (:,k)
473 : endwhere
474 : end do
475 :
476 : !----------------------------------------------------------------------
477 : ! ... NAT (20% of total HNO3 logic)
478 : ! ... using total H2O and gas-phase HNO3 after STS calc
479 : !----------------------------------------------------------------------
480 : ! NOTE: NAT only sees 20% of the total HNO3 (Wegner et al., JGR, 2013a)
481 0 : hno3_avail(:,:) = hno3_total(:,:)*twenty_percent
482 0 : call nat_sat_temp( ncol, hno3_avail, h2o_avail, press, tsat_nat, mask_sts)
483 :
484 0 : do k = sad_topp,pver
485 0 : do i = 1,ncol
486 0 : if( .not. mask_lbs(i,k) ) then
487 0 : mask_nat(i,k) = temp(i,k) <= tsat_nat(i,k)
488 : else
489 0 : mask_nat(i,k) = .false.
490 : end if
491 : end do
492 : end do
493 :
494 0 : do k = sad_topp,pver
495 0 : where( mask_nat(:,k) )
496 0 : h2o_avail (:,k) = h2o_gas (:,k)
497 : hno3_avail(:,k) = hno3_total (:,k)*twenty_percent
498 : endwhere
499 0 : if( any(mask_nat(:,k)) ) then
500 0 : where( mask_nat(:,k) )
501 : z_val(:) = hno3_avail(:,k) == 0._r8
502 : elsewhere
503 : z_val(:) = .false.
504 : endwhere
505 0 : if( any( z_val(:) ) ) then
506 0 : write(iulog,*) 'sad_nat_calc: After CHEM Sulf_SAD_CALC_1 has zero hno3_avail at lchnk,k = ',lchnk,k
507 : end if
508 : end if
509 : end do
510 :
511 : call nat_sad_calc( ncol, press, temp, h2o_avail, hno3_avail, m, &
512 : hno3_gas_nat, hno3_cond_nat, &
513 0 : sad_nat, radius_nat, mask_nat )
514 :
515 :
516 : ! NOTE: Add in gas-phase from STS with gas-phase from NAT
517 0 : do k = sad_topp,pver
518 0 : where( mask_nat(:,k) )
519 0 : sad_strat (:,k,2) = sad_nat (:,k)
520 : radius_strat(:,k,2) = radius_nat (:,k)
521 : hno3_gas (:,k) = hno3_gas_sulf (:,k) + hno3_gas_nat (:,k)
522 : hno3_cond (:,k,2) = hno3_cond_nat (:,k)
523 : endwhere
524 : end do
525 :
526 : ! NOTE: If NAT does not form (in STS region), need to add the 20% of the total HNO3 back to gas-phase
527 0 : do k = sad_topp,pver
528 0 : do i = 1,ncol
529 0 : if ( .not. mask_nat(i,k) ) then
530 0 : if ( mask_sts(i,k) ) then
531 0 : hno3_gas (i,k) = hno3_gas_sulf(i,k) + hno3_total(i,k)*twenty_percent
532 : end if
533 : end if
534 : end do
535 : end do
536 :
537 : end if sts_nat_sad
538 :
539 :
540 0 : call outfld( 'H2SO4M_C', h2so4m(:ncol,:), ncol, lchnk )
541 :
542 0 : end subroutine sad_strat_calc
543 :
544 0 : subroutine nat_sat_temp( ncol, hno3_total, h2o_avail, press, tsat_nat, mask )
545 :
546 : implicit none
547 :
548 : !----------------------------------------------------------------------
549 : ! ... dummy arguments
550 : !----------------------------------------------------------------------
551 : integer, intent(in) :: ncol
552 : real(r8), intent(in) :: press(ncol,pver)
553 : real(r8), intent(in) :: h2o_avail(ncol,pver)
554 : real(r8), intent(in) :: hno3_total(ncol,pver)
555 : real(r8), intent(out) :: tsat_nat(ncol,pver)
556 : logical, intent(in) :: mask(ncol,pver)
557 :
558 : !----------------------------------------------------------------------
559 : ! ... local variables
560 : !----------------------------------------------------------------------
561 : real(r8), parameter :: ssrNAT = 10.0_r8
562 : real(r8), parameter :: ssrNATi = .1_r8
563 : real(r8), parameter :: aa = -2.7836_r8, &
564 : bb = -0.00088_r8, &
565 : cc = 38.9855_r8, &
566 : dd = -11397.0_r8, &
567 : ee = 0.009179_r8
568 : integer :: k, i
569 : real(r8) :: bbb ! temporary variable
570 : real(r8) :: ccc ! temporary variable
571 : real(r8) :: wrk ! temporary variable
572 : real(r8) :: tmp ! temporary variable
573 : real(r8) :: phno3 ! hno3 partial pressure
574 : real(r8) :: ph2o ! h2o partial pressure
575 :
576 0 : tsat_nat(:,:) = 0._r8
577 :
578 : !----------------------------------------------------------------------
579 : ! ... Derive HNO3 and H2O partial pressure (torr)
580 : ! where: 0.7501 = 760/1013.
581 : !----------------------------------------------------------------------
582 0 : do k = sad_topp,pver
583 0 : do i = 1,ncol
584 0 : if( mask(i,k) ) then
585 :
586 0 : bbb = press(i,k) * .7501_r8
587 0 : phno3 = hno3_total(i,k) * bbb
588 0 : ph2o = h2o_avail(i,k) * bbb
589 :
590 0 : if( phno3 > 0._r8 ) then
591 : !----------------------------------------------------------------------
592 : ! ... Calculate NAT Saturation Threshold Temperature
593 : ! Hanson and Mauersberger: GRL, Vol.15, 8, p855-858, 1988.
594 : ! Substitute m(T) and b(T) into Equation (1). Rearrange and
595 : ! solve quadratic eqation.
596 : !----------------------------------------------------------------------
597 0 : tmp = log10( ph2o )
598 0 : wrk = 1._r8 / (ee + bb*tmp)
599 0 : bbb = (aa*tmp - log10( phno3*ssrNATi ) + cc) * wrk
600 0 : ccc = dd *wrk
601 0 : tsat_nat(i,k) = .5_r8 * (-bbb + sqrt( bbb*bbb - 4._r8*ccc ))
602 : endif
603 : endif
604 : enddo
605 : end do
606 :
607 0 : end subroutine nat_sat_temp
608 :
609 0 : subroutine calc_radius_lbs( ncol, mask, sad_sage, radius_lbs )
610 :
611 : implicit none
612 :
613 : !----------------------------------------------------------------------
614 : ! ... dummy arguments
615 : !----------------------------------------------------------------------
616 : integer, intent(in) :: ncol
617 : real(r8), intent(in) :: sad_sage(ncol,pver)
618 : real(r8), intent(out) :: radius_lbs(ncol,pver)
619 : logical, intent(in) :: mask(ncol,pver)
620 :
621 : !----------------------------------------------------------------------
622 : ! ... local variables
623 : !----------------------------------------------------------------------
624 : integer :: k
625 0 : real(r8) :: lbs_vol_dens(ncol) ! Vol Density (cm3 aer / cm3 air)
626 :
627 : !----------------------------------------------------------------------
628 : ! ... parameters
629 : !----------------------------------------------------------------------
630 : real(r8), parameter :: lbs_part_dens = 10._r8, sigma_lbs = 1.6_r8
631 :
632 : !----------------------------------------------------------------------
633 : ! ... calculate the volume density (cm3 aerosol / cm3 air)
634 : ! calculate the mean radius for binary soln
635 : !----------------------------------------------------------------------
636 0 : do k = sad_topp,pver
637 0 : where( mask(:,k) )
638 : lbs_vol_dens(:) = ((sad_sage(:,k)**1.5_r8)/3._r8)/sqrt( four_pi*lbs_part_dens ) &
639 : *exp( 1.5_r8*(log( sigma_lbs ))**2 )
640 : radius_lbs(:,k) = (3._r8*lbs_vol_dens(:)/(four_pi*lbs_part_dens))**one_thrd &
641 : *exp( -1.5_r8*(log( sigma_lbs ))**2 )
642 : endwhere
643 : end do
644 :
645 0 : end subroutine calc_radius_lbs
646 :
647 0 : subroutine ice_sad_calc( ncol, press, temp, m, h2o_avail, &
648 0 : sad_ice, radius_ice, mask )
649 :
650 : implicit none
651 :
652 : !----------------------------------------------------------------------
653 : ! ... dummy arguments
654 : !----------------------------------------------------------------------
655 : integer, intent(in) :: ncol
656 : real(r8), intent(in) :: press (ncol,pver)
657 : real(r8), intent(in) :: temp (pcols,pver)
658 : real(r8), intent(in) :: m (ncol,pver)
659 : real(r8), intent(in) :: h2o_avail (ncol,pver)
660 : real(r8), intent(out) :: sad_ice (ncol,pver)
661 : real(r8), intent(out) :: radius_ice(ncol,pver)
662 : logical, intent(in) :: mask (ncol,pver)
663 :
664 : !----------------------------------------------------------------------
665 : ! ... local variables
666 : !----------------------------------------------------------------------
667 : real(r8), parameter :: &
668 : avo_num = 6.02214e23_r8, &
669 : aconst = -2663.5_r8, &
670 : bconst = 12.537_r8, &
671 : ice_mass_dens = 1._r8, &
672 : ice_part_dens = 1.e-1_r8, &
673 : mwh2o = 18._r8, &
674 : sigma_ice = 1.6_r8, &
675 : ice_dens_aer = ice_mass_dens / (mwh2o/avo_num), &
676 : ice_dens_aeri = 1._r8/ice_dens_aer
677 :
678 : integer :: k
679 0 : real(r8) :: h2o_cond_ice(ncol) ! Condensed phase H2O (from CAM)
680 0 : real(r8) :: voldens_ice (ncol) ! Volume Density, um3 cm-3
681 :
682 0 : do k = sad_topp,pver
683 0 : where( mask(:,k) )
684 : !----------------------------------------------------------------------
685 : ! .... Convert condensed phase to molecules cm-3 units
686 : !----------------------------------------------------------------------
687 : h2o_cond_ice(:) = h2o_avail(:,k) * m(:,k)
688 : !----------------------------------------------------------------------
689 : ! .... ICE volume density .....
690 : !----------------------------------------------------------------------
691 : voldens_ice(:) = h2o_cond_ice(:)*ice_dens_aeri
692 : !----------------------------------------------------------------------
693 : ! .... Calculate the SAD from log normal distribution .....
694 : !----------------------------------------------------------------------
695 : sad_ice(:,k) = (four_pi*ice_part_dens)**one_thrd &
696 : *(3._r8*voldens_ice(:))**two_thrd &
697 : *exp( -(log( sigma_ice ))**2 )
698 : !----------------------------------------------------------------------
699 : ! .... Calculate the radius from log normal distribution .....
700 : !----------------------------------------------------------------------
701 : radius_ice(:,k) = (3._r8*h2o_cond_ice(:) &
702 : /(ice_dens_aer*four_pi*ice_part_dens))**one_thrd &
703 : *exp( -1.5_r8*(log( sigma_ice ))**2 )
704 : endwhere
705 : end do
706 :
707 0 : end subroutine ice_sad_calc
708 :
709 0 : subroutine sulfate_sad_calc( ncol, press, temp, h2o_avail, hno3_avail, hcl_avail, &
710 0 : sad_sage, m, hno3_gas, hno3_cond, &
711 0 : hcl_gas, hcl_cond, sad_sulfate, &
712 0 : radius_sulfate, mask, lchnk, flag, h2so4m, is_chem)
713 : implicit none
714 :
715 : !----------------------------------------------------------------------
716 : ! ... dummy arguments
717 : !----------------------------------------------------------------------
718 : integer, intent(in) :: ncol
719 : integer, intent(in) :: lchnk, flag
720 : real(r8), intent(in) :: temp (pcols,pver)
721 : real(r8), intent(in) :: press (ncol,pver)
722 : real(r8), intent(in) :: m (ncol,pver)
723 : real(r8), intent(in) :: h2o_avail (ncol,pver)
724 : real(r8), intent(in) :: hno3_avail (ncol,pver)
725 : real(r8), intent(in) :: hcl_avail (ncol,pver)
726 : real(r8), intent(in) :: sad_sage (ncol,pver)
727 : real(r8), intent(out) :: hno3_gas (ncol,pver) ! Gas-phase HNO3, mole fraction
728 : real(r8), intent(out) :: hno3_cond (ncol,pver) ! Condensed phase HNO3, mole fraction
729 : real(r8), intent(out) :: hcl_gas (ncol,pver) ! Gas-phase HCL, mole fraction
730 : real(r8), intent(out) :: hcl_cond (ncol,pver) ! Condensed phase HCL, mole fraction
731 : real(r8), intent(out) :: sad_sulfate(ncol,pver)
732 : real(r8), intent(out) :: radius_sulfate(ncol,pver)
733 : real(r8), intent(inout) :: h2so4m (ncol,pver) ! mass per volume, micro grams m-3
734 : logical, intent(in) :: is_chem ! chemistry calc switch
735 : logical, intent(in) :: mask (ncol,pver)
736 :
737 : !----------------------------------------------------------------------
738 : ! ... local variables
739 : !----------------------------------------------------------------------
740 : real(r8), parameter :: t_limit = 200._r8
741 : real(r8), parameter :: avo_num = 6.02214e23_r8
742 : real(r8), parameter :: mwh2so4 = 98.076_r8
743 : real(r8), parameter :: sigma_sulfate = 1.6_r8
744 : real(r8), parameter :: sulfate_part_dens = 10._r8
745 :
746 : integer :: i, k
747 : integer :: cnt1, cnt2
748 : real(r8) :: ratio
749 : real(r8) :: vals(2)
750 0 : real(r8) :: h2so4_aer_dens (ncol,pver) ! grams cm-3 solution
751 0 : real(r8) :: h2so4_cond (ncol,pver) ! Condensed H2SO4 (moles cm-3 air)
752 0 : real(r8) :: sulfate_vol_dens(ncol,pver) ! Volume Density, cm3 aerosol cm-3 air
753 0 : real(r8) :: wtf (ncol,pver) ! weight fraction of H2SO4 in ternary soln
754 0 : real(r8) :: wts (ncol,pver) ! weight percent of ternary solution
755 :
756 : ! Carslaw HCl solubility
757 0 : real(r8) :: wts0 (ncol,pver) ! weight percent of H2SO4 is LBS
758 0 : real(r8) :: wtn (ncol,pver) ! weight percent of HNO3 in STS
759 : real(r8) :: ch2so4 (ncol,pver) ! Total H2SO4 (moles / cm3 of air)
760 0 : real(r8) :: molh2so4 (ncol,pver) ! Equil molality of H2SO4 in STS
761 0 : real(r8) :: molhno3 (ncol,pver) ! Equil molality of HNO3 in STS
762 0 : real(r8) :: AD (ncol,pver) ! air density (molecules cm-3)
763 0 : real(r8) :: xmf (ncol,pver) !
764 0 : real(r8) :: hhcl (ncol,pver) ! henry's solubility of hcl in binary
765 0 : real(r8) :: phcl0 (ncol,pver) ! partial pressure of hcl (hPa)
766 0 : real(r8) :: h2so4vmr (ncol,pver) ! atmospheric mole fraction of H2SO4
767 0 : real(r8) :: nsul (ncol,pver) ! moles / m3- H2SO4 pure liquid
768 0 : real(r8) :: mcl (ncol,pver) ! molality of hcl in ?
769 0 : real(r8) :: wtcl (ncol,pver) !
770 0 : real(r8) :: phcl (ncol,pver) ! partial pressure of hcl (over aerosol)
771 0 : real(r8) :: parthcl (ncol,pver) ! fraction of HCl in gas-phase
772 : !
773 : real(r8) :: packer (ncol*pver)
774 : logical :: do_equil ! local mask
775 0 : logical :: mask_lt (ncol,pver) ! local temperature mask
776 : logical :: maskx (ncol,pver)
777 0 : logical :: converged (ncol,pver) ! EQUIL convergence test
778 :
779 0 : hcl_gas (:,:) = 0.0_r8
780 0 : hcl_cond(:,:) = 0.0_r8
781 0 : parthcl (:,:) = 0.0_r8
782 0 : phcl0 (:,:) = 0.0_r8
783 :
784 0 : do k = sad_topp,pver
785 0 : mask_lt(:,k) = mask(:,k)
786 : end do
787 :
788 : !----------------------------------------------------------------------
789 : ! ... derive H2SO4 (micro grams / m3) from SAGEII SAD
790 : !----------------------------------------------------------------------
791 : call sad2h2so4( h2o_avail, press, sad_sage, temp, sulfate_vol_dens, &
792 0 : h2so4_aer_dens, h2so4m, mask, ncol )
793 :
794 : !----------------------------------------------------------------------
795 : ! ... limit h2so4m
796 : !----------------------------------------------------------------------
797 :
798 0 : do k = sad_topp,pver
799 0 : do i = 1,ncol
800 0 : if( mask(i,k) ) then
801 0 : if( h2so4m(i,k) <= 0._r8 ) then
802 0 : h2so4m(i,k) = 1.e-2_r8
803 : end if
804 : end if
805 : end do
806 : end do
807 :
808 0 : if( is_chem ) then
809 : else
810 0 : do k = sad_topp,pver
811 0 : where( mask(:ncol,k) )
812 0 : mask_lt(:ncol,k) = temp(:ncol,k) < t_limit
813 : end where
814 : end do
815 : end if
816 :
817 0 : if( is_chem ) then
818 : do_equil = .true.
819 : else
820 0 : do_equil = any( mask_lt(:,sad_topp:pver) )
821 : end if
822 :
823 : !----------------------------------------------------------------------
824 : ! .... Calculate the ternary soln volume density
825 : !----------------------------------------------------------------------
826 0 : if( do_equil ) then
827 :
828 : call equil( temp, h2so4m, hno3_avail, h2o_avail, press, &
829 : hno3_cond, h2so4_cond, wts, wtn, wts0, molh2so4, molhno3, mask_lt, ncol, &
830 0 : lchnk, flag, is_chem, converged )
831 :
832 0 : do k = sad_topp,pver
833 :
834 0 : where( ( mask_lt(:ncol,k) ) .AND. ( converged(:ncol,k) ) )
835 :
836 : !----------------------------------------------------------------------
837 : ! .... convert h2o, hno3 from moles cm-3 air to molecules cm-3 air
838 : !----------------------------------------------------------------------
839 : hno3_cond(:ncol,k) = min( hno3_cond(:ncol,k),hno3_avail(:ncol,k) )
840 : hno3_gas(:ncol,k) = hno3_avail(:ncol,k) - hno3_cond(:ncol,k)
841 :
842 : !----------------------------------------------------------------------
843 : ! .... Derive ternary volume density (cm3 aer / cm3 air)
844 : !----------------------------------------------------------------------
845 : wtf(:ncol,k) = .01_r8* wts(:ncol,k)
846 : sulfate_vol_dens(:ncol,k) = h2so4_cond(:ncol,k)*mwh2so4/(wtf(:ncol,k)*h2so4_aer_dens(:ncol,k))
847 :
848 : ! Carslaw solubility
849 : !----------------------------------------------------------------------
850 : ! .... Partition HCl (gas/condensed) *** Carslaw
851 : !----------------------------------------------------------------------
852 : ! THE SOLUBILITY OF HCL
853 : ! HHCl (MOL/KG/ATM) taken form Shi et al., JGR 2001
854 : !
855 :
856 : ! .... Convert weight % to weight fraction
857 : wtn(:ncol,k) = wtn(:ncol,k) * 0.01_r8
858 : wts0(:ncol,k) = wts0(:ncol,k) * 0.01_r8
859 :
860 : ! .... Derive xmf (mole fraction H2SO4 in LBS )
861 : xmf(:ncol,k) = (wts0(:ncol,k)*100.0_r8)/((wts0(:ncol,k)*100.0_r8)+ &
862 : (100.0_r8-(wts0(:ncol,k)*100._r8))*98.0_r8/18.0_r8)
863 :
864 : ! .... Derive hhcl (henry's solubility of hcl in binary)
865 : hhcl(:ncol,k) = (0.094_r8-0.61_r8*xmf(:ncol,k)+1.2_r8*xmf(:ncol,k)**2.0_r8) &
866 0 : *exp(-8.68_r8+(8515.0_r8-10718.0_r8*xmf(:ncol,k)**(0.7_r8))/temp(:ncol,k))
867 :
868 : ! .... Derive phcl0 (partial pressure of hcl( hPa))
869 : phcl0(:ncol,k) = hcl_avail(:ncol,k)*press(:ncol,k) / 1013.26_r8
870 :
871 : ! .... Derive H2SO4 vmr (h2so4_cond = mole / cm-3)
872 : AD(:ncol,k) = (6.022098e23_r8 * press(:ncol,k) / 1013.26_r8) &
873 : / (temp(:ncol,k)*8.2058e-2_r8*1000.0_r8)
874 : h2so4vmr(:ncol,k) = (h2so4_cond(:ncol,k)*6.022098e23_r8) / AD(:ncol,k)
875 :
876 : ! .... Derive nsul (moles / m3 H2SO4 pure liquid )
877 : nsul(:ncol,k) = h2so4vmr(:ncol,k) * press(:ncol,k) * 100.0_r8 / 8.314_r8 / temp(:ncol,k)
878 :
879 : ! .... Derive mcl (molality of hcl)
880 : mcl(:ncol,k) = (1.0_r8/8.314e-5_r8/temp(:ncol,k)*phcl0(:ncol,k))/(nsul(:ncol,k)/molh2so4(:ncol,k) + &
881 : 1.0_r8/(8.314e-5_r8)/temp(:ncol,k)/hhcl(:ncol,k))
882 :
883 : ! .... Derive wtcl ( )
884 : wtcl(:ncol,k) = mcl(:ncol,k)*36.5_r8/(1000.0_r8 + 98.12_r8*molh2so4(:ncol,k) + 63.03_r8*molhno3(:ncol,k))
885 :
886 : ! .... Derive phcl (partial pressure over the aerosol)
887 : phcl(:ncol,k) = mcl(:ncol,k)/hhcl(:ncol,k)
888 :
889 : ! .... Derive parhcl (fraction of HCl in gas-phase)
890 : where(phcl0(:ncol,k)>0._r8)
891 : parthcl(:ncol,k) = 1.0_r8 - (phcl0(:ncol,k) - phcl(:ncol,k)) / phcl0(:ncol,k)
892 : elsewhere
893 : parthcl(:ncol,k) = 0._r8
894 : endwhere
895 :
896 : ! .... Partition HCl (gas/condensed)
897 : hcl_gas (:ncol,k) = hcl_avail(:ncol,k) * parthcl(:ncol,k)
898 : hcl_cond(:ncol,k) = hcl_avail(:ncol,k) - hcl_gas(:ncol,k)
899 :
900 : elsewhere
901 : hno3_cond(:ncol,k) = 0.0_r8
902 : hno3_gas(:ncol,k) = hno3_avail(:ncol,k)
903 : hcl_cond (:ncol,k) = 0.0_r8
904 : hcl_gas (:ncol,k) = hcl_avail(:ncol,k)
905 :
906 : endwhere
907 : end do
908 :
909 : end if
910 :
911 0 : do k = sad_topp,pver
912 0 : where( mask(:,k) )
913 : !----------------------------------------------------------------------
914 : ! .... Calculate the SAD (assuming ternary solution)
915 : !----------------------------------------------------------------------
916 : sad_sulfate(:,k) = (four_pi*sulfate_part_dens)**one_thrd &
917 : *(3._r8*sulfate_vol_dens(:,k))**two_thrd &
918 : *exp( -(log( sigma_sulfate ))**2 )
919 :
920 : !----------------------------------------------------------------------
921 : ! .... Calculate the radius (assuming ternary solution) (in cm?)
922 : !----------------------------------------------------------------------
923 : radius_sulfate(:,k) = (3._r8*sulfate_vol_dens(:,k) &
924 : /(four_pi*sulfate_part_dens))**one_thrd &
925 : *exp( -1.5_r8*(log( sigma_sulfate ))**2 )
926 : endwhere
927 :
928 : end do
929 :
930 0 : end subroutine sulfate_sad_calc
931 :
932 :
933 0 : subroutine nat_sad_calc( ncol, press, temp, h2o_avail, hno3_avail, m, &
934 0 : hno3_gas, hno3_cond, sad_nat, radius_nat, mask )
935 :
936 : implicit none
937 :
938 : !----------------------------------------------------------------------
939 : ! ... dummy arguments
940 : !----------------------------------------------------------------------
941 : integer, intent(in) :: ncol
942 : real(r8), intent(in) :: press (ncol,pver)
943 : real(r8), intent(in) :: m (ncol,pver)
944 : real(r8), intent(in) :: temp (pcols,pver)
945 : real(r8), intent(in) :: h2o_avail (ncol,pver)
946 : real(r8), intent(in) :: hno3_avail(ncol,pver)
947 : real(r8), intent(out) :: hno3_cond (ncol,pver) ! HNO3 in condensed phase (mole fraction)
948 : real(r8), intent(out) :: hno3_gas (ncol,pver) ! HNO3 in gas-phase (mole fraction)
949 : real(r8), intent(out) :: sad_nat (ncol,pver)
950 : real(r8), intent(out) :: radius_nat(ncol,pver)
951 : logical, intent(in) :: mask(ncol,pver) ! grid mask
952 :
953 : !----------------------------------------------------------------------
954 : ! ... local variables
955 : !----------------------------------------------------------------------
956 : integer :: k, i
957 0 : real(r8) :: nat_dens_condphase(ncol, pver) ! Condensed phase NAT, molec cm-3
958 0 : real(r8) :: voldens_nat (ncol, pver) ! Volume Density, um3 cm-3
959 0 : real(r8) :: hno3_cond_total (ncol, pver) ! Total Condensed phase HNO3
960 :
961 : !----------------------------------------------------------------------
962 : ! ... parameters
963 : !----------------------------------------------------------------------
964 : real(r8), parameter :: avo_num = 6.02214e23_r8, &
965 : nat_mass_dens = 1.6_r8, &
966 : nat_part_dens = 5.0e-4_r8, &
967 : mwnat = 117._r8, &
968 : sigma_nat = 1.6_r8, &
969 : nat_dens_aer = nat_mass_dens / (mwnat/avo_num), &
970 : nat_dens_aeri = 1._r8/nat_dens_aer
971 :
972 : !----------------------------------------------------------------------
973 : ! ... Derive HNO3 paritioning (call A. Tabazedeh routine for NAT)
974 : !----------------------------------------------------------------------
975 : call nat_cond( ncol, press, temp, h2o_avail, hno3_avail, &
976 0 : hno3_gas, hno3_cond_total, mask )
977 :
978 0 : do k = sad_topp,pver
979 0 : do i = 1,ncol
980 0 : masked : if( mask(i,k) ) then
981 :
982 : !----------------------------------------------------------------------
983 : ! .... Set Condensed phase for return arguments
984 : !----------------------------------------------------------------------
985 0 : hno3_cond(i,k) = hno3_cond_total(i,k)
986 :
987 : !----------------------------------------------------------------------
988 : ! .... Calculated Condensed Phase NAT (i.e. HNO3) in
989 : ! molecules cm-3 of air units.
990 : !----------------------------------------------------------------------
991 0 : nat_dens_condphase(i,k) = hno3_cond_total(i,k) * m(i,k)
992 :
993 : !----------------------------------------------------------------------
994 : ! .... Calculate the Volume Density
995 : !----------------------------------------------------------------------
996 0 : voldens_nat(i,k) = nat_dens_condphase(i,k) * nat_dens_aeri
997 :
998 : !----------------------------------------------------------------------
999 : ! .... Calculate the SAD from log normal distribution
1000 : ! .... Assuming sigma and nad_part_dens (# particles per cm3 of air)
1001 : !----------------------------------------------------------------------
1002 : sad_nat(i,k) = (four_pi*nat_part_dens)**(one_thrd) &
1003 : *(3._r8*voldens_nat(i,k))**(two_thrd) &
1004 0 : *exp( -(log( sigma_nat )**2 ))
1005 :
1006 : !----------------------------------------------------------------------
1007 : ! .... Calculate the radius of NAT from log normal distribution
1008 : ! .... Assuming sigma and nat_part_dens (# particles per cm3
1009 : ! .... of air)
1010 : !----------------------------------------------------------------------
1011 : radius_nat(i,k) = (3._r8*nat_dens_condphase(i,k) &
1012 : /(nat_dens_aer*four_pi*nat_part_dens))**(one_thrd) &
1013 0 : *exp( -1.5_r8*(log( sigma_nat ))**2 )
1014 :
1015 : end if masked
1016 : end do
1017 : end do
1018 :
1019 0 : end subroutine nat_sad_calc
1020 :
1021 0 : subroutine nat_cond( ncol, press, temp, h2o_avail, hno3_avail, &
1022 0 : hno3_gas, hno3_cond, mask )
1023 :
1024 : implicit none
1025 :
1026 : !----------------------------------------------------------------------
1027 : ! ... dummy arguments
1028 : !----------------------------------------------------------------------
1029 : integer, intent(in) :: ncol
1030 : real(r8), intent(in) :: press(ncol,pver)
1031 : real(r8), intent(in) :: temp(pcols,pver)
1032 : real(r8), intent(in) :: h2o_avail(ncol,pver)
1033 : real(r8), intent(in) :: hno3_avail(ncol,pver)
1034 : real(r8), intent(out) :: hno3_gas(ncol,pver)
1035 : real(r8), intent(out) :: hno3_cond(ncol,pver)
1036 : logical, intent(in) :: mask(ncol,pver)
1037 :
1038 : !----------------------------------------------------------------------
1039 : ! ... local variables
1040 : !----------------------------------------------------------------------
1041 : real(r8), parameter :: aa = -2.7836_r8, &
1042 : bb = -0.00088_r8, &
1043 : cc = 38.9855_r8, &
1044 : dd = -11397.0_r8, &
1045 : ee = 0.009179_r8
1046 :
1047 : integer :: i, k
1048 : real(r8) :: bt ! temporary variable
1049 : real(r8) :: mt ! temporary variable
1050 : real(r8) :: t ! temporary variable
1051 : real(r8) :: logPhno3 ! temporary variable
1052 : real(r8) :: phno3 ! hno3 partial pressure
1053 : real(r8) :: ph2o ! h2o partial pressure
1054 : real(r8) :: phno3_eq ! partial pressure above NAT
1055 : real(r8) :: wrk
1056 :
1057 0 : do k = sad_topp,pver
1058 0 : do i = 1,ncol
1059 : !----------------------------------------------------------------------
1060 : ! .... Derive HNO3 and H2O partial pressure (torr)
1061 : ! where: 0.7501 = 760/1013.
1062 : !----------------------------------------------------------------------
1063 0 : if( mask(i,k) ) then
1064 0 : wrk = press(i,k) * .7501_r8
1065 0 : phno3 = hno3_avail(i,k) * wrk
1066 0 : ph2o = h2o_avail(i,k) * wrk
1067 : !----------------------------------------------------------------------
1068 : ! Calculating the temperature coefficients for the variation of HNO3
1069 : ! and H2O vapor pressure (torr) over a trihydrate solution of HNO3/H2O
1070 : ! The coefficients are taken from Hanson and Mauersberger:
1071 : ! GRL, Vol.15, 8, p855-858, 1988.
1072 : !----------------------------------------------------------------------
1073 0 : t = temp(i,k)
1074 0 : bt = cc + dd/t + ee*t
1075 0 : mt = aa + bb*t
1076 :
1077 0 : logphno3 = mt*log10( ph2o ) + bt
1078 0 : phno3_eq = 10._r8**logphno3
1079 :
1080 0 : if( phno3 > phno3_eq ) then
1081 0 : wrk = 1._r8 / wrk
1082 0 : hno3_cond(i,k) = (phno3 - phno3_eq) * wrk
1083 0 : hno3_gas(i,k) = phno3_eq * wrk
1084 : else
1085 0 : hno3_cond(i,k) = 0._r8
1086 0 : hno3_gas(i,k) = hno3_avail(i,k)
1087 : end if
1088 : end if
1089 : end do
1090 : end do
1091 :
1092 0 : end subroutine nat_cond
1093 :
1094 0 : subroutine sad2h2so4( h2o, press, sad_sage, temp, lbs_vol_dens, &
1095 0 : h2so4_aer_dens, h2so4m, mask, ncol )
1096 :
1097 : implicit none
1098 :
1099 : !----------------------------------------------------------------------
1100 : ! ... dummy arguments
1101 : !----------------------------------------------------------------------
1102 : integer, intent(in) :: ncol ! columns in chunk
1103 : real(r8), intent(in) :: h2o(ncol,pver) ! h2o mole fraction
1104 : real(r8), intent(in) :: sad_sage(ncol,pver) ! sad from SAGEII cm2 aer, cm-3 air
1105 : real(r8), intent(in) :: press(ncol,pver) ! pressure (hPa)
1106 : real(r8), intent(in) :: temp(pcols,pver) ! temperature (K)
1107 : real(r8), intent(inout) :: h2so4m(ncol,pver) ! microgram/m**3 of air,
1108 : real(r8), intent(out) :: h2so4_aer_dens(ncol,pver) ! units: grams / cm3-aerosol
1109 : real(r8), intent(out) :: lbs_vol_dens(ncol,pver) ! cm3 aer / cm3 air
1110 : logical, intent(in) :: mask(ncol,pver) ! activation mask
1111 :
1112 : !----------------------------------------------------------------------
1113 : ! ... local variables
1114 : !----------------------------------------------------------------------
1115 : real(r8), parameter :: lbs_part_dens = 10._r8
1116 : real(r8), parameter :: sigma_lbs = 1.6_r8
1117 : real(r8), parameter :: t_floor = 180._r8
1118 :
1119 : integer :: i, k, l
1120 : real(r8) :: wts0
1121 : real(r8) :: p ! pressure, torr
1122 : real(r8) :: tr ! inverse temperature
1123 : real(r8) :: c(6)
1124 :
1125 : !----------------------------------------------------------------------
1126 : ! ... Calculate the volume density (cm3 aerosol / cm3 air)
1127 : !----------------------------------------------------------------------
1128 0 : do k = sad_topp,pver
1129 0 : do i = 1,ncol
1130 0 : if( mask(i,k) ) then
1131 : lbs_vol_dens(i,k) = ((sad_sage(i,k)**1.5_r8)/3._r8)/sqrt( four_pi*lbs_part_dens ) &
1132 0 : *exp( 1.5_r8*(log( sigma_lbs ))**2 )
1133 : !----------------------------------------------------------------------
1134 : ! ... calculate Molality from Tabazadeh EQUIL routine (binary soln)
1135 : !----------------------------------------------------------------------
1136 : ! ... DEK, added a minimum to temperature
1137 : !----------------------------------------------------------------------
1138 0 : p = h2o(i,k) * press(i,k) * .7501_r8
1139 0 : tr = 1._r8 / max( t_floor,temp(i,k) )
1140 :
1141 0 : do l = 1,6
1142 0 : c(l) = exp( a(1,l) + tr*(a(2,l) + tr*(a(3,l) + tr*(a(4,l) + tr*a(5,l)))) )
1143 : end do
1144 : !----------------------------------------------------------------------
1145 : ! ... H2SO4/H2O pure weight percent and molality (moles gram-1)
1146 : !----------------------------------------------------------------------
1147 0 : wts0 = max( 0._r8,c(1) + p*(-c(2) + p*(c(3) + p*(-c(4) + p*(c(5) - p*c(6))))) )
1148 : !----------------------------------------------------------------------
1149 : ! ... derive weight fraction for density routine
1150 : !----------------------------------------------------------------------
1151 0 : wts0 = .01_r8 *wts0
1152 : !----------------------------------------------------------------------
1153 : ! ... calculate the binary solution density, grams / cm3-aerosol
1154 : !----------------------------------------------------------------------
1155 0 : h2so4_aer_dens(i,k) = max( 0._r8,density( temp(i,k), wts0 ) )
1156 : !----------------------------------------------------------------------
1157 : ! ... calculate the H2SO4 micrograms m-3 abundance for binary soln
1158 : !----------------------------------------------------------------------
1159 0 : h2so4m(i,k) = lbs_vol_dens(i,k)*h2so4_aer_dens(i,k)*wts0*1.e12_r8
1160 : end if
1161 : end do
1162 : end do
1163 :
1164 0 : end subroutine sad2h2so4
1165 :
1166 : !======================================================================
1167 : !
1168 : ! ROUTINE
1169 : ! EQUIL
1170 : !
1171 : ! Date...
1172 : ! 7 October 1999
1173 : !
1174 : ! Programmed by...
1175 : ! A. Tabazadeh
1176 : !
1177 : ! DESCRIPTION
1178 : ! Ternary solution routine
1179 : !
1180 : ! ARGUMENTS
1181 : !
1182 : !.... INPUT:
1183 : !
1184 : ! H2SO4m = microgram/m**3 of air
1185 : ! HNO3r = mole fraction
1186 : ! H2Or = mole fraction
1187 : ! PTOTAL = hPa
1188 : !
1189 : !.... Output
1190 : !
1191 : ! Cwater = Total moles of liguid water / cm**3 of air
1192 : ! hno3_cond = HNO3 Condensed phase (mole fraction)
1193 : ! CH2SO4 = Total H2SO4 moles / cm**3 of air
1194 : ! WTS = Weight percent of H2SO4 in the ternary aerosol
1195 : !
1196 : !======================================================================
1197 0 : subroutine equil( temper, h2so4m, hno3_avail, h2o_avail, press, &
1198 0 : hno3_cond, ch2so4, wts, wtn, wts0, molh2so4, molhno3, mask, ncol, &
1199 0 : lchnk, flag, is_chem, converged)
1200 : !----------------------------------------------------------------------
1201 : ! Written by Azadeh Tabazadeh (1993)
1202 : ! (modified from EQUISOLV -- M. Z. Jacobson -- see below)
1203 : ! NASA Ames Research Center , Tel. (415) 604 - 1096
1204 : !
1205 : ! This program solves the equilibrium composition for the ternary
1206 : ! system of H2SO4/HNO3/H2O under typical stratospheric conditions.
1207 : ! The formulation of this work is described by Tabazadeh, A.,
1208 : ! Turco, R. P., and Jacobson, M. Z. (1994), "A model for studying
1209 : ! the composition and chemical effects of stratospheric aerosols,"
1210 : ! J. Geophys. Res., 99, 12,897, 1994. *
1211 : !
1212 : ! The solution mechanism for the equilibrium equations is des-
1213 : ! cribed by Jacobson, M. Z., Turco, R. P., and Tabazadeh, A.
1214 : ! (1994), "Simulating Equilibrium within aerosols and non-equil-
1215 : ! ibrium between gases and aerosols," J. Geophys. Res., in review.
1216 : ! The mechanism is also codified in the fortran program, EQUISOLV,
1217 : ! by M.Z. Jacobson (1991-3). EQUISOLV solves any number of
1218 : ! gas / liquid / solid / ionic equilibrium equations simultan-
1219 : ! eously and includes treatment of the water equations and act-
1220 : ! ivity coefficients. The activity coeffients currently in
1221 : ! EQUISOLV are valid for tropospheric temperatures. The acitiv-
1222 : ! ities listed here are valid for stratospheric temperatures only.
1223 : !
1224 : ! DEFINING PARAMETERS
1225 : !
1226 : ! *NOTE* Solver parameters including, F, Z, QN, QD, and deltaX
1227 : ! are described in Jacobson et al.
1228 : !
1229 : ! PTOTAL = Total atmospheric pressure in mb
1230 : ! H2SO4m = Total mass of H2SO4 (microgram/m**3 of air)
1231 : ! HNO3r = HNO3 mixing ratio
1232 : ! H2Or = H2O mixing ratio
1233 : ! P = Partial pressure of water in units of torr
1234 : ! pures = molality for a pure H2SO4/H2O system
1235 : ! puren = molality for a pure HNO3/H2O sytem
1236 : ! WTS0 = weight percent of H2SO4 in a pure H2SO4/H2O system
1237 : ! WTN0 = weight percent of HNO3 in a pure HNO3/H2O system
1238 : ! WTS = weight percent of H2SO4 in the ternary aerosol
1239 : ! WTN = weight percent of HNO3 in the ternary aerosol
1240 : ! PHNO3 = HNO3 vapor pressure over the ternary system in atm
1241 : ! HNO3 = HNO3 vapor concentration over the ternary system (#/cm3)
1242 : ! CH2SO4 = Total H2SO4 moles / cm**3 of air
1243 : ! CHNO3 = Total HNO3 moles / cm**3 of air
1244 : ! CHplus = Total H+ moles / cm**3 0f air
1245 : ! CPHNO3 = Total moles of HNO3 gas / cm**3 of air
1246 : ! CNO3 = Total moles of NO3- / cm**3 0f air
1247 : ! Cwater = Total moles of liguid water / cm**3 of air
1248 : ! KS = Solubility constant for dissolution of HNO3 in
1249 : ! water ( HNO3(gas) === H+(aq) + NO3- (aq) )
1250 : ! nm = HNO3 molality at the STREN of the ternary solution
1251 : ! sm = H2SO4 molality at the STREN of the ternary solution
1252 : ! molHNO3 = Equilibrium molality of HNO3 in the ternary solution
1253 : ! molH2SO4= Equilibrium molality of H2SO4 in the ternary solution
1254 : ! STREN = ionic strenght for the ternary solutin, which in
1255 : ! this case is = 3 * molH2SO4 + molHNO3
1256 : ! acts = Pure mean binary activity coefficient for the H2SO4/
1257 : ! H2O system evaluated at the STREN of the ternary system
1258 : ! actn = Pure mean binary activity coefficient for the HNO3/
1259 : ! H2O system evaluated at the STREN of the ternary system
1260 : ! ymix = Mixed binary activity coefficient for the HNO3/H2O in
1261 : ! the ternary solution
1262 : !----------------------------------------------------------------------
1263 :
1264 : use cam_abortutils, only : endrun
1265 :
1266 : implicit none
1267 :
1268 : !----------------------------------------------------------------------
1269 : ! ... dummy arguments
1270 : !----------------------------------------------------------------------
1271 : integer, intent(in) :: lchnk
1272 : integer, intent(in) :: flag
1273 : integer, intent(in) :: ncol ! columns in chunk
1274 : real(r8), intent(in) :: h2so4m(ncol,pver)
1275 : real(r8), intent(in) :: hno3_avail(ncol,pver)
1276 : real(r8), intent(in) :: h2o_avail(ncol,pver)
1277 : real(r8), intent(in) :: press(ncol,pver)
1278 : real(r8), intent(in) :: temper(pcols,pver)
1279 : real(r8), intent(out) :: hno3_cond(ncol,pver)
1280 : real(r8), intent(out) :: ch2so4(ncol,pver)
1281 : real(r8), intent(out) :: wts(ncol,pver)
1282 : real(r8), intent(out) :: wtn(ncol,pver)
1283 : real(r8), intent(out) :: wts0(ncol,pver)
1284 : real(r8), intent(out) :: molh2so4(ncol,pver)
1285 : real(r8), intent(out) :: molhno3(ncol,pver)
1286 : logical, intent(in) :: is_chem
1287 : logical, intent(in) :: mask(ncol,pver) ! activation mask
1288 : logical, intent(out) :: converged(ncol,pver)
1289 : !----------------------------------------------------------------------
1290 : ! ... local variables
1291 : !----------------------------------------------------------------------
1292 : ! integer, parameter :: itermax = 50
1293 : integer, parameter :: itermax = 100
1294 : real(r8), parameter :: con_lim = .00005_r8
1295 : real(r8), parameter :: t0 = 298.15_r8
1296 : real(r8), parameter :: ks0 = 2.45e6_r8
1297 : real(r8), parameter :: lower_delx = 1.e-10_r8
1298 : real(r8), parameter :: upper_delx = .98_r8
1299 : real(r8), parameter :: con_crit_chem = 5.e-5_r8
1300 :
1301 : integer :: i, iter, k, l, nstep
1302 : real(r8) :: reduction_factor
1303 : real(r8) :: p
1304 : real(r8) :: tr
1305 : real(r8) :: wtn0
1306 : real(r8) :: pures
1307 : real(r8) :: puren
1308 : real(r8) :: chno3
1309 : real(r8) :: chplus
1310 : real(r8) :: cno3
1311 : real(r8) :: wrk
1312 : real(r8) :: z, num, den
1313 : real(r8) :: deltax
1314 : real(r8) :: chplusnew
1315 : real(r8) :: cno3new
1316 : real(r8) :: stren
1317 : real(r8) :: sm
1318 : real(r8) :: actn
1319 : real(r8) :: acts
1320 : real(r8) :: nm
1321 : real(r8) :: ks
1322 : real(r8) :: lnks
1323 : real(r8) :: lnks0
1324 : real(r8) :: mixyln
1325 : real(r8) :: wrk_h2so4
1326 : real(r8) :: cphno3new
1327 : real(r8) :: con_val
1328 : real(r8) :: t, t1, t2, f, f1, f2, ymix, hplus, wtotal, ratio
1329 : real(r8) :: con_crit
1330 0 : real(r8) :: h2o_cond(ncol,pver)
1331 : real(r8) :: fratio(0:itermax)
1332 : real(r8) :: delx(0:itermax)
1333 : real(r8) :: delz(0:itermax)
1334 : real(r8) :: c(12)
1335 : real(r8) :: d(13:22)
1336 : logical :: interval_set
1337 : logical :: positive
1338 :
1339 0 : converged(:,:) = .false.
1340 :
1341 0 : lnks0 = log( ks0 )
1342 : if( is_chem ) then
1343 : con_crit = con_crit_chem
1344 : else
1345 : con_crit = con_crit_chem
1346 : end if
1347 0 : Level_loop : do k = sad_topp,pver
1348 0 : Column_loop : do i = 1,ncol
1349 0 : if( mask(i,k) ) then
1350 0 : p = h2o_avail(i,k) * press(i,k) * .7501_r8
1351 : !----------------------------------------------------------------------
1352 : ! Calculating the molality for pure binary systems of H2SO4/H2O
1353 : ! and HNO3/H2O at a given temperature and water vapor pressure
1354 : ! profile (relative humiditiy). Water activities were used to
1355 : ! calculate the molalities as described in Tabazadeh et al. (1994).
1356 : !----------------------------------------------------------------------
1357 0 : t = max( 180._r8,temper(i,k) )
1358 0 : tr = 1._r8/t
1359 0 : do l = 1,12
1360 0 : c(l) = exp( a(1,l) + tr*(a(2,l) + tr*(a(3,l) + tr*(a(4,l) + tr*a(5,l)))) )
1361 : end do
1362 : !----------------------------------------------------------------------
1363 : ! ... H2SO4/H2O pure weight percent and molality
1364 : !----------------------------------------------------------------------
1365 0 : wts0(i,k) = max( 0.01_r8,c(1) + p*(-c(2) + p*(c(3) + p*(-c(4) + p*(c(5) - p*c(6))))) )
1366 0 : pures = (wts0(i,k) * 1000._r8)/(100._r8 - wts0(i,k))
1367 0 : pures = pures / 98._r8
1368 : !----------------------------------------------------------------------
1369 : ! ... HNO3/H2O pure weight percent and molality
1370 : !----------------------------------------------------------------------
1371 0 : puren = max( 0._r8,c(7) + p*(-c(8) + p*(c(9) + p*(-c(10) + p*(c(11) - p*c(12))))) )
1372 : ! wtn0 = (puren * 6300._r8) /(puren * 63._r8 + 1000._r8)
1373 : !----------------------------------------------------------------------
1374 : ! The solving scheme is described both in Jacobson et al. and Tabazadeh
1375 : ! et al.. Assumptions:
1376 : ! (1) H2SO4 is present only in the aqueous-phase
1377 : ! (2) H2SO4 and HNO3 in solution are fully dissocated into H+
1378 : ! SO42- and NO3-
1379 : ! (3) PHNO3 + NO3- = constant
1380 : !----------------------------------------------------------------------
1381 0 : ch2so4(i,k) = (h2so4m(i,k)*1.e-12_r8) / 98._r8
1382 0 : if( pures > 0._r8 ) then
1383 0 : wrk_h2so4 = (1000._r8*ch2so4(i,k))/(pures*18._r8)
1384 : else
1385 : wrk_h2so4 = 0._r8
1386 : end if
1387 0 : chno3 = 1.2029e-5_r8 * press(i,k) * tr * hno3_avail(i,k)
1388 0 : do l = 13,22
1389 0 : d(l) = b(1,l) + t*(b(2,l) + t*(b(3,l) + t*(b(4,l) + t*b(5,l))))
1390 : end do
1391 : !----------------------------------------------------------------------
1392 : ! Note that KS depends only on the temperature
1393 : !----------------------------------------------------------------------
1394 0 : t1 = (t - t0)/(t*t0)
1395 0 : t2 = t0/t - 1._r8 - log( t0/t )
1396 0 : lnks = lnks0 - 8792.3984_r8 * t1 - 16.8439_r8 * t2
1397 0 : ks = exp( lnks )
1398 :
1399 0 : converged(i,k) = .false.
1400 : !----------------------------------------------------------------------
1401 : ! Setting up initial guesses for the equations above. Note that
1402 : ! for the initial choices the mass and the charge must be conserved.
1403 : !----------------------------------------------------------------------
1404 0 : delx(0) = .5_r8
1405 0 : z = .5_r8
1406 0 : delz(0) = .5_r8
1407 0 : fratio(0) = 0._r8
1408 0 : reduction_factor = .1_r8
1409 0 : interval_set = .false.
1410 0 : Iter_loop : do iter = 1,itermax
1411 : !----------------------------------------------------------------------
1412 : ! Cwater is the water equation as described in Tabazadeh et
1413 : ! al. and Jacobson et al.
1414 : !----------------------------------------------------------------------
1415 0 : cno3new = chno3 * delx(iter-1)
1416 0 : cphno3new = chno3 * (1._r8 - delx(iter-1))
1417 0 : if( puren > 0._r8 ) then
1418 0 : t1 = (1000._r8*cno3new)/(puren*18._r8)
1419 : else
1420 : t1 = 0._r8
1421 : end if
1422 0 : h2o_cond(i,k) = t1 + wrk_h2so4
1423 0 : if( h2o_cond(i,k) > 0._r8 ) then
1424 0 : wrk = 1.e3_r8 / (18._r8 * h2o_cond(i,k))
1425 0 : molhno3(i,k) = cno3new * wrk
1426 0 : molh2so4(i,k) = ch2so4(i,k) * wrk
1427 : else
1428 0 : molhno3(i,k) = 0._r8
1429 0 : molh2so4(i,k) = 0._r8
1430 : end if
1431 0 : stren = molhno3(i,k) + 3._r8 * molh2so4(i,k)
1432 : !----------------------------------------------------------------------
1433 : ! (1) Calculate the activity of H2SO4 at a given STREN
1434 : !----------------------------------------------------------------------
1435 0 : sm = stren/3._r8
1436 0 : acts = d(13) + sm*(d(14) + sm*(d(15) + sm*d(16)))
1437 : !----------------------------------------------------------------------
1438 : ! (2) Calculate the activity for HNO3 at a given STREN
1439 : !----------------------------------------------------------------------
1440 0 : nm = stren
1441 0 : actn = d(17) + nm*(d(18) + nm*(d(19) + nm*(d(20) + nm*(d(21) + nm*d(22)))))
1442 : !----------------------------------------------------------------------
1443 : ! (3) Calculate the mixed activity coefficient for HNO3 at STREN
1444 : ! as described by Tabazadeh et al.
1445 : !----------------------------------------------------------------------
1446 0 : f1 = 2._r8 * (molh2so4(i,k) + molhno3(i,k)) * actn
1447 0 : f2 = 2.25_r8 * molh2so4(i,k) * acts
1448 :
1449 0 : if (stren > 0._r8) then
1450 0 : mixyln = (f1 + f2) / (2._r8 * stren)
1451 : else
1452 0 : mixyln = 0._r8
1453 : end if
1454 0 : ymix = exp( mixyln )
1455 0 : hplus = 2._r8 * molh2so4(i,k) + molhno3(i,k)
1456 0 : num = ymix**2 * hplus * molhno3(i,k)
1457 0 : den = 1000._r8 * cphno3new * .0820578_r8 * t * ks
1458 0 : if( chno3 == 0._r8 ) then
1459 0 : converged(i,k) = .true.
1460 0 : exit Iter_loop
1461 : end if
1462 : !----------------------------------------------------------------------
1463 : ! the denominator
1464 : ! Calculate the ratio F, check convergence
1465 : !----------------------------------------------------------------------
1466 : !----------------------------------------------------------------------
1467 : ! Calculate the ratio F and reset the deltaX (see Jacobson et al.)
1468 : !----------------------------------------------------------------------
1469 : !!DEK
1470 : ! When the numerator is zero, it can drive the denominator
1471 : ! to 0, which resulted in a NaN for f and also the fraction
1472 : ! ratio. Assume that in this case, the limit of f would
1473 : ! really approach 1, not infinity and thus converge the
1474 : ! solution.
1475 0 : if ((num .eq. 0._r8) .and. (den .eq. 0._r8)) then
1476 0 : f = 1._r8
1477 : else
1478 0 : f = num / den
1479 : end if
1480 0 : fratio(iter) = abs( f ) - 1._r8
1481 0 : con_val = abs( f - 1._r8 )
1482 0 : if( con_val <= con_lim ) then
1483 0 : converged(i,k) = .true.
1484 0 : exit Iter_loop
1485 : end if
1486 : !----------------------------------------------------------------------
1487 : ! non-convergence; setup next iterate
1488 : !----------------------------------------------------------------------
1489 0 : if( interval_set ) then
1490 0 : z = reduction_factor * z
1491 0 : delz(iter) = z
1492 0 : if( f > 1._r8 ) then
1493 0 : deltax = -z
1494 : else
1495 : deltax = z
1496 : end if
1497 0 : delx(iter) = delx(iter-1) + deltax
1498 : else
1499 0 : if( iter == 1 ) then
1500 0 : if( fratio(iter) >= 1._r8 ) then
1501 : positive = .false.
1502 : else
1503 0 : positive = .true.
1504 : end if
1505 : end if
1506 0 : if( fratio(iter)*fratio(iter-1) < 0._r8 ) then
1507 0 : interval_set = .true.
1508 0 : reduction_factor = .5_r8
1509 0 : delx(iter) = .5_r8*(delx(iter-1) + delx(iter-2))
1510 0 : z = .5_r8*abs( delx(iter-1) - delx(iter-2) )
1511 : else
1512 0 : if( .not. positive ) then
1513 0 : delx(iter) = reduction_factor * delx(iter-1)
1514 : else
1515 0 : delx(iter) = reduction_factor + delx(iter-1)
1516 0 : if( delx(iter) > upper_delx ) then
1517 0 : delx(iter) = .5_r8
1518 0 : interval_set = .true.
1519 0 : reduction_factor = .5_r8
1520 : end if
1521 : end if
1522 : end if
1523 : end if
1524 : end do Iter_loop
1525 :
1526 0 : wtotal = molhno3(i,k) * 63._r8 + molh2so4(i,k) * 98._r8 + 1000._r8
1527 0 : wts(i,k) = (molh2so4(i,k) * 9800._r8) / wtotal
1528 0 : wtn(i,k) = (molhno3(i,k) *6300._r8)/ wtotal
1529 0 : if( cno3new /= 0._r8 .or. cphno3new /= 0._r8 ) then
1530 0 : ratio = max( 0._r8,min( 1._r8,cno3new/(cphno3new + cno3new) ) )
1531 0 : hno3_cond(i,k) = ratio*hno3_avail(i,k)
1532 : else
1533 0 : hno3_cond(i,k) = 0._r8
1534 : end if
1535 0 : if( .not. converged(i,k) ) then
1536 0 : write(iulog,*) 'equil: Failed to converge @ is_chem,flag,lchnk,i,k,f = ',is_chem,flag,lchnk,i,k,f
1537 0 : write(iulog,*) ' wts0,pures,puren,chno3,ch2so4 = ',wts0(i,k),pures,puren,chno3,ch2so4(i,k)
1538 0 : write(iulog,*) ' stren,mixyln,ymix,hplus,num,den = ',stren,mixyln,ymix,hplus,num,den
1539 0 : write(iulog,*) ' h2o_avail,hno3_avail,p,t = ',h2o_avail(i,k),hno3_avail(i,k),press(i,k),temper(i,k)
1540 0 : write(iulog,*) ' molhno3,molh2so4,h2o_cond,hno3_cond = ', &
1541 0 : molhno3(i,k),molh2so4(i,k),h2o_cond(i,k),hno3_cond(i,k)
1542 0 : if( con_val > .05_r8 ) then
1543 0 : write(iulog,*) ' '
1544 0 : write(iulog,*) 'equil; diagnostics at lchnk, flag, i, k, iter = ',lchnk,flag,i,k,iter
1545 0 : write(iulog,*) 'equil; fratio'
1546 0 : write(iulog,'(5(1pg15.7))') fratio(0:iter-1)
1547 0 : write(iulog,*) ' '
1548 0 : write(iulog,*) 'equil; delx'
1549 0 : write(iulog,'(5(1pg15.7))') delx(0:iter-1)
1550 0 : write(iulog,*) ' '
1551 0 : write(iulog,*) 'equil; delz'
1552 0 : write(iulog,'(5(1pg15.7))') delz(0:iter-1)
1553 0 : write(iulog,*) ' '
1554 0 : else if( iter > 50 ) then
1555 0 : write(iulog,*) 'equil: Iterations are beyond 50, number of iter = ',iter
1556 0 : write(iulog,*) 'equil: converged @ is_chem,flag,lchnk,i,k = '
1557 0 : write(iulog,*) is_chem,flag,lchnk,i,k
1558 0 : write(iulog,*) 'equil: converged @ f, num, den = '
1559 0 : write(iulog,*) f, num, den
1560 0 : write(iulog,*) ' h2o_avail,hno3_avail,p,t = '
1561 0 : write(iulog,*) h2o_avail(i,k),hno3_avail(i,k),press(i,k),temper(i,k)
1562 0 : write(iulog,*) ' molhno3(i,k),molh2so4(i,k),h2o_cond,hno3_cond = '
1563 0 : write(iulog,*) molhno3(i,k),molh2so4(i,k),h2o_cond(i,k),hno3_cond(i,k)
1564 : end if
1565 : end if
1566 : end if
1567 : end do Column_loop
1568 : end do Level_loop
1569 :
1570 0 : end subroutine equil
1571 :
1572 : !======================================================================
1573 : !
1574 : !
1575 : ! ROUTINE
1576 : ! DENSITY
1577 : !
1578 : ! Date...
1579 : ! 7 October 1999
1580 : !
1581 : ! Programmed by...
1582 : ! A. Tabazadeh
1583 : !
1584 : ! DESCRIPTION
1585 : ! Calculates the density (g cm-3) of a binary sulfate solution.
1586 : !
1587 : ! ARGUMENTS
1588 : ! INPUT
1589 : ! T Temperature
1590 : ! w Weight fraction
1591 : !
1592 : ! OUTPUT
1593 : ! den Density of the Binary Solution (g cm-3)
1594 : !
1595 : !======================================================================
1596 :
1597 0 : function density( temp, w )
1598 :
1599 : implicit none
1600 :
1601 : !----------------------------------------------------------------------
1602 : ! ... Dummy arguments
1603 : !----------------------------------------------------------------------
1604 : real(r8), intent(in) :: temp, w
1605 :
1606 : !----------------------------------------------------------------------
1607 : ! ... Function declarations
1608 : !----------------------------------------------------------------------
1609 : real(r8) :: density
1610 :
1611 : !----------------------------------------------------------------------
1612 : ! ... Local variables
1613 : !----------------------------------------------------------------------
1614 : real(r8), parameter :: a9 = -268.2616e4_r8, a10 = 576.4288e3_r8
1615 :
1616 : real(r8) :: a0, a1, a2, a3, a4, a5, a6, a7 ,a8
1617 : real(r8) :: c1, c2, c3, c4
1618 :
1619 : !----------------------------------------------------------------------
1620 : ! ... Temperature variables
1621 : !----------------------------------------------------------------------
1622 0 : c1 = temp - 273.15_r8
1623 0 : c2 = c1**2
1624 0 : c3 = c1*c2
1625 0 : c4 = c1*c3
1626 : !----------------------------------------------------------------------
1627 : ! Polynomial Coefficients
1628 : !----------------------------------------------------------------------
1629 0 : a0 = 999.8426_r8 + 334.5402e-4_r8*c1 - 569.1304e-5_r8*c2
1630 0 : a1 = 547.2659_r8 - 530.0445e-2_r8*c1 + 118.7671e-4_r8*c2 + 599.0008e-6_r8*c3
1631 0 : a2 = 526.295e1_r8 + 372.0445e-1_r8*c1 + 120.1909e-3_r8*c2 - 414.8594e-5_r8*c3 + 119.7973e-7_r8*c4
1632 0 : a3 = -621.3958e2_r8 - 287.7670_r8*c1 - 406.4638e-3_r8*c2 + 111.9488e-4_r8*c3 + 360.7768e-7_r8*c4
1633 0 : a4 = 409.0293e3_r8 + 127.0854e1_r8*c1 + 326.9710e-3_r8*c2 - 137.7435e-4_r8*c3 - 263.3585e-7_r8*c4
1634 0 : a5 = -159.6989e4_r8 - 306.2836e1_r8*c1 + 136.6499e-3_r8*c2 + 637.3031e-5_r8*c3
1635 0 : a6 = 385.7411e4_r8 + 408.3717e1_r8*c1 - 192.7785e-3_r8*c2
1636 0 : a7 = -580.8064e4_r8 - 284.4401e1_r8*c1
1637 0 : a8 = 530.1976e4_r8 + 809.1053_r8*c1
1638 : !----------------------------------------------------------------------
1639 : ! ... Summation
1640 : !----------------------------------------------------------------------
1641 0 : density = .001_r8*(a0 + w*(a1 + w*(a2 + w*(a3 + w*(a4 + w*(a5 + w*(a6 + w*(a7 + w*(a8 + w*(a9 + w*a10))))))))))
1642 :
1643 0 : end function density
1644 :
1645 : end module mo_sad
|
