conv_dif_u_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_DIF_U_s                                            C
!  Purpose: Determine convection diffusion terms for U_s momentum eqs  C
!  The off-diagonal coefficients calculated here must be posotive.     C
!  The center coefficient and the source vector are negative.          C
!  See source_u_s                                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 24-DEC-96  C
!  Reviewer:                                          Date:            C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE conv_dif_u_s(A_M, B_M)

! 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_x_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)

! 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_x_eq(m)) THEN

             IF(def_cor)THEN
! USE DEFERRED CORRECTION TO SOLVE U_S
                CALL store_a_u_s0 (a_m(1,-3,m), m)
                IF (discretize(3) > 1) CALL store_a_u_sdc (m, b_m)

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

            CALL dif_u_is (mu_s(1,m), a_m, m)
          ENDIF
        ENDIF
      ENDDO

      RETURN
      END SUBROUTINE conv_dif_u_s

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

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

      USE cutcell, only: cut_u_treatment_at
      USE cutcell, only: theta_ue, theta_ue_bar
      USE cutcell, only: theta_u_ne, theta_u_nw
      USE cutcell, only: theta_u_te, theta_u_tw
      USE cutcell, only: alpha_ue_c, alpha_un_c, alpha_ut_c

      USE fldvar, only: u_s, v_s, w_s

      USE fun_avg, only: avg_x_e, avg_x
      USE functions, only: ip_of
      USE geometry, only: do_k
      USE indices, only: i_of, ip1

      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, I, IP, IPJK
! for cartesian grid
      DOUBLE PRECISION :: AW, HW, VELW
!---------------------------------------------------------------------//

!!!$omp parallel do private(IJK,I,IP,IPJK)
      DO ijk = ijkstart3, ijkend3
         i = i_of(ijk)
         ip = ip1(i)
         ipjk = ip_of(ijk)

         IF(cut_u_treatment_at(ijk)) THEN

! East face (i+1, j, k)
            u(ijk) = (theta_ue_bar(ijk) * u_s(ijk,m) + &
                      theta_ue(ijk) * u_s(ipjk,m))
            CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
               alpha_ue_c(ijk), aw, hw, velw)
            u(ijk) = u(ijk) * aw

! North face (i+1/2, j+1/2, k)
            v(ijk) = (theta_u_nw(ijk) * v_s(ijk,m) + &
                      theta_u_ne(ijk) * v_s(ipjk,m))
            CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
               alpha_un_c(ijk), aw, hw, velw)
            v(ijk) = v(ijk) * aw

! Top face (i+1/2, j, k+1/2)
            IF (do_k) THEN
               ww(ijk) = (theta_u_tw(ijk) * w_s(ijk,m) + &
                          theta_u_te(ijk) * w_s(ipjk,m))
               CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
                  alpha_ut_c(ijk), aw, hw, velw)
               ww(ijk) = ww(ijk) * aw
            ENDIF

         ELSE
            u(ijk) = avg_x_e(u_s(ijk,m),u_s(ipjk,m),ip)
            v(ijk) = avg_x(v_s(ijk,m),v_s(ipjk,m),i)
            IF (do_k) ww(ijk) = avg_x(w_s(ijk,m),w_s(ipjk,m),i)
         ENDIF
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE get_ucell_svterms


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

! Modules
!---------------------------------------------------------------------//
      USE cutcell, only: cut_u_treatment_at
      USE cutcell, only: theta_ue, theta_ue_bar
      USE cutcell, only: theta_u_ne, theta_u_nw
      USE cutcell, only: theta_u_te, theta_u_tw
      USE cutcell, only: alpha_ue_c, alpha_un_c, alpha_ut_c
      USE functions, only: ip_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
! ijk index
      INTEGER, INTENT(IN) :: ijk

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: imjk, ijmk, ijkm
      INTEGER :: ipjk, ipjmk, ipjkm

! for cartesian grid
      DOUBLE PRECISION :: AW, HW, VELW
