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