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