conv_dif_v_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_DIF_V_s                                            C
!  Purpose: Determine convection diffusion terms for V_s momentum eqs  C
!           The off-diagonal coefficients calculated here must be      C
!           positive. The center coefficient and the source vector     C
!           are negative;                                              C
!           See source_v_s                                             C
!                                                                      C
!  Author: M. Syamlal                                 Date: 24-DEC-96  C
!  Reviewer:                                          Date:            C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE conv_dif_v_s(A_M, B_M, IER)

! Modules
!---------------------------------------------------------------------//
      USE param, only: dimension_3, dimension_m
      USE run, only: kt_type_enum
      USE run, only: ghd_2007
      USE run, only: def_cor
      USE run, only: momentum_y_eq
      USE run, only: discretize
      USE physprop, only: mmax
      USE visc_s, only: mu_s
      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)
! Error index
      INTEGER, INTENT(INOUT) :: IER

! Local Variables
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER :: M
!---------------------------------------------------------------------//

      DO m = 1, mmax
        IF(kt_type_enum /= ghd_2007 .OR. &
           (kt_type_enum == ghd_2007 .AND. m==mmax)) THEN



          IF  (momentum_y_eq(m)) THEN

             IF (def_cor) THEN
! USE DEFERRED CORRECTION TO SOLVE V_S
                CALL store_a_v_s0 (a_m(1,-3,m), m)
                IF (discretize(4) > 1)CALL store_a_v_sdc (m, b_m)

             ELSE
! DO NOT USE DEFERRED CORRECTION TO SOLVE FOR V_S
                IF (discretize(4) == 0) THEN         ! 0 & 1 => FOUP
                   CALL store_a_v_s0 (a_m(1,-3,m), m)
                ELSE
                   CALL store_a_v_s1 (a_m(1,-3,m), m)
                ENDIF
             ENDIF

             CALL dif_v_is (mu_s(1,m), a_m, m, ier)
          ENDIF
        ENDIF
      ENDDO

      RETURN
      END SUBROUTINE conv_dif_v_s

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: Calculate the components of velocity on the east, north,   C
!  and top face of a v-momentum cell                                   C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE get_vcell_svterms(U, V, WW, M)

! Modules
!---------------------------------------------------------------------//
      USE compar, only: ijkstart3, ijkend3

      USE cutcell, only: cut_v_treatment_at
      USE cutcell, only: theta_vn, theta_vn_bar
      USE cutcell, only: theta_v_ne, theta_v_se
      USE cutcell, only: theta_v_nt, theta_v_st
      USE cutcell, only: alpha_ve_c, alpha_vn_c, alpha_vt_c

      USE fldvar, only: u_s, v_s, w_s

      USE fun_avg, only: avg_y_n, avg_y
      USE functions, only: jp_of
      USE geometry, only: do_k
      USE indices, only: j_of

      USE param, only: dimension_3
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! velocity components
      DOUBLE PRECISION, INTENT(OUT) :: U(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: V(dimension_3)
      DOUBLE PRECISION, INTENT(OUT) :: WW(dimension_3)

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: IJK, J, IJPK
! for cartesian grid
      DOUBLE PRECISION :: AW, HW, VELW
!---------------------------------------------------------------------//

!!!$omp parallel do private(IJK,J,IJPK)
      DO ijk = ijkstart3, ijkend3
         j = j_of(ijk)
         ijpk = jp_of(ijk)

         IF(cut_v_treatment_at(ijk)) THEN

! East face (i+1/2, j+1/2, k)
            u(ijk) = (theta_v_se(ijk) * u_s(ijk,m) + &
                      theta_v_ne(ijk) * u_s(ijpk,m))
            CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM', &
               alpha_ve_c(ijk), aw, hw, velw)
            u(ijk) = u(ijk) * aw

! North face (i, j+1, k)
            v(ijk) = (theta_vn_bar(ijk) * v_s(ijk,m) + &
                      theta_vn(ijk) * v_s(ijpk,m))
            CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
               alpha_vn_c(ijk), aw, hw, velw)
            v(ijk) = v(ijk) * aw

