File: N:\mfix\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 fldvar
22           USE physprop
23           USE geometry
24           USE run
25           USE indices
26           USE compar
27           USE sendrecv
28           USE fun_avg
29           USE functions
30           IMPLICIT NONE
31     !-----------------------------------------------
32     ! Dummy Arguments
33     !-----------------------------------------------
34     ! Septadiagonal matrix A_m
35           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
36     ! Vector b_m
37           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
38     !-----------------------------------------------
39     ! Local Variables
40     !-----------------------------------------------
41     ! Indices
42           INTEGER :: IJK, IJKN, IJMK
43     ! Phase index
44           INTEGER :: M
45     !-----------------------------------------------
46     
47           DO M = 1, MMAX
48              IF (DRAG_TYPE_ENUM == GHD_2007 .AND. M /= MMAX) CYCLE
49              IF (MOMENTUM_Y_EQ(M)) THEN
50     
51     !!$omp     parallel do private(IJK,IJKN,IJMK)
52                 DO IJK = ijkstart3, ijkend3
53                    IF (ABS(A_M(IJK,0,M)) < SMALL_NUMBER) THEN
54                       A_M(IJK,east,M) = ZERO
55                       A_M(IJK,west,M) = ZERO
56                       A_M(IJK,north,M) = ZERO
57                       A_M(IJK,south,M) = ZERO
58                       A_M(IJK,top,M) = ZERO
59                       A_M(IJK,bottom,M) = ZERO
60                       A_M(IJK,0,M) = -ONE
61                       IF (B_M(IJK,M) < ZERO) THEN
62                          IJKN = NORTH_OF(IJK)
63                          IF (ROP_S(IJKN,M)*AXZ(IJK) > SMALL_NUMBER) THEN
64                             B_M(IJK,M) = SQRT((-B_M(IJK,M)/(ROP_S(IJKN,M)*&
65                                AVG_Y_N(ONE,ZERO)*AXZ(IJK))))
66                          ELSE
67                             B_M(IJK,M) = ZERO
68                          ENDIF
69                       ELSEIF (B_M(IJK,M) > ZERO) THEN
70                          IJMK = JM_OF(IJK)
71                          IF (ROP_S(IJK,M)*AXZ(IJK) > SMALL_NUMBER) THEN
72                             B_M(IJK,M) = SQRT(B_M(IJK,M)/(ROP_S(IJK,M)*&
73                                AVG_Y_N(ZERO,ONE)*AXZ(IJK)))
74                          ELSE
75                             B_M(IJK,M) = ZERO
76                          ENDIF
77                       ENDIF
78                    ENDIF    ! end if (abs(a_m(ijk,0,m))<small_number)
79                 ENDDO    ! end do loop (ijk=ijkstart3,ijkend3)
80     
81              ENDIF   ! end if (momentum_y_eq(m))
82           ENDDO   ! end do loop (m=1,mmax)
83     
84           RETURN
85           END SUBROUTINE ADJUST_A_V_S
86