calc_mflux.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_MFLUX                                              C
!  Purpose: Calculate the convection fluxes. Master routine.           C
!                                                                      C
!  Author: M. Syamlal                                 Date: 31-MAY-05  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_mflux()

! Modules
!---------------------------------------------------------------------//
      USE fldvar, only: u_g, v_g, w_g
      USE fldvar, only: u_s, v_s, w_s
      USE mflux, only: rop_ge, rop_gn, rop_gt
      USE mflux, only: rop_se, rop_sn, rop_st
      USE mflux, only: flux_ge, flux_gn, flux_gt
      USE mflux, only: flux_se, flux_sn, flux_st
      USE mflux, only: flux_gse, flux_gsn, flux_gst
      USE mflux, only: flux_sse, flux_ssn, flux_sst
      USE physprop, only: mmax
      USE run, only: added_mass, m_am
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
! solids phase index
      INTEGER ::  M
!---------------------------------------------------------------------//

      IF(.NOT.added_mass) THEN
         CALL calc_mflux0 (u_g, v_g, w_g, rop_ge, rop_gn, rop_gt, &
                           flux_ge, flux_gn, flux_gt)
         DO m = 1, mmax
           CALL calc_mflux0 (u_s(1,m), v_s(1,m), w_s(1,m), &
                             rop_se(1,m), rop_sn(1,m), rop_st(1,m), &
                             flux_se(1,m), flux_sn(1,m), flux_st(1,m))
         ENDDO

      ELSE
! New fluxes are defined for gas based on added mass
         CALL calc_mflux_am (u_g, v_g, w_g, rop_ge, rop_gn, rop_gt, &
                           rop_se(1,m_am), rop_sn(1,m_am), rop_st(1,m_am), &
                           flux_ge, flux_gn, flux_gt, m_am)

         CALL calc_mflux0 (u_g, v_g, w_g, rop_ge, rop_gn, rop_gt, &
                           flux_gse, flux_gsn, flux_gst)
         DO m = 1, mmax
! New fluxes are defined for M = M_am where virtual mass force is added.
           IF(m==m_am) THEN
              CALL calc_mflux_am (u_s(1,m), v_s(1,m), w_s(1,m), &
                             rop_se(1,m), rop_sn(1,m), rop_st(1,m), &
                             rop_ge, rop_gn, rop_gt, &
                             flux_se(1,m), flux_sn(1,m), flux_st(1,m),&
                             m_am)
              CALL calc_mflux0 (u_s(1,m), v_s(1,m), w_s(1,m), &
                             rop_se(1,m), rop_sn(1,m), rop_st(1,m), &
                             flux_sse, flux_ssn, flux_sst)
           ELSE
              CALL calc_mflux0 (u_s(1,m), v_s(1,m), w_s(1,m), &
                             rop_se(1,m), rop_sn(1,m), rop_st(1,m), &
                             flux_se(1,m), flux_sn(1,m), flux_st(1,m))
           ENDIF
         ENDDO
      ENDIF

      RETURN
      END SUBROUTINE calc_mflux

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine:: CALC_MFLUX0                                            C
!  Purpose: Calculate the convection fluxes.                           C
!                                                                      C
!  Author: M. Syamlal                                 Date: 31-MAY-05  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_mflux0(U, V, W, ROP_E, ROP_N, ROP_T,&
                             flux_e, flux_n, flux_t)