!---------------------------------------------------------------------//

! indices
      ipjk = ip_of(ijk)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ipjmk = ip_of(ijmk)
      ipjkm = ip_of(ijkm)

! First calculate the fluxes at the faces
      IF(cut_u_treatment_at(ijk)) THEN
! east face: i+1, j, k
         flux_e = (theta_ue_bar(ijk) * flux_se(ijk,m) + &
                   theta_ue(ijk) * flux_se(ipjk,m))
         CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
            alpha_ue_c(ijk), aw, hw, velw)
         flux_e = flux_e * aw
! west face: i-1, j, k
         flux_w = (theta_ue_bar(imjk) * flux_se(imjk,m) + &
                   theta_ue(imjk) * flux_se(ijk,m))
         CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
            alpha_ue_c(imjk), aw, hw, velw)
         flux_w = flux_w * aw


! north face: i+1/2, j+1/2, k
         flux_n = (theta_u_nw(ijk) * flux_sn(ijk,m) + &
                   theta_u_ne(ijk) * flux_sn(ipjk,m))
         CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM', &
            alpha_un_c(ijk), aw, hw, velw)
         flux_n = flux_n * aw
! south face: i+1/2, j-1/2, k
         flux_s = (theta_u_nw(ijmk) * flux_sn(ijmk,m) + &
                   theta_u_ne(ijmk) * flux_sn(ipjmk,m))
         CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
            alpha_un_c(ijmk), aw, hw, velw)
         flux_s = flux_s * aw

         IF (do_k) THEN
! top face: i+1/2, j, k+1/2
            flux_t = (theta_u_tw(ijk) * flux_st(ijk,m) + &
                     theta_u_te(ijk) * flux_st(ipjk,m))
            CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM', &
               alpha_ut_c(ijk), aw, hw, velw)
            flux_t = flux_t * aw
! bottom face: i+1/2, j, k-1/2
            flux_b = (theta_u_tw(ijkm) * flux_st(ijkm,m) + &
                      theta_u_te(ijkm) * flux_st(ipjkm,m))
            CALL get_interpolation_terms_s(ijk, m, 'U_MOMENTUM',&
               alpha_ut_c(ijkm), aw, hw, velw)
            flux_b = flux_b * aw
         ENDIF
      ELSE
         flux_e = half * (flux_se(ijk,m)  + flux_se(ipjk,m))
         flux_w = half * (flux_se(imjk,m) + flux_se(ijk,m))
         flux_n = half * (flux_sn(ijk,m)  + flux_sn(ipjk,m))
         flux_s = half * (flux_sn(ijmk,m) + flux_sn(ipjmk,m))

         IF (do_k) THEN
            flux_t = half * (flux_st(ijk,m)  + flux_st(ipjk,m))
            flux_b = half * (flux_st(ijkm,m) + flux_st(ipjkm,m))
         ENDIF
      ENDIF   ! end if/else cut_u_treatment_at
      RETURN
      END SUBROUTINE get_ucell_scflux_terms


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

! Modules
!---------------------------------------------------------------------//
      USE cutcell, only: cut_u_treatment_at
      USE cutcell, only: oneodx_e_u, oneody_n_u, oneodz_t_u

      USE functions, only: wall_at
      USE functions, only: east_of, north_of, top_of
      USE functions, only: south_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, ody_n, odz_t
      USE geometry, only: do_k
      USE geometry, only: ox_e
      USE geometry, only: ayz_u, axz_u, axy_u

      USE indices, only: i_of, j_of, k_of
      USE indices, only: ip1, jm1, km1

      USE visc_s, only: mu_s
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! diffusion through faces of given ijk u-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, ip, jm, km
      INTEGER :: ijkc, ijke, ijkn, ijkne, ijks, ijkse
      INTEGER :: ijkt, ijkte, ijkb, ijkbe

