File: N:\mfix\model\des\dif_phi_source_des.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: SOURCE_phi                                              !
4     !  Author: M. Syamlal                                 Date: 30-APR-97  !
5     !                                                                      !
6     !  Purpose: Determine source terms for phi eq. The terms               !
7     !     appear in the center coefficient and RHS vector. The center      !
8     !     coefficient and source vector are negative. The off-diagonal     !
9     !     coefficients are positive. S_p must be positive.                 !
10     !                                                                      !
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
12           SUBROUTINE DIF_PHI_SOURCE_DES(PHI, M, A_M, B_M, lDT)
13     
14     !-----------------------------------------------
15     ! Modules
16     !-----------------------------------------------
17           USE param
18           USE param1
19           USE parallel
20           USE scales
21           USE physprop
22           USE fldvar
23           USE rxns
24           USE run
25           USE toleranc
26           USE geometry
27           USE indices
28           USE compar
29           USE fun_avg
30           USE functions
31     
32           IMPLICIT NONE
33     !-----------------------------------------------
34     ! Dummy Arguments
35     !-----------------------------------------------
36     ! Source term on RHS
37           DOUBLE PRECISION :: S_C
38     ! Phi
39           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
40     ! phase index
41           INTEGER, INTENT(IN) :: M
42     ! Septadiagonal matrix A_m
43           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
44     ! Vector b_m
45           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
46     ! Diffusion time-step
47           DOUBLE PRECISION, INTENT(IN) :: lDT
48     !-----------------------------------------------
49     ! Local variables
50     !-----------------------------------------------
51     ! Indices
52           INTEGER :: IJK
53     
54           DOUBLE PRECISION :: APO
55           DOUBLE PRECISION :: lOoDT
56     !-----------------------------------------------
57     
58     
59           lOoDT = 1.0d0/lDT
60     
61           DO IJK = IJKSTART3, IJKEND3
62     
63              IF (FLUID_AT(IJK)) THEN
64     
65                 APO = VOL(IJK)*lOoDT
66     
67     ! Collect the terms
68                 A_M(IJK,0,M) = -(A_M(IJK,east,M)+A_M(IJK,west,M)+&
69                                  A_M(IJK,north,M)+A_M(IJK,south,M)+&
70                                  A_M(IJK,top,M)+A_M(IJK,bottom,M)+ APO)
71     
72                 S_C = APO*PHI(IJK)
73     
74                 IF(B_M(IJK,M) < S_C .OR. PHI(IJK) == ZERO) THEN
75                    B_M(IJK,M) = B_M(IJK,M) - S_C
76                 ELSE
77                    A_M(IJK,0,M) = A_M(IJK,0,M) - B_M(IJK,M)/PHI(IJK)
78                    B_M(IJK,M) = -S_C
79                 ENDIF
80     
81              ENDIF
82           ENDDO
83     
84     
85           RETURN
86           END SUBROUTINE DIF_PHI_SOURCE_DES
87