source_rop_g.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOURCE_ROP_g                                            C
!  Purpose: Determine source terms for continuity equation.            C
!                                                                      C
!  Notes: The off-diagonal coefficients are positive. The center       C
!         coefficient and the source vector are negative.              C
!         See conv_rop_g for additional details.                       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 source_rop_g(A_M, B_M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE fldvar
      USE rxns
      USE run
      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)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! DEL dot V
      DOUBLE PRECISION :: DEL_V
! Mass source
      DOUBLE PRECISION :: Src
! Indices
      INTEGER :: I, J, K, IJK, IMJK, IJMK, IJKM
! error message
      CHARACTER(LEN=80) :: LINE

!-----------------------------------------------

!!$omp  parallel do private( I, J, K, IJK, IMJK, IJMK, IJKM,  DEL_V, &
!!$omp&  Src, LINE) &
!!$omp&  schedule(static)
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk) .AND. phase_4_p_g(ijk)/=0) THEN
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)

            del_v = u_g(ijk)*ayz(ijk) - u_g(imjk)*ayz(imjk) + &
                    v_g(ijk)*axz(ijk) - v_g(ijmk)*axz(ijmk) + &
                    w_g(ijk)*axy(ijk) - w_g(ijkm)*axy(ijkm)

            IF (rop_g(ijk) > zero) THEN
               src = vol(ijk)*zmax((-sum_r_g(ijk)))/rop_g(ijk)
            ELSE
               src = zero
            ENDIF

            a_m(ijk,0,0) = -(a_m(ijk,east,0)+a_m(ijk,west,0)+a_m(ijk,north,0)+&
                             a_m(ijk,south,0)+a_m(ijk,top,0)+a_m(ijk,bottom,0)+&
                             vol(ijk)*odt+zmax(del_v)+src)
            b_m(ijk,0) = -(rop_go(ijk)*vol(ijk)*odt+&
                           zmax((-del_v))*rop_g(ijk)+&
                           zmax(sum_r_g(ijk))*vol(ijk))

            IF (abs(a_m(ijk,0,0)) < small_number) THEN
               IF (abs(b_m(ijk,0)) < small_number) THEN
                  a_m(ijk,0,0) = -one
                  b_m(ijk,0) = zero
               ELSE
!!$omp             critical
                  WRITE (line, '(A,I6,A,I1,A,G12.5)') 'Error: At IJK = ', ijk, &
                     ' M = ', 0, ' A = 0 and b = ', b_m(ijk,0)
                  CALL write_error ('SOURCE_ROP_g', line, 1)
!!$omp             end critical
               ENDIF
            ENDIF
         ELSE
! set the value of rop_g in all wall and flow boundary cells to what is
! known for that cell
            a_m(ijk,east,0) = zero
            a_m(ijk,west,0) = zero
            a_m(ijk,north,0) = zero
            a_m(ijk,south,0) = zero
            a_m(ijk,top,0) = zero
            a_m(ijk,bottom,0) = zero
            a_m(ijk,0,0) = -one
            b_m(ijk,0) = -rop_g(ijk)
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE source_rop_g