conv_pp_g.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_Pp_g                                               C
!  Purpose: Determine convection terms for Pressure correction         C
!           equation.                                                  C
!                                                                      C
!  Notes: The off-diagonal coefficients calculated here must be        C
!         positive. The center coefficient and the source vector are   C
!         negative.                                                    C
!         Multiplication with factors d_e, d_n, and d_t are carried    C
!         out in source_pp_g.  Constant pressure boundaries are        C
!         handled by holding the Pp_g at the boundaries zero.  For     C
!         specified mass flow boundaries (part of) a's are calculated  C
!         here since b is calculated from a's in source_pp_g.  After   C
!         calculating b, a's are multiplied by d and at the flow       C
!         boundaries are set to zero.                                  C
!                                                                      C
!  Author: M. Syamlal                                 Date: 20-JUN-96  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number: 1                                                  C
!  Purpose: to use face densities calculated in CONV_ROP               C
!  Author: M. Syamlal                                 Date: 1-JUN-05   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE conv_pp_g(A_M, B_M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE fldvar
      USE run
      USE parallel
      USE physprop
      USE geometry
      USE indices
      USE pgcor
      USE compar
      USE mflux
      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 :: IJK, IPJK, IJPK, IJKP
      INTEGER :: IMJK, IJMK, IJKM
      INTEGER :: M
! local value of A_m
      DOUBLE PRECISION :: am
!-----------------------------------------------


! Calculate convection fluxes through each of the faces

!$omp  parallel do default(none) private( IJK, IPJK, IJPK, IJKP, M, AM, IMJK, IJMK, IJKM) &
!$omp              shared(ijkstart3, ijkend3, rop_ge, rop_gn, rop_gt, axz, axy, ayz, do_k, a_m)
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
            ipjk = ip_of(ijk)
            ijpk = jp_of(ijk)
            ijkp = kp_of(ijk)

! East face (i+1/2, j, k)
            am = rop_ge(ijk)*ayz(ijk)
            a_m(ijk,east,0) = am
            a_m(ipjk,west,0) = am

! North face (i, j+1/2, k)
            am = rop_gn(ijk)*axz(ijk)
            a_m(ijk,north,0) = am
            a_m(ijpk,south,0) = am

! Top face (i, j, k+1/2)
            IF (do_k) THEN
               am = rop_gt(ijk)*axy(ijk)
               a_m(ijk,top,0) = am
               a_m(ijkp,bottom,0) = am
            ENDIF

! West face (i-1/2, j, k)
            imjk = im_of(ijk)
            IF (.NOT.fluid_at(imjk)) THEN
               am = rop_ge(imjk)*ayz(imjk)
               a_m(ijk,west,0) = am
            ENDIF

! South face (i, j-1/2, k)
            ijmk = jm_of(ijk)
            IF (.NOT.fluid_at(ijmk)) THEN
               am = rop_gn(ijmk)*axz(ijmk)
               a_m(ijk,south,0) = am
            ENDIF

! Bottom face (i, j, k-1/2)
            IF (do_k) THEN
               ijkm = km_of(ijk)
               IF (.NOT.fluid_at(ijkm)) THEN
                  am = rop_gt(ijkm)*axy(ijkm)
                  a_m(ijk,bottom,0) = am
               ENDIF
            ENDIF
         ENDIF   ! end if (fluid_at(ijk))
      ENDDO   ! end do ijk=ijkstart3,ijkend3


! solids phase terms are needed for those solids phases that do not
! become close packed. these terms are used to essentially make the
! matrix equation for gas pressure correction a mixture pressure
! correction equation by adding solids phases that do not become
! close packed
      DO m = 1, mmax
         IF (.NOT.close_packed(m)) THEN
            DO ijk = ijkstart3, ijkend3
               IF (fluid_at(ijk)) THEN
                  ipjk = ip_of(ijk)
                  ijpk = jp_of(ijk)
                  ijkp = kp_of(ijk)

! East face (i+1/2, j, k)
                  am = rop_se(ijk,m)*ayz(ijk)
                  a_m(ijk,east,m) = am
                  a_m(ipjk,west,m) = am

! North face (i, j+1/2, k)
                  am = rop_sn(ijk,m)*axz(ijk)
                  a_m(ijk,north,m) = am
                  a_m(ijpk,south,m) = am

! Top face (i, j, k+1/2)
                  IF (do_k) THEN
                     am = rop_st(ijk,m)*axy(ijk)
                     a_m(ijk,top,m) = am
                     a_m(ijkp,bottom,m) = am
                  ENDIF

! West face (i-1/2, j, k)
                  imjk = im_of(ijk)
                  IF (.NOT.fluid_at(imjk)) THEN
                     am = rop_se(imjk,m)*ayz(imjk)
                     a_m(ijk,west,m) = am
                  ENDIF

! South face (i, j-1/2, k)
                  ijmk = jm_of(ijk)
                  IF (.NOT.fluid_at(ijmk)) THEN
                     am = rop_sn(ijmk,m)*axz(ijmk)
                     a_m(ijk,south,m) = am
                  ENDIF

! Bottom face (i, j, k-1/2)
                  IF (do_k) THEN
                     ijkm = km_of(ijk)
                     IF (.NOT.fluid_at(ijkm)) THEN
                        am = rop_st(ijkm,m)*axy(ijkm)
                        a_m(ijk,bottom,m) = am
                     ENDIF
                  ENDIF

               ENDIF   ! end if(fluid_at(ijk))
            ENDDO   ! end do ijk=ijkstart3,ijkend3
         ENDIF   ! end if(.not.close_packed)
      ENDDO   ! end do (m=1,mmax)

      RETURN
      END SUBROUTINE conv_pp_g