! Modules
!---------------------------------------------------------------------//
      USE compar, only: ijkstart3, ijkend3
      USE functions, only: fluid_at
      USE functions, only: im_of, jm_of, km_of
      USE geometry, only: do_k
      USE geometry, only: ayz, axz, axy
      USE param, only: dimension_3
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Velocity components
      DOUBLE PRECISION, INTENT(IN) :: U(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: V(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: W(dimension_3)
! Face value of density (for calculating convective fluxes)
      DOUBLE PRECISION, INTENT(IN) :: ROP_E(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: ROP_N(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: ROP_T(dimension_3)
! Convective mass fluxes
      DOUBLE PRECISION, INTENT(OUT) :: Flux_E(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: Flux_N(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: Flux_T(dimension_3)

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK, IMJK, IJMK, IJKM
!---------------------------------------------------------------------//

!$omp  parallel do default(none) private( IJK, IMJK, IJMK, IJKM) &
!$omp              shared(ijkstart3, ijkend3, flux_e, flux_n, flux_t,&
!$omp                     rop_e, rop_n, rop_t, axz, axy, ayz, u, v, w, do_k)

      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN

            imjk = im_of(ijk)
            ijmk = jm_of(ijk)

! East face (i+1/2, j, k)
            flux_e(ijk) = rop_e(ijk)*ayz(ijk)*u(ijk)
! West face (i-1/2, j, k)
            IF (.NOT.fluid_at(imjk)) THEN
               flux_e(imjk) = rop_e(imjk)*ayz(imjk)*u(imjk)
            ENDIF

! North face (i, j+1/2, k)
            flux_n(ijk) = rop_n(ijk)*axz(ijk)*v(ijk)
! South face (i, j-1/2, k)
            IF (.NOT.fluid_at(ijmk)) THEN
              flux_n(ijmk) = rop_n(ijmk)*axz(ijmk)*v(ijmk)
            ENDIF

            IF (do_k) THEN
               ijkm = km_of(ijk)
! Top face (i, j, k+1/2)
               flux_t(ijk) = rop_t(ijk)*axy(ijk)*w(ijk)
! Bottom face (i, j, k-1/2)
               IF (.NOT.fluid_at(ijkm)) THEN
                 flux_t(ijkm) = rop_t(ijkm)*axy(ijkm)*w(ijkm)
               ENDIF
            ENDIF   ! end if do_k
         ENDIF   ! end if fluid_at
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE calc_mflux0

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!                                                                      C
!  Subroutine: CALC_MFLUX_AM                                           C
!  Purpose: Calculate the convection fluxes for case with added_mass.  C
!                                                                      C
!  Author: S. Benyahia                                 Date:           C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_mflux_am(U, V, W, ROP_E, ROP_N, ROP_T, &
                             ropa_e, ropa_n, ropa_t, flux_e, &
                             flux_n, flux_t, m_am)

! Modules
!---------------------------------------------------------------------//
      USE compar, only: ijkstart3, ijkend3
      USE functions, only: fluid_at
      USE functions, only: im_of, jm_of, km_of
      USE fldvar, only: ro_s
      USE geometry, only: do_k
      USE geometry, only: ayz, axz, axy
      USE param, only: dimension_3
      USE param1, only: one
      USE physprop, only: cv
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Velocity components
      DOUBLE PRECISION, INTENT(IN) :: U(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: V(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: W(dimension_3)
! Face value of density (for calculating convective fluxes)
      DOUBLE PRECISION, INTENT(IN) :: ROP_E(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: ROP_N(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: ROP_T(dimension_3)
! Face value of density of added mass phase
      DOUBLE PRECISION, INTENT(IN) :: ROPa_E(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: ROPa_N(dimension_3)
      DOUBLE PRECISION, INTENT(IN) :: ROPa_T(dimension_3)
! Convective mass fluxes
      DOUBLE PRECISION, INTENT(OUT) :: Flux_E(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: Flux_N(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: Flux_T(dimension_3)
! phase index for phase of added_mass
      INTEGER, INTENT(IN) :: M_AM

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK, IMJK, IJMK, IJKM
!---------------------------------------------------------------------//

!!!$omp  parallel do private( IJK, IMJK, IJMK, IJKM) &
!!!$omp&  schedule(static)
      DO ijk = ijkstart3, ijkend3

         IF (fluid_at(ijk)) THEN

            imjk = im_of(ijk)
            ijmk = jm_of(ijk)

! East face (i+1/2, j, k)
            flux_e(ijk) = rop_e(ijk)* &
                          (one + cv*ropa_e(ijk)/ro_s(ijk,m_am))*&
                          ayz(ijk)*u(ijk)
! West face (i-1/2, j, k)
            IF (.NOT.fluid_at(imjk)) THEN
               flux_e(imjk) = rop_e(imjk)*&
                              (one + cv*ropa_e(imjk)/ro_s(ijk,m_am))*&
                              ayz(imjk)*u(imjk)
            ENDIF

! North face (i, j+1/2, k)
            flux_n(ijk) = rop_n(ijk)*&
                          (one + cv*ropa_n(ijk)/ro_s(ijk,m_am))*&
                          axz(ijk)*v(ijk)
! South face (i, j-1/2, k)
            IF (.NOT.fluid_at(ijmk)) THEN
              flux_n(ijmk) = rop_n(ijmk)*&
                             (one + cv*ropa_n(ijmk)/ro_s(ijk,m_am))*&
                             axz(ijmk)*v(ijmk)
            ENDIF

            IF (do_k) THEN
               ijkm = km_of(ijk)

! Top face (i, j, k+1/2)
               flux_t(ijk) = rop_t(ijk)*&
                             (one + cv*ropa_t(ijk)/ro_s(ijk,m_am))*&
                             axy(ijk)*w(ijk)
! Bottom face (i, j, k-1/2)
               IF (.NOT.fluid_at(ijkm)) THEN
                 flux_t(ijkm) = rop_t(ijkm)*&
                                (one + cv*ropa_t(ijkm)/ro_s(ijk,m_am))*&
                                axy(ijkm)*w(ijkm)
               ENDIF
            ENDIF   ! end if do_k
         ENDIF   ! end if fluid_at
      ENDDO    ! end do ijk

      RETURN
      END SUBROUTINE calc_mflux_am