Line data Source code
1 :
2 : !======================================================================
3 : !
4 : ! ROUTINE
5 : ! strat_aer_settling
6 : !
7 : ! Date...
8 : ! 8 November 1999
9 : !
10 : ! Programmed by...
11 : ! Douglas E. Kinnison
12 : !
13 : ! Modified for WACCM2 on 3 September 2004 - Removed ICE
14 : ! Modified for WACCM3 on 8 November 2004 - added NAT 7um mode.
15 : !
16 : ! DESCRIPTION
17 : ! This routine vertically redistributes condensed phase HNO3.
18 : !
19 : ! For each aerosol type the terminal velocity is calculated. This
20 : ! quantity is dependent on: 1) the mass density of the aerosol;
21 : ! 2) the radius; 3) the dynamic viscosity; 4) shape; and the Cunningham
22 : ! correction factor for spherical particles. See Fuchs, The Mechanics
23 : ! of Aerosols Oxford, Pergmann Press, pp 27-31, 1964 and Kasten,
24 : ! Falling Speed of Aerosol Particles, J. of Appl. Met., 7, 944-947,
25 : ! 1968 for details. For aerosol with a radius of 3 microns (e.g., NAT)
26 : ! and 10 microns (e.g., ICE) the terminal velocity (cm sec-1) is 0.2
27 : ! and 1.7 respectively.
28 : !
29 : ! The flux of condensed phase HNO3 is then derived using the
30 : ! following equation:
31 : ! Flux (molec cm-2 s-1) = V * C * exp (8*ln^2 sigma).
32 : !
33 : ! where: V is terminal velocity (cm sec-1)
34 : ! C is condensed phase conc (molecules cm-3)
35 : ! sigma is the width of the log normal distribution.
36 : !
37 : ! The approach of settling the entire aerosol size distribution is
38 : ! based on work of Considine et al., JGR, 1999.
39 : !
40 : ! The routine is a straighforward approach, starting at the top zone
41 : ! (highest altitude) and progressing towards the surface. The gross
42 : ! settling of condensed phase HNO3 is simply the quanity:
43 : ! Flux * dt/dz.
44 : !
45 : !
46 : ! ARGUMENTS
47 : !
48 : ! All of the following components are grid center quanitities.
49 : !
50 : ! INPUT:
51 : ! ad Air Density (molecules cm-3)
52 : ! press Pressure (Pascals)
53 : ! timestep Gross chemistry timestep (in seconds)
54 : ! temp Temperature (K)
55 : ! hno3_cond Condensed phase HNO3 (mole fraction)
56 : ! radius_nat Mean radius of NAT (cm)
57 : ! zstar log pressure altitude coordinate (km)
58 : !
59 : ! OUTPUT:
60 : ! hno3_cond vertically modified
61 : !
62 : !======================================================================
63 :
64 : module mo_aero_settling
65 :
66 : use shr_kind_mod, only : r8 => shr_kind_r8
67 :
68 : private
69 : public :: strat_aer_settling
70 : public :: strat_aer_settl_init
71 :
72 : contains
73 :
74 0 : subroutine strat_aer_settl_init
75 :
76 : use cam_history, only : addfld
77 :
78 : implicit none
79 :
80 0 : call addfld( 'VEL_NAT1', (/ 'lev' /), 'I', 'cm/s', 'small nat settling velocity' )
81 0 : call addfld( 'VEL_NAT2', (/ 'lev' /), 'I', 'cm/s', 'large nat settling velocity' )
82 :
83 0 : end subroutine strat_aer_settl_init
84 :
85 :
86 0 : subroutine strat_aer_settling( ad, press, timestep, zstar, temp, &
87 0 : hno3_cond, radius_nat, ncol, lchnk, aero_ndx )
88 :
89 0 : use ppgrid, only : pcols, pver
90 : use physconst, only : gravit
91 : use cam_history, only : outfld
92 :
93 : !-----------------------------------------------------------------------
94 : ! ... dummy arguments
95 : !-----------------------------------------------------------------------
96 : integer, intent(in) :: ncol ! columns in chunk
97 : integer, intent(in) :: lchnk ! chunk id
98 : integer, intent(in) :: aero_ndx ! aerosol index
99 : real(r8), intent(in) :: timestep ! model time step (s)
100 : real(r8), intent(in) :: ad(ncol,pver) ! Air density (molecules cm-3)
101 : real(r8), intent(in) :: radius_nat(ncol,pver) ! Mean radius of NAT (cm)
102 : real(r8), intent(in) :: zstar(ncol,pver) ! altitude (km)
103 : real(r8), intent(in) :: press(pcols,pver) ! Pressure (Pa)
104 : real(r8), intent(in) :: temp(pcols,pver) ! temperature (K)
105 : real(r8), intent(inout) :: hno3_cond(ncol,pver) ! Condensed Phase HNO3 (VMR)
106 :
107 : !-----------------------------------------------------------------------
108 : ! ... local variables
109 : !-----------------------------------------------------------------------
110 : real(r8), parameter :: avo_num = 6.022e23_r8, & ! molecules/mole
111 : MW_air = 28.8_r8, & ! grams/mole air
112 : nat_dens = 1.6_r8, & ! g/cm^3
113 : shape_fac_nat = 1.0_r8, & ! TBD
114 : sigma_nat = 1.6_r8, & ! Width of distribution
115 : av_const = 2.117265e4_r8, & ! (8*8.31448*1000 / PI)
116 : km2cm = 1.e5_r8, & ! km to cm
117 : m2cm = 1.e2_r8, & ! m to cm
118 : c1 = 2._r8/9._r8
119 :
120 : integer :: i, k, kp1
121 : real(r8) :: gravity ! gravity cm/s^2
122 : real(r8) :: Cc_nat ! Cunningham Correction Factor
123 : real(r8) :: dt_dz ! dt / dz, sec cm-1
124 : real(r8) :: flux_nat ! aerosol flux, molec cm-2 sec-1
125 : real(r8) :: mean_vel ! mean velocity, cm sec-1
126 : real(r8) :: mfp ! Mean Free Path
127 : real(r8) :: vel_nat ! Terminal velocity (cm sec-1)
128 : real(r8) :: rad_nat ! wrk radius nat (cm)
129 : real(r8) :: visc ! Dynamic viscosity of air
130 : real(r8) :: depos ! molecules/cm**3 deposited
131 : real(r8) :: atm_dens, atm_densa ! total atm density and inverse
132 : real(r8) :: cond_hno3 ! wrk variables
133 : real(r8) :: const_nat ! wrk variables
134 : real(r8) :: t ! working temperatue
135 0 : real(r8) :: velnat(ncol,pver) ! holding variable for output
136 0 : logical :: lon_mask(ncol) ! longitude logic mask
137 :
138 0 : gravity = gravit*m2cm ! (cm/s^2)
139 :
140 0 : const_nat = exp( 8._r8*(log(sigma_nat))**2 )
141 0 : do k = 1,pver
142 0 : velnat(:,k) = 0._r8
143 : end do
144 :
145 :
146 : !----------------------------------------------------------------------
147 : ! ... derive Aerosol Settling (explicit approach)
148 : !----------------------------------------------------------------------
149 : Level_loop : &
150 0 : do k = 2,pver-1
151 : !----------------------------------------------------------------------
152 : ! ... operate between 2.0hPa and 300hPa and only where nat exist
153 : !----------------------------------------------------------------------
154 0 : kp1 = k + 1
155 0 : lon_mask(:) = press(:ncol,k) >= 2.e2_r8 .and. press(:ncol,k) <= 300.e2_r8 &
156 0 : .and. (radius_nat(:,k) > 0._r8 )
157 0 : if( any( lon_mask(:) ) ) then
158 : Column_loop : &
159 0 : do i = 1,ncol
160 0 : if( lon_mask(i) ) then
161 0 : t = temp(i,k)
162 : !----------------------------------------------------------------------
163 : ! ... General Setup for NAT
164 : ! Calculate the settling of the NAT Aerosol
165 : ! NOTE: Index "k" is the box that is being calculated
166 : ! Index "k-1" is flux from above into the k box
167 : ! A positive NET {flux*dt/dz} adds to the "k" box
168 : !----------------------------------------------------------------------
169 : ! ... mean Molecular Velocity, cm sec-1
170 : !----------------------------------------------------------------------
171 0 : mean_vel = sqrt( av_const*t/MW_air )*100._r8
172 : !----------------------------------------------------------------------
173 : ! ... dynamic Viscosity, g cm-1 sec-1
174 : !----------------------------------------------------------------------
175 0 : visc = (t*1.458e-6_r8)**1.5_r8 /(t + 110.4_r8)*10000._r8
176 : !----------------------------------------------------------------------
177 : ! ... mean Free Path, cm
178 : !----------------------------------------------------------------------
179 0 : atm_dens = ad(i,k)
180 0 : atm_densa = 1._r8/atm_dens
181 0 : mfp = visc* avo_num / (.499_r8*atm_dens*MW_air*mean_vel)
182 : !----------------------------------------------------------------------
183 : ! ... dt / dz, sec cm-1; NOTE: zstar is in km
184 : !----------------------------------------------------------------------
185 0 : dt_dz = timestep / ((zstar(i,k-1) - zstar(i,k))*km2cm)
186 : !----------------------------------------------------------------------
187 : ! ... calculate NAT Aerosol Settling
188 : !----------------------------------------------------------------------
189 0 : rad_nat = radius_nat(i,k)
190 0 : cond_hno3 = hno3_cond(i,k)*atm_dens
191 : !----------------------------------------------------------------------
192 : ! ... Cunningham Correction Factor, Unitless
193 : !----------------------------------------------------------------------
194 0 : Cc_nat = 1._r8 + (mfp/rad_nat)*(1.246_r8 + .42_r8*exp( -.87_r8*rad_nat/mfp ))
195 : !----------------------------------------------------------------------
196 : ! ... terminal Velocity of Aerosol, cm sec-1
197 : !----------------------------------------------------------------------
198 0 : vel_nat = c1*rad_nat**2 * nat_dens*gravity*Cc_nat/visc*shape_fac_nat
199 0 : velnat(i,k) = vel_nat
200 : !----------------------------------------------------------------------
201 : ! ... aerosol Flux, Cond-phase molecules cm-2 sec-1
202 : !----------------------------------------------------------------------
203 0 : flux_nat = cond_hno3*vel_nat* const_nat
204 : !----------------------------------------------------------------------
205 : ! ... calculate NAT Aerosol Settling (i.e., HNO3 redistribution)
206 : !----------------------------------------------------------------------
207 0 : depos = min( cond_hno3,flux_nat*dt_dz )
208 : !----------------------------------------------------------------------
209 : ! ... modify the HNO3_cond in level "k"
210 : !----------------------------------------------------------------------
211 0 : hno3_cond(i,k) = (cond_hno3 - depos)*atm_densa
212 : !----------------------------------------------------------------------
213 : ! ... modify the HNO3_cond in level "k+1"
214 : !----------------------------------------------------------------------
215 0 : hno3_cond(i,kp1) = hno3_cond(i,kp1) + depos/ad(i,kp1)
216 : end if
217 : end do Column_loop
218 : end if
219 : end do Level_loop
220 :
221 : !----------------------------------------------------------------------
222 : ! ... output nat velocity
223 : !----------------------------------------------------------------------
224 0 : if( aero_ndx == 1 ) then
225 0 : call outfld( 'VEL_NAT1', velnat, ncol, lchnk )
226 0 : else if( aero_ndx == 2 ) then
227 0 : call outfld( 'VEL_NAT2', velnat, ncol, lchnk )
228 : end if
229 :
230 0 : end subroutine strat_aer_settling
231 :
232 : end module mo_aero_settling
|
