dif_phi_source_des.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: SOURCE_phi                                              !
!  Author: M. Syamlal                                 Date: 30-APR-97  !
!                                                                      !
!  Purpose: Determine source terms for phi eq. The terms               !
!     appear in the center coefficient and RHS vector. The center      !
!     coefficient and source vector are negative. The off-diagonal     !
!     coefficients are positive. S_p must be positive.                 !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE dif_phi_source_des(PHI, M, A_M, B_M, lDT)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE scales
      USE physprop
      USE fldvar
      USE rxns
      USE run
      USE toleranc
      USE geometry
      USE indices
      USE compar
      USE fun_avg
      USE functions

      IMPLICIT NONE
!-----------------------------------------------
! Dummy Arguments
!-----------------------------------------------
! Source term on RHS
      DOUBLE PRECISION :: S_C
! Phi
      DOUBLE PRECISION, INTENT(IN) :: Phi(dimension_3)
! phase index
      INTEGER, INTENT(IN) :: M
! 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)
! Diffusion time-step
      DOUBLE PRECISION, INTENT(IN) :: lDT
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: IJK

      DOUBLE PRECISION :: APO
      DOUBLE PRECISION :: lOoDT
!-----------------------------------------------


      loodt = 1.0d0/ldt

      DO ijk = ijkstart3, ijkend3

         IF (fluid_at(ijk)) THEN

            apo = vol(ijk)*loodt

! Collect the terms
            a_m(ijk,0,m) = -(a_m(ijk,east,m)+a_m(ijk,west,m)+&
                             a_m(ijk,north,m)+a_m(ijk,south,m)+&
                             a_m(ijk,top,m)+a_m(ijk,bottom,m)+ apo)

            s_c = apo*phi(ijk)

            IF(b_m(ijk,m) < s_c .OR. phi(ijk) == zero) THEN
               b_m(ijk,m) = b_m(ijk,m) - s_c
            ELSE
               a_m(ijk,0,m) = a_m(ijk,0,m) - b_m(ijk,m)/phi(ijk)
               b_m(ijk,m) = -s_c
            ENDIF

         ENDIF
      ENDDO


      RETURN
      END SUBROUTINE dif_phi_source_des