conv_rop_g.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_ROP_g                                              C
!  Purpose: Determine convection terms for gas continuity              C
!           equation.  Master routine.                                 C
!                                                                      C
!  Author: M. Syamlal                                 Date: 19-MAR-97  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE conv_rop_g(A_M, B_M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE fldvar
      USE run
      USE compar
      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)
!-----------------------------------------------

      IF (discretize(1) == 0) THEN               ! 0 & 1 => first order upwinding
         CALL conv_rop_g0 (a_m)
      ELSE
         CALL conv_rop_g1 (a_m)
      ENDIF


      RETURN
      END SUBROUTINE conv_rop_g


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_ROP_g0                                             C
!  Purpose: Determine convection terms for gas continuity equation     C
!                                                                      C
!  Notes: The off-diagonal coefficients calculated here must be        C
!         positive. The center coefficient and the source vector       C
!         are negative.  -- First order upwinding                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 2-JUL-96   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE conv_rop_g0(A_M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE fldvar
      USE run
      USE parallel
      USE physprop
      USE geometry
      USE indices
      USE pgcor
      USE compar
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
      INTEGER :: IMJK, IJMK, IJKM
!-----------------------------------------------


! Calculate convection fluxes through each of the faces

!!$omp  parallel do private( J, K, IJK, IPJK, IJPK, IJKP, &
!!$omp&  IMJK, IJMK, IJKM) &
!!$omp&  schedule(static)
      DO ijk = ijkstart3, ijkend3
         IF (phase_4_p_g(ijk) /= 0) THEN
! if phase_4_p_g is assigned the index of a solids phase then this
! section is entered. currently, phase_4_p_g is always assigned to
! the gas phase index 0 and so this section will never be entered.

! it was previously possible for phase_4_p_g to be assigned the index of
! the gas phase in some cells, and the index of a solids phase in other
! cells if that solids phase has close_packed=F and was in higher
! concentrations than the gas. in such a case this branch would become
! activated, while the corresponding solids continuity would become
! skipped. moreover, that solids phase's volume fraction would be
! corrected rather than the gas phases void fraction in calc_vol_fr.
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
            ipjk = ip_of(ijk)
            ijpk = jp_of(ijk)
            ijkp = kp_of(ijk)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)

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

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

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

! Modify west (i-1/2,j,k), south (i j-1/2,k) and bottom (i,j,k-1/2)
! faces if the neighboring west, south, bottom cells have
! phase_4_p_g of 0
            IF(phase_4_p_g(imjk)==0) a_m(ijk,west,0)=zmax(u_g(imjk))*ayz(imjk)
            IF(phase_4_p_g(ijmk)==0) a_m(ijk,south,0)=zmax(v_g(ijmk))*axz(ijmk)
            IF (do_k) THEN
               IF(phase_4_p_g(ijkm)==0) a_m(ijk,bottom,0)=zmax(w_g(ijkm))*axy(ijkm)
            ENDIF

         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE conv_rop_g0


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_ROP_g1                                             C
!  Purpose: Determine convection terms for gas continuity equation.    C
!                                                                      C
!  Notes: The off-diagonal coefficients calculated here must be        C
!         positive. The center coefficient and the source vector       C
!         are negative.  -- Higher order schemes                       C
!                                                                      C
!  Author: M. Syamlal                                 Date:19-MAR-97   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE conv_rop_g1(A_M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE compar
      USE fldvar
      USE functions
      USE geometry
      USE indices
      USE parallel
      USE param
      USE param1
      USE pgcor
      USE physprop
      USE run
      USE xsi
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
      INTEGER :: IMJK, IJMK, IJKM
! loezos
      INTEGER :: incr
! temporary use of global arrays:
! xsi_array: convection weighting factors
      DOUBLE PRECISION :: XSI_e(dimension_3), XSI_n(dimension_3),&
                          XSI_t(dimension_3)
!-----------------------------------------------


! loezos
      incr=0

! Calculate convection factors
      CALL calc_xsi (discretize(1), rop_g, u_g, v_g, w_g, xsi_e, &
         xsi_n, xsi_t,incr)

! Calculate convection fluxes through each of the faces

!!$omp  parallel do private( J, K, IJK, IPJK, IJPK, IJKP,  &
!!$omp&  IMJK, IJMK, IJKM) &
!!$omp&  schedule(static)
      DO ijk = ijkstart3, ijkend3
         IF (phase_4_p_g(ijk) /= 0) THEN
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
            ipjk = ip_of(ijk)
            ijpk = jp_of(ijk)
            ijkp = kp_of(ijk)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)

! East face (i+1/2, j, k)
            a_m(ijk,east,0) = -xsi_e(ijk)*u_g(ijk)*ayz(ijk)
            a_m(ipjk,west,0) = (one - xsi_e(ijk))*u_g(ijk)*ayz(ijk)

! North face (i, j+1/2, k)
            a_m(ijk,north,0) = -xsi_n(ijk)*v_g(ijk)*axz(ijk)
            a_m(ijpk,south,0) = (one - xsi_n(ijk))*v_g(ijk)*axz(ijk)

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

! Modify west (i-1/2,j,k), south (i j-1/2,k) and bottom (i,j,k-1/2)
! faces if the neighboring west, south, bottom cells have
! phase_4_p_g of 0
            IF (phase_4_p_g(imjk) == 0) a_m(ijk,west,0) = (one - xsi_e(imjk))*u_g(&
               imjk)*ayz(imjk)
            IF (phase_4_p_g(ijmk) == 0) a_m(ijk,south,0) = (one - xsi_n(ijmk))*v_g(&
               ijmk)*axz(ijmk)
            IF (do_k) THEN
               IF (phase_4_p_g(ijkm) == 0) a_m(ijk,bottom,0) = (one - xsi_t(ijkm))*&
                  w_g(ijkm)*axy(ijkm)
            ENDIF

         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE conv_rop_g1