! Top face (i, j+1/2, k+1/2)
            IF (do_k) THEN
               ww(ijk) = (theta_v_nt(ijk) * w_s(ijk,m) + &
                          theta_v_st(ijk) * w_s(ijpk,m))
               CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
                  alpha_vt_c(ijk), aw, hw, velw)
               ww(ijk) = ww(ijk) * aw
            ENDIF

         ELSE
            u(ijk) = avg_y(u_s(ijk,m),u_s(ijpk,m),j)
            v(ijk) = avg_y_n(v_s(ijk,m),v_s(ijpk,m))
            IF (do_k) ww(ijk) = avg_y(w_s(ijk,m),w_s(ijpk,m),j)
         ENDIF
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE get_vcell_svterms



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: Calculate the convective fluxes through the faces of a     C
!  v-momentum cell. Note the fluxes are calculated at all faces        C
!  regardless of flow_at_n of condition of the west, south, or         C
!  bottom face.                                                        C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE get_vcell_scflux_terms(FLUX_E, FLUX_W, FLUX_N, &
         flux_s, flux_t, flux_b, ijk, m)

! Modules
!---------------------------------------------------------------------//
      USE cutcell, only: cut_v_treatment_at
      USE cutcell, only: theta_vn, theta_vn_bar
      USE cutcell, only: theta_v_ne, theta_v_se
      USE cutcell, only: theta_v_nt, theta_v_st
      USE cutcell, only: alpha_ve_c, alpha_vn_c, alpha_vt_c

      USE functions, only: jp_of, im_of, jm_of, km_of

      USE geometry, only: do_k

      USE mflux, only: flux_se, flux_sn, flux_st

      USE param1, only: half
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! fluxes through faces of given ijk u-momentum cell
      DOUBLE PRECISION, INTENT(OUT) :: flux_e, flux_w
      DOUBLE PRECISION, INTENT(OUT) :: flux_n, flux_s
      DOUBLE PRECISION, INTENT(OUT) :: flux_t, flux_b
! indices
      INTEGER, INTENT(IN) :: ijk

! Local variables
!---------------------------------------------------------------------//
      INTEGER :: imjk, ijmk, ijkm
      INTEGER :: ijpk, imjpk, ijpkm

! for cartesian grid
      DOUBLE PRECISION :: AW, HW, VELW
!---------------------------------------------------------------------//
! indices
      ijpk = jp_of(ijk)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ijpkm = jp_of(ijkm)
      imjpk = jp_of(imjk)

      IF(cut_v_treatment_at(ijk)) THEN
! East face (i+1/2, j+1/2, k)
         flux_e = (theta_v_se(ijk) * flux_se(ijk,m) + &
                 theta_v_ne(ijk) * flux_se(ijpk,m))
         CALL get_interpolation_terms_s(ijk, m, 'V_MOMENTUM',&
            alpha_ve_c(ijk), aw, hw, velw)
! West face (i-1/2, j+1/2, k)
         flux_e = flux_e * aw
         flux_w = (theta_v_se(imjk) * flux_se(imjk,m) + &
                 theta_v_ne(imjk) * flux_se(imjpk,m))
         CALL get_interpolation_terms_s(ijk, m, 'V_MOMENTUM',&
            alpha_ve_c(imjk), aw, hw, velw)
         flux_w = flux_w * aw

! North face (i, j+1, k)
         flux_n = (theta_vn_bar(ijk) * flux_sn(ijk,m) + &
                 theta_vn(ijk) * flux_sn(ijpk,m))
         CALL get_interpolation_terms_s(ijk, m, 'V_MOMENTUM',&
            alpha_vn_c(ijk), aw, hw, velw)
         flux_n = flux_n * aw
! South face (i, j, k)
         flux_s = (theta_vn_bar(ijmk) * flux_sn(ijmk,m) + &
                 theta_vn(ijmk) * flux_sn(ijk,m))
         CALL get_interpolation_terms_s(ijk, m, 'V_MOMENTUM',&
            alpha_vn_c(ijmk), aw, hw, velw)
         flux_s = flux_s * aw

         IF (do_k) THEN
! Top face (i, j+1/2, k+1/2)
            flux_t = (theta_v_nt(ijk) * flux_st(ijk,m) + &
                    theta_v_st(ijk) * flux_st(ijpk,m))
            CALL get_interpolation_terms_s(ijk, m, 'V_MOMENTUM',&
               alpha_vt_c(ijk), aw, hw, velw)
            flux_t = flux_t * aw
