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 sets up mapping arrays for coagulation. It only computes varaibles that
6 : !! are independent of the model state. The calculation of factors needed for coagulation
7 : !! that depend on state are calculated in <i>setupckern</i>.
8 : !!
9 : !! @author Eric Jensen
10 : !! @ version Oct-1995
11 1536 : subroutine setupcoag(carma, rc)
12 :
13 : ! types
14 : use carma_precision_mod
15 : use carma_enums_mod
16 : use carma_constants_mod
17 : use carma_types_mod
18 : use carma_mod
19 :
20 : implicit none
21 :
22 : type(carma_type), intent(inout) :: carma !! the CARMA object
23 : integer, intent(inout) :: rc !! return code, negative indicates failure
24 :
25 : ! Local declarations
26 : integer :: ielem, isolto, icompto, igto, ig, iepart
27 : integer :: icompfrom, ic, iecore
28 : integer :: isolfrom
29 : integer :: igrp, jg, i, j , ipair
30 : real(kind=f) :: rmsum
31 : integer :: ibin
32 : real(kind=f) :: rmkbin
33 : integer :: kb, ncg
34 : real(kind=f) :: rmk
35 : logical :: fill_bot ! used for filling <icoag>
36 : integer :: irow, icol
37 : logical :: isCoag
38 : integer :: igtest
39 : real(kind=f) :: pkernl, pkernu
40 :
41 :
42 : ! NOTE: Moved this section from from setupckern.f, since it is not dependent on the
43 : ! model's state.
44 : !
45 : ! Fill <icoag>, maintaining diagonal symmetry
46 : ! -------------------------------------------
47 : ! Fill bottom of matrix if non-zero term(s) in upper half;
48 : ! also check for non-zero, non-matching, non-diagonal terms.
49 1536 : fill_bot = .true.
50 3072 : do irow = 2, NGROUP
51 4608 : do icol = 1, irow-1
52 3072 : if( icoag(irow,icol) .ne. 0 )then
53 0 : fill_bot = .false.
54 0 : if( icoag(icol,irow) .ne. 0 .and. &
55 : icoag(icol,irow) .ne. icoag(irow,icol) )then
56 0 : if (do_print) write(LUNOPRT, *) 'setupcoag::ERROR bad icoag array'
57 0 : rc = -1
58 0 : return
59 : endif
60 : endif
61 : enddo
62 : enddo
63 :
64 3072 : do ig = 2, NGROUP
65 4608 : do jg = 1, ig-1
66 1536 : if( fill_bot )then
67 1536 : irow = ig
68 1536 : icol = jg
69 : else
70 0 : irow = jg
71 0 : icol = ig
72 : endif
73 3072 : icoag(irow,icol) = icoag(icol,irow)
74 : enddo
75 : enddo
76 :
77 : ! Initialize <icoagelem> with zeros
78 12288 : do ielem = 1,NELEM
79 33792 : do ig = 1,NGROUP
80 21504 : icoagelem(ielem,ig) = 0
81 32256 : icoagelem_cm(ielem,ig) = 0
82 : enddo
83 : enddo
84 :
85 : ! For each element <ielem> and each group <ig>, determine which element in <ig>
86 : ! contributes to production in <ielem>: <icoagelem(ielem,ig)>.
87 : ! If no elements in <ig> are transfered into element <ielem> during coagulation,
88 : ! then set <icoagelem(ielem,ig)> to 0.
89 12288 : do ielem = 1,NELEM
90 10752 : isolto = isolelem(ielem) ! target solute type
91 10752 : icompto = icomp(ielem) ! target element compound
92 10752 : igto = igelem(ielem) ! target group
93 :
94 33792 : do ig = 1, NGROUP ! source group
95 : ! source particle number concentration element
96 21504 : iepart = ienconc(ig)
97 :
98 : ! source element compound
99 21504 : icompfrom = icomp(iepart)
100 :
101 : ! Check to see if the target group is produced by coagulation of any
102 : ! group with the source group.
103 21504 : isCoag = .FALSE.
104 :
105 64512 : do igtest = 1, NGROUP
106 64512 : if (icoag(ig, igtest) .eq. igto .or. icoag(igtest, ig) .eq. igto) then
107 29184 : isCoag = .TRUE.
108 : endif
109 : end do
110 :
111 : ! Only find the source production element if the group igto can
112 : ! be produced by coagulation from group ig.
113 32256 : if (isCoag) then
114 :
115 : ! If <ig> only has no cores, then the only way to make particles
116 : ! would be if the one element <iepart> is the same type as the
117 : ! source.
118 19968 : if( ncore(ig) .eq. 0 ) then
119 :
120 10752 : if( icompfrom .eq. icompto )then
121 3072 : icoagelem(ielem,ig) = iepart
122 : endif
123 : else
124 :
125 : ! Search the elements in the group to see if one has the same
126 : ! type as the source.
127 :
128 : ! First check the particle number concentration element of the group.
129 : !
130 : ! NOTE: No matter what else happens, you need to adjust the total
131 : ! particle mass.
132 9216 : if( icompfrom .eq. icompto )then
133 1536 : icoagelem(ielem,ig) = iepart
134 : else
135 :
136 : ! Now check the other cores for a match.
137 46080 : do ic = 1,ncore(ig)
138 38400 : iecore = icorelem(ic,ig) ! absolute element number of core
139 38400 : icompfrom = icomp(iecore) ! source element compound
140 :
141 46080 : if( icompfrom .eq. icompto ) then
142 :
143 : ! For core second moment elements, we need additional pairs of source
144 : ! elements c to account for core moment production due to products
145 : ! of source particle core mass.
146 7680 : if( itype(ielem) .eq. I_CORE2MOM )then
147 0 : icoagelem_cm(ielem,ig) = iecore
148 0 : icoagelem(ielem,ig) = imomelem(ig)
149 : else
150 7680 : icoagelem(ielem,ig) = iecore
151 : endif
152 : endif
153 : enddo
154 : endif
155 : endif
156 :
157 : ! If <ielem> is a core mass type and <ig> is a pure CN group and the
158 : ! solutes don't match, then set <icoagelem> to zero to make sure no
159 : ! coag production occurs.
160 19968 : if( itype(ielem) .eq. I_COREMASS .and. &
161 0 : itype(ienconc(ig)).eq. I_INVOLATILE &
162 39936 : .and. ncore(ig) .eq. 0 ) then
163 0 : isolfrom = isolelem(ienconc(ig))
164 0 : if( isolfrom .ne. isolto ) then
165 0 : icoagelem(ielem,ig) = 0
166 : endif
167 : endif
168 :
169 : ! If there is a source and this is a multi-component group,
170 : ! then we need to make sure that the particle concentration
171 : ! of the group also gets updated, since this keeps track of
172 : ! the total mass.
173 19968 : if (icoagelem(ielem,ig) .ne. 0) then
174 12288 : if (ncore(igto) .ne. 0 .and. ielem .ne. ienconc(igto)) then
175 7680 : icoagelem(ienconc(igto), ig) = iepart
176 : endif
177 : endif
178 :
179 : endif
180 : enddo ! end of (ig = 1, NGROUP)
181 : enddo ! end of (ielem = 1,NELEM)
182 :
183 :
184 : ! Coagulation won't work properly if any of the elements are produced by
185 : ! items that come later in the element list than themselves. Report an
186 : ! error if that is the case.
187 12288 : do ielem = 1, NELEM
188 33792 : do ig = 1, NGROUP
189 32256 : if (icoagelem(ielem, ig) .gt. ielem) then
190 0 : if (do_print) write(LUNOPRT, '(a,i3,a,i3,a)') &
191 0 : 'setupcoag::ERROR For coagulation, element (', &
192 0 : icoagelem(ielem,ig), ') must come before (', ielem, &
193 0 : ') in the element list.'
194 0 : rc = -1
195 0 : return
196 : endif
197 : enddo
198 : enddo
199 :
200 :
201 : ! Calculate lower bin <kbin> which coagulated particle goes into
202 : ! and make sure it is less than <NBIN>+1
203 : !
204 : ! Colliding particles come from group <ig>, bin <i> and group <jg>, bin <j>
205 : ! Resulting particle lands in group <igrp>, between <ibin> and <ibin> + 1
206 4608 : do igrp = 1, NGROUP
207 10752 : do ig = 1, NGROUP
208 21504 : do jg = 1, NGROUP
209 264192 : do i = 1, NBIN
210 5173248 : do j = 1, NBIN
211 :
212 4915200 : rmsum = rmass(i,ig) + rmass(j,jg)
213 :
214 98304000 : do ibin = 1, NBIN-1
215 98304000 : if( rmsum .ge. rmass(ibin,igrp) .and. rmsum .lt. rmass(ibin+1,igrp) ) then
216 3592704 : kbin(igrp,ig,jg,i,j) = ibin
217 : endif
218 : enddo
219 :
220 4915200 : ibin = NBIN
221 5160960 : if( rmsum .ge. rmass(ibin,igrp) ) kbin(igrp,ig,jg,i,j) = NBIN
222 : enddo
223 : enddo
224 : enddo
225 : enddo
226 : enddo
227 :
228 : ! Calculate partial loss fraction
229 : !
230 : ! This fraction is needed because when a particle in bin <i> collides
231 : ! with a particle in bin <j> resulting in a particle whose mass falls
232 : ! between <i> and <i>+1, only partial loss occurs from bin <i>.
233 : !
234 : ! Since different particle groups have different radius grids, this
235 : ! fraction is a function of the colliding groups and the resulting group.
236 4608 : do igrp = 1, NGROUP
237 10752 : do ig = 1, NGROUP
238 21504 : do jg = 1, NGROUP
239 :
240 18432 : if( igrp .eq. icoag(ig,jg) ) then
241 :
242 129024 : do i = 1, NBIN
243 2586624 : do j = 1,NBIN
244 2457600 : volx(igrp,ig,jg,i,j) = 1._f
245 :
246 2580480 : if(kbin(igrp,ig,jg,i,j).eq.i) then
247 :
248 1208832 : ibin = kbin(igrp,ig,jg,i,j)
249 1208832 : rmkbin = rmass(ibin,igrp)
250 0 : volx(igrp,ig,jg,i,j) = 1._f - &
251 2417664 : (rmrat(igrp)*rmkbin-rmass(i,ig)-rmass(j,jg)) &
252 : /(rmrat(igrp)*rmkbin-rmkbin)* &
253 3626496 : rmass(i,ig)/(rmass(i,ig) + rmass(j,jg))
254 : endif
255 : enddo
256 : enddo
257 : endif
258 : enddo
259 : enddo
260 : enddo
261 :
262 : ! Calculate mapping functions that specify sets of quadruples
263 : ! (group pairs and bin pairs) that contribute to production
264 : ! in each bin. Mass transfer from <ig,i> to <igrp,ibin> occurs due to
265 : ! collisions between particles in <ig,i> and particles in <jg,j>.
266 : ! 2 sets of quadruples must be generated:
267 : ! low: k = ibin and (k != i or ig != igrp) and icoag(ig,jg) = igrp
268 : ! up: k+1 = ibin and icoag(ig,jg) = igrp
269 : !
270 : ! npair#(igrp,ibin) is the number of pairs in each set (# = l,u)
271 : ! i#, j#, ig#, and jg# are the bin pairs and group pairs in each
272 : ! set (# = low, up)
273 4608 : do igrp = 1, NGROUP
274 66048 : do ibin = 1, NBIN
275 :
276 61440 : npairl(igrp,ibin) = 0
277 61440 : npairu(igrp,ibin) = 0
278 :
279 187392 : do ig = 1, NGROUP
280 430080 : do jg = 1, NGROUP
281 5283840 : do i = 1, NBIN
282 103464960 : do j = 1, NBIN
283 98304000 : kb = kbin(igrp,ig,jg,i,j)
284 98304000 : ncg = icoag(ig,jg)
285 :
286 98304000 : if( kb+1.eq.ibin .and. ncg.eq.igrp ) then
287 2276352 : npairu(igrp,ibin) = npairu(igrp,ibin) + 1
288 2276352 : iup(igrp,ibin,npairu(igrp,ibin)) = i
289 2276352 : jup(igrp,ibin,npairu(igrp,ibin)) = j
290 2276352 : igup(igrp,ibin,npairu(igrp,ibin)) = ig
291 2276352 : jgup(igrp,ibin,npairu(igrp,ibin)) = jg
292 : endif
293 :
294 103219200 : if( kb.eq.ibin .and. ncg.eq.igrp .and. (i.ne.ibin .or. ig.ne.igrp) ) then
295 1279488 : npairl(igrp,ibin) = npairl(igrp,ibin) + 1
296 1279488 : ilow(igrp,ibin,npairl(igrp,ibin)) = i
297 1279488 : jlow(igrp,ibin,npairl(igrp,ibin)) = j
298 1279488 : iglow(igrp,ibin,npairl(igrp,ibin)) = ig
299 1279488 : jglow(igrp,ibin,npairl(igrp,ibin)) = jg
300 : endif
301 : enddo
302 : enddo
303 : enddo
304 : enddo
305 : enddo
306 : enddo
307 :
308 :
309 : ! NOTE: Split ckernel out of pkernel, so that it can be made independent of model state.
310 : ! It also reduces the size of the tables and should improve the intialization time.
311 :
312 : ! Calculate variables needed in routine coagp.f
313 4608 : do igrp = 1, NGROUP
314 10752 : do jg = 1, NGROUP
315 21504 : do ig = 1, NGROUP
316 :
317 18432 : if( igrp .eq. icoag(ig,jg) ) then
318 :
319 129024 : do j = 1, NBIN
320 2586624 : do i = 1, NBIN
321 :
322 2457600 : ibin = kbin(igrp,ig,jg,i,j)
323 2457600 : rmk = rmass(ibin,igrp)
324 2457600 : rmsum = rmass(i,ig) + rmass(j,jg)
325 :
326 2457600 : pkernl = (rmrat(igrp)*rmk - rmsum) / (rmrat(igrp)*rmk - rmk)
327 :
328 2457600 : pkernu = (rmsum - rmk) / (rmrat(igrp)*rmk - rmk)
329 :
330 2457600 : if( ibin .eq. NBIN )then
331 181248 : pkernl = rmsum / rmass(ibin,igrp)
332 181248 : pkernu = 0._f
333 : endif
334 :
335 2457600 : pkernel(i,j,ig,jg,igrp,1) = pkernu * rmass(i,ig)/rmsum
336 2457600 : pkernel(i,j,ig,jg,igrp,2) = pkernl * rmass(i,ig)/rmsum
337 2457600 : pkernel(i,j,ig,jg,igrp,3) = pkernu * rmk*rmrat(igrp)/rmsum
338 2457600 : pkernel(i,j,ig,jg,igrp,4) = pkernl * rmk/rmsum
339 2457600 : pkernel(i,j,ig,jg,igrp,5) = pkernu * ( rmk*rmrat(igrp)/rmsum )**2
340 2580480 : pkernel(i,j,ig,jg,igrp,6) = pkernl * ( rmk/rmsum )**2
341 : enddo
342 : enddo
343 : endif
344 : enddo
345 : enddo
346 : enddo
347 :
348 : ! Do some extra debugging reports (normally commented)
349 1536 : if (do_print_init) then
350 2 : write(LUNOPRT,*) ' '
351 2 : write(LUNOPRT,*) 'Coagulation group mapping:'
352 6 : do ig = 1, NGROUP
353 14 : do jg = 1, NGROUP
354 12 : write(LUNOPRT,*) 'ig jg icoag = ', ig, jg, icoag(ig,jg)
355 : enddo
356 : enddo
357 2 : write(LUNOPRT,*) ' '
358 2 : write(LUNOPRT,*) 'Coagulation element mapping:'
359 16 : do ielem = 1, NELEM
360 44 : do ig = 1, NGROUP
361 28 : write(LUNOPRT,*) 'ielem ig icoagelem icomp(ielem) = ', &
362 70 : ielem, ig, icoagelem(ielem,ig), icomp(ielem)
363 : enddo
364 : enddo
365 2 : write(LUNOPRT,*) ' '
366 2 : write(LUNOPRT,*) 'Coagulation bin mapping arrays'
367 6 : do igrp = 1, NGROUP
368 18 : do ibin = 1,3
369 12 : write(LUNOPRT,*) 'igrp, ibin = ',igrp, ibin
370 112 : do ipair = 1,npairl(igrp,ibin)
371 100 : write(LUNOPRT,*) 'low:np,ig,jg,i,j ', &
372 100 : ipair,iglow(igrp,ibin,ipair), &
373 100 : jglow(igrp,ibin,ipair), ilow(igrp,ibin,ipair), &
374 212 : jlow(igrp,ibin,ipair)
375 : enddo
376 136 : do ipair = 1,npairu(igrp,ibin)
377 120 : write(LUNOPRT,*) 'up:np,ig,jg,i,j ', &
378 120 : ipair,igup(igrp,ibin,ipair), &
379 120 : jgup(igrp,ibin,ipair), iup(igrp,ibin,ipair), &
380 252 : jup(igrp,ibin,ipair)
381 : enddo
382 : enddo
383 : enddo
384 : endif
385 :
386 : ! Return to caller with coagulation mapping arrays defined
387 : return
388 1536 : end
|
