File: C:\NOAA\NEMS_11731\src\chem\gocart\src\Components\GOCART_GridComp\CARMA_GridComp\setupcoag.F90

1            subroutine setupcoag ( carma, rc )
2     
3     !     types
4           use carma_types_mod
5     
6           implicit none
7     
8           integer, intent(out) :: rc
9      
10     !  @(#) setupcoag.f  Jensen  Oct-1995
11     !  This routine sets up mapping arrays and precomputed
12     !  factors for coagulation.
13     !
14     !  This routine requires that <ckernel> has been defined.
15     !  (setupckern.f must be called before this)
16     !
17     !
18     !  Argument list input:
19     !    None.
20     !
21     !  Argument list output:
22     !    None.
23     !
24     !
25     !  Include global constants and variables
26     !
27     !      include 'globaer.h'
28     !
29     !
30     !      dimension kbin(NGROUP,NGROUP,NGROUP,NBIN,NBIN)
31     integer, allocatable :: kbin(:,:,:,:,:)
32     integer :: ios
33     
34     !--
35      
36     integer :: ielem, isolto, icompto, igto, ig, iepart 
37     integer :: icompfrom, ic, iecore  
38     integer :: isolfrom
39     integer :: igrp, jg, i, j , ipair
40     real(kind=f) :: rmsum
41     integer :: ibin
42     real(kind=f) :: rmkbin
43     integer :: kb, ncg 
44     real(kind=f) :: rmk 
45     integer :: k
46     real(kind=f) :: pkernl, pkernu
47     integer :: ix, iy
48     
49     !--
50     
51     #include "carma_globaer.h"
52     
53     !--
54     !  Announce entry to this routine
55     !
56     !      if( DEBUG ) write(LUNOPRT,'(/,a)') 'Enter setupcoag'
57     #ifdef DEBUG
58       write(*,*) '+ setupcoag'
59     #endif
60     
61     rc = 0
62     
63     allocate( kbin( carma%ngroup, carma%ngroup, carma%ngroup, &
64                     carma%nbin, carma%nbin ), stat=ios)
65     if(ios /= 0) then
66      rc = 1
67      return
68     endif
69     
70     !-----------------------------------------------------------------
71     !
72     !
73     !  For each element <ielem> and each group <ig>, determine which element in <ig>
74     !  contributes to production  in <ielem>: <icoagelem(ielem,ig)>.
75     !  If no elements in <ig> are transfered into element <ielem> during coagulation,
76     !  then set <icoagelem(ielem,ig)> to 0.
77     !
78           do ielem = 1,NELEM
79     
80             isolto = isolelem(ielem)           ! target solute type
81             icompto = icomp(ielem)             ! target element compound
82             igto = igelem(ielem)               ! target group
83     
84             do ig = 1, NGROUP                 ! source group
85     
86               iepart = ienconc(ig)            ! source particle number concentration element
87     !
88     !
89     !  If <ielem> is particle number concentration type or <ig> is pure group, then
90     !  the source element is the particle number concentration element of group <ig>.
91     !
92               if( ( itype(ielem) .eq. I_INVOLATILE .or. &
93                     itype(ielem) .eq. I_VOLATILE ) .or. &
94                     ncore(ig) .eq. 0 ) then
95                 icoagelem(ielem,ig) = iepart
96     !
97     !
98     !  Otherwise, use element compound names to determine which source element matches
99     !  target core mass element.
100     !
101               else
102                 icompfrom = icomp(iepart)       ! source element compound
103                 if( icompfrom .eq. icompto )then
104                   icoagelem(ielem,ig) = iepart
105                 else
106                   do ic = 1,ncore(ig)
107                     iecore = icorelem(ic,ig)       ! absolute element number of core
108                     icompfrom = icomp(iecore)    ! source element compound
109                     if( icompfrom .eq. icompto ) &
110                       icoagelem(ielem,ig) = iecore
111                   enddo
112                 endif
113               endif
114     !
115     !
116     !  Otherwise, use solute types to determine which source element matches
117     !  target core mass element.
118     !
119     !          else
120     !            isolfrom = isolelem(ienconc(ig))   ! source solute type
121     !            if( isolfrom .eq. isolto ) then
122     !              icoagelem(ielem,ig) = ienconc(ig)
123     !            else
124     !              do ic = 1,ncore(ig)
125     !                iecore = icorelem(ic,ig)       ! absolute element number of core
126     !                isolfrom = isolelem(iecore)    ! source solute type
127     !                if( isolfrom .eq. isolto )
128     !     $            icoagelem(ielem,ig) = iecore
129     !              enddo
130     !            endif
131     !          endif
132     !
133     !
134     !  If <ielem> is a core mass type and <ig> is a pure CN group and the
135     !  solutes don't match, then set <icoagelem> to zero to make sure no
136     !  coag production occurs.
137     !
138               if( itype(ielem) .eq. I_COREMASS .and. &
139                   itype(ienconc(ig)).eq. I_INVOLATILE &
140                   .and. ncore(ig) .eq. 0 ) then
141                 isolfrom = isolelem(ienconc(ig))
142                 if( isolfrom .ne. isolto ) then
143                   icoagelem(ielem,ig) = 0
144                 endif
145               endif
146     
147             enddo          ! end of (ig = 1, NGROUP)
148     
149           enddo            ! end of (ielem = 1,NELEM)
150     
151     !
152     !
153     !  Calculate lower bin <kbin> which coagulated particle goes into
154     !  and make sure it is less than <NBIN>+1
155     !
156     !  Colliding particles come from group <ig>, bin <i> and group <jg>, bin <j>
157     !  Resulting particle lands in group <igrp>, between <ibin> and <ibin> + 1
158     !
159           do igrp = 1, NGROUP
160             do ig = 1, NGROUP
161               do jg = 1, NGROUP
162                 do i = 1, NBIN
163                   do j = 1, NBIN
164     
165                     rmsum = rmass(i,ig) + rmass(j,jg)
166     
167                     do ibin = 1, NBIN-1
168                       if( rmsum .ge. rmass(ibin,igrp) .and. &
169                           rmsum .lt. rmass(ibin+1,igrp) ) then
170                         kbin(igrp,ig,jg,i,j) = ibin
171                       endif
172                     enddo
173     
174                     if( rmsum .ge. rmass(NBIN,igrp) ) &
175                              kbin(igrp,ig,jg,i,j) = NBIN
176     
177                   enddo
178                 enddo
179               enddo
180             enddo
181           enddo
182     !
183     !
184     !  Calculate partial loss fraction
185     !
186     !  This fraction is needed because when a particle in bin <i> collides
187     !  with a particle in bin <j> resulting in a particle whose mass falls
188     !  between <i> and <i>+1, only partial loss occurs from bin <i>.
189     !
190     !  Since different particle groups have different radius grids, this
191     !  fraction is a function of the colliding groups and the resulting group.
192     !
193           do igrp = 1, NGROUP
194             do ig = 1, NGROUP
195               do jg = 1, NGROUP
196     
197                 if( igrp .eq. icoag(ig,jg) ) then
198     
199                   do i = 1, NBIN
200                     do j = 1,NBIN
201     
202                       volx(igrp,ig,jg,i,j) = 1.
203      
204                       if(kbin(igrp,ig,jg,i,j).eq.i) then
205      
206                         rmkbin = rmass(kbin(igrp,ig,jg,i,j),igrp)
207                         volx(igrp,ig,jg,i,j) = 1. -                       &
208                             (rmrat(igrp)*rmkbin-rmass(i,ig)-rmass(j,jg))  &
209                             /(rmrat(igrp)*rmkbin-rmkbin)*                 &
210                             rmass(i,ig)/(rmass(i,ig) + rmass(j,jg))
211     
212                       endif
213     
214                     enddo
215                   enddo
216                 endif
217               enddo
218             enddo
219           enddo
220     !
221     !
222     !  Calculate mapping functions that specify sets of quadruples
223     !  (group pairs and bin pairs) that contribute to production
224     !  in each bin. Mass transfer from <ig,i> to <igrp,ibin> occurs due to
225     !  collisions between particles in <ig,i> and particles in <jg,j>.
226     !  2 sets of quadruples must be generated:
227     !     low: k = ibin and (k != i or ig != igrp)  and  icoag(ig,jg) = igrp
228     !      up: k+1 = ibin        and  icoag(ig,jg) = igrp
229     !
230     !  npair#(igrp,ibin) is the number of pairs in each set (# = l,u)
231     !  i#, j#, ig#, and jg# are the bin pairs and group pairs in each
232     !  set (# = low, up)
233     !
234           do igrp = 1, NGROUP
235             do ibin = 1, NBIN
236     
237               npairl(igrp,ibin) = 0
238               npairu(igrp,ibin) = 0
239     
240               do ig = 1, NGROUP
241               do jg = 1, NGROUP
242                 do i = 1, NBIN
243                 do j = 1, NBIN
244     
245                   kb = kbin(igrp,ig,jg,i,j)
246                   ncg = icoag(ig,jg)
247     
248                   if( kb+1.eq.ibin .and. ncg.eq.igrp ) then
249                     npairu(igrp,ibin) = npairu(igrp,ibin) + 1
250                     iup(igrp,ibin,npairu(igrp,ibin)) = i
251                     jup(igrp,ibin,npairu(igrp,ibin)) = j
252                     igup(igrp,ibin,npairu(igrp,ibin)) = ig
253                     jgup(igrp,ibin,npairu(igrp,ibin)) = jg
254                   endif
255     
256                   if( kb.eq.ibin .and. ncg.eq.igrp .and. &
257                     (i.ne.ibin .or. ig.ne.igrp) ) then
258                     npairl(igrp,ibin) = npairl(igrp,ibin) + 1
259                     ilow(igrp,ibin,npairl(igrp,ibin)) = i
260                     jlow(igrp,ibin,npairl(igrp,ibin)) = j
261                     iglow(igrp,ibin,npairl(igrp,ibin)) = ig
262                     jglow(igrp,ibin,npairl(igrp,ibin)) = jg
263                   endif
264     
265                 enddo
266                 enddo
267               enddo
268               enddo
269             enddo
270           enddo
271     !
272     !
273     !  Do some extra debugging reports  (normally commented)
274     !
275     #ifdef DEBUG
276           write(LUNOPRT,*) ' '
277           write(LUNOPRT,*) 'Coagulation group mapping:'
278           do ig = 1, NGROUP
279             do jg = 1, NGROUP
280               write(LUNOPRT,*) 'ig jg icoag = ', ig, jg, icoag(ig,jg)
281             enddo
282           enddo
283           write(LUNOPRT,*) ' '
284           write(LUNOPRT,*) 'Coagulation element mapping:'
285           do ielem = 1, NELEM
286             do ig = 1, NGROUP
287               write(LUNOPRT,*) 'ielem ig icoagelem icomp(ielem) = ', &
288                ielem, ig, icoagelem(ielem,ig), icomp(ielem)
289             enddo
290           enddo
291           write(LUNOPRT,*) ' '
292           write(LUNOPRT,*) 'Coagulation bin mapping arrays'
293           do igrp = 1, NGROUP
294             do ibin = 1,3
295               write(LUNOPRT,*) 'igrp, ibin = ',igrp, ibin
296               do ipair = 1,npairl(igrp,ibin)
297                 write(LUNOPRT,*) 'low:np,ig,jg,i,j ', &
298                     ipair,iglow(igrp,ibin,ipair), &
299                 jglow(igrp,ibin,ipair), ilow(igrp,ibin,ipair), &
300                       jlow(igrp,ibin,ipair)
301               enddo
302               do ipair = 1,npairu(igrp,ibin)
303                 write(LUNOPRT,*) 'up:np,ig,jg,i,j ', &
304                     ipair,igup(igrp,ibin,ipair), &
305                 jgup(igrp,ibin,ipair), iup(igrp,ibin,ipair), &
306                      jup(igrp,ibin,ipair)
307               enddo
308             enddo
309           enddo
310     stop
311     #endif
312     !
313     !
314     !  Calculate variables needed in routine coagp.f
315     !
316           do ix = 1, NX
317           do iy = 1, NY
318     
319           ckernel => carma%ckernel(ix,iy)%data5d
320           pkernel => carma%pkernel(ix,iy)%data7d
321     
322           do igrp = 1, NGROUP
323             do ig = 1, NGROUP
324             do jg = 1, NGROUP
325     
326               if( igrp .eq. icoag(ig,jg) ) then
327     
328                 do i = 1, NBIN
329                 do j = 1, NBIN
330     
331                   rmk = rmass(kbin(igrp,ig,jg,i,j),igrp)
332                   rmsum = rmass(i,ig) + rmass(j,jg)
333     
334                   do k = 1, NZ
335     
336                     pkernl = ckernel(k,i,j,ig,jg)* &
337                              (rmrat(igrp)*rmk - rmsum) / &
338                              (rmrat(igrp)*rmk - rmk)
339     
340                     pkernu = ckernel(k,i,j,ig,jg)* &
341                              (rmsum - rmk) / &
342                              (rmrat(igrp)*rmk - rmk)
343     
344                     if( kbin(igrp,ig,jg,i,j) .eq. NBIN )then
345                       pkernl = ckernel(k,i,j,ig,jg)* &
346                                rmsum / rmass(NBIN,igrp)
347                       pkernu = 0.
348                     endif
349       
350                     pkernel(k,i,j,ig,jg,igrp,1) = pkernu * &
351                                                   rmass(i,ig)/rmsum
352                     pkernel(k,i,j,ig,jg,igrp,2) = pkernl * &
353                                                   rmass(i,ig)/rmsum
354                     pkernel(k,i,j,ig,jg,igrp,3) = pkernu * &
355                                                   rmk*rmrat(igrp)/rmsum
356                     pkernel(k,i,j,ig,jg,igrp,4) = pkernl * &
357                                                   rmk/rmsum
358                     pkernel(k,i,j,ig,jg,igrp,5) = pkernu * &
359                                            ( rmk*rmrat(igrp)/rmsum )**2
360                     pkernel(k,i,j,ig,jg,igrp,6) = pkernl * &
361                                                   ( rmk/rmsum )**2
362     
363                   enddo  ! k
364     
365                 enddo
366                 enddo
367     
368               endif
369     
370             enddo
371             enddo
372           enddo
373     
374           enddo  ! iy
375           enddo  ! ix
376     
377     !
378     !
379     !  Return to caller with coagulation mapping arrays defined
380     !
381           deallocate(kbin, stat = ios)
382           if(ios /= 0) rc = 1
383     
384           return
385           end
386