! Bottom face (i, j+1/2, k-1/2)
            flux_b = (theta_v_nt(ijkm) * flux_st(ijkm,m) + &
                    theta_v_st(ijkm) * flux_st(ijpkm,m))
            CALL get_interpolation_terms_s(ijk, m, 'V_MOMENTUM',&
               alpha_vt_c(ijkm), aw, hw, velw)
            flux_b = flux_b * aw
         ENDIF
      ELSE
         flux_e = half * (flux_se(ijk,m)  + flux_se(ijpk,m))
         flux_w = half * (flux_se(imjk,m) + flux_se(imjpk,m))
         flux_n = half * (flux_sn(ijk,m)  + flux_sn(ijpk,m))
         flux_s = half * (flux_sn(ijmk,m) + flux_sn(ijk,m))

         IF (do_k) THEN
            flux_t = half * (flux_st(ijk,m)  + flux_st(ijpk,m))
            flux_b = half * (flux_st(ijkm,m) + flux_st(ijpkm,m))
         ENDIF
      ENDIF   ! end if/else cut_v_treatmeant_at

      RETURN
      END SUBROUTINE get_vcell_scflux_terms


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: Calculate the components of diffusive flux through the     C
!  the faces of a v-momentum cell. Note the fluxes are calculated at   C
!  all faces regardless of flow_at_n of condition of the west, south   C
!  or bottom face.                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE get_vcell_sdiff_terms(D_FE, D_FW, D_FN, D_FS, &
         d_ft, d_fb, ijk, m)

! Modules
!---------------------------------------------------------------------//
      USE cutcell, only: cut_v_treatment_at
      USE cutcell, only: oneodx_e_v, oneody_n_v, oneodz_t_v

      USE functions, only: wall_at
      USE functions, only: east_of, north_of, top_of
      USE functions, only: west_of, bottom_of
      USE functions, only: im_of, jm_of, km_of

      USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h

      USE geometry, only: odx_e, ody, odz_t
      USE geometry, only: do_k
      USE geometry, only: ox
      USE geometry, only: ayz_v, axz_v, axy_v

      USE indices, only: i_of, j_of, k_of
      USE indices, only: jp1, im1, km1

      USE visc_s, only: mu_s
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! diffusion through faces of given ijk v-momentum cell
      DOUBLE PRECISION, INTENT(OUT) :: d_fe, d_fw
      DOUBLE PRECISION, INTENT(OUT) :: d_fn, d_fs
      DOUBLE PRECISION, INTENT(OUT) :: d_ft, d_fb
! ijk index
      INTEGER, INTENT(IN) :: ijk

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: imjk, ijmk, ijkm
      INTEGER :: i, j, k, jp, im, km
      INTEGER :: ijkc, ijkn, ijke, ijkne, ijkw, ijknw
      INTEGER :: ijkt, ijktn, ijkb, ijkbn
! length terms
      DOUBLE PRECISION :: C_AE, C_AW, C_AN, C_AS, C_AT, C_AB
!---------------------------------------------------------------------//

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

      i = i_of(ijk)
      j = j_of(ijk)
      k = k_of(ijk)
      im = im1(i)
      jp = jp1(j)
      km = km1(k)

      ijkn = north_of(ijk)
      IF (wall_at(ijk)) THEN
         ijkc = ijkn
      ELSE
         ijkc = ijk
      ENDIF
      ijke = east_of(ijk)
      ijkne = east_of(ijkn)
      ijkw = west_of(ijk)
      ijknw = north_of(ijkw)

      IF(cut_v_treatment_at(ijk)) THEN
         c_ae = oneodx_e_v(ijk)
         c_aw = oneodx_e_v(imjk)
         c_an = oneody_n_v(ijk)
         c_as = oneody_n_v(ijmk)
         c_at = oneodz_t_v(ijk)
         c_ab = oneodz_t_v(ijkm)
      ELSE
         c_ae = odx_e(i)
         c_aw = odx_e(im)
         c_an = ody(jp)
         c_as = ody(j)
         c_at = odz_t(k)
         c_ab = odz_t(km)
      ENDIF

! East face (i+1/2, j+1/2, k)
      d_fe = avg_y_h(avg_x_h(mu_s(ijkc,m),mu_s(ijke,m),i),&
                     avg_x_h(mu_s(ijkn,m),mu_s(ijkne,m),i),j)*&
             c_ae*ayz_v(ijk)