! 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)
      ip = ip1(i)
      jm = jm1(j)
      km = km1(k)

      ijke = east_of(ijk)
      IF (wall_at(ijk)) THEN
         ijkc = ijke
      ELSE
         ijkc = ijk
      ENDIF
      ijkn = north_of(ijk)
      ijkne = east_of(ijkn)
      ijks = south_of(ijk)
      ijkse = east_of(ijks)

      IF(cut_u_treatment_at(ijk)) THEN
         c_ae = oneodx_e_u(ijk)
         c_aw = oneodx_e_u(imjk)
         c_an = oneody_n_u(ijk)
         c_as = oneody_n_u(ijmk)
         c_at = oneodz_t_u(ijk)
         c_ab = oneodz_t_u(ijkm)
      ELSE
         c_ae = odx(ip)
         c_aw = odx(i)
         c_an = ody_n(j)
         c_as = ody_n(jm)
         c_at = odz_t(k)
         c_ab = odz_t(km)
      ENDIF

! East face (i+1, j, k)
      d_fe = mu_s(ijke,m)*c_ae*ayz_u(ijk)
! West face (i-1, j, k)
      d_fw = mu_s(ijkc,m)*c_aw*ayz_u(imjk)


! North face (i+1/2, j+1/2, k)
      d_fn = avg_x_h(avg_y_h(mu_s(ijkc,m),mu_s(ijkn,m),j),&
                     avg_y_h(mu_s(ijke,m),mu_s(ijkne,m),j),i)*&
             c_an*axz_u(ijk)
! South face (i+1/2, j-1/2, k)
      d_fs = avg_x_h(avg_y_h(mu_s(ijks,m),mu_s(ijkc,m),jm),&
                     avg_y_h(mu_s(ijkse,m),mu_s(ijke,m),jm),i)*&
             c_as*axz_u(ijmk)

      IF (do_k) THEN
         ijkt = top_of(ijk)
         ijkte = east_of(ijkt)
         ijkb = bottom_of(ijk)
         ijkbe = east_of(ijkb)

! Top face (i+1/2, j, k+1/2)
         d_ft = avg_x_h(avg_z_h(mu_s(ijkc,m),mu_s(ijkt,m),k),&
                        avg_z_h(mu_s(ijke,m),mu_s(ijkte,m),k),i)*&
                ox_e(i)*c_at*axy_u(ijk)
! Bottom face (i+1/2, j, k-1/2)
         d_fb = avg_x_h(avg_z_h(mu_s(ijkb,m),mu_s(ijkc,m),km),&
                        avg_z_h(mu_s(ijkbe,m),mu_s(ijke,m),km),i)*&
                ox_e(i)*c_ab*axy_u(ijkm)
      ENDIF
      RETURN
      END SUBROUTINE get_ucell_sdiff_terms


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_U_s0                                            C
!  Purpose: Determine convection diffusion terms for U_s momentum eqs. C
!  The off-diagonal coefficients calculated here must be positive.     C
!  The center coefficient and the source vector are negative. See      C
!  source_u_s.                                                         C
!  Implement FOUP discretization.                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 29-APR-96  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_u_s0(A_U_S, M)

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

      USE functions, only: flow_at_e
      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
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! Septadiagonal matrix A_U_s
      DOUBLE PRECISION, INTENT(INOUT) :: A_U_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, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,         &
!$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)  &
!$omp     shared(ijkstart3, ijkend3, do_k, a_u_s, m)

      DO ijk = ijkstart3, ijkend3

         IF (flow_at_e(ijk)) THEN

! Calculate convection-diffusion fluxes through each of the faces
            CALL get_ucell_scflux_terms(flux_e, flux_w, flux_n, &
               flux_s, flux_t, flux_b, ijk, m)
            CALL get_ucell_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, j, k)
            IF (flux_e >= zero) THEN
               a_u_s(ijk,east,m) = d_fe
               a_u_s(ipjk,west,m) = d_fe + flux_e
            ELSE
               a_u_s(ijk,east,m) = d_fe - flux_e
               a_u_s(ipjk,west,m) = d_fe
            ENDIF
