adjust_a_u_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: ADJUST_A_U_s                                            C
!  Purpose: Handle the special case of the center coefficient in       C
!           U_s momentum eq. becoming zero.                            C
!                                                                      C
!  Author: M. Syamlal                                 Date:  2-AUG-96  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE adjust_a_u_s(A_M, B_M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE fldvar
      USE physprop
      USE geometry
      USE run
      USE indices
      USE compar
      USE sendrecv
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy Arguments
!-----------------------------------------------
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! Indices
      INTEGER :: I, IP, IJK, IJKE, IMJK
! Phase index
      INTEGER :: M
!-----------------------------------------------

      DO m = 1, mmax
         IF (drag_type_enum == ghd_2007 .AND. m /= mmax) cycle
         IF (momentum_x_eq(m)) THEN

!!$omp     parallel do private(I, IP, IJK, IJKE, IMJK )

            DO ijk = ijkstart3, ijkend3
               IF (abs(a_m(ijk,0,m)) < small_number) THEN
                  a_m(ijk,east,m) = zero
                  a_m(ijk,west,m) = zero
                  a_m(ijk,north,m) = zero
                  a_m(ijk,south,m) = zero
                  a_m(ijk,top,m) = zero
                  a_m(ijk,bottom,m) = zero
                  a_m(ijk,0,m) = -one
                  IF (b_m(ijk,m) < zero) THEN
                     ijke = east_of(ijk)
                     ip = ip1(i_of(ijk))
                     IF (rop_s(ijke,m)*ayz_u(ijk) > small_number) THEN
                        b_m(ijk,m) = sqrt((-b_m(ijk,m)/(rop_s(ijke,m)*&
                           avg_x_e(one,zero,ip)*ayz_u(ijk))))
                     ELSE
                        b_m(ijk,m) = zero
                     ENDIF
                  ELSEIF (b_m(ijk,m) > zero) THEN
                     i = i_of(ijk)
                     imjk = im_of(ijk)
                     IF (rop_s(ijk,m)*ayz_u(imjk) > small_number) THEN
                        b_m(ijk,m) = sqrt(b_m(ijk,m)/(rop_s(ijk,m)*&
                           avg_x_e(zero,one,i)*ayz_u(imjk)))
                     ELSE
                        b_m(ijk,m) = zero
                     ENDIF
                  ENDIF
               ENDIF    ! end if (abs(a_m(ijk,0,m))<small_number)
            ENDDO    ! end do loop (ijk=ijkstart3,ijkend3)

         ENDIF   ! end if (momentum_x_eq(m))
      ENDDO   ! end do loop (m=1,mmax)

      RETURN
      END SUBROUTINE adjust_a_u_s