! West face (i-1/2, j+1/2, k)
      d_fw = avg_y_h(avg_x_h(mu_s(ijkw,m),mu_s(ijkc,m),im),&
                     avg_x_h(mu_s(ijknw,m),mu_s(ijkn,m),im),j)*&
             c_aw*ayz_v(imjk)

! North face (i, j+1, k)
      d_fn = mu_s(ijkn,m)*c_an*axz_v(ijk)
! South face (i, j, k)
      d_fs = mu_s(ijkc,m)*c_as*axz_v(ijmk)

      IF (do_k) THEN
         ijkt = top_of(ijk)
         ijktn = north_of(ijkt)
         ijkb = bottom_of(ijk)
         ijkbn = north_of(ijkb)

! Top face (i, j+1/2, k+1/2)
         d_ft = avg_y_h(avg_z_h(mu_s(ijkc,m),mu_s(ijkt,m),k),&
                        avg_z_h(mu_s(ijkn,m),mu_s(ijktn,m),k),j)*&
                ox(i)*c_at*axy_v(ijk)
! Bottom face (i, j+1/2, k-1/2)
         d_fb = avg_y_h(avg_z_h(mu_s(ijkb,m),mu_s(ijkc,m),km),&
                        avg_z_h(mu_s(ijkbn,m),mu_s(ijkn,m),km),j)*&
                ox(i)*c_ab*axy_v(ijkm)
      ENDIF   ! end if (do_k)

      RETURN
      END SUBROUTINE get_vcell_sdiff_terms


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_V_s0                                            C
!  Purpose: Determine convection diffusion terms for V_s momentum eqs  C
!           The off-diagonal coefficients calculated here must be      C
!           positive. The center coefficient and the source vector     C
!           are negative. FOUP                                         C
!           See source_v_s                                             C
!                                                                      C
!  Author: M. Syamlal                                 Date: 7-JUN-96   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number: 1                                                  C
!  Purpose: To incorporate Cartesian grid modifications                C
!  Author: Jeff Dietiker                              Date: 01-Jul-09  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE store_a_v_s0(A_V_S, M)

! Modules
!---------------------------------------------------------------------//
      USE compar, only: ijkstart3, ijkend3

      USE functions, only: flow_at_n
      USE functions, only: ip_of, jp_of, kp_of
      USE functions, only: im_of, jm_of, km_of

      USE geometry, only: do_k

      USE param
      USE param1, only: zero
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER, INTENT(IN) :: M
! Septadiagonal matrix A_V_s
      DOUBLE PRECISION, INTENT(INOUT) :: A_V_s(dimension_3, -3:3, m:m)

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK
      INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
! Face mass flux
      DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
      DOUBLE PRECISION :: flux_t, flux_b
! Diffusion parameter
      DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb

!---------------------------------------------------------------------//

!$omp     parallel do default(none)                              &
!$omp     private(IJK, IMJK, IJMK, IPJK, IJPK, IJKM, IJKP,       &
!$omp             D_fe, d_fw, d_fn, d_fs, d_ft, d_fb,            &
!$omp             flux_e, flux_w, flux_n, flux_s, flux_t,        &
!$omp             flux_b)                                        &
!$omp     shared(ijkstart3, ijkend3, do_k, a_v_s, m)
      DO ijk = ijkstart3, ijkend3

         IF (flow_at_n(ijk)) THEN

! Calculate convection-diffusion fluxes through each of the faces
            CALL get_vcell_scflux_terms(flux_e, flux_w, flux_n, &
               flux_s, flux_t, flux_b, ijk, m)
            CALL get_vcell_sdiff_terms(d_fe, d_fw, d_fn, d_fs, &
               d_ft, d_fb, ijk, m)
            ipjk = ip_of(ijk)
            ijpk = jp_of(ijk)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)

! East face (i+1/2, j+1/2, k)
            IF (flux_e >= zero) THEN
               a_v_s(ijk,east,m) = d_fe
               a_v_s(ipjk,west,m) = d_fe + flux_e
            ELSE
               a_v_s(ijk,east,m) = d_fe - flux_e
               a_v_s(ipjk,west,m) = d_fe
            ENDIF
