File: N:\mfix\model\source_rop_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_ROP_g                                            C
4     !  Purpose: Determine source terms for continuity equation.            C
5     !                                                                      C
6     !  Notes: The off-diagonal coefficients are positive. The center       C
7     !         coefficient and the source vector are negative.              C
8     !         See conv_rop_g for additional details.                       C
9     !                                                                      C
10     !  Author: M. Syamlal                                 Date: 2 -JUL-96  C
11     !  Reviewer:                                          Date:            C
12     !                                                                      C
13     !  Literature/Document References:                                     C
14     !  Variables referenced:                                               C
15     !  Variables modified:                                                 C
16     !  Local variables:                                                    C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     
20           SUBROUTINE SOURCE_ROP_G(A_M, B_M)
21     
22     !-----------------------------------------------
23     ! Modules
24     !-----------------------------------------------
25           USE param
26           USE param1
27           USE parallel
28           USE fldvar
29           USE rxns
30           USE run
31           USE geometry
32           USE indices
33           USE pgcor
34           USE compar
35           USE functions
36           IMPLICIT NONE
37     !-----------------------------------------------
38     ! Dummy arguments
39     !-----------------------------------------------
40     ! Septadiagonal matrix A_m
41           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
42     ! Vector b_m
43           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
44     !-----------------------------------------------
45     ! Local variables
46     !-----------------------------------------------
47     ! DEL dot V
48           DOUBLE PRECISION :: DEL_V
49     ! Mass source
50           DOUBLE PRECISION :: Src
51     ! Indices
52           INTEGER :: I, J, K, IJK, IMJK, IJMK, IJKM
53     ! error message
54           CHARACTER(LEN=80) :: LINE
55     
56     !-----------------------------------------------
57     
58     !!$omp  parallel do private( I, J, K, IJK, IMJK, IJMK, IJKM,  DEL_V, &
59     !!$omp&  Src, LINE) &
60     !!$omp&  schedule(static)
61           DO IJK = ijkstart3, ijkend3
62              IF (FLUID_AT(IJK) .AND. PHASE_4_P_G(IJK)/=0) THEN
63                 I = I_OF(IJK)
64                 J = J_OF(IJK)
65                 K = K_OF(IJK)
66                 IMJK = IM_OF(IJK)
67                 IJMK = JM_OF(IJK)
68                 IJKM = KM_OF(IJK)
69     
70                 DEL_V = U_G(IJK)*AYZ(IJK) - U_G(IMJK)*AYZ(IMJK) + &
71                         V_G(IJK)*AXZ(IJK) - V_G(IJMK)*AXZ(IJMK) + &
72                         W_G(IJK)*AXY(IJK) - W_G(IJKM)*AXY(IJKM)
73     
74                 IF (ROP_G(IJK) > ZERO) THEN
75                    SRC = VOL(IJK)*ZMAX((-SUM_R_G(IJK)))/ROP_G(IJK)
76                 ELSE
77                    SRC = ZERO
78                 ENDIF
79     
80                 A_M(IJK,0,0) = -(A_M(IJK,east,0)+A_M(IJK,west,0)+A_M(IJK,north,0)+&
81                                  A_M(IJK,south,0)+A_M(IJK,top,0)+A_M(IJK,bottom,0)+&
82                                  VOL(IJK)*ODT+ZMAX(DEL_V)+SRC)
83                 B_M(IJK,0) = -(ROP_GO(IJK)*VOL(IJK)*ODT+&
84                                ZMAX((-DEL_V))*ROP_G(IJK)+&
85                                ZMAX(SUM_R_G(IJK))*VOL(IJK))
86     
87                 IF (ABS(A_M(IJK,0,0)) < SMALL_NUMBER) THEN
88                    IF (ABS(B_M(IJK,0)) < SMALL_NUMBER) THEN
89                       A_M(IJK,0,0) = -ONE
90                       B_M(IJK,0) = ZERO
91                    ELSE
92     !!$omp             critical
93                       WRITE (LINE, '(A,I6,A,I1,A,G12.5)') 'Error: At IJK = ', IJK, &
94                          ' M = ', 0, ' A = 0 and b = ', B_M(IJK,0)
95                       CALL WRITE_ERROR ('SOURCE_ROP_g', LINE, 1)
96     !!$omp             end critical
97                    ENDIF
98                 ENDIF
99              ELSE
100     ! set the value of rop_g in all wall and flow boundary cells to what is
101     ! known for that cell
102                 A_M(IJK,east,0) = ZERO
103                 A_M(IJK,west,0) = ZERO
104                 A_M(IJK,north,0) = ZERO
105                 A_M(IJK,south,0) = ZERO
106                 A_M(IJK,top,0) = ZERO
107                 A_M(IJK,bottom,0) = ZERO
108                 A_M(IJK,0,0) = -ONE
109                 B_M(IJK,0) = -ROP_G(IJK)
110              ENDIF
111           ENDDO
112     
113           RETURN
114           END SUBROUTINE SOURCE_ROP_G
115     
116     
117