File: /nfs/home/0/users/jenkins/mfix.git/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 matrix
29           USE fldvar
30           USE rxns
31           USE run
32           USE geometry
33           USE indices
34           USE pgcor
35           USE compar
36           USE functions
37           IMPLICIT NONE
38     !-----------------------------------------------
39     ! Dummy arguments
40     !-----------------------------------------------
41     ! Septadiagonal matrix A_m
42           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
43     ! Vector b_m
44           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
45     !-----------------------------------------------
46     ! Local variables
47     !-----------------------------------------------
48     ! DEL dot V
49           DOUBLE PRECISION :: DEL_V
50     ! Mass source
51           DOUBLE PRECISION :: Src
52     ! Indices
53           INTEGER :: I, J, K, IJK, IMJK, IJMK, IJKM
54     ! error message
55           CHARACTER(LEN=80) :: LINE
56     
57     !-----------------------------------------------
58     
59     !!$omp  parallel do private( I, J, K, IJK, IMJK, IJMK, IJKM,  DEL_V, &
60     !!$omp&  Src, LINE) &
61     !!$omp&  schedule(static)
62           DO IJK = ijkstart3, ijkend3
63              IF (FLUID_AT(IJK) .AND. PHASE_4_P_G(IJK)/=0) THEN
64                 I = I_OF(IJK)
65                 J = J_OF(IJK)
66                 K = K_OF(IJK)
67                 IMJK = IM_OF(IJK)
68                 IJMK = JM_OF(IJK)
69                 IJKM = KM_OF(IJK)
70     
71                 DEL_V = U_G(IJK)*AYZ(IJK) - U_G(IMJK)*AYZ(IMJK) + &
72                         V_G(IJK)*AXZ(IJK) - V_G(IJMK)*AXZ(IJMK) + &
73                         W_G(IJK)*AXY(IJK) - W_G(IJKM)*AXY(IJKM)
74     
75                 IF (ROP_G(IJK) > ZERO) THEN
76                    SRC = VOL(IJK)*ZMAX((-SUM_R_G(IJK)))/ROP_G(IJK)
77                 ELSE
78                    SRC = ZERO
79                 ENDIF
80     
81                 A_M(IJK,0,0) = -(A_M(IJK,E,0)+A_M(IJK,W,0)+A_M(IJK,N,0)+&
82                                  A_M(IJK,S,0)+A_M(IJK,T,0)+A_M(IJK,B,0)+&
83                                  VOL(IJK)*ODT+ZMAX(DEL_V)+SRC)
84                 B_M(IJK,0) = -(ROP_GO(IJK)*VOL(IJK)*ODT+&
85                                ZMAX((-DEL_V))*ROP_G(IJK)+&
86                                ZMAX(SUM_R_G(IJK))*VOL(IJK))
87     
88                 IF (ABS(A_M(IJK,0,0)) < SMALL_NUMBER) THEN
89                    IF (ABS(B_M(IJK,0)) < SMALL_NUMBER) THEN
90                       A_M(IJK,0,0) = -ONE
91                       B_M(IJK,0) = ZERO
92                    ELSE
93     !!$omp             critical
94                       WRITE (LINE, '(A,I6,A,I1,A,G12.5)') 'Error: At IJK = ', IJK, &
95                          ' M = ', 0, ' A = 0 and b = ', B_M(IJK,0)
96                       CALL WRITE_ERROR ('SOURCE_ROP_g', LINE, 1)
97     !!$omp             end critical
98                    ENDIF
99                 ENDIF
100              ELSE
101     ! set the value of rop_g in all wall and flow boundary cells to what is
102     ! known for that cell
103                 A_M(IJK,E,0) = ZERO
104                 A_M(IJK,W,0) = ZERO
105                 A_M(IJK,N,0) = ZERO
106                 A_M(IJK,S,0) = ZERO
107                 A_M(IJK,T,0) = ZERO
108                 A_M(IJK,B,0) = ZERO
109                 A_M(IJK,0,0) = -ONE
110                 B_M(IJK,0) = -ROP_G(IJK)
111              ENDIF
112           ENDDO
113     
114           RETURN
115           END SUBROUTINE SOURCE_ROP_G
116     
117     
118