! West face (i-1/2, j+1/2, k)
            IF (.NOT.flow_at_n(imjk)) THEN
               IF (flux_w >= zero) THEN
                  a_v_s(ijk,west,m) = d_fw + flux_w
               ELSE
                  a_v_s(ijk,west,m) = d_fw
               ENDIF
            ENDIF


! North face (i, j+1, k)
            IF (flux_n >= zero) THEN
               a_v_s(ijk,north,m) = d_fn
               a_v_s(ijpk,south,m) = d_fn + flux_n
            ELSE
               a_v_s(ijk,north,m) = d_fn - flux_n
               a_v_s(ijpk,south,m) = d_fn
            ENDIF
! South face (i, j, k)
            IF (.NOT.flow_at_n(ijmk)) THEN
               IF (flux_s >= zero) THEN
                  a_v_s(ijk,south,m) = d_fs + flux_s
               ELSE
                  a_v_s(ijk,south,m) = d_fs
               ENDIF
            ENDIF


            IF (do_k) THEN
               ijkp = kp_of(ijk)
               ijkm = km_of(ijk)

! Top face (i, j+1/2, k+1/2)
               IF (flux_t >= zero) THEN
                  a_v_s(ijk,top,m) = d_ft
                  a_v_s(ijkp,bottom,m) = d_ft + flux_t
               ELSE
                  a_v_s(ijk,top,m) = d_ft - flux_t
                  a_v_s(ijkp,bottom,m) = d_ft
               ENDIF
! Bottom face (i, j+1/2, k-1/2)
               IF (.NOT.flow_at_n(ijkm)) THEN
                  IF (flux_b >= zero) THEN
                     a_v_s(ijk,bottom,m) = d_fb + flux_b
                  ELSE
                     a_v_s(ijk,bottom,m) = d_fb
                  ENDIF
               ENDIF
            ENDIF   ! end if (do_k)

         ENDIF   ! end if (flow_at_n)
      ENDDO   ! end do ijk
!$omp end parallel do

      RETURN
      END SUBROUTINE store_a_v_s0


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_V_sdc                                           C
!  Purpose: To use deferred correction method to solve the v-momentum  C
!           equation. This method combines first order upwind and a    C
!           user specified higher order method                         C
!                                                                      C
!  Author: C. GUENTHER                                 Date: 8-APR-99  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number: 1                                                  C
!  Purpose: To incorporate Cartesian grid modifications                C
!  Author: Jeff Dietiker                              Date: 01-Jul-09  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE store_a_v_sdc(M, B_M)

! Modules
!---------------------------------------------------------------------//
      USE compar, only: ijkstart3, ijkend3

      USE discretization, only: fpfoi_of

      USE fldvar, only: v_s

      USE function3, only: funijk3
      USE functions, only: flow_at_n
      USE functions, only: ip_of, jp_of, kp_of
      USE functions, only: im_of, jm_of, km_of

      USE geometry, only: do_k

      USE indices, only: i_of, j_of, k_of

      USE param, only: dimension_3, dimension_m, dimension_4
      USE param1, only: zero

      USE run, only: discretize, fpfoi
      USE sendrecv3, only: send_recv3
      USE xsi, only: calc_xsi

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER, INTENT(IN) :: M
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: I, J, K, IJK
      INTEGER :: IPJK, IMJK, IJMK, IJPK, IJKM, IJKP
      INTEGER :: IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
      INTEGER :: IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
! indication for shear
      INTEGER :: incr
! deferred corrction contribution form high order method
      DOUBLE PRECISION :: MOM_HO
! low order approximation
      DOUBLE PRECISION :: MOM_LO
! convection factor at each face
      DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
      DOUBLE PRECISION :: flux_t, flux_b
! deferred correction contributions from each face
      DOUBLE PRECISION :: EAST_DC
      DOUBLE PRECISION :: WEST_DC
      DOUBLE PRECISION :: NORTH_DC
      DOUBLE PRECISION :: SOUTH_DC
      DOUBLE PRECISION :: TOP_DC
      DOUBLE PRECISION :: BOTTOM_DC

! temporary use of global arrays:
! array1 (locally u)  - the x directional velocity
      DOUBLE PRECISION :: U(dimension_3)
! array2 (locally v)  - the y directional velocity
      DOUBLE PRECISION :: V(dimension_3)
