File: RELATIVE:/../../../mfix.git/model/adjust_a_v_s.f

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