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