! array3 (locally ww) - the z directional velocity
      DOUBLE PRECISION :: WW(dimension_3)
!---------------------------------------------------------------------//
      DOUBLE PRECISION :: TMP4(dimension_4)
      DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t

      CALL get_vcell_svterms(u, v, ww, m)

! Send recv the third ghost layer
      IF (fpfoi) THEN
         Do ijk = ijkstart3, ijkend3
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
            ijk4 = funijk3(i,j,k)
            tmp4(ijk4) = v_s(ijk,m)
         ENDDO
         CALL send_recv3(tmp4)
      ENDIF

! shear indicator:
      incr=2
      CALL calc_xsi (discretize(4), v_s(1,m), u, v, ww, xsi_e, xsi_n, &
                     xsi_t, incr)

!!!$omp      parallel do                                             &
!!!$omp&     private( I, J, K, IJK, IPJK, IMJK, IJPK, IJMK,          &
!!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_b, flux_t, &
!!!$omp&             MOM_HO, MOM_LO, EAST_DC, WEST_DC, NORTH_DC,     &
!!!$omp&             SOUTH_DC, TOP_DC, BOTTOM_DC)
      DO ijk = ijkstart3, ijkend3

         IF (flow_at_n(ijk)) THEN

! Calculate convection fluxes through each of the faces
            CALL get_vcell_scflux_terms(flux_e, flux_w, flux_n, &
               flux_s, flux_t, flux_b, ijk, m)

            ipjk = ip_of(ijk)
            imjk = im_of(ijk)
            ijpk = jp_of(ijk)
            ijmk = jm_of(ijk)
            ijkp = kp_of(ijk)
            ijkm = km_of(ijk)

! Third Ghost layer information
            ippp  = ip_of(ip_of(ipjk))
            ippp4 = funijk3(i_of(ippp), j_of(ippp), k_of(ippp))
            immm  = im_of(im_of(imjk))
            immm4 = funijk3(i_of(immm), j_of(immm), k_of(immm))
            jppp  = jp_of(jp_of(ijpk))
            jppp4 = funijk3(i_of(jppp), j_of(jppp), k_of(jppp))
            jmmm  = jm_of(jm_of(ijmk))
            jmmm4 = funijk3(i_of(jmmm), j_of(jmmm), k_of(jmmm))
            kppp  = kp_of(kp_of(ijkp))
            kppp4 = funijk3(i_of(kppp), j_of(kppp), k_of(kppp))
            kmmm  = km_of(km_of(ijkm))
            kmmm4 = funijk3(i_of(kmmm), j_of(kmmm), k_of(kmmm))


! East face (i+1/2, j+1/2, k)
            IF(u(ijk) >= zero)THEN
               mom_lo = v_s(ijk,m)
               IF (fpfoi) mom_ho = fpfoi_of(v_s(ipjk,m), v_s(ijk,m), &
                                   v_s(imjk,m), v_s(im_of(imjk),m))
            ELSE
                mom_lo = v_s(ipjk,m)
                IF (fpfoi) mom_ho = fpfoi_of(v_s(ijk,m), v_s(ipjk,m), &
                                    v_s(ip_of(ipjk),m), tmp4(ippp4))
            ENDIF
            IF (.NOT. fpfoi) mom_ho = xsi_e(ijk)*v_s(ipjk,m)+ &
                                      (1.0-xsi_e(ijk))*v_s(ijk,m)
            east_dc = flux_e*(mom_lo-mom_ho)

! West face (i-1/2, j+1/2, k)
            IF(u(imjk) >= zero)THEN
              mom_lo = v_s(imjk,m)
              IF (fpfoi) mom_ho = fpfoi_of(v_s(ijk,m), v_s(imjk,m), &
                                  v_s(im_of(imjk),m), tmp4(immm4))
            ELSE
               mom_lo = v_s(ijk,m)
               IF (fpfoi) mom_ho = fpfoi_of(v_s(imjk,m), v_s(ijk,m), &
                                   v_s(ipjk,m), v_s(ip_of(ipjk),m))
            ENDIF
            IF (.NOT. fpfoi) mom_ho = xsi_e(imjk)*v_s(ijk,m)+ &
                                      (1.0-xsi_e(imjk))*v_s(imjk,m)
            west_dc = flux_w*(mom_lo-mom_ho)


