File: /nfs/home/0/users/jenkins/mfix.git/model/qmomk/qmomk_collision_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module:  qmomk_collisions                                           C
4     !  Purpose: Collision operators for the kinetic equation               C
5     !  Contains the following subroutines and functions:                   C
6     !      collisions_istantaneous, collisions_bgk,                        C
7     !      compute_collision_time, collisions_boltzmann_one_specie,        C
8     !      collisions_boltzmann_two_specie,                                C
9     !      solve_boltzmann_collisions_one_specie,                          C
10     !      solve_boltzmann_collisions_two_specie, radial_g0                C
11     !                                                                      C
12     !  Author: A. Passalacqua                           Date:              C
13     !  Reviewer:                                        Date:              C
14     !                                                                      C
15     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16     MODULE qmomk_collision
17     
18     
19       USE qmomk_parameters
20       USE qmomk_quadrature
21     
22       IMPLICIT NONE
23     
24       PRIVATE
25     
26       PUBLIC :: COLLISIONS_ISTANTANEOUS
27       PUBLIC :: COLLISIONS_BGK
28       PUBLIC :: COMPUTE_COLLISION_TIME
29       !PUBLIC :: GLOBAL_COLLISION_TIME
30       PUBLIC :: SOLVE_BOLTZMANN_COLLISIONS_ONE_SPECIE
31       PUBLIC :: SOLVE_BOLTZMANN_COLLISIONS_TWO_SPECIES
32     
33       PUBLIC :: RADIAL_G0
34     
35     CONTAINS
36     
37     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
38     !                                                                      C
39     !  Instantaneous collisions                                            C
40     !                                                                      C
41     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
42     
43       SUBROUTINE COLLISIONS_ISTANTANEOUS(M, dp)
44     
45         IMPLICIT NONE
46     
47         DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: M
48         DOUBLE PRECISION, INTENT(IN) :: dp
49     
50         DOUBLE PRECISION :: m0, m1, m2, m3, m11, m22, m33
51     
52         DOUBLE PRECISION :: up1, up2, up3
53         DOUBLE PRECISION :: sig1, sig2, sig3, seq
54     
55         DOUBLE PRECISION :: d11, d12, d13, d22, d23, d33
56         DOUBLE PRECISION :: d111, d112, d113, d122, d123, d133, d222, d223, d233, d333
57     
58         m0 = M(1)
59         m1 = M(2)
60         m2 = M(3)
61         m3 = M(4)
62         m11 = M(5)
63         m22 = M(8)
64         m33 = M(10)
65     
66         up1 = m1/m0
67         up2 = m2/m0
68         up3 = m3/m0
69         sig1 = max(0.D0, m11/m0 - up1**2)
70         sig2 = max(0.D0, m22/m0 - up2**2)
71         sig3 = max(0.D0, m33/m0 - up3**2)
72         seq = ( sig1 + sig2 + sig3 )/3.D0
73     
74         d11 = m0*(seq + up1**2)
75         d12 = m0*up1*up2
76         d13 = m0*up1*up3
77         d22 = m0*(seq + up2**2)
78         d23 = m0*up2*up3
79         d33 = m0*(seq + up3**2)
80         d111 = m1*(3.D0*seq + up1**2)
81         d112 = m2*(seq + up1**2)
82         d113 = m3*(seq + up1**2)
83         d122 = m1*(seq + up2**2)
84         d123 = m0*up1*up2*up3
85         d133 = m1*(seq + up3**2)
86         d222 = m2*(3.D0*seq + up2**2)
87         d223 = m3*(seq + up2**2)
88         d233 = m2*(seq + up3**2)
89         d333 = m3*(3.D0*seq + up3**2)
90     
91         M(5)  =  d11
92         M(6)  =  d12
93         M(7)  =  d13
94         M(8)  =  d22
95         M(9)  =  d23
96         M(10) =  d33
97         M(11) =  d111
98         M(12) =  d112
99         M(13) =  d113
100         M(14) =  d122
101         M(15) =  d123
102         M(16) =  d133
103         M(17) =  d222
104         M(18) =  d223
105         M(19) =  d233
106         M(20) =  d333
107     
108       END SUBROUTINE COLLISIONS_ISTANTANEOUS
109     
110     
111     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
112     !                                                                      C
113     !  Bathnagar-Gross-Krook (BGK) collision operator                      C
114     !                                                                      C
115     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
116     
117       SUBROUTINE COLLISIONS_BGK(M, Dt, tcol, dp, e)
118     
119         USE param1, only: small_number
120         USE constant, only: Pi, EP_STAR
121     
122         IMPLICIT NONE
123     
124         DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: M
125         DOUBLE PRECISION, INTENT(IN) :: Dt, tcol, dp, e
126     
127         DOUBLE PRECISION :: m0, m1, m2, m3, m11, m12, m13, m22, m23, m33
128         DOUBLE PRECISION :: m111, m112, m113, m122, m123, m133, m222
129         DOUBLE PRECISION :: m223, m233, m333
130     
131         DOUBLE PRECISION :: up1, up2, up3, b, a1, b1, gamma, om
132         DOUBLE PRECISION :: sig1, sig2, sig3, sig11, sig12, sig13, sig22, sig23
133         DOUBLE PRECISION :: sig33, seq, T
134     
135         DOUBLE PRECISION :: d11, d12, d13, d22, d23, d33
136         DOUBLE PRECISION :: d111, d112, d113, d122, d123, d133, d222, d223, d233, d333
137     
138         DOUBLE PRECISION :: kap, ALPHA_MAX
139     
140         m0   = M(1)
141         m1   = M(2)
142         m2   = M(3)
143         m3   = M(4)
144         m11  = M(5)
145         m12  = M(6)
146         m13  = M(7)
147         m22  = M(8)
148         m23  = M(9)
149         m33  = M(10)
150         m111 = M(11)
151         m112 = M(12)
152         m113 = M(13)
153         m122 = M(14)
154         m123 = M(15)
155         m133 = M(16)
156         m222 = M(17)
157         m223 = M(18)
158         m233 = M(19)
159         m333 = M(20)
160     
161         up1 = m1/m0
162         up2 = m2/m0
163         up3 = m3/m0
164         sig1 = MAX(0.D0, m11/m0 - up1**2)
165         sig2 = MAX(0.D0, m22/m0 - up2**2)
166         sig3 = MAX(0.D0, m33/m0 - up3**2)
167         T = ( sig1 + sig2 + sig3 )/3.D0
168     
169         !A.P. This parameter is hardcoded:
170         !     b = 0 for BGK, b = -1/2 for ES-BGK
171         b = 0.D0
172         gamma = 1.D0 - b
173     
174         om = (1 + e)/2.D0
175         a1 = gamma*(om**2)
176         b1 = a1 - 2.D0*gamma*om + 1.D0
177     
178         sig11 = a1*T + b1*sig1
179         sig22 = a1*T + b1*sig2
180         sig33 = a1*T + b1*sig3
181     
182         sig12 = b1*(m12/m0 - up1*up2)
183         sig13 = b1*(m13/m0 - up1*up3)
184         sig23 = b1*(m23/m0 - up2*up3)
185     
186         seq = ( sig11 + sig22 + sig33 )/3.D0
187     
188         d11 = m0*(sig11 + up1**2)
189         d12 = m0*(sig12 + up1*up2)
190         d13 = m0*(sig13 + up1*up3)
191         d22 = m0*(sig22 + up2**2)
192         d23 = m0*(sig23 + up2*up3)
193         d33 = m0*(sig33 + up3**2)
194     
195         d111 = -(-3.D0*d11 + 2.D0*m0*up1**2)*up1
196         d112 = -(-d11*up2 - (-2.D0*d12 + 2.D0*up1*up2*m0)*up1)
197     
198         d113 = -(-d11*up3+(-2.D0*d13+2.D0*up1*up3*m0)*up1)
199         d122 = -(-2.D0*d12*up2+(2.D0*(up2**2)*m0-d22)*up1)
200         d123 = -(-d12*up3-d13*up2+(2.D0*up2*up3*m0-d23)*up1)
201         d133 = -(-2.D0*d12*up3+(2.D0*(up3**2)*m0-d33)*up1)
202         d222 = -(-3.D0*d22+2.D0*(up2**2)*m0)*up2
203         d223 = -(-d22*up2+(-2.D0*d23+2.D0*up2*up3*m0)*up2)
204         d233 = -(-2.D0*d23*up3+(2.D0*(up3**2)*m0-d33)*up2)
205         d333 = -(-3.D0*d33+2.D0*(up3**2)*m0)*up3
206     
207     ! JEC: Probable BUG: alpha_max was not defined. include definition as found in
208     ! other places
209         ALPHA_MAX = 1. - EP_STAR
210     
211         IF (m0 > ALPHA_MAX - SMALL_NUMBER) THEN
212           kap = 0.
213         ELSE
214           kap = EXP(-12.0*RADIAL_G0(m0)*m0*Dt*sqrt(T/pi)/(dp*gamma))
215         END IF
216     
217         M(5)  =  d11  - kap*( d11 - m11 )
218         M(6)  =  d12  - kap*( d12 - m12 )
219         M(7)  =  d13  - kap*( d13 - m13 )
220         M(8)  =  d22  - kap*( d22 - m22 )
221         M(9)  =  d23  - kap*( d23 - m23 )
222         M(10) =  d33  - kap*( d33 - m33 )
223         M(11) =  d111 - kap*( d111 - m111 )
224         M(12) =  d112 - kap*( d112 - m112 )
225         M(13) =  d113 - kap*( d113 - m113 )
226         M(14) =  d122 - kap*( d122 - m122 )
227         M(15) =  d123 - kap*( d123 - m123 )
228         M(16) =  d133 - kap*( d133 - m133 )
229         M(17) =  d222 - kap*( d222 - m222 )
230         M(18) =  d223 - kap*( d223 - m223 )
231         M(19) =  d233 - kap*( d233 - m233 )
232         M(20) =  d333 - kap*( d333 - m333 )
233       END SUBROUTINE COLLISIONS_BGK
234     
235     
236     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
237     !                                                                      C
238     !  Collision time                                                      C
239     !                                                                      C
240     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
241     
242       SUBROUTINE COMPUTE_COLLISION_TIME (M, dp, theta, tcol, dt)
243     
244         USE param1, only: small_number
245         USE constant, only: Pi, EP_STAR
246     ! JEC: should ep_star be replaced by ep_star_array?
247     
248         IMPLICIT NONE
249     
250         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NMOM) :: M
251         DOUBLE PRECISION, INTENT(IN) :: theta, dp, dt
252         DOUBLE PRECISION, INTENT(OUT) :: tcol
253     
254         DOUBLE PRECISION :: ALPHA_MAX
255     
256         !A.P.: Assuming m0 = M(1) = solids volume fraction due to the moments scaling
257         ALPHA_MAX = 1. - EP_STAR
258     
259         IF (M(1) > ALPHA_MAX - SMALL_NUMBER) THEN
260           tcol = dt ! If this condition is satisfied, moments are set to
261                     ! equilibrium. No need to integrate on more than one time step
262         ELSE
263           tcol = dp*SQRT(pi/(theta+SMALL_NUMBER))/(12.0*M(1)*RADIAL_G0(M(1))+SMALL_NUMBER)
264         END IF
265       END SUBROUTINE COMPUTE_COLLISION_TIME
266     
267       !     Global collision time based on the mean volume fraction (1 specie)
268       !SUBROUTINE GLOBAL_COLLISION_TIME (alpha, dp, theta, tcol, dt)
269       !  DOUBLE PRECISION, INTENT(IN) :: alpha
270       !  DOUBLE PRECISION, INTENT(IN) :: theta, dp, dt
271       !  DOUBLE PRECISION, INTENT(OUT) :: tcol
272     
273       !  DOUBLE PRECISION :: ALPHA_MAX
274     
275       !  ALPHA_MAX = 1. - EP_STAR
276     
277       ! IF (alpha > ALPHA_MAX - SMALL_NUMBER) THEN
278       !    tcol = dt ! If this condition is satisfied, moments are set to
279                     ! equilibrium. No need to integrate on more than one time step
280       !  ELSE
281       !    tcol = dp*SQRT(pi/(theta+SMALL_NUMBER))/(12.0*alpha*RADIAL_G0(M(1))+SMALL_NUMBER)
282       !  END IF
283       !END SUBROUTINE
284     
285     
286     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
287     !                                                                      C
288     !  Boltzmann collision operator for the monodispersed case             C
289     !                                                                      C
290     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
291     
292       SUBROUTINE COLLISIONS_BOLTZMANN_ONE_SPECIE (N, U, V, W, dp, e, Coll)
293     
294         USE constant, only: Pi
295     
296         IMPLICIT NONE
297     
298         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: N, U, V, W
299         DOUBLE PRECISION, INTENT(IN) :: dp, e
300         DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NMOM) :: Coll
301     
302         DOUBLE PRECISION, DIMENSION(3) :: a1, a, vr, v1r, v2r, g1, g2
303     
304         DOUBLE PRECISION :: C200, C110, C101, C020, C011, C002, C300, C210, C201
305         DOUBLE PRECISION :: C120, C111, C102, C030, C021, C012, C003
306         DOUBLE PRECISION :: d, g, gsqr, o, ni, ui, vi, wi, nj, uj, vj, wj, cross
307         DOUBLE PRECISION :: denom, L11, L12, L13, L21, L22, L23, L31, L32, L33
308     
309         INTEGER :: i, j
310     
311         a1 = (/0.0462,0.0971,0.8235/) ! chosen at random (cannot be colinear with vr)
312         a = 1/((a1(1))**2 + (a1(2))**2 + (a1(3))**2) * a1;
313     
314         C200 = 0.
315         C110 = 0.
316         C101 = 0.
317         C020 = 0.
318         C011 = 0.
319         C002 = 0.
320         C300 = 0.
321         C210 = 0.
322         C201 = 0.
323         C120 = 0.
324         C111 = 0.
325         C102 = 0.
326         C030 = 0.
327         C021 = 0.
328         C012 = 0.
329         C003 = 0.
330     
331         o = (1+e)/2
332         DO i= 1, QMOMK_NN
333            ni = N(i)
334            ui = U(i)
335            vi = V(i)
336            wi = W(i)
337            v1r(1) = ui
338            v1r(2) = vi
339            v1r(3) = wi
340            DO j = 1, 8
341               IF (i == j) THEN
342                  CYCLE ! vi - vj = 0 always => nothing to add to source terms
343               END IF
344               nj = N(j)
345               uj = U(j)
346               vj = V(j)
347               wj = W(j)
348     
349               v2r(1) = uj
350               v2r(2) = vj
351               v2r(3) = wj
352     
353               vr(1) = v1r(1) - v2r(1)
354               vr(2) = v1r(2) - v2r(2)
355               vr(3) = v1r(3) - v2r(3)
356     
357               g = SQRT( vr(1)**2 + vr(2)**2 + vr(3)**2 )
358     
359               g1(1) = vr(1)/g
360               g1(2) = vr(2)/g
361               g1(3) = vr(3)/g
362     
363               cross = a(1)*g1(1) + a(2)*g1(2) + a(3)*g1(3)
364     
365               g2(1) = a(1) - cross*g1(1)
366               g2(2) = a(2) - cross*g1(2)
367               g2(3) = a(3) - cross*g1(3)
368     
369               denom = SQRT( g2(1)**2 + g2(2)**2 + g2(3)**2 )
370     
371               L11 = g2(1)/denom
372               L12 = g2(2)/denom
373               L13 = g2(3)/denom
374               L21 = ( g1(3)*a(2)-g1(2)*a(3) )/denom
375               L22 = ( g1(1)*a(3)-g1(3)*a(1) )/denom
376               L23 = ( g1(2)*a(1)-g1(1)*a(2) )/denom
377               L31 = g1(1)
378               L32 = g1(2)
379               L33 = g1(3)
380     
381               gsqr = g**2
382     
383               C200 = C200 + ni*nj*(pi*gsqr*o*(4*g*o*(L31**2)+6*uj*L31-6*ui*L31+g*o*(L21**2)+g*o*(L11**2)))/6
384     
385               C110 = C110 + ni*nj*(pi*gsqr*o*(4*g*o*L31*L32+3*uj*L32-3*ui*L32+3*vj*L31-3*vi*L31+g*o*L21*L22 + &
386                    g*o*L11*L12))/6
387     
388               C101 = C101 + ni*nj*(pi*gsqr*o*(4*g*o*L31*L33+3*uj*L33-3*ui*L33+3*wj*L31-3*wi*L31+g*o*L21*L23 + &
389                    g*o*L11*L13))/6
390     
391               C020 = C020 + ni*nj*(pi*gsqr*o*(4*g*o*(L32**2)+6*vj*L32-6*vi*L32+g*o*(L22**2)+g*o*(L12**2)))/6
392     
393               C011 = C011 + ni*nj*(pi*gsqr*o*(4*g*o*L32*L33+3*vj*L33-3*vi*L33+3*wj*L32-3*wi*L32+g*o*L22*L23 + &
394                    g*o*L12*L13))/6
395     
396               C002 = C002 + ni*nj*(pi*gsqr*o*(4*g*o*(L33**2)+6*wj*L33-6*wi*L33+g*o*(L23**2)+g*o*(L13**2)))/6
397     
398               C300 = C300 + ni*nj*(pi*gsqr*o*(uj+ui)*(4*g*o*(L31**2)+6*uj*L31-6*ui*L31+g*o*(L21**2)+g*o*(L11**2)))/4
399     
400               C210 = C210 + ni*nj*(pi*gsqr*o*(8*g*o*uj*L31*L32+8*g*o*ui*L31*L32+6*(uj**2)*L32-6*(ui**2)*L32 + &
401                    4*g*o*vj*(L31**2)+4*g*o*vi*(L31**2)+12*uj*vj*L31-12*ui*vi*L31+2*g*o*uj*L21*L22 + &
402                    2*g*o*ui*L21*L22+g*o*vj*(L21**2)+g*o*vi*(L21**2)+2*g*o*uj*L11*L12+2*g*o*ui*L11*L12 + &
403                    g*o*vj*(L11**2)+g*o*vi*(L11**2)))/12
404     
405               C201 = C201 + ni*nj*(((8*pi*(g**3)*(o**2)*uj+8*pi*(g**3)*(o**2)*ui)*L31+6*pi*gsqr*o*(uj**2)-6*pi*gsqr*o*(ui**2))*L33 + &
406                    (4*pi*(g**3)*(o**2)*wj+4*pi*(g**3)*(o**2)*wi)*(L31**2)+(12*pi*(g**2)*o*uj*wj-12*pi*(g**2)*o*ui*wi)*L31 + &
407                    (2*pi*(g**3)*(o**2)*uj+2*pi*(g**3)*(o**2)*ui)*L21*L23+(pi*(g**3)*(o**2)*wj+pi*(g**3)*(o**2)*wi)*(L21**2) + &
408                    (2*pi*(g**3)*(o**2)*uj+2*pi*(g**3)*(o**2)*ui)*L11*L13+(pi*(g**3)*(o**2)*wj+pi*(g**3)*(o**2)*wi)*(L11**2))/12
409     
410               C120 = C120 + ni*nj*(pi*gsqr*o*(4*g*o*uj*(L32**2)+4*g*o*ui*(L32**2)+8*g*o*vj*L31*L32+8*g*o*vi*L31*L32 + &
411                    12*uj*vj*L32-12*ui*vi*L32+6*(vj**2)*L31-6*(vi**2)*L31+g*o*uj*(L22**2)+g*o*ui*(L22**2) + &
412                    2*g*o*vj*L21*L22+2*g*o*vi*L21*L22+g*o*uj*(L12**2)+g*o*ui*(L12**2)+2*g*o*vj*L11*L12 + &
413                    2*g*o*vi*L11*L12))/12
414     
415               C111 = C111 + ni*nj*(pi*gsqr*o*(4*g*o*uj*L32*L33+4*g*o*ui*L32*L33+4*g*o*vj*L31*L33+4*g*o*vi*L31*L33 + &
416                    6*uj*vj*L33-6*ui*vi*L33+4*g*o*wj*L31*L32+4*g*o*wi*L31*L32+6*uj*wj*L32-6*ui*wi*L32 + &
417                    6*vj*wj*L31-6*vi*wi*L31+g*o*uj*L22*L23+g*o*ui*L22*L23+g*o*vj*L21*L23+g*o*vi*L21*L23 + &
418                    g*o*wj*L21*L22+g*o*wi*L21*L22+g*o*uj*L12*L13+g*o*ui*L12*L13+g*o*vj*L11*L13+g*o*vi*L11*L13 + &
419                    g*o*wj*L11*L12+g*o*wi*L11*L12))/12
420     
421               C102 = C102 + ni*nj*(gsqr*o*pi*(4*g*o*vj*(L33**2)+4*g*o*vi*(L33**2)+8*g*o*wj*L32*L33+8*g*o*wi*L32*L33 + &
422                    12*vj*wj*L33-12*vi*wi*L33+6*(wj**2)*L32-6*(wi**2)*L32+g*o*vj*(L23**2)+g*o*vi*(L23**2)+2*g*o*wj*L22*L23 + &
423                    2*g*o*wi*L22*L23+g*o*vj*(L13**2)+g*o*vi*(L13**2)+2*g*o*wj*L12*L13+2*g*o*wi*L12*L13))/12
424     
425               C030 = C030 + ni*nj*(pi*gsqr*o*(vj+vi)*(4*g*o*(L32**2)+6*vj*L32-6*vi*L32+g*o*(L22**2)+g*o*(L12**2)))/4
426     
427               C021 = C021 + ni*nj*(pi*gsqr*o*(8*g*o*vj*L32*L33+8*g*o*vi*L32*L33+6*(vj**2)*L33-6*(vi**2)*L33 + &
428                    4*g*o*wj*(L32**2)+4*g*o*wi*(L32**2)+12*vj*wj*L32-12*vi*wi*L32+2*g*o*vj*L22*L23+2*g*o*vi*L22*L23 + &
429                    g*o*wj*(L22**2)+g*o*wi*(L22**2)+2*g*o*vj*L12*L13+2*g*o*vi*L12*L13+g*o*wj*(L12**2)+g*o*wi*(L12**2)))/12
430     
431               C012 = C012 + ni*nj*(gsqr*o*pi*(4*g*o*vj*(L33**2)+4*g*o*vi*(L33**2)+8*g*o*wj*L32*L33+8*g*o*wi*L32*L33 + &
432                    12*vj*wj*L33-12*vi*wi*L33+6*(wj**2)*L32-6*(wi**2)*L32+g*o*vj*(L23**2)+g*o*vi*(L23**2)+2*g*o*wj*L22*L23 + &
433                    2*g*o*wi*L22*L23+g*o*vj*(L13**2)+g*o*vi*(L13**2)+2*g*o*wj*L12*L13+2*g*o*wi*L12*L13))/12
434     
435               C003 = C003 + ni*nj*(pi*gsqr*o*(wj+wi)*(4*g*o*(L33**2)+6*wj*L33-6*wi*L33+g*o*(L23**2)+g*o*(L13**2)))/4
436            END DO
437         END DO
438     
439         d = dp*dp/2 ! Scaling
440         Coll(1) = 0
441         Coll(2) = 0
442         Coll(3) = 0
443         Coll(4) = 0
444         Coll(5) = d*C200
445         Coll(6) = d*C110
446         Coll(7) = d*C101
447         Coll(8) = d*C020
448         Coll(9) = d*C011
449         Coll(10) = d*C002
450         Coll(11) = d*C300
451         Coll(12) = d*C210
452         Coll(13) = d*C201
453         Coll(14) = d*C120
454         Coll(15) = d*C111
455         Coll(16) = d*C102
456         Coll(17) = d*C030
457         Coll(18) = d*C021
458         Coll(19) = d*C012
459         Coll(20) = d*C003
460     
461       END SUBROUTINE COLLISIONS_BOLTZMANN_ONE_SPECIE
462     
463     
464     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
465     !                                                                      C
466     !  Boltzmann collision operator for bidisperse case                    C
467     !                                                                      C
468     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
469     
470       SUBROUTINE COLLISIONS_BOLTZMANN_TWO_SPECIES (N1, U1, V1, W1, N2, U2, V2, W2, m1, m2, dp1, dp2, e, Coll)
471     
472         USE constant, only: Pi
473     
474         IMPLICIT NONE
475     
476         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: N1, U1, V1, W1, N2, U2, V2, W2
477         DOUBLE PRECISION, INTENT(IN) :: m1, m2, dp1, dp2, e
478         DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NMOM) :: Coll
479     
480         DOUBLE PRECISION, DIMENSION(3) :: a1, a, vr, v1r, v2r, g1, g2
481     
482         DOUBLE PRECISION :: C100, C010, C001,  C200, C110, C101, C020, C011
483         DOUBLE PRECISION :: C002, C300, C210, C201, C120, C111, C102, C030, C021, C012, C003
484     
485         DOUBLE PRECISION :: dp, d, g, gsqr, o, mu,  ni, ui, vi, wi, nj, uj, vj, wj
486         DOUBLE PRECISION :: cross, denom, L11, L12, L13, L21, L22, L23, L31, L32, L33
487     
488         INTEGER :: i, j
489     
490         a1 = (/0.0462,0.0971,0.8235/) ! chosen at random (cannot be colinear with vr)
491         a = 1/((a1(1))**2 + (a1(2))**2 + (a1(3))**2) * a1;
492     
493         C100 = 0.D0
494         C010 = 0.D0
495         C001 = 0.D0
496         C200 = 0.D0
497         C110 = 0.D0
498         C101 = 0.D0
499         C020 = 0.D0
500         C011 = 0.D0
501         C002 = 0.D0
502         C300 = 0.D0
503         C210 = 0.D0
504         C201 = 0.D0
505         C120 = 0.D0
506         C111 = 0.D0
507         C102 = 0.D0
508         C030 = 0.D0
509         C021 = 0.D0
510         C012 = 0.D0
511         C003 = 0.D0
512     
513         ! Reduced mass
514         mu = m1*m2/(m1+m2)
515         o = -(1+e)*mu/m1
516         ! Collision diameter
517         dp = (dp1 + dp2)/2
518         ! Scaling factor
519         d = (dp**2)/m2
520     
521         DO i = 1, QMOMK_NN
522            ui = U1(i)
523            vi = V1(i)
524            wi = W1(i)
525            ni = N1(i)
526     
527            v1r(1) = ui
528            v1r(2) = vi
529            v1r(3) = wi
530     
531            DO j = 1, QMOMK_NN
532               IF (i == j) THEN
533                  CYCLE ! vi - vj = 0 always => nothing to add to source term
534               END IF
535     
536               nj = N2(j)
537               uj = U2(j)
538               vj = V2(j)
539               wj = W2(j)
540     
541               v2r(1) = uj
542               v2r(2) = vj
543               v2r(3) = wj
544     
545               vr(1) = v1r(1) - v2r(1)
546               vr(2) = v1r(2) - v2r(2)
547               vr(3) = v1r(3) - v2r(3)
548     
549               g = SQRT(vr(1)*vr(1) + vr(2)*vr(2) + vr(3)*vr(3))
550     
551               IF (g == 0) THEN ! Avoiding division by zero in the transformation matrix
552                  CYCLE
553               END IF
554     
555               g1(1) = vr(1)/g
556               g1(2) = vr(2)/g
557               g1(3) = vr(3)/g
558     
559               cross = a(1)*g1(1) + a(2)*g1(2) + a(3)*g1(3)
560     
561               g2(1) = a(1) - cross*g1(1)
562               g2(2) = a(2) - cross*g1(2)
563               g2(3) = a(3) - cross*g1(3)
564     
565               denom = SQRT( g2(1)*g2(1) + g2(2)*g2(2) + g2(3)*g2(3))
566     
567               L11 = g2(1)/denom
568               L12 = g2(2)/denom
569               L13 = g2(3)/denom
570               L21 = ( g1(3)*a(2)-g1(2)*a(3) )/denom
571               L22 = ( g1(1)*a(3)-g1(3)*a(1) )/denom
572               L23 = ( g1(2)*a(1)-g1(1)*a(2) )/denom
573               L31 = g1(1)
574               L32 = g1(2)
575               L33 = g1(3)
576     
577               gsqr = g**2
578     
579               C100 = C100 + ni*nj*(pi*gsqr*o*L31)/2
580               C010 = C010 + ni*nj*(pi*gsqr*o*L32)/2
581               C001 = C001 + ni*nj*(pi*gsqr*o*L33)/2
582     
583               C200 = C200 + ni*nj*(4*pi*(g**3)*(o**2)*(L31**2)+12*pi*gsqr*o*ui*L31 &
584                    + pi*(g**3)*(o**2)*(L21**2)+pi*(g**3)*(o**2)*(L11**2))/12
585               C110 = C110 + ni*nj*((4*pi*(g**3)*(o**2)*L31+6*pi*gsqr*o*ui)*L32+6*pi*gsqr*o*vi*L31+pi*(g**3)*(o**2)*L21*L22 + &
586                    pi*(g**3)*(o**2)*L11*L12)/12
587     
588               C101 = C101 + ni*nj*((4*pi*(g**3)*(o**2)*L31+6*pi*gsqr*o*ui)*L33+6*pi*gsqr*o*wi*L31+pi*(g**3)*(o**2)*L21*L23 + &
589                    pi*(g**3)*(o**2)*L11*L13)/12
590     
591               C020 = C020 + ni*nj*(4*pi*(g**3)*(o**2)*(L32**2)+12*pi*gsqr*o*vi*L32 &
592                    + pi*(g**3)*(o**2)*(L22**2)+pi*(g**3)*(o**2)*(L12**2))/12
593     
594               C011 = C011 + ni*nj*((4*pi*(g**3)*(o**2)*L32+6*pi*gsqr*o*vi)*L33+6*pi*gsqr*o*wi*L32+pi*(g**3)*(o**2)*L22*L23 + &
595                    pi*(g**3)*(o**2)*L12*L13)/12
596     
597               C002 = C002 + ni*nj*(4*pi*(g**3)*(o**2)*(L33**2)+12*pi*gsqr*o*wi*L33 &
598                    + pi*(g**3)*(o**2)*(L23**2)+pi*(g**3)*(o**2)*(L13**2))/12
599     
600               C300 = C300 + ni*nj*(2*pi*(g**4)*(o**3)*(L31**3)+8*pi*(g**3)*(o**2)*ui*(L31**2) &
601                    + (pi*(g**4)*(o**3)*(L21**2)+pi*(g**4)*(o**3)*(L11**2) + 12*pi*gsqr*o*(ui**2))*L31 &
602                    + 2*pi*(g**3)*(o**2)*ui*(L21**2)+2*pi*(g**3)*(o**2)*ui*(L11**2))/8
603     
604               C210 = C210 + ni*nj*((6*pi*(g**4)*(o**3)*(L31**2)+16*pi*(g**3)*(o**2)*ui*L31 &
605                    + pi*(g**4)*(o**3)*(L21**2)+pi*(g**4)*(o**3)*(L11**2) + 12*pi*gsqr*o*(ui**2))*L32 &
606                    + 8*pi*(g**3)*(o**2)*vi*(L31**2)+(2*pi*(g**4)*(o**3)*L21*L22+2*pi*(g**4)*(o**3)*L11*L12 &
607                    + 24*pi*(g**2)*o*ui*vi)*L31+4*pi*(g**3)*(o**2)*ui*L21*L22+2*pi*(g**3)*(o**2)*vi*(L21**2) &
608                    + 4*pi*(g**3)*(o**2)*ui*L11*L12+2*pi*(g**3)*(o**2)*vi*(L11**2))/(24)
609     
610               C201 = C201 + ni*nj*((6*pi*(g**4)*(o**3)*(L31**2)+16*pi*(g**3)*(o**2)*ui*L31 &
611                    + pi*(g**4)*(o**3)*(L21**2)+pi*(g**4)*(o**3)*(L11**2) + 12*pi*gsqr*o*(ui**2))*L33 &
612                    + 8*pi*(g**3)*(o**2)*wi*(L31**2)+(2*pi*(g**4)*(o**3)*L21*L23+2*pi*(g**4)*(o**3)*L11*L13 &
613                    + 24*pi*gsqr*o*ui*wi)*L31+4*pi*(g**3)*(o**2)*ui*L21*L23+2*pi*(g**3)*(o**2)*wi*(L21**2) &
614                    + 4*pi*(g**3)*(o**2)*ui*L11*L13+2*pi*(g**3)*(o**2)*wi*(L11**2))/(24)
615     
616               C120 = C120 + ni*nj*((6*pi*(g**4)*(o**3)*L31+8*pi*(g**3)*(o**2)*ui)*(L32**2)+(16*pi*(g**3)*(o**2)*vi*L31 + &
617                    2*pi*(g**4)*(o**3)*L21*L22+2*pi*(g**4)*(o**3)*L11*L12+24*pi*gsqr*o*ui*vi)*L32+(pi*(g**4)*(o**3)*(L22**2) + &
618                    pi*(g**4)*(o**3)*(L12**2)+12*pi*gsqr*o*(vi**2))*L31+2*pi*(g**3)*(o**2)*ui*(L22**2)+4*pi*(g**3)*(o**2)*vi*L21*L22 + &
619                    2*pi*(g**3)*(o**2)*ui*(L12**2)+4*pi*(g**3)*(o**2)*vi*L11*L12)/(24)
620     
621               C111 = C111 + ni*nj*(((6*pi*(g**4)*(o**3)*L31+8*pi*(g**3)*(o**2)*ui)*L32 &
622                    + 8*pi*(g**3)*(o**2)*vi*L31+pi*(g**4)*(o**3)*L21*L22 + pi*(g**4)*(o**3)*L11*L12 &
623                    + 12*pi*gsqr*o*ui*vi)*L33+(8*pi*(g**3)*(o**2)*wi*L31+pi*(g**4)*(o**3)*L21*L23 + &
624                    pi*(g**4)*(o**3)*L11*L13+12*pi*gsqr*o*ui*wi)*L32+(pi*(g**4)*(o**3)*L22*L23+pi*(g**4)*(o**3)*L12*L13+ &
625                    12*pi*gsqr*o*vi*wi)*L31+(2*pi*(g**3)*(o**2)*ui*L22+2*pi*(g**3)*(o**2)*vi*L21)*L23 + &
626                    2*pi*(g**3)*(o**2)*wi*L21*L22+(2*pi*(g**3)*(o**2)*ui*L12+2*pi*(g**3)*(o**2)*vi*L11)*L13 + &
627                    2*pi*(g**3)*(o**2)*wi*L11*L12)/(24)
628     
629               C102 = C102 + ni*nj*((6*pi*(g**4)*(o**3)*L31+8*pi*(g**3)*(o**2)*ui)*(L33**2)+(16*pi*(g**3)*(o**2)*wi*L31 + &
630                    2*pi*(g**4)*(o**3)*L21*L23+2*pi*(g**4)*(o**3)*L11*L13+24*pi*(g**2)*o*ui*wi)*L33 + &
631                    (pi*(g**4)*(o**3)*(L23**2)+pi*(g**4)*(o**3)*(L13**2)+12*pi*(g**2)*o*(wi**2))*L31+2*pi*(g**3)*(o**2)*ui*(L23**2)+ &
632                    4*pi*(g**3)*(o**2)*wi*L21*L23+2*pi*(g**3)*(o**2)*ui*(L13**2)+4*pi*(g**3)*(o**2)*wi*L11*L13)/(24)
633     
634               C030 = C030 + ni*nj*(2*pi*(g**4)*(o**3)*(L32**3)+8*pi*(g**3)*(o**2)*vi*(L32**2) &
635                    + (pi*(g**4)*(o**3)*(L22**2)+pi*(g**4)*(o**3)*(L12**2) + 12*pi*(g**2)*o*(vi**2))*L32 &
636                    + 2*pi*(g**3)*(o**2)*vi*(L22**2)+2*pi*(g**3)*(o**2)*vi*(L12**2))/8
637     
638               C021 = C021 + ni*nj*((6*pi*(g**4)*(o**3)*(L32**2)+16*pi*(g**3)*(o**2)*vi*L32 &
639                    + pi*(g**4)*(o**3)*(L22**2)+pi*(g**4)*(o**3)*(L12**2) + 12*pi*gsqr*o*(vi**2))*L33 &
640                    + 8*pi*(g**3)*(o**2)*wi*(L32**2)+(2*pi*(g**4)*(o**3)*L22*L23+2*pi*(g**4)*(o**3)*L12*L13 &
641                    + 24*pi*(g**2)*o*vi*wi)*L32+4*pi*(g**3)*(o**2)*vi*L22*L23+2*pi*(g**3)*(o**2)*wi*(L22**2) &
642                    + 4*pi*(g**3)*(o**2)*vi*L12*L13+2*pi*(g**3)*(o**2)*wi*(L12**2))/(24)
643     
644               C012 = C012 + ni*nj*((6*pi*(g**4)*(o**3)*L32+8*pi*(g**3)*(o**2)*vi)*(L33**2)+(16*pi*(g**3)*(o**2)*wi*L32 + &
645                    2*pi*(g**4)*(o**3)*L22*L23+2*pi*(g**4)*(o**3)*L12*L13+24*pi*gsqr*o*vi*wi)*L33+(pi*(g**4)*(o**3)*(L23**2) + &
646                    pi*(g**4)*(o**3)*(L13**2)+12*pi*gsqr*o*(wi**2))*L32+2*pi*(g**3)*(o**2)*vi*(L23**2)+4*pi*(g**3)*(o**2)*wi*L22*L23 + &
647                    2*pi*(g**3)*(o**2)*vi*(L13**2)+4*pi*(g**3)*(o**2)*wi*L12*L13)/(24)
648     
649               C003 = C003 + ni*nj*(2*pi*(g**4)*(o**3)*(L33**3)+8*pi*(g**3)*(o**2)*wi*(L33**2) &
650                    + (pi*(g**4)*(o**3)*(L23**2)+pi*(g**4)*(o**3)*(L13**2) + 12*pi*gsqr*o*(wi**2))*L33 &
651                    + 2*pi*(g**3)*(o**2)*wi*(L23**2)+2*pi*(g**3)*(o**2)*wi*(L13**2))/8;
652     
653            END DO
654         END DO
655     
656         Coll(1) = 0;
657         Coll(2) = d*C100;
658         Coll(3) = d*C010;
659         Coll(4) = d*C001;
660         Coll(5) = d*C200;
661         Coll(6) = d*C110;
662         Coll(7) = d*C101;
663         Coll(8) = d*C020;
664         Coll(9) = d*C011;
665         Coll(10) = d*C002;
666         Coll(11) = d*C300;
667         Coll(12) = d*C210;
668         Coll(13) = d*C201;
669         Coll(14) = d*C120;
670         Coll(15) = d*C111;
671         Coll(16) = d*C102;
672         Coll(17) = d*C030;
673         Coll(18) = d*C021;
674         Coll(19) = d*C012;
675         Coll(20) = d*C003;
676     
677       END SUBROUTINE COLLISIONS_BOLTZMANN_TWO_SPECIES
678     
679     
680     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
681     !                                                                      C
682     !   Calculates the rate of change in the moments due to collisions     C
683     !   (1 specie)                                                         C
684     !                                                                      C
685     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
686     
687       SUBROUTINE SOLVE_BOLTZMANN_COLLISIONS_ONE_SPECIE(M, N, U, V, W, dt, &
688          e, dp, order)
689     
690         IMPLICIT NONE
691     
692         DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: M
693         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: N, U, V, W
694         DOUBLE PRECISION, INTENT(IN) ::  dt, e, dp
695         INTEGER, INTENT(IN) :: order
696     
697         DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: Coll, Mtmp
698         DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Ntmp, Utmp, Vtmp, Wtmp
699         DOUBLE PRECISION :: h
700     
701     
702         IF (order == 0) THEN      ! Euler method - First order
703            h = dt
704            CALL COLLISIONS_BOLTZMANN_ONE_SPECIE (N, U, V, W, dp, e, Coll)
705            Mtmp = M + h*Coll
706         ELSE IF (order == 1) THEN ! Runge - Kutta second order
707            h = dt ;
708            CALL COLLISIONS_BOLTZMANN_ONE_SPECIE (N, U, V, W, dp, e, Coll)
709            Mtmp = M + 0.5*h*Coll
710            CALL EIGHT_NODE_3D (Mtmp, Ntmp, Utmp, Vtmp, Wtmp) ! Projection step
711            CALL COLLISIONS_BOLTZMANN_ONE_SPECIE (Ntmp, Utmp, Vtmp, Wtmp, dp, e, Coll)
712            Mtmp = M + h*Coll
713         ELSE
714            PRINT *,'QMOMK: Wrong order in Boltzmann collisions'
715            STOP
716         END IF
717         M= Mtmp
718       END SUBROUTINE SOLVE_BOLTZMANN_COLLISIONS_ONE_SPECIE
719     
720     
721     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
722     !                                                                      C
723     !   Calculates the rate of change in the moments due to collisions     C
724     !   (2 specie)                                                         C
725     !                                                                      C
726     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
727     
728       SUBROUTINE SOLVE_BOLTZMANN_COLLISIONS_TWO_SPECIES(M1, N1, U1, V1, &
729          W1, M2, N2, U2, V2, W2, dt, m_1, m_2, dp1, dp2, e11, e12, order)
730     
731         IMPLICIT NONE
732     
733         DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: M1
734         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NMOM) ::  M2
735         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: N1, U1, V1, W1, N2, U2, V2, W2
736         DOUBLE PRECISION, INTENT(IN) ::  dt, dp1, dp2, m_1, m_2, e11, e12
737         INTEGER, INTENT(IN) :: order
738     
739         DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: Coll, Colltmp, M1tmp
740         DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: N1tmp, U1tmp, V1tmp, W1tmp
741         DOUBLE PRECISION :: h
742     
743         IF (order == 0) THEN      ! Euler method - First order
744            h = dt ;
745            !     Collisions inside the same specie
746            CALL COLLISIONS_BOLTZMANN_TWO_SPECIES (N1, U1, V1, W1, N1, U1, V1, W1, m_1, m_1, dp1, dp1, e11, Colltmp)
747            Coll = Colltmp
748            !     Collisions between the two species
749            CALL COLLISIONS_BOLTZMANN_TWO_SPECIES (N1, U1, V1, W1, N2, U2, V2, W2, m_1, m_2, dp1, dp2, e12, Colltmp)
750            Coll = Coll + Colltmp
751            M1tmp = M1 + h*Coll
752         ELSE IF (order == 1) THEN     ! Runge-Kutta second order method
753            h = dt ;
754            !     Collisions inside the same specie
755            CALL COLLISIONS_BOLTZMANN_TWO_SPECIES (N1, U1, V1, W1, N1, U1, V1, W1, m_1, m_1, dp1, dp1, e11, Colltmp)
756            Coll = Colltmp
757            !     Collisions between the two species
758            CALL COLLISIONS_BOLTZMANN_TWO_SPECIES (N1, U1, V1, W1, N2, U2, V2, W2, m_1, m_2, dp1, dp2, e12, Colltmp)
759            Coll = Coll + Colltmp
760            M1tmp = M1 + 0.5*h*Coll ;
761            CALL EIGHT_NODE_3D (M1tmp, N1tmp, U1tmp, V1tmp, W1tmp) ! Projection step
762            !     Collisions inside the same specie
763     
764            CALL COLLISIONS_BOLTZMANN_TWO_SPECIES (N1tmp, U1tmp, V1tmp, W1tmp, &
765                 N1tmp, U1tmp, V1tmp, W1tmp, m_1, m_1, dp1, dp1, e11, Colltmp)
766     
767            Coll = Colltmp
768            !     Collisions between the two species
769            CALL COLLISIONS_BOLTZMANN_TWO_SPECIES (N1tmp, U1tmp, V1tmp, W1tmp, N2, U2, V2, W2, m_1, m_2, dp1, dp2, e12, Colltmp)
770            Coll = Coll + Colltmp
771            M1tmp = M1 + h*Coll;
772         ELSE
773            PRINT *,'QMOMK: Wrong order in Boltzmann collisions'
774            STOP
775         END IF
776         M1 = M1tmp;
777       END SUBROUTINE SOLVE_BOLTZMANN_COLLISIONS_TWO_SPECIES
778     
779     
780     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
781     !                                                                      C
782     !                                                                     C
783     !                                                                      C
784     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
785     
786       DOUBLE PRECISION FUNCTION RADIAL_G0(ALPHA)
787     
788         USE param1, only: small_number, large_number
789         USE constant, only: Pi, EP_STAR
790     
791         IMPLICIT NONE
792     
793         DOUBLE PRECISION, INTENT(IN) :: ALPHA
794     
795         DOUBLE PRECISION ALPHA_MAX
796     
797         ALPHA_MAX = 1. - EP_STAR
798     
799         RADIAL_G0 = 1.
800         IF (ALPHA < ALPHA_MAX - SMALL_NUMBER) THEN
801           RADIAL_G0 = 1./(1.0-ALPHA) + 3.0*ALPHA/(2.*(1.-ALPHA)**2) + (ALPHA**2)/(2*(1.0-ALPHA)**3)
802         ELSE
803           RADIAL_G0 = LARGE_NUMBER
804         END IF
805     
806       END FUNCTION RADIAL_G0
807     
808     END MODULE qmomk_collision
809