Line data Source code
1 : #include "carma_globaer.h"
2 :
3 : !! This routine calculates coagulation production terms <coagpe>.
4 : !! See [Jacobson, et al., Atmos. Env., 28, 1327, 1994] for details
5 : !! on the coagulation algorithm.
6 : !!
7 : !! @author Eric Jensen
8 : !! @version Oct-1995
9 147087360 : subroutine coagp(carma, cstate, ibin, ielem, rc)
10 :
11 : ! types
12 : use carma_precision_mod
13 : use carma_enums_mod
14 : use carma_constants_mod
15 : use carma_types_mod
16 : use carmastate_mod
17 : use carma_mod
18 :
19 : implicit none
20 :
21 : type(carma_type), intent(in) :: carma !! the carma object
22 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
23 : integer, intent(in) :: ibin !! bin index
24 : integer, intent(in) :: ielem !! element index
25 : integer, intent(inout) :: rc !! return code, negative indicates failure
26 :
27 : ! Local Variables
28 : integer :: iz
29 : integer :: igrp
30 : integer :: i_pkern
31 : integer :: iquad
32 : integer :: ig
33 : integer :: jg
34 : integer :: i
35 : integer :: j
36 : integer :: iefrom
37 : integer :: iefrom_cm
38 : integer :: je
39 : integer :: je_cm
40 : integer :: ic
41 : integer :: iecore
42 : real(kind=f) :: totmass
43 : real(kind=f) :: rmasscore
44 : real(kind=f) :: fracmass
45 : real(kind=f) :: elemass
46 : real(kind=f) :: rmi
47 : real(kind=f) :: rmj
48 :
49 :
50 : ! Definition of i,j,k,n used in comments: colision between i and j bins
51 : ! yields particle between bins k and k+1. Production in bin n is calculated.
52 :
53 :
54 : ! Determine group that particles are produced in
55 147087360 : igrp = igelem(ielem)
56 :
57 : ! Particle number production
58 : !
59 : ! Coagulation between particle in group <ig> bin <i> with particle in
60 : ! group <jg> bin <j> results in particle with mass between bins k and k+1.
61 : ! First, loop over group-bin quads <ig,i,jg,j> resulting in production in
62 : ! bin <ibin> = k+1. The set of quads <igup,jgup,iup,jup> is
63 : ! defined in setupcoag.
64 :
65 7592859648 : do iquad = 1, npairu(igrp,ibin)
66 :
67 7445772288 : ig = igup(igrp,ibin,iquad) ! source group
68 7445772288 : jg = jgup(igrp,ibin,iquad) ! source group
69 7445772288 : i = iup(igrp,ibin,iquad) ! source bin
70 7445772288 : j = jup(igrp,ibin,iquad) ! source bin
71 :
72 7445772288 : iefrom = icoagelem(ielem,ig) ! source element for <i> particle
73 :
74 7445772288 : if( if_sec_mom(igrp) )then
75 0 : iefrom_cm = icoagelem_cm(ielem,ig) ! core mass moment source element
76 : endif
77 :
78 : ! If <iefrom> = 0 then there is no contribution to production
79 7592859648 : if( iefrom .ne. 0 ) then
80 :
81 5449586688 : je = ienconc(jg) ! source element for <j> particle
82 :
83 5449586688 : if( if_sec_mom(igrp) )then
84 0 : je_cm = icoagelem_cm(ielem,jg) ! core mass moment source element
85 : endif
86 :
87 : ! If ielem is core mass type and <ig> is a CN type and <ig> is different
88 : ! from <igrp>, then we must multiply production by mass
89 : ! per particle (<rmi>) of element <icoagelem>. (this is <rmass> for all source
90 : ! elements except particle number concentration in a multicomponent CN group).
91 >38692*10^7 : do iz = 1, NZ
92 :
93 : ! Bypass calculation if few source particles present
94 >38147*10^7 : if( pconmax(iz,ig) .gt. FEW_PC .and. &
95 5449586688 : pconmax(iz,jg) .gt. FEW_PC )then
96 :
97 >10029*10^7 : rmi = 1._f
98 >10029*10^7 : i_pkern = 1
99 :
100 >10029*10^7 : if( itype(ielem) .eq. I_COREMASS .or. &
101 : itype(ielem) .eq. I_VOLCORE )then ! core mass
102 :
103 72892956410 : i_pkern = 3 ! Use different kernel for core mass prod.
104 :
105 >14578*10^7 : if( ( itype(ienconc(ig)) .eq. I_INVOLATILE .or. &
106 : itype(ienconc(ig)) .eq. I_VOLATILE ) &
107 >14578*10^7 : .and. ig .ne. igrp ) then
108 :
109 : ! CN source and ig different from igrp
110 :
111 0 : if( ncore(ig) .eq. 0 )then ! No cores in source group
112 :
113 0 : if(icomp(ienconc(ig)) .eq. icomp(ielem)) then
114 0 : rmi = rmass(i,ig)
115 : else
116 : rmi = 0._f
117 : endif
118 :
119 0 : elseif( itype(iefrom) .eq. I_INVOLATILE .or. &
120 : itype(iefrom) .eq. I_VOLATILE ) then
121 :
122 : ! Source element is number concentration elem of mixed CN group
123 0 : totmass = pc(iz,i,iefrom) * rmass(i,ig)
124 0 : rmasscore = pc(iz,i,icorelem(1,ig))
125 :
126 0 : do ic = 2,ncore(ig)
127 0 : iecore = icorelem(ic,ig)
128 0 : rmasscore = rmasscore + pc(iz,i,iecore)
129 : enddo
130 :
131 0 : fracmass = 1._f - rmasscore/totmass
132 0 : elemass = fracmass * rmass(i,ig)
133 0 : rmi = elemass
134 : endif
135 : endif ! ig is a CCN and not igrp
136 :
137 27403646724 : elseif( itype(ielem) .eq. I_CORE2MOM )then ! core mass^2
138 :
139 0 : i_pkern = 5 ! Use different kernel for core mass^2 production
140 0 : rmj = 1._f
141 :
142 0 : if( itype(ienconc(ig)) .eq. I_INVOLATILE ) then
143 0 : rmi = rmass(i,ig)
144 0 : rmj = rmass(j,jg)
145 : endif
146 :
147 : endif ! itype(ielem) is a coremass or core2mom
148 :
149 : ! For each spatial grid point, sum up coagulation production
150 : ! contributions from each quad.
151 >10029*10^7 : if( itype(ielem) .ne. I_CORE2MOM )then
152 0 : coagpe(iz,ibin,ielem) = coagpe(iz,ibin,ielem) + &
153 0 : pc(iz,i,iefrom)*pcl(iz,j,je)*rmi * &
154 0 : ckernel(iz,i,j,ig,jg) * &
155 >10029*10^7 : pkernel(i,j,ig,jg,igrp,i_pkern)
156 : else
157 0 : coagpe(iz,ibin,ielem) = coagpe(iz,ibin,ielem) + &
158 0 : ( pc(iz,i,iefrom)*pcl(iz,j,je)*rmi**2 + &
159 0 : pc(iz,i,iefrom_cm)*rmi* &
160 0 : pcl(iz,j,je_cm)*rmj ) * &
161 0 : ckernel(iz,i,j,ig,jg) * &
162 0 : pkernel(i,j,ig,jg,igrp,i_pkern)
163 : endif
164 : endif ! end of ( pconmax .gt. FEW_PC )
165 : enddo ! iz = 1, NZ
166 : endif ! iefrom .ne. 0
167 : enddo ! iquad
168 :
169 : ! Next, loop over group-bin quads for production in bin <ibin> = k from
170 : ! bin <i> due to collision between bins <i> and <j>.
171 : ! Production will only occur if either k != <i> or igrp != <ig>
172 4400013312 : do iquad = 1, npairl(igrp,ibin)
173 :
174 4252925952 : ig = iglow(igrp,ibin,iquad)
175 4252925952 : jg = jglow(igrp,ibin,iquad)
176 4252925952 : i = ilow(igrp,ibin,iquad)
177 4252925952 : j = jlow(igrp,ibin,iquad)
178 :
179 4252925952 : iefrom = icoagelem(ielem,ig) ! source element for <i> particle
180 :
181 4252925952 : if( if_sec_mom(igrp) )then
182 0 : iefrom_cm = icoagelem_cm(ielem,ig) ! core mass moment source element
183 : endif
184 :
185 4400013312 : if( iefrom .ne. 0 ) then
186 :
187 2151677952 : je = ienconc(jg) ! source element for <j> particle
188 :
189 2151677952 : if( if_sec_mom(igrp) )then
190 0 : je_cm = icoagelem_cm(ielem,jg) ! core mass moment source element
191 : endif
192 :
193 >15276*10^7 : do iz = 1, NZ
194 :
195 : ! Bypass calculation if few particles present
196 >15061*10^7 : if( pconmax(iz,ig) .gt. FEW_PC .and. &
197 2151677952 : pconmax(iz,jg) .gt. FEW_PC )then
198 :
199 43199642230 : rmi = 1._f
200 43199642230 : i_pkern = 2
201 :
202 43199642230 : if( itype(ielem) .eq. I_COREMASS .or. &
203 : itype(ielem) .eq. I_VOLCORE )then ! core mass
204 :
205 28031558375 : i_pkern = 4 ! Use different kernel for core mass production
206 :
207 56063116750 : if( ( itype(ienconc(ig)) .eq. I_INVOLATILE .or. &
208 : itype(ienconc(ig)) .eq. I_VOLATILE ) &
209 56063116750 : .and. ig .ne. igrp ) then
210 :
211 : ! CN source and ig different from igrp
212 :
213 0 : if( ncore(ig) .eq. 0 )then ! No cores in source group
214 0 : rmi = rmass(i,ig)
215 :
216 0 : elseif( itype(iefrom) .eq. I_INVOLATILE .or. &
217 : itype(iefrom) .eq. I_VOLATILE ) then
218 :
219 : ! Source element is number concentration elem of mixed CN group
220 :
221 0 : totmass = pc(iz,i,iefrom) * rmass(i,ig)
222 0 : rmasscore = pc(iz,i,icorelem(1,ig))
223 0 : do ic = 2,ncore(ig)
224 0 : iecore = icorelem(ic,ig)
225 0 : rmasscore = rmasscore + pc(iz,i,iecore)
226 : enddo
227 0 : fracmass = 1._f - rmasscore/totmass
228 0 : elemass = fracmass * rmass(i,ig)
229 0 : rmi = elemass
230 :
231 : endif ! pure CN group or CN group w/ cores
232 :
233 : endif ! src group is CN and different from the target group
234 :
235 15168083855 : elseif( itype(ielem) .eq. I_CORE2MOM )then ! core mass^2
236 :
237 0 : i_pkern = 6 ! Use different kernel for core mass^2 production
238 0 : rmj = 1._f
239 0 : if( itype(ienconc(ig)) .eq. I_INVOLATILE ) then
240 0 : rmi = rmass(i,ig)
241 0 : rmj = rmass(j,jg)
242 : endif
243 : endif ! itype(ielem)
244 :
245 43199642230 : if( itype(ielem) .ne. I_CORE2MOM )then
246 0 : coagpe(iz,ibin,ielem) = coagpe(iz,ibin,ielem) + &
247 0 : pc(iz,i,iefrom)*pcl(iz,j,je)*rmi * &
248 0 : ckernel(iz,i,j,ig,jg) * &
249 43199642230 : pkernel(i,j,ig,jg,igrp,i_pkern)
250 : else
251 0 : coagpe(iz,ibin,ielem) = coagpe(iz,ibin,ielem) + &
252 0 : ( pc(iz,i,iefrom)*pcl(iz,j,je)*rmi**2 + &
253 0 : pc(iz,i,iefrom_cm)*rmi* &
254 0 : pcl(iz,j,je_cm)*rmj ) * &
255 0 : ckernel(iz,i,j,ig,jg) * &
256 0 : pkernel(i,j,ig,jg,igrp,i_pkern)
257 : endif
258 : endif ! end of ( pconmax .gt. FEW_PC )
259 : enddo ! iz = 1, NZ
260 : endif ! end of (iefrom .ne. 0)
261 : enddo ! iquad
262 :
263 : ! Return to caller with coagulation production terms evaluated.
264 :
265 147087360 : return
266 147087360 : end
|