! North face (i, j+1, k)
            IF(v(ijk) >= zero)THEN
               mom_lo = v_s(ijk,m)
               IF (fpfoi) mom_ho = fpfoi_of(v_s(ijpk,m), v_s(ijk,m), &
                                   v_s(ijmk,m), v_s(jm_of(ijmk),m))
            ELSE
               mom_lo = v_s(ijpk,m)
               IF (fpfoi) mom_ho = fpfoi_of(v_s(ijk,m), v_s(ijpk,m), &
                                   v_s(jp_of(ijpk),m), tmp4(jppp4))
            ENDIF
            IF (.NOT. fpfoi) mom_ho = xsi_n(ijk)*v_s(ijpk,m)+ &
                                      (1.0-xsi_n(ijk))*v_s(ijk,m)
            north_dc = flux_n*(mom_lo-mom_ho)

! South face (i, j, k)
            IF(v(ijmk) >= zero)THEN
               mom_lo = v_s(ijmk,m)
               IF (fpfoi) mom_ho = fpfoi_of(v_s(ijk,m), v_s(ijmk,m), &
                                   v_s(jm_of(ijmk),m), tmp4(jmmm4))
            ELSE
               mom_lo = v_s(ijk,m)
               IF (fpfoi) mom_ho = fpfoi_of(v_s(ijmk,m), v_s(ijk,m), &
                                   v_s(ijpk,m), v_s(jp_of(ijpk),m))
            ENDIF
            IF (.NOT. fpfoi) mom_ho = xsi_n(ijmk)*v_s(ijk,m)+ &
                                      (1.0-xsi_n(ijmk))*v_s(ijmk,m)
            south_dc = flux_s*(mom_lo-mom_ho)


            IF (do_k) THEN

! Top face (i, j+1/2, k+1/2)
               IF(ww(ijk) >= zero)THEN
                  mom_lo = v_s(ijk,m)
                  IF (fpfoi) mom_ho = fpfoi_of(v_s(ijkp,m), v_s(ijk,m), &
                                      v_s(ijkm,m), v_s(km_of(ijkm),m))
               ELSE
                  mom_lo = v_s(ijkp,m)
                  IF (fpfoi) mom_ho = fpfoi_of(v_s(ijk,m), v_s(ijkp,m), &
                                      v_s(kp_of(ijkp),m), tmp4(kppp4))
               ENDIF
               IF (.NOT. fpfoi) mom_ho = xsi_t(ijk)*v_s(ijkp,m)+ &
                                         (1.0-xsi_t(ijk))*v_s(ijk,m)
               top_dc = flux_t*(mom_lo-mom_ho)

! Bottom face (i, j+1/2, k-1/2)
               IF(ww(ijk) >= zero)THEN
                  mom_lo = v_s(ijkm,m)
                  IF (fpfoi) mom_ho = fpfoi_of(v_s(ijk,m), v_s(ijkm,m), &
                                      v_s(km_of(ijkm),m), tmp4(kmmm4))
               ELSE
                  mom_lo = v_s(ijk,m)
                  IF (fpfoi) mom_ho = fpfoi_of(v_s(ijkm,m), v_s(ijk,m), &
                                      v_s(ijkp,m), v_s(kp_of(ijkp),m))
               ENDIF
               IF (.NOT. fpfoi) mom_ho = xsi_t(ijkm)*v_s(ijk,m)+ &
                                         (1.0-xsi_t(ijkm))*v_s(ijkm,m)
               bottom_dc = flux_b*(mom_lo-mom_ho)
            ELSE
               top_dc = zero
               bottom_dc = zero
            ENDIF   ! end if do_k


! CONTRIBUTION DUE TO DEFERRED CORRECTION
            b_m(ijk,m) = b_m(ijk,m)+west_dc-east_dc+south_dc-north_dc+&
                         bottom_dc-top_dc

         ENDIF   ! end if flow_at_n
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE store_a_v_sdc

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_V_s1                                            C
!  Purpose: Determine convection diffusion terms for V_s momentum eqs  C
!           The off-diagonal coefficients calculated here must be      C
!           positive. The center coefficient and the source vector     C
!           are negative. Higher order                                 C
!  See source_v_s                                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 20-MAR-97  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number: 1                                                  C
!  Purpose: To incorporate Cartesian grid modifications                C
!  Author: Jeff Dietiker                              Date: 01-Jul-09  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE store_a_v_s1(A_V_S, M)

