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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module: qmomk_quadrature                                            C
4     !  Purpose: Routines to perform the quadrature calculations            C
5     !  Contains the following subroutines:                                 C
6     !     quadrature_bounded, moments_twenty_eight_nodes,                  C
7     !     check_moments_twenty, eight_node_3d, bind_theta                  C
8     !                                                                      C
9     !  Author: A. Passalacqua                      Date: 18-Jun-2008       C
10     !  Reviewer:                                Date: dd-mmm-yy            C
11     !                                                                      C
12     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
13     MODULE qmomk_quadrature
14     
15       USE qmomk_tools
16       USE qmomk_parameters
17     
18       IMPLICIT NONE
19     
20       PRIVATE
21     
22       PUBLIC :: MOMENTS_TWENTY_EIGHT_NODES
23       PUBLIC :: EIGHT_NODE_3D
24       PUBLIC :: BIND_THETA
25     
26     CONTAINS
27     
28     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
29     !                                                                      C
30     !  Two-nodes quadrature                                                C
31     !                                                                      C
32     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
33     
34       SUBROUTINE QUADRATURE_BOUNDED (m, w, u)
35     
36         IMPLICIT NONE
37     
38         DOUBLE PRECISION, INTENT(IN), DIMENSION(4) :: m
39         DOUBLE PRECISION, INTENT(OUT), DIMENSION(2) :: w
40         DOUBLE PRECISION, INTENT(OUT), DIMENSION(2) :: u
41     
42         DOUBLE PRECISION :: eps = 1.D-8
43         DOUBLE PRECISION :: smax = 1.D6
44     
45         DOUBLE PRECISION :: rho = 0.D0
46         DOUBLE PRECISION :: um, sigma, x, q, s
47     
48         rho = m(1)
49         um = m(2)/rho
50     
51         sigma = (m(3)*rho - m(2)*m(2))/(rho**2)
52     
53         IF (sigma <= 0.D0) THEN
54            PRINT *,'QMOMK: Negative variance in quadrature'
55            sigma = 0.D0
56            x = 0
57         ELSE
58            sigma = SQRT(sigma)
59            q = m(4)/rho - um**3 - 3.D0*um*(sigma**2)
60            s = MIN(smax, ABS(q)/(sigma**3))
61            x = 0.D0
62            IF ((s**2) < eps) THEN
63               x = 0.D0
64            ELSE
65               x = 0.5 * SIGN(1.D0, q)/SQRT(1.D0 + 4.D0/(s**2))
66            END IF
67         END IF
68     
69         w(1) = rho*(0.5 + x)
70         u(1) = um - SQRT((0.5-x)/(0.5+x))*sigma
71     
72         w(2) = rho*(0.5-x)
73         u(2) = um + SQRT((0.5+x)/(0.5-x))*sigma
74     
75         IF (ABS(x) >= 0.5D0) THEN
76            PRINT *,'QMOMK: Warning: x > 0.5'
77            PRINT *,x
78         END IF
79       END SUBROUTINE QUADRATURE_BOUNDED
80     
81     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
82     !                                                                      C
83     !  Calculation of the twenty moments from weights and abscissas        C
84     !                                                                      C
85     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
86     
87       SUBROUTINE MOMENTS_TWENTY_EIGHT_NODES (n, u, v, w, mom)
88     
89         USE param1, only: small_number
90         IMPLICIT NONE
91     
92         INTEGER :: i
93     
94         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: n
95         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: u
96         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: v
97         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: w
98         DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NMOM) :: mom
99     
100         DOUBLE PRECISION, dimension(QMOMK_NMOM) :: M
101     
102         M=0.D0
103     
104         DO i = 1, QMOMK_NN
105            M(1) = M(1) + n(i)
106            M(2) = M(2) + n(i)*u(i)
107            M(3) = M(3) + n(i)*v(i)
108            M(4) = M(4) + n(i)*w(i)
109            M(5) = M(5) + n(i)*(u(i)**2)
110            M(6) = M(6) + n(i)*u(i)*v(i)
111            M(7) = M(7) + n(i)*u(i)*w(i)
112            M(8) = M(8) + n(i)*(v(i)**2)
113            M(9) = M(9) + n(i)*v(i)*w(i)
114            M(10)= M(10)+ n(i)*w(i)**2
115            M(11)= M(11)+ n(i)*(u(i)**3)
116            M(12)= M(12)+ n(i)*(u(i)**2)*v(i)
117            M(13)= M(13)+ n(i)*(u(i)**2)*w(i)
118            M(14)= M(14)+ n(i)*u(i)*(v(i)**2)
119            M(15)= M(15)+ n(i)*u(i)*v(i)*w(i)
120            M(16)= M(16)+ n(i)*u(i)*(w(i)**2)
121            M(17)= M(17)+ n(i)*(v(i)**3)
122            M(18)= M(18)+ n(i)*(v(i)**2)*w(i)
123            M(19)= M(19)+ n(i)*v(i)*(w(i)**2)
124            M(20)= M(20)+ n(i)*(w(i)**3)
125         END DO
126     
127         IF (M(1) < SMALL_NUMBER) THEN
128           PRINT *,'QMOMK: Zero order moment is very small!'
129         END IF
130     
131         mom = M
132       END SUBROUTINE MOMENTS_TWENTY_EIGHT_NODES
133     
134     
135     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
136     !                                                                      C
137     !  Moments check                                                       C
138     !                                                                      C
139     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
140     
141       SUBROUTINE CHECK_MOMENTS_TWENTY(mom, n, u, v, w, Check)
142     
143         IMPLICIT NONE
144     
145         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NMOM) :: mom
146         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: n
147         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: u
148         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: v
149         DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: w
150         DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NMOM) :: Check
151     
152         DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: F
153     
154         INTEGER :: i
155     
156         F = 0.D0
157     
158         F(1) = - mom(1)
159         F(2) = - mom(2)
160         F(3) = - mom(3)
161         F(4) = - mom(4)
162         F(5) = - mom(5)
163         F(6) = - mom(6)
164         F(7) = - mom(7)
165         F(8) = - mom(8)
166         F(9) = - mom(9)
167         F(10)= - mom(10)
168         F(11)= - mom(11)
169         F(12)= - mom(12)
170         F(13)= - mom(13)
171         F(14)= - mom(14)
172         F(15)= - mom(15)
173         F(16)= - mom(16)
174         F(17)= - mom(17)
175         F(18)= - mom(18)
176         F(19)= - mom(19)
177         F(20)= - mom(20)
178     
179         DO i = 1, QMOMK_NN
180            F(1) = F(1) + n(i)
181            F(2) = F(2) + n(i)*u(i)
182            F(3) = F(3) + n(i)*v(i)
183            F(4) = F(4) + n(i)*w(i)
184            F(5) = F(5) + n(i)*(u(i)**2)
185            F(6) = F(6) + n(i)*u(i)*v(i)
186            F(7) = F(7) + n(i)*u(i)*w(i)
187            F(8) = F(8) + n(i)*(v(i)**2)
188            F(9) = F(9) + n(i)*v(i)*w(i)
189            F(10)= F(10)+ n(i)*(w(i)**2)
190            F(11)= F(11)+ n(i)*(u(i)**3)
191            F(12)= F(12)+ n(i)*(u(i)**2)*v(i)
192            F(13)= F(13)+ n(i)*(u(i)**2)*w(i)
193            F(14)= F(14)+ n(i)*u(i)*(v(i)**2)
194            F(15)= F(15)+ n(i)*u(i)*v(i)*w(i)
195            F(16)= F(16)+ n(i)*u(i)*(w(i)**2)
196            F(17)= F(17)+ n(i)*(v(i)**3)
197            F(18)= F(18)+ n(i)*(v(i)**2)*w(i)
198            F(19)= F(19)+ n(i)*v(i)*(w(i)**2)
199            F(20)= F(20)+ n(i)*(w(i)**3)
200         END DO
201     
202         DO i = 1, QMOMK_NMOM
203            F(i) = ABS(F(i))
204         END DO
205     
206         Check = F
207       END SUBROUTINE CHECK_MOMENTS_TWENTY
208     
209     
210     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
211     !                                                                      C
212     !  Calculation of weights and abscissas                                C
213     !                                                                      C
214     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
215     
216       SUBROUTINE EIGHT_NODE_3D (mom, n, u, v, w)
217     
218         IMPLICIT NONE
219     
220         DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: mom
221         DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NN) :: n
222         DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NN) :: u
223         DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NN) :: v
224         DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NN) :: w
225     
226         DOUBLE PRECISION :: um, vm, wm, s1, s2, s3, s12, s13, s23
227         DOUBLE PRECISION :: q000, q101, q110, q011, q200, q020, q002, q210, q300
228         DOUBLE PRECISION :: q201, q120, q111, q102, q030, q021, q012, q003
229     
230         DOUBLE PRECISION :: b11, b12, b13, b21, b22, b23, b31, b32, b33, Fmax
231     
232         DOUBLE PRECISION :: m000, m100, m010, m001, m200, m020, m002, m300, m030, m003, m111
233     
234         DOUBLE PRECISION, DIMENSION (QMOMK_NCOV, QMOMK_NCOV) :: sig
235         DOUBLE PRECISION, DIMENSION (QMOMK_NCOV, QMOMK_NCOV) :: chol
236         DOUBLE PRECISION, DIMENSION (QMOMK_NCOV, QMOMK_NCOV) :: scal, invScale
237         DOUBLE PRECISION, DIMENSION (QMOMK_NCOV, QMOMK_NCOV) :: Umat, Vmat
238     
239         DOUBLE PRECISION, DIMENSION (4) :: alpha1, alpha2, alpha3
240         DOUBLE PRECISION, DIMENSION (2) :: aq1, aq2, aq3, wq1, wq2, wq3
241     
242         DOUBLE PRECISION :: a, b, c, u1, u2, v1, v2, w1, w2, rho123
243     
244         DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: nstar, ones
245     
246         DOUBLE PRECISION, DIMENSION(32) :: x
247     
248         DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: F
249     
250         q000 = mom(1)
251         um   = mom(2)/q000
252         vm   = mom(3)/q000
253         wm   = mom(4)/q000
254         q200 = mom(5)/q000
255         q110 = mom(6)/q000
256         q101 = mom(7)/q000
257         q020 = mom(8)/q000
258         q011 = mom(9)/q000
259         q002 = mom(10)/q000
260         q300 = mom(11)/q000
261         q210 = mom(12)/q000
262         q201 = mom(13)/q000
263         q120 = mom(14)/q000
264         q111 = mom(15)/q000
265         q102 = mom(16)/q000
266         q030 = mom(17)/q000
267         q021 = mom(18)/q000
268         q012 = mom(19)/q000
269         q003 = mom(20)/q000
270     
271         sig = 0.D0
272     
273         sig(1,1) = q200 - um**2
274         sig(2,2) = q020 - vm**2
275         sig(3,3) = q002 - wm**2
276     
277         IF ( sig(1,1) < 0.D0 ) THEN ! check that u variance is nonnegative
278            PRINT *,'QMOMK: Negative u variance',sig
279            sig(1,1) = 0.D0
280            STOP
281         END IF
282     
283         IF ( sig(2,2) < 0.D0 ) THEN ! check that v variance is nonnegative
284            PRINT *,'QMOMK: Negative v variance',sig
285            sig(2,2) = 0.D0
286            STOP
287         END IF
288     
289         IF ( sig(3,3) < 0.D0 ) THEN ! check that v variance is nonnegative
290            PRINT *,'QMOMK: Negative w variance',sig
291            sig(3,3) = 0.D0
292            STOP
293         END IF
294     
295         s1 = SQRT(sig(1,1))
296         s2 = SQRT(sig(2,2))
297         s3 = SQRT(sig(3,3))
298     
299         sig(1,2) = q110 - um*vm
300         sig(2,1) = sig(1,2)
301         s12 = sig(1,2)
302     
303         sig(1,3) = q101 - um*wm
304         sig(3,1) = sig(1,3)
305         s13 = sig(1,3)
306     
307         sig(2,3) = q011 - vm*wm
308         sig(3,2) = sig(2,3)
309         s23 = sig(2,3)
310     
311         IF ( s1*s2 < ABS(s12) ) THEN ! check that covariance matrix is nonnegative
312            PRINT *,'QMOMK: unphysical uv correlation', sig(1,2), s1, s2
313            sig(1,2) = SIGN(1.D0, s12)*s1*s2
314            sig(2,1) = sig(1,2)
315     
316            !q110 = 0.99*q110
317            !mom(6) = q000*q110
318            !sig(1,2) = q110 - um*vm
319            !sig(2,1) = sig(1,2)
320            PRINT *,'Sigma = ',sig
321            PRINT *,'M = ', mom
322            PRINT *,'N = ', n
323            PRINT *,'U = ', u
324            PRINT *,'V = ', v
325            PRINT *,'W = ', w
326            STOP
327         END IF
328     
329         IF ( s1*s3 < ABS(s13) ) THEN ! check that covariance matrix is nonnegative
330            sig(1,3) = SIGN(1.D0, s13)*s1*s3
331            sig(3,1) = sig(1,3)
332            PRINT *,'QMOMK: unphysical uw correlation',sig(1,3)
333            STOP
334         END IF
335     
336         IF ( s2*s3 < ABS(s23) ) THEN ! check that covariance matrix is nonnegative
337            sig(2,3) = SIGN(1.D0, s23)*s2*s3
338            sig(3,2) = sig(2,3)
339            PRINT *,'QMOMK: unphysical vw correlation'
340            STOP
341         END IF
342     
343         CALL CHOLESKY3(sig, chol) ! Bottom cholesky decomposition
344         CALL DIAG3(s1, s2, s3, scal)
345         CALL INV3(scal, invScale)
346         CALL MULTIPLYMATRIX3(chol, invScale, Umat)
347         CALL INV3(Umat, Vmat)
348     
349         b11 = Vmat(1,1)
350         b12 = Vmat(1,2)
351         b13 = Vmat(1,3)
352     
353         b21 = Vmat(2,1)
354         b22 = Vmat(2,2)
355         b23 = Vmat(2,3)
356     
357         b31 = Vmat(3,1)
358         b32 = Vmat(3,2)
359         b33 = Vmat(3,3)
360     
361         m000 = 1
362         m100 = 0
363         m010 = 0
364         m001 = 0
365     
366         m200 = 2.D0*b11*b12*q110 - 2.D0*b11*um*b12*vm + 2.D0*b11*b13*q101 - &
367              2.D0*b11*um*b13*wm + 2.D0*b12*b13*q011 - 2.D0*b12*vm*b13*wm - &
368              (b11**2)*(um**2) - (b12**2)*(vm**2) - (b13**2)*(wm**2) + &
369              (b11**2)*q200 + (b12**2)*q020 + (b13**2)*q002
370     
371         m020 = - (b21**2)*(um**2) - (b22**2)*(vm**2) - (b23**2)*(wm**2) + &
372              2.D0*b21*b23*q101 - 2.D0*b21*um*b23*wm + 2.D0*b22*b23*q011 - &
373              2.D0*b22*vm*b23*wm - 2.D0*b21*um*b22*vm + 2.D0*b21*b22*q110 + &
374              (b21**2)*q200 + (b22**2)*q020 + (b23**2)*q002
375     
376         m002 = - (b33**2)*(wm**2) - (b31**2)*(um**2) - (b32**2)*(vm**2) - &
377              2.D0*b31*um*b32*vm + 2.D0*b31*b32*q110 + (b31**2)*q200 + &
378              (b32**2)*q020 + (b33**2)*q002 + 2.D0*b31*b33*q101 - 2.D0*b31*um*b33*wm + &
379              2.D0*b32*b33*q011 - 2.D0*b32*vm*b33*wm
380     
381         m300 = - 3.D0*(b13**3)*wm*q002 + 3.D0*b11*(b12**2)*q120 + 3.D0*(b12**2)*b13*q021 + &
382              3.D0*(b11**2)*b12*q210 + 3.D0*(b11**2)*b13*q201 + 3.D0*b11*(b13**2)*q102 + &
383              3.D0*b12*(b13**2)*q012 + 6.D0*b12*vm*(b13**2)*(wm**2) - 3.D0*(b12**3)*vm*q020 - &
384              3.D0*(b11**3)*um*q200 - 6.D0*b11*b12*b13*wm*q110 - 6.D0*b11*b12*vm*b13*q101 + &
385              (b11**3)*q300 + 2.D0*(b11**3)*(um**3) - 6.D0*b11*um*b12*b13*q011 + 2.D0*(b12**3)*(vm**3) + &
386              2.D0*(b13**3)*(wm**3) + (b13**3)*q003 + (b12**3)*q030 - 3.D0*b11*um*(b12**2)*q020 - &
387              3.D0*b11*um*(b13**2)*q002 - 3.D0*(b12**2)*b13*wm*q020 + 6.D0*(b11**2)*(um**2)*b12*vm + &
388              6.D0*(b12**2)*(vm**2)*b13*wm + 6.D0*(b11**2)*(um**2)*b13*wm + 6.D0*b11*um*(b12**2)*(vm**2) + &
389              6.D0*b11*um*(b13**2)*(wm**2) - 6.D0*b11*(b12**2)*vm*q110 - 6.D0*b11*(b13**2)*wm*q101 - &
390              6.D0*(b12**2)*vm*b13*q011 + 12.D0*b11*um*b12*vm*b13*wm - 6.D0*b12*(b13**2)*wm*q011 - &
391              6.D0*(b11**2)*um*b13*q101 + 6.D0*b11*b12*b13*q111 - 3.D0*b12*vm*(b13**2)*q002 - &
392              6.D0*(b11**2)*um*b12*q110 - 3.D0*(b11**2)*b13*wm*q200 - 3.D0*(b11**2)*b12*vm*q200
393     
394         m030 = 6.D0*(b21**2)*(um**2)*b22*vm + 6.D0*(b22**2)*(vm**2)*b23*wm - 3.D0*(b21**3)*um*q200 + &
395              2.D0*(b22**3)*(vm**3) + 3.D0*(b21**2)*b23*q201 + 2.D0*(b21**3)*(um**3) - 3.D0*(b23**3)*wm*q002 - &
396              3.D0*(b21**2)*b22*vm*q200 + 6.D0*b21*um*(b22**2)*(vm**2) + 3.D0*b21*(b22**2)*q120 + &
397              3.D0*(b21**2)*b22*q210 + 6.D0*(b21**2)*(um**2)*b23*wm - 3.D0*b21*um*(b22**2)*q020 + &
398              3.D0*b21*(b23**2)*q102 + 3.D0*(b22**2)*b23*q021 + 6.D0*b22*vm*(b23**2)*(wm**2) - &
399              3.D0*(b21**2)*b23*wm*q200 + (b22**3)*q030 + (b21**3)*q300 - 3.D0*(b22**3)*vm*q020 + &
400              6.D0*b21*b22*b23*q111 - 6.D0*(b21**2)*um*b22*q110 + (b23**3)*q003 - 3.D0*b22*vm*(b23**2)*q002 - &
401              6.D0*(b21**2)*um*b23*q101 + 3.D0*b22*(b23**2)*q012 - 6.D0*(b22**2)*vm*b23*q011 - &
402              3.D0*b21*um*(b23**2)*q002 - 3.D0*(b22**2)*b23*wm*q020 - 6.D0*b22*(b23**2)*wm*q011 - &
403              6.D0*b21*b22*b23*wm*q110 - 6.D0*b21*um*b22*b23*q011 + 6.D0*b21*um*(b23**2)*(wm**2) + &
404              2.D0*(b23**3)*(wm**3) - 6.D0*b21*(b23**2)*wm*q101 - 6.D0*b21*b22*vm*b23*q101 + &
405              12.D0*b21*um*b22*vm*b23*wm - 6.D0*b21*(b22**2)*vm*q110
406     
407         m003 = (b31**3)*q300 + 6.D0*b31*um*(b33**2)*(wm**2) + (b32**3)*q030 + (b33**3)*q003 - &
408              3.D0*b31*um*(b33**2)*q002  + 6.D0*b32*vm*(b33**2)*(wm**2) + 6.D0*b31*um*(b32**2)*(vm**2) + &
409              2.D0*(b31**3)*(um**3) + 2.D0*(b32**3)*(vm**3) + 3.D0*(b31**2)*b33*q201 - &
410              3.D0*b32*vm*(b33**2)*q002 - 6.D0*b31*(b32**2)*vm*q110 - 3.D0*(b31**2)*b32*vm*q200 - &
411              3.D0*(b31**2)*b33*wm*q200 - 6.D0*b31*(b33**2)*wm*q101 - 3.D0*(b33**3)*wm*q002 + &
412              6.D0*(b31**2)*(um**2)*b33*wm + 3.D0*b32*(b33**2)*q012 + 6.D0*b31*b32*b33*q111 - &
413              6.D0*(b31**2)*um*b33*q101 - 6.D0*b32*(b33**2)*wm*q011 + 6.D0*(b31**2)*(um**2)*b32*vm - &
414              3.D0*b31*um*(b32**2)*q020 - 6.D0*b31*b32*b33*wm*q110 + 3.D0*b31*(b33**2)*q102 - &
415              3.D0*(b32**2)*b33*wm*q020 + 12.D0*b31*um*b32*vm*b33*wm + 3.D0*(b31**2)*b32*q210 - &
416              6.D0*(b32**2)*vm*b33*q011 + 3.D0*b31*(b32**2)*q120 + 6.D0*(b32**2)*(vm**2)*b33*wm + &
417              3.D0*(b32**2)*b33*q021 - 3.D0*(b31**3)*um*q200 - 3.D0*(b32**3)*vm*q020 - &
418              6.D0*b31*um*b32*b33*q011 + 2.D0*(b33**3)*(wm**3) - 6.D0*b31*b32*vm*b33*q101 - &
419              6.D0*(b31**2)*um*b32*q110
420     
421         m111 = - b12*b23*wm*b31*q110 + 2.D0*b11*(um**2)*b22*vm*b31 - b12*b21*b33*wm*q110 - b12*vm*b23*b33*q002 + &
422              2.D0*b12*vm*b23*wm*b31*um - b12*vm*b21*b31*q200 - b13*b23*b31*um*q002 - b12*b22*b31*um*q020 + &
423              b12*b22*b31*q120 - b11*b21*b33*wm*q200 + 2.D0*b12*(vm**2)*b23*wm*b32 - b12*b23*wm*b32*q020 - &
424              b11*b21*b32*vm*q200 - b12*b21*um*b32*q020 - 2.D0*b12*b22*vm*b33*q011 - b11*um*b23*b33*q002 - &
425              2.D0*b12*b22*vm*b31*q110 - 3.D0*b11*b21*b31*um*q200 + b11*b21*b32*q210 + b11*b21*b31*q300 - &
426              2.D0*b12*b23*b33*wm*q011 + b13*b23*b33*q003 + b12*b23*b33*q012 + b12*b22*b33*q021 - &
427              b12*b23*b31*um*q011 - b12*b22*b33*wm*q020 + b13*b23*b31*q102 - 2.D0*b12*b21*b31*um*q110 + &
428              b12*b21*b32*q120 + b12*b23*b31*q111 - 2.D0*b12*b21*b32*vm*q110 + 2.D0*b12*(vm**2)*b21*um*b32 - &
429              3.D0*b12*b22*b32*vm*q020 + 2.D0*b12*(vm**2)*b22*b31*um + 2.D0*b12*vm*b21*(um**2)*b31 - &
430              2.D0*b12*b23*b32*vm*q011 - b12*b21*um*b33*q011 + 2.D0*b12*(vm**2)*b22*b33*wm + b12*b23*b32*q021 + &
431              b11*b21*b33*q201 - 3.D0*b13*b23*b33*wm*q002 - 2.D0*b11*b21*um*b32*q110 + 2.D0*b12*vm*b23*(wm**2)*b33 - &
432              2.D0*b13*b21*b33*wm*q101 + 2.D0*b13*(wm**2)*b22*vm*b33 - b13*b21*um*b32*q011 - b12*vm*b21*b33*q101 - &
433              2.D0*b11*b21*um*b33*q101 - b13*b22*vm*b33*q002 - b13*b21*um*b33*q002 + 2.D0*b11*(um**3)*b21*b31 + &
434              2.D0*b11*b21*(um**2)*b32*vm + b13*b21*b31*q201 + b13*b22*b33*q012 + b13*b22*b31*q111 + &
435              2.D0*b13*b23*(wm**2)*b31*um - b13*b23*b32*vm*q002 + b11*b23*b31*q201 - b11*b22*vm*b31*q200 + &
436              2.D0*b13*wm*b21*(um**2)*b31 - b13*b22*vm*b31*q101 + 2.D0*b11*b21*(um**2)*b33*wm + b13*b23*b32*q012 - &
437              2.D0*b11*b22*b31*um*q110 - 2.D0*b13*b22*b32*vm*q011 + b12*b21*b31*q210 - b13*wm*b22*b32*q020 - &
438              2.D0*b13*b23*wm*b31*q101 + b11*b22*b33*q111 + 2.D0*b13*(wm**2)*b21*um*b33 - b13*b21*b32*vm*q101 + &
439              b12*b22*b32*q030 + 2.D0*b13*b23*(wm**2)*b32*vm - b11*b22*b33*wm*q110 - b12*vm*b23*b31*q101 - &
440              b13*wm*b21*b32*q110 + 2.D0*b11*um*b23*wm*b32*vm + 2.D0*b11*um*b22*(vm**2)*b32 + 2*b13*(wm**3)*b23*b33 + &
441              2.D0*b12*(vm**3)*b22*b32 - 2.D0*b11*b23*b31*um*q101 - b11*b22*vm*b33*q101 - b11*um*b22*b32*q020 + &
442              2.D0*b13*wm*b22*(vm**2)*b32 + b13*b21*b33*q102 + b13*b21*b32*q111 - 2.D0*b13*b21*b31*um*q101 - &
443              b13*b22*b31*um*q011 + b11*b22*b32*q120 - b13*wm*b22*b31*q110 + 2.D0*b11*um*b23*(wm**2)*b33 + &
444              2.D0*b11*(um**2)*b23*wm*b31 + b11*b23*b32*q111 + b13*b22*b32*q021 + 2.D0*b12*vm*b21*um*b33*wm + &
445              2.D0*b13*wm*b22*vm*b31*um - 2.D0*b13*b23*wm*b32*q011 + b12*b21*b33*q111 - 2.D0*b11*b23*b33*wm*q101 + &
446              2.D0*b11*um*b22*vm*b33*wm - b11*um*b23*b32*q011 - b11*um*b22*b33*q011 - 2.D0*b11*b22*b32*vm*q110 + &
447              b11*b22*b31*q210 - b13*wm*b21*b31*q200 - 2.D0*b13*b22*b33*wm*q011 - b11*b23*wm*b32*q110 + &
448              b11*b23*b33*q102 - b11*b23*wm*b31*q200 + 2.D0*b13*wm*b21*um*b32*vm - b11*b23*b32*vm*q101
449     
450         alpha1(1) = m000
451         alpha1(2) = m100
452         alpha1(3) = m200
453         alpha1(4) = m300
454     
455         CALL QUADRATURE_BOUNDED(alpha1, wq1, aq1)
456     
457         alpha2(1) = m000
458         alpha2(2) = m010
459         alpha2(3) = m020
460         alpha2(4) = m030
461     
462         call QUADRATURE_BOUNDED(alpha2, wq2, aq2)
463     
464         alpha3(1) = m000
465         alpha3(2) = m001
466         alpha3(3) = m002
467         alpha3(4) = m003
468     
469         call QUADRATURE_BOUNDED(alpha3, wq3, aq3)
470     
471         a = wq1(1)
472         b = wq2(1)
473         c = wq3(1)
474     
475         u1 = aq1(1)
476         u2 = aq1(2)
477     
478         v1 = aq2(1)
479         v2 = aq2(2)
480     
481         w1 = aq3(1)
482         w2 = aq3(2)
483     
484         rho123 = SQRT(a*(1-a)*b*(1-b)*c*(1-c))*m111/SQRT(m200*m020*m002)
485     
486         IF (rho123 > 0.D0) THEN
487            rho123 = MIN(rho123 , a*b*c , (1-a)*(1-b)*c , (1-a)*b*(1-c) , a*(1-b)*(1-c))
488         ELSE
489            rho123 = -MIN(-rho123 , (1-a)*b*c , a*(1-b)*c , a*b*(1-c) , (1-a)*(1-b)*(1-c))
490         END IF
491     
492         IF (a > 1.) THEN
493            PRINT *,a
494         END IF
495         IF (b > 1.) THEN
496            PRINT *,b
497         END IF
498         IF (c > 1.) THEN
499            PRINT *,c
500         END IF
501     
502         nstar(1) = a*b*c                    - rho123
503         nstar(2) = (1.D0-a)*b*c             + rho123
504         nstar(3) = a*(1.D0-b)*c             + rho123
505         nstar(4) = (1.D0-a)*(1.D0-b)*c      - rho123
506         nstar(5) = a*b*(1.D0-c)             + rho123
507         nstar(6) = (1.D0-a)*b*(1.D0-c)      - rho123
508         nstar(7) = a*(1.D0-b)*(1.D0-c)      - rho123
509         nstar(8) = (1.D0-a)*(1.D0-b)*(1.D0-c)   + rho123
510     
511         IF (MINVAL(nstar) < -1.D-16) THEN
512            PRINT *,'QMOMK: Negative weight in quadrature'
513            PRINT *,nstar
514            nstar = ABS(nstar)
515            STOP
516         END IF
517     
518         x(1:8) = nstar
519         x(9:16) = (/u1, u2, u1, u2, u1, u2, u1, u2/)
520         x(17:24) = (/v1, v1, v2, v2, v1, v1, v2, v2/)
521         x(25:32) = (/w1, w1, w1, w1, w2, w2, w2, w2/)
522     
523         ones = (/1.,1.,1.,1.,1.,1.,1.,1./)
524     
525         n = x(1:8)*q000
526         u = Umat(1,1)*x(9:16) + Umat(1,2)*x(17:24) + Umat(1,3)*x(25:32) + um*ones(1:8)
527         v = Umat(2,1)*x(9:16) + Umat(2,2)*x(17:24) + Umat(2,3)*x(25:32) + vm*ones(1:8)
528         w = Umat(3,1)*x(9:16) + Umat(3,2)*x(17:24) + Umat(3,3)*x(25:32) + wm*ones(1:8)
529     
530         CALL CHECK_MOMENTS_TWENTY(mom, n, u, v, w, F)
531         Fmax = MAXVAL(F(1:10))
532     
533         IF (Fmax > 1.D-10) THEN
534            PRINT *,'QMOMK: First 10 moments not satisfied.'
535            PRINT *,'F(1:10) = ',F(1:10)
536            PRINT *,'M(1:10) = ',mom(1:10)
537            STOP
538         ENDIF
539       END SUBROUTINE EIGHT_NODE_3D
540     
541     
542     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
543     !                                                                      C
544     !  Temperature limiter                                                 C
545     !                                                                      C
546     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
547     
548       SUBROUTINE BIND_THETA (mom, tau_p)
549     
550         DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: mom
551         DOUBLE PRECISION, INTENT(IN) :: tau_p
552     
553         DOUBLE PRECISION :: sigma2
554     
555         sigma2 = MIN(MINIMUM_THETA/tau_p, MAXIMUM_SIGMA)
556     
557         ! Second order moments
558         mom(5) = mom(5) + 2.0*sigma2*mom(1)
559         mom(8) = mom(8) + 2.0*sigma2*mom(1)
560         mom(10) = mom(10) + 2.0*sigma2*mom(1)
561     
562         ! Third order moments
563         mom(11) = mom(11) + 6.0*sigma2*mom(2)
564         mom(12) = mom(12) + 2.0*sigma2*mom(3)
565         mom(13) = mom(13) + 2.0*sigma2*mom(4)
566         mom(14) = mom(14) + 2.0*sigma2*mom(2)
567     
568         mom(16) = mom(16) + 2.0*sigma2*mom(2)
569         mom(17) = mom(17) + 6.0*sigma2*mom(3)
570         mom(18) = mom(18) + 2.0*sigma2*mom(4)
571     
572         mom(19) = mom(19) + 2.0*sigma2*mom(3)
573         mom(20) = mom(20) + 6.0*sigma2*mom(4)
574     
575       END SUBROUTINE BIND_THETA
576     
577     END MODULE qmomk_quadrature
578