calc_e.f

Go to the documentation of this file.

! TO DO:
! Check the formulation based on MCp
!
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_e_e                                                C
!  Purpose: Calculate coefficients linking velocity correction to      C
!           solids volume fraction correction -- East                  C
!                                                                      C
!  Author: M. Syamlal                                 Date: 1-JUL-96   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE calc_e_e(A_M, MCP, E_E)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE fldvar
      USE geometry
      USE indices
      USE physprop
      USE run
      USE constant
      USE compar
      USE sendrecv
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Index of close packed solids phase. note mcp is the lowest index of
! those solids phases that have close_packed=t and of the solids
! phase that is used for the solids correction equation (when mmax=1).
      INTEGER, INTENT(IN) :: Mcp
! Coefficient for solids correction equation. These coefficients are
! initialized to zero in the subroutine time_march before the time loop.
! Note that the E_e coefficients are only used when mcp is assigned
! (when any solids phase has close_packed=T).
      DOUBLE PRECISION, INTENT(INOUT) :: e_e(dimension_3)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: IJK
!-----------------------------------------------

! returning if mcp is not defined.
      IF (mcp == undefined_i) RETURN

! returning if the x-momentum equation of this phase is not to be
! solved.
      IF (.NOT.momentum_x_eq(mcp)) RETURN

!!$omp parallel do private(ijk)
      DO ijk = ijkstart3, ijkend3
         IF (sip_at_e(ijk) .OR. mflow_at_e(ijk)) THEN
            e_e(ijk) = zero
         ELSE
            IF (a_m(ijk,0,mcp) /= zero) THEN
! calculating the correction coefficient
               e_e(ijk) = ayz(ijk)/(-a_m(ijk,0,mcp))
            ELSE
               e_e(ijk) = large_number
            ENDIF
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_e_e

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_e_n                                                C
!  Purpose: Calculate coefficients linking velocity correction to      C
!           solids volume fraction correction -- North                 C
!                                                                      C
!  Author: M. Syamlal                                 Date: 1-JUL-96   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE calc_e_n(A_M, MCP, E_N)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE fldvar
      USE geometry
      USE indices
      USE physprop
      USE run
      USE constant
      USE compar
      USE sendrecv
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Index of close packed solids phase. note mcp is the lowest index of
! those solids phases that have close_packed=t and of the solids
! phase that is used for the solids correction equation (when mmax=1).
      INTEGER, INTENT(IN) :: Mcp
! Coefficient for solids correction equation. These coefficients are
! initialized to zero in the subroutine time_march before the time loop.
! Note that the E_n coefficients are only used when mcp is assigned
! (when the solids phase has close_packed=T)
      DOUBLE PRECISION, INTENT(INOUT) :: e_n(dimension_3)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: IJK
!-----------------------------------------------

! returning if mcp is not defined.
      IF (mcp == undefined_i) RETURN
! returning if the y-momentum equation of this phase is not to be
! solved.
      IF (.NOT.momentum_y_eq(mcp)) RETURN

!!$omp parallel do private(IJK)
      DO ijk = ijkstart3, ijkend3
         IF (sip_at_n(ijk) .OR. mflow_at_n(ijk)) THEN
            e_n(ijk) = zero
         ELSE
            IF ((-a_m(ijk,0,mcp)) /= zero) THEN
! calculating the correction coefficient
               e_n(ijk) = axz(ijk)/(-a_m(ijk,0,mcp))
            ELSE
               e_n(ijk) = large_number
            ENDIF
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_e_n

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_e_t                                                C
!  Purpose: Calculate coefficients linking velocity correction to      C
!           solids volume fraction correction -- Top                   C
!                                                                      C
!  Author: M. Syamlal                                 Date: 1-JUL-96   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE calc_e_t(A_M, MCP, E_T)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE fldvar
      USE geometry
      USE indices
      USE physprop
      USE run
      USE constant
      USE compar
      USE sendrecv
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Index of close packed solids phase. note mcp is the lowest index of
! those solids phases that have close_packed=t and of the solids
! phase that is used for the solids correction equation (when mmax=1).
      INTEGER, INTENT(IN) :: Mcp
! Coefficient for solids correction equation. These coefficients are
! initialized to zero in the subroutine time_march before the time loop.
! Note that the E_t coefficients are only used when mcp is assigned
! (when the solids phase has close_packed=T)
      DOUBLE PRECISION, INTENT(INOUT) :: e_t(dimension_3)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: IJK
!-----------------------------------------------

! returning if mcp is not defined.
      IF (mcp == undefined_i) RETURN
! returning if the z-momentum equation of this phase is not to be
! solved
      IF (.NOT.momentum_z_eq(mcp)) RETURN

!!$omp parallel do private(IJK)
      DO ijk = ijkstart3, ijkend3
         IF (sip_at_t(ijk) .OR. mflow_at_t(ijk)) THEN
            e_t(ijk) = zero
         ELSE
            IF ((-a_m(ijk,0,mcp)) /= zero) THEN
! calculating the correction coefficient
               e_t(ijk) = axy(ijk)/(-a_m(ijk,0,mcp))
            ELSE
               e_t(ijk) = large_number
            ENDIF
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_e_t