! Modules
!---------------------------------------------------------------------//
      USE compar, only: ijkstart3, ijkend3
      USE fldvar, only: v_s

      USE functions, only: flow_at_n
      USE functions, only: ip_of, jp_of, kp_of
      USE functions, only: im_of, jm_of, km_of

      USE geometry, only: do_k

      USE param
      USE param1, only: one

      USE run, only: discretize
      USE xsi, only: calc_xsi
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! Septadiagonal matrix A_V_s
      DOUBLE PRECISION, INTENT(INOUT) :: A_V_s(dimension_3, -3:3, m:m)

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK, IPJK, IMJK, IJPK, IJMK, IJKP, IJKM
! indicator for shear
      INTEGER :: incr
! Diffusion parameter
      DOUBLE PRECISION :: d_fe, d_fw, d_fn, d_fs, d_ft, d_fb
! Face mass flux
      DOUBLE PRECISION :: Flux_e, flux_w, flux_n, flux_s
      DOUBLE PRECISION :: flux_t, flux_b

! temporary use of global arrays:
! array1 (locally u)  - the x directional velocity
      DOUBLE PRECISION :: U(dimension_3)
! array2 (locally v)  - the y directional velocity
      DOUBLE PRECISION :: V(dimension_3)
! array3 (locally ww) - the z directional velocity
      DOUBLE PRECISION :: WW(dimension_3)
!---------------------------------------------------------------------//
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TMP4
      DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t

      CALL get_vcell_svterms(u, v, ww,m)

! shear indicator:
      incr=2
      CALL calc_xsi (discretize(4), v_s(1,m), u, v, ww, xsi_e, xsi_n, &
                     xsi_t, incr)

!!!$omp      parallel do                                             &
!!!$omp&     private(IJK, IMJK, IJMK, IPJK, IJPK, IJKM, IJKP,        &
!!!$omp&             d_fe, d_fw, d_fn, d_fs, d_ft, d_fb,             &
!!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b)
      DO ijk = ijkstart3, ijkend3

         IF (flow_at_n(ijk)) THEN

! Calculate convection-diffusion fluxes through each of the faces
            CALL get_vcell_scflux_terms(flux_e, flux_w, flux_n, &
               flux_s, flux_t, flux_b, ijk, m)
            CALL get_vcell_sdiff_terms(d_fe, d_fw, d_fn, d_fs, &
               d_ft, d_fb, ijk, m)

            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ipjk = ip_of(ijk)
            ijpk = jp_of(ijk)

! East face (i+1/2, j+1/2, k)
            a_v_s(ijk,east,m) = d_fe - xsi_e(ijk)*flux_e
            a_v_s(ipjk,west,m) = d_fe + (one - xsi_e(ijk))*flux_e
! West face (i-1/2, j+1/2, k)
            IF (.NOT.flow_at_n(imjk)) THEN
               a_v_s(ijk,west,m) = d_fw + (one - xsi_e(imjk))*flux_w
            ENDIF


! North face (i, j+1, k)
            a_v_s(ijk,north,m) = d_fn - xsi_n(ijk)*flux_n
            a_v_s(ijpk,south,m) = d_fn + (one - xsi_n(ijk))*flux_n
! South face (i, j, k)
            IF (.NOT.flow_at_n(ijmk)) THEN
               a_v_s(ijk,south,m) = d_fs + (one - xsi_n(ijmk))*flux_s
            ENDIF


            IF (do_k) THEN
               ijkp = kp_of(ijk)
               ijkm = km_of(ijk)
! Top face (i, j+1/2, k+1/2)
               a_v_s(ijk,top,m) = d_ft - xsi_t(ijk)*flux_t
               a_v_s(ijkp,bottom,m) = d_ft + (one - xsi_t(ijk))*flux_t
! Bottom face (i, j+1/2, k-1/2)
               IF (.NOT.flow_at_n(ijkm)) THEN
                  a_v_s(ijk,bottom,m) = d_fb + (one - xsi_t(ijkm))*flux_b
               ENDIF
            ENDIF   ! end if (do_k)

         ENDIF   ! end if flow_at_n
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE store_a_v_s1