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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: ADJUST_A_U_s                                            C
4     !  Purpose: Handle the special case of the center coefficient in       C
5     !           U_s momentum eq. becoming zero.                            C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date:  2-AUG-96  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Literature/Document References:                                     C
11     !                                                                      C
12     !                                                                      C
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14     
15           SUBROUTINE ADJUST_A_U_S(A_M, B_M)
16     
17     !-----------------------------------------------
18     ! Modules
19     !-----------------------------------------------
20           USE param
21           USE param1
22           USE parallel
23           USE matrix
24           USE fldvar
25           USE physprop
26           USE geometry
27           USE run
28           USE indices
29           USE compar
30           USE sendrecv
31           USE fun_avg
32           USE functions
33           IMPLICIT NONE
34     !-----------------------------------------------
35     ! Dummy Arguments
36     !-----------------------------------------------
37     ! Septadiagonal matrix A_m
38           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
39     ! Vector b_m
40           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
41     !-----------------------------------------------
42     ! Local Variables
43     !-----------------------------------------------
44     ! Indices
45           INTEGER :: I, IP, IJK, IJKE, IMJK
46     ! Phase index
47           INTEGER :: M
48     !-----------------------------------------------
49     
50           DO M = 1, MMAX
51              IF (DRAG_TYPE_ENUM == GHD_2007 .AND. M /= MMAX) CYCLE
52              IF (MOMENTUM_X_EQ(M)) THEN
53     
54     !!$omp     parallel do private(I, IP, IJK, IJKE, IMJK )
55     
56                 DO IJK = ijkstart3, ijkend3
57                    IF (ABS(A_M(IJK,0,M)) < SMALL_NUMBER) THEN
58                       A_M(IJK,E,M) = ZERO
59                       A_M(IJK,W,M) = ZERO
60                       A_M(IJK,N,M) = ZERO
61                       A_M(IJK,S,M) = ZERO
62                       A_M(IJK,T,M) = ZERO
63                       A_M(IJK,B,M) = ZERO
64                       A_M(IJK,0,M) = -ONE
65                       IF (B_M(IJK,M) < ZERO) THEN
66                          IJKE = EAST_OF(IJK)
67                          IP = IP1(I_OF(IJK))
68                          IF (ROP_S(IJKE,M)*AYZ_U(IJK) > SMALL_NUMBER) THEN
69                             B_M(IJK,M) = SQRT((-B_M(IJK,M)/(ROP_S(IJKE,M)*&
70                                AVG_X_E(ONE,ZERO,IP)*AYZ_U(IJK))))
71                          ELSE
72                             B_M(IJK,M) = ZERO
73                          ENDIF
74                       ELSEIF (B_M(IJK,M) > ZERO) THEN
75                          I = I_OF(IJK)
76                          IMJK = IM_OF(IJK)
77                          IF (ROP_S(IJK,M)*AYZ_U(IMJK) > SMALL_NUMBER) THEN
78                             B_M(IJK,M) = SQRT(B_M(IJK,M)/(ROP_S(IJK,M)*&
79                                AVG_X_E(ZERO,ONE,I)*AYZ_U(IMJK)))
80                          ELSE
81                             B_M(IJK,M) = ZERO
82                          ENDIF
83                       ENDIF
84                    ENDIF    ! end if (abs(a_m(ijk,0,m))<small_number)
85                 ENDDO    ! end do loop (ijk=ijkstart3,ijkend3)
86     
87              ENDIF   ! end if (momentum_x_eq(m))
88           ENDDO   ! end do loop (m=1,mmax)
89     
90           RETURN
91           END SUBROUTINE ADJUST_A_U_S
92