! West face (i, j, k)
            IF (.NOT.flow_at_e(imjk)) THEN
               IF (flux_w >= zero) THEN
                  a_u_s(ijk,west,m) = d_fw + flux_w
               ELSE
                  a_u_s(ijk,west,m) = d_fw
               ENDIF
            ENDIF


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


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

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

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


      RETURN
      END SUBROUTINE store_a_u_s0



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_U_sdc                                           C
!  Purpose: Use deferred correction method to solve the u-momentum     C
!  equation. This method combines first order upwind and a user        C
!  specified higher order method                                       C
!                                                                      C
!  Author: C. GUENTHER                                Date: 8-APR-99   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_u_sdc(M, B_M)

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

      USE discretization, only: fpfoi_of

      USE fldvar, only: u_s

      USE function3, only: funijk3
      USE functions, only: flow_at_e
      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
!---------------------------------------------------------------------//
! 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 :: IMJK, IPJK, 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 correction contribution from 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 contribution 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_ucell_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) = u_s(ijk,m)
         ENDDO
         CALL send_recv3(tmp4)
      ENDIF

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

!!!$omp      parallel do                                              &
!!!$omp&     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,         &
!!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b,  &
!!!$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_e(ijk)) THEN

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

            ipjk = ip_of(ijk)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)
            ijpk = jp_of(ijk)
            ijkp = kp_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, j, k)
            IF(u(ijk) >= zero)THEN
               mom_lo = u_s(ijk,m)
               IF (fpfoi) mom_ho = fpfoi_of(u_s(ipjk,m), u_s(ijk,m), &
                                   u_s(imjk,m), u_s(im_of(imjk),m))
            ELSE
               mom_lo = u_s(ipjk,m)
               IF (fpfoi) mom_ho = fpfoi_of(u_s(ijk,m), u_s(ipjk,m), &
                                   u_s(ip_of(ipjk),m), tmp4(ippp4))
            ENDIF

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


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

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


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

! Bottom face (i+1/2, j, k-1/2)
               IF(ww(ijk) >= zero)THEN
                  mom_lo = u_s(ijkm,m)
                  IF (fpfoi) mom_ho = fpfoi_of(u_s(ijk,m), u_s(ijkm,m), &
                                      u_s(km_of(ijkm),m), tmp4(kmmm4))
               ELSE
                  mom_lo = u_s(ijk,m)
                  IF (fpfoi) mom_ho = fpfoi_of(u_s(ijkm,m), u_s(ijk,m), &
                                      u_s(ijkp,m), u_s(kp_of(ijkp),m))
               ENDIF
               IF (.NOT. fpfoi) mom_ho = xsi_t(ijkm)*u_s(ijk,m)+ &
                                         (1.0-xsi_t(ijkm))*u_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_e
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE store_a_u_sdc

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_U_s1                                            C
!  Purpose: Determine convection diffusion terms for U_s momentum eqs  C
!  The off-diagonal coefficients calculated here must be positive.     C
!  The center coefficient and the source vector are negative.          C
!  Implements higher order discretization.                             C
!  See source_u_s                                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 20-MAR-97  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_u_s1(A_U_S, M)
! Modules
!---------------------------------------------------------------------//
      USE compar, only: ijkstart3, ijkend3
      USE fldvar, only: u_s

      USE functions, only: flow_at_e
      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
!---------------------------------------------------------------------//
! Solids phase
      INTEGER, INTENT(IN) :: M
! Septadiagonal matrix A_U_s
      DOUBLE PRECISION, INTENT(INOUT) :: A_U_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_ucell_svterms(u, v, ww, m)

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

!!!$omp      parallel do                                                 &
!!!$omp&     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,            &
!!!$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_e(ijk)) THEN

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

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

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


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


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

         ENDIF   ! end if flow_at_e
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE store_a_u_s1