Line data Source code
1 : ! Include shortname defintions, so that the F77 code does not have to be modified to
2 : ! reference the CARMA structure.
3 : #include "carma_globaer.h"
4 :
5 : !! This routine evaluates particle loss rates due to nucleation <rnuclg>:
6 : !! heterogeneous deposition nucleation only. The parameters are adjusted
7 : !! for mesospheric conditions, based upon the recommendations of Keesee.
8 : !!
9 : !! Based on expressions from ...
10 : !! Keesee [JGR,1989]
11 : !! Pruppacher and Klett [2000]
12 : !! Rapp and Thomas [JASTP, 2006]
13 : !! Trainer et al. [2008]
14 : !!
15 : !! The loss rates for all particle elements in a particle group are equal.
16 : !!
17 : !! To avoid nucleation into an evaporating bin, this subroutine must
18 : !! be called after growp, which evaluates evaporation loss rates <evaplg>.
19 : !!
20 : !! @author Eric Jensen, Chuck Bardeen
21 : !! @version Oct-2000, Jan-2010
22 1418703438 : subroutine hetnucl(carma, cstate, iz, rc)
23 :
24 : ! types
25 : use carma_precision_mod
26 : use carma_enums_mod
27 : use carma_constants_mod
28 : use carma_types_mod
29 : use carmastate_mod
30 : use carma_mod
31 :
32 : implicit none
33 :
34 : type(carma_type), intent(in) :: carma !! the carma object
35 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
36 : integer, intent(in) :: iz !! z index
37 : integer, intent(inout) :: rc !! return code, negative indicates failure
38 :
39 : ! Local declarations
40 : integer :: igas ! gas index
41 : integer :: igroup ! group index
42 : integer :: ibin ! bin index
43 : integer :: iepart ! element for condensing group index
44 : integer :: inuc ! nucleating element index
45 : integer :: ienucto ! index of target nucleation element
46 : integer :: ignucto ! index of target nucleation group
47 : real(kind=f) :: rmw
48 : real(kind=f) :: R_H2O
49 : real(kind=f) :: rnh2o
50 : real(kind=f) :: rlogs
51 : real(kind=f) :: ag
52 : real(kind=f) :: contang
53 : real(kind=f) :: xh
54 : real(kind=f) :: phih
55 : real(kind=f) :: rath
56 : real(kind=f) :: fv3h
57 : real(kind=f) :: fv4h
58 : real(kind=f) :: fh
59 : real(kind=f) :: delfg
60 : real(kind=f) :: expon
61 :
62 : ! Heterogeneous nucleation factors
63 : real(kind=f), parameter :: gdes = 2.9e-13_f
64 : real(kind=f), parameter :: gsd = 2.9e-14_f
65 : real(kind=f), parameter :: zeld = 0.1_f
66 : real(kind=f), parameter :: vibfreq = 1.e13_f
67 : real(kind=f), parameter :: diflen = 0.1e-7_f
68 : real(kind=f) :: rmiv
69 :
70 1418703438 : rmiv = 0.95_f
71 :
72 : ! rmiv - Eq. 2, Trainer et al. [2008]
73 : ! rmiv = 0.94_f - (6005._f * exp(-0.065_f * max(150._f, t(iz))))
74 : ! rmiv = max(0._f, 0.94_f - (6005._f * exp(-0.065_f * t(iz))))
75 :
76 : ! Loop over particle groups.
77 4256110314 : do igroup = 1, NGROUP
78 :
79 2837406876 : igas = inucgas(igroup) ! condensing gas
80 :
81 4256110314 : if (igas .ne. 0) then
82 :
83 1418703438 : iepart = ienconc(igroup) ! particle number density element
84 :
85 1418703438 : rmw = gwtmol(igas) / AVG
86 1418703438 : R_H2O = RGAS / gwtmol(igas)
87 1418703438 : rnh2o = gc(iz,igas) * R_H2O / BK
88 :
89 : ! Calculate nucleation loss rates. Do not allow nucleation into
90 : ! an evaporating bin.
91 : !
92 : ! <ienucto> is index of target nucleation element;
93 : ! <ignucto> is index of target nucleation group.
94 2837406876 : do inuc = 1, nnuc2elem(iepart)
95 :
96 1418703438 : ienucto = inuc2elem(inuc,iepart)
97 :
98 1418703438 : if (ienucto .ne. 0) then
99 1418703438 : ignucto = igelem(ienucto)
100 : else
101 : ignucto = 0
102 : endif
103 :
104 : ! Only compute nucleation rate for heterogenous nucleation
105 2837406876 : if (inucproc(iepart,ienucto) .eq. I_HETNUC) then
106 :
107 : ! Loop over particle bins. Loop from largest to smallest for
108 : ! evaluation of index of smallest bin nucleated during time step <inucstep>.
109 0 : do ibin = NBIN, 1, -1
110 :
111 : ! Bypass calculation if few particles are present
112 0 : if (pconmax(iz,igroup) .gt. FEW_PC) then
113 :
114 : ! Only proceed if ice supersaturated
115 : !
116 : ! NOTE: We are only trying to model PMC partcles, so turn of nucleation
117 : ! where the CAM microphysics takes over (~1 mb = 1000 dyne).
118 0 : if ((p(iz) .lt. 1.e3_f) .and. (supsati(iz,igas) .gt. 0._f)) then
119 0 : rlogs = log(supsati(iz,igas) + 1._f)
120 :
121 : ! Critical ice germ radius formed in the sulfate solution
122 : !
123 : ! Eq. 2, Rapp & Thomas [2006]
124 0 : ag = 2._f * gwtmol(igas) * surfctia(iz) / rgas / t(iz) / RHO_I / rlogs
125 :
126 : ! Heterogeneous nucleation geometric factor
127 : !
128 : ! Eq. 9-22, Pruppacher & Klett [2000]
129 0 : contang = acos(rmiv)
130 0 : xh = r(ibin,igroup) / ag
131 0 : phih = sqrt(1._f - 2._f * rmiv * xh + xh**2 )
132 0 : rath = (xh-rmiv) / phih
133 0 : fv3h = xh**3 * (2._f - 3._f * rath + rath**3 )
134 0 : fv4h = 3._f * rmiv * xh**2 * (rath - 1._f)
135 :
136 0 : if (abs(rath) .gt. 1._f - 1.e-8_f) fv3h = 0._f
137 0 : if (abs(rath) .gt. 1._f - 1.e-10_f) fv4h = 0._f
138 :
139 0 : fh = 0.5_f * (1._f + ((1._f - rmiv * xh) / phih)**3 + fv3h + fv4h)
140 :
141 : ! Gibbs free energy of ice germ formation in the ice/sulfate solution
142 : !
143 : ! Eq. 3, Rapp & Thomas [2006]
144 0 : delfg = 4._f * PI * ag**2 * surfctia(iz) - 4._f * PI * RHO_I * ag**3 *BK * t(iz) * rlogs / 3._f / rmw
145 :
146 : ! Ice nucleation rate in a 0.2 micron aerosol (/sec)
147 0 : expon = (2._f * gdes - gsd - fh*delfg) / BK / t(iz)
148 :
149 : ! NOTE: Excessive nucleation makes it difficult for the substepping to find a
150 : ! stable solution, so put a cap on really large nucleation values that can be produced.
151 0 : rnuclg(ibin,igroup,ignucto) = min(1e10_f, zeld * BK * t(iz) * diflen * ag * sin(contang) * &
152 0 : 4._f * PI * r(ibin,igroup)**2 * rnh2o**2 / (fh * rmw * vibfreq) * exp(expon))
153 : endif
154 : endif ! pconmax(ixyz,igroup) .gt. FEW_PC
155 : enddo ! ibin = 1,NBIN
156 : endif ! inucproc(iepart,ienucto) .eq. I_DROPACT
157 : enddo ! inuc = 1,nnuc2elem(iepart)
158 : endif ! (igas = inucgas(igroup) .ne. 0)
159 : enddo ! igroup = 1,NGROUP
160 :
161 : ! Return to caller with particle loss rates due to nucleation evaluated.
162 1418703438 : return
163 1418703438 : end
|
