conv_dif_w_g.f

Go to the documentation of this file.

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

! Modules
!---------------------------------------------------------------------//
      USE param, only: dimension_3, dimension_m
      USE run, only: momentum_z_eq
      USE run, only: discretize
      USE run, only: def_cor
      USE visc_g, only: epmu_gt
      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)
!---------------------------------------------------------------------//

      IF (.NOT.momentum_z_eq(0)) RETURN

      IF (def_cor) THEN
! USE DEFERRED CORRECTION TO SOLVE W_G
         CALL store_a_w_g0 (a_m)
         IF (discretize(5) > 1) CALL store_a_w_gdc(b_m(1,0))

      ELSE
! DO NOT USE DEFERRED CORRECTOIN TO SOLVE FOR W_G
         IF (discretize(5) == 0) THEN               ! 0 & 1 => FOUP
            CALL store_a_w_g0 (a_m)
         ELSE
            CALL store_a_w_g1 (a_m)
         ENDIF
      ENDIF

      CALL dif_w_is(epmu_gt, a_m, 0)

      RETURN
      END SUBROUTINE conv_dif_w_g


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

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

      USE cutcell, only: cut_w_treatment_at
      USE cutcell, only: theta_wt, theta_wt_bar
      USE cutcell, only: theta_w_be, theta_w_te
      USE cutcell, only: theta_w_bn, theta_w_tn
      USE cutcell, only: alpha_we_c, alpha_wn_c, alpha_wt_c

      USE fldvar, only: u_g, v_g, w_g

      USE fun_avg, only: avg_z_t, avg_z
      USE functions, only: kp_of
      USE indices, only: k_of

      USE param, only: dimension_3
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
      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, K, IJKP
! for cartesian grid
      DOUBLE PRECISION :: AW, HW, VELW
!---------------------------------------------------------------------//


!!!$omp parallel do private(IJK,K,IJKT,IJKP)
      DO ijk = ijkstart3, ijkend3
         k = k_of(ijk)
         ijkp = kp_of(ijk)

         IF(cut_w_treatment_at(ijk)) THEN

! East face (i+1/2, j, k+1/2)
            u(ijk) = (theta_w_be(ijk) * u_g(ijk) + &
                      theta_w_te(ijk) * u_g(ijkp))
            CALL get_interpolation_terms_g(ijk,'W_MOMENTUM', &
               alpha_we_c(ijk), aw, hw, velw)
            u(ijk) = u(ijk) * aw

! North face (i, j+1/2, k+1/2)
            v(ijk) = (theta_w_bn(ijk) * v_g(ijk) + &
                      theta_w_tn(ijk) * v_g(ijkp))
            CALL get_interpolation_terms_g(ijk,'W_MOMENTUM', &
               alpha_wn_c(ijk), aw, hw, velw)
            v(ijk) = v(ijk) * aw

! Top face (i, j, k+1)
            ww(ijk) = (theta_wt_bar(ijk) * w_g(ijk) + &
                       theta_wt(ijk) * w_g(ijkp))
            CALL get_interpolation_terms_g(ijk,'W_MOMENTUM', &
               alpha_wt_c(ijk), aw, hw, velw)
            ww(ijk) = ww(ijk) * aw

         ELSE
            u(ijk) = avg_z(u_g(ijk),u_g(ijkp),k)
            v(ijk) = avg_z(v_g(ijk),v_g(ijkp),k)
            ww(ijk) = avg_z_t(w_g(ijk),w_g(ijkp))
         ENDIF   ! end if/else cut_w_treatment_at
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE get_wcell_gvterms


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

! Modules
!---------------------------------------------------------------------//
      USE cutcell, only: cut_w_treatment_at
      USE cutcell, only: theta_wt, theta_wt_bar
      USE cutcell, only: theta_w_be, theta_w_te
      USE cutcell, only: theta_w_bn, theta_w_tn
      USE cutcell, only: alpha_we_c, alpha_wn_c, alpha_wt_c

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

      USE mflux, only: flux_ge, flux_gn, flux_gt

      USE param1, only: half
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! 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 :: ijkp, imjkp, ijmkp

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

      ijkp = kp_of(ijk)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      imjkp = kp_of(imjk)
      ijmkp = kp_of(ijmk)

      IF(cut_w_treatment_at(ijk)) THEN
! East face (i+1/2, j, k+1/2)
         flux_e = (theta_w_be(ijk) * flux_ge(ijk) + &
                 theta_w_te(ijk) * flux_ge(ijkp))
         CALL get_interpolation_terms_g(ijk,'W_MOMENTUM',&
            alpha_we_c(ijk), aw, hw, velw)
         flux_e = flux_e * aw
! West face (i-1/2, j, k+1/2)
         flux_w = (theta_w_be(imjk) * flux_ge(imjk) + &
                   theta_w_te(imjk) * flux_ge(imjkp))
         CALL get_interpolation_terms_g(ijk,'W_MOMENTUM',&
            alpha_we_c(imjk), aw, hw, velw)
         flux_w = flux_w * aw


! North face (i, j+1/2, k+1/2)
         flux_n = (theta_w_bn(ijk) * flux_gn(ijk) + &
                 theta_w_tn(ijk) * flux_gn(ijkp))
         CALL get_interpolation_terms_g(ijk,'W_MOMENTUM',&
            alpha_wn_c(ijk), aw, hw, velw)
         flux_n = flux_n * aw
! South face (i, j-1/2, k+1/2)
         flux_s = (theta_w_bn(ijmk) * flux_gn(ijmk) + &
                   theta_w_tn(ijmk) * flux_gn(ijmkp))
         CALL get_interpolation_terms_g(ijk,'W_MOMENTUM',&
            alpha_wn_c(ijmk), aw, hw, velw)
         flux_s = flux_s * aw


! Top face (i, j, k+1)
         flux_t = (theta_wt_bar(ijk) * flux_gt(ijk) + &
                   theta_wt(ijk) * flux_gt(ijkp))
         CALL get_interpolation_terms_g(ijk,'W_MOMENTUM',&
            alpha_wt_c(ijk),aw,hw,velw)
         flux_t = flux_t * aw
! Bottom face (i, j, k)
         flux_b = (theta_wt_bar(ijkm) * flux_gt(ijkm) + &
                   theta_wt(ijkm) * flux_gt(ijk))
         CALL get_interpolation_terms_g(ijk,'W_MOMENTUM',&
            alpha_wt_c(ijkm), aw, hw, velw)
         flux_b = flux_b * aw
      ELSE
         flux_e = half * (flux_ge(ijk) + flux_ge(ijkp))
         flux_w = half * (flux_ge(imjk) + flux_ge(imjkp))
         flux_n = half * (flux_gn(ijk) + flux_gn(ijkp))
         flux_s = half * (flux_gn(ijmk) + flux_gn(ijmkp))
         flux_t = half * (flux_gt(ijk) + flux_gt(ijkp))
         flux_b = half * (flux_gt(ijkm) + flux_gt(ijk))
      ENDIF   ! end if/else cut_w_treatment_at

      RETURN
      END SUBROUTINE get_wcell_gcflux_terms


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

! Modules
!---------------------------------------------------------------------//
      USE cutcell, only: cut_w_treatment_at
      USE cutcell, only: oneodx_e_w, oneody_n_w, oneodz_t_w

      USE fldvar, only: epg_jfac

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

      USE geometry, only: odx_e, ody_n, odz
      USE geometry, only: ox
      USE geometry, only: ayz_w, axz_w, axy_w

      USE indices, only: i_of, j_of, k_of
      USE indices, only: kp1, im1, jm1

      USE visc_g, only: epmu_gt, df_gw
      USE param
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! diffusion through faces of given ijk w-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 idnex
      INTEGER, INTENT(IN) :: ijk

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: imjk, ijmk, ijkm
      INTEGER :: i, j, k, kp, im, jm
      INTEGER :: ijkc, ijkt, ijke, ijkte, ijkw, ijkwt
      INTEGER :: ijkn, ijktn, ijks, ijkst
! length terms
      DOUBLE PRECISION :: C_AE, C_AW, C_AN, C_AS, C_AT, C_AB
! average voidage
      DOUBLE PRECISION :: EPGA
!---------------------------------------------------------------------//

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

      i = i_of(ijk)
      j = j_of(ijk)
      k = k_of(ijk)
      kp = kp1(k)
      im = im1(i)
      jm = jm1(j)

      ijkt = top_of(ijk)
      IF (wall_at(ijk)) THEN
         ijkc = ijkt
      ELSE
         ijkc = ijk
      ENDIF
      ijke = east_of(ijk)
      ijkte = east_of(ijkt)
      ijkn = north_of(ijk)
      ijktn = north_of(ijkt)
      ijkw = west_of(ijk)
      ijkwt = top_of(ijkw)
      ijks = south_of(ijk)
      ijkst = top_of(ijks)

      IF(cut_w_treatment_at(ijk)) THEN
         c_ae = oneodx_e_w(ijk)
         c_aw = oneodx_e_w(imjk)
         c_an = oneody_n_w(ijk)
         c_as = oneody_n_w(ijmk)
         c_at = oneodz_t_w(ijk)
         c_ab = oneodz_t_w(ijkm)
      ELSE
         c_ae = odx_e(i)
         c_aw = odx_e(im)
         c_an = ody_n(j)
         c_as = ody_n(jm)
         c_at = odz(kp)
         c_ab = odz(k)
      ENDIF

! East face (i+1/2, j, k+1/2)
      d_fe = avg_z_h(avg_x_h(epmu_gt(ijkc),epmu_gt(ijke),i),&
                     avg_x_h(epmu_gt(ijkt),epmu_gt(ijkte),i),k)*&
               c_ae*ayz_w(ijk)
! West face (i-1/2, j, k+1/2)
      d_fw = avg_z_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijkc),im),&
                     avg_x_h(epmu_gt(ijkwt),epmu_gt(ijkt),im),k)*&
                c_aw*ayz_w(imjk)

! North face (i, j+1/2, k+1/2)
      d_fn = avg_z_h(avg_y_h(epmu_gt(ijkc),epmu_gt(ijkn),j),&
                     avg_y_h(epmu_gt(ijkt),epmu_gt(ijktn),j),k)*&
                c_an*axz_w(ijk)
! South face (i, j-1/2, k+1/2)
      d_fs = avg_z_h(avg_y_h(epmu_gt(ijks),epmu_gt(ijkc),jm),&
                     avg_y_h(epmu_gt(ijkst),epmu_gt(ijkt),jm),k)*&
                c_as*axz_w(ijmk)

! Top face (i, j, k+1)
      d_ft = epmu_gt(ijkt)*ox(i)*c_at*axy_w(ijk)
! Bottom face (i, j, k)
      d_fb = epmu_gt(ijk)*ox(i)*c_ab*axy_w(ijkm)

      df_gw(ijk,east) = d_fe
      df_gw(ijk,west) = d_fw
      df_gw(ijk,north) = d_fn
      df_gw(ijk,south) = d_fs
      df_gw(ijk,top) = d_ft
      df_gw(ijk,bottom) = d_fb

! if jackson, implement jackson style governing equations: multiply by
! the void fraction otherwise multiply by 1
      epga = avg_z(epg_jfac(ijkc), epg_jfac(ijkt), k)
      d_fe = epga*d_fe
      d_fw = epga*d_fw
      d_fn = epga*d_fn
      d_fs = epga*d_fs
      d_ft = epga*d_ft
      d_fb = epga*d_fb

      RETURN

    CONTAINS

      include 'fun_avg.inc'

    END SUBROUTINE get_wcell_gdiff_terms

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_W_g0                                            C
!  Purpose: Determine convection diffusion terms for W_g 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_w_g.                                                         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_w_g0(A_W_G)

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

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

      USE param
      USE param1, only: zero

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Septadiagonal matrix A_U_g
      DOUBLE PRECISION, INTENT(INOUT) :: A_W_g(dimension_3, -3:3, 0:dimension_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, a_w_g)

      DO ijk = ijkstart3, ijkend3

         IF (flow_at_t(ijk)) THEN

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

            CALL get_wcell_gdiff_terms(d_fe, d_fw, d_fn, d_fs, &
               d_ft, d_fb, ijk)

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

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


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


! Top face (i, j, k+1)
            IF (flux_t >= zero) THEN
               a_w_g(ijk,top,0) = d_ft
               a_w_g(ijkp,bottom,0) = d_ft + flux_t
            ELSE
               a_w_g(ijk,top,0) = d_ft - flux_t
               a_w_g(ijkp,bottom,0) = d_ft
            ENDIF
! Bottom face (i, j, k)
            IF (.NOT.flow_at_t(ijkm)) THEN
               IF (flux_b >= zero) THEN
                  a_w_g(ijk,bottom,0) = d_fb + flux_b
               ELSE
                  a_w_g(ijk,bottom,0) = d_fb
               ENDIF
            ENDIF
         ENDIF   ! end if (flow_at_t)
      ENDDO   ! end do (ijk)
!$omp end parallel do

      RETURN
      END SUBROUTINE store_a_w_g0


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_W_GDC                                           C
!  Purpose: Use deferred correction method to solve the w-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_w_gdc(B_M)

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

      USE discretization, only: fpfoi_of

      USE fldvar, only: w_g

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

      USE indices, only: i_of, j_of, k_of

      USE param, only: dimension_3, 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
!---------------------------------------------------------------------//
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3)

! 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_wcell_gvterms(u, v, ww)

! 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) = w_g(ijk)
         ENDDO
         CALL send_recv3(tmp4)
      ENDIF

! shear indicator:
      incr=0
      CALL calc_xsi (discretize(5), w_g, 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_t(ijk)) THEN

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

            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, k+1/2)
            IF(u(ijk) >= zero)THEN
               mom_lo = w_g(ijk)
               IF (fpfoi) mom_ho = fpfoi_of(w_g(ipjk), w_g(ijk), &
                                   w_g(imjk), w_g(im_of(imjk)))
            ELSE
               mom_lo = w_g(ipjk)
               IF (fpfoi) mom_ho = fpfoi_of(w_g(ijk), w_g(ipjk), &
                                   w_g(ip_of(ipjk)), tmp4(ippp4))
            ENDIF
            IF (.NOT. fpfoi) mom_ho = xsi_e(ijk)*w_g(ipjk)+ &
                                      (1.0-xsi_e(ijk))*w_g(ijk)
            east_dc = flux_e*(mom_lo-mom_ho)


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


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


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


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


! Bottom face (i, j, k)
            IF(ww(ijk) >= zero)THEN
               mom_lo = w_g(ijkm)
               IF (fpfoi) mom_ho = fpfoi_of(w_g(ijk), w_g(ijkm), &
                                   w_g(km_of(ijkm)), tmp4(kmmm4))
            ELSE
               mom_lo = w_g(ijk)
               IF (fpfoi) mom_ho = fpfoi_of(w_g(ijkm), w_g(ijk), &
                                   w_g(ijkp), w_g(kp_of(ijkp)))
            ENDIF
            IF (.NOT. fpfoi) mom_ho = xsi_t(ijkm)*w_g(ijk)+ &
                                      (1.0-xsi_t(ijkm))*w_g(ijkm)
            bottom_dc = flux_b*(mom_lo-mom_ho)


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

         ENDIF   ! end if flow_at_t
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE store_a_w_gdc

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_W_g1                                            C
!  Purpose: Determine convection diffusion terms for W_g 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_w_g                                                      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_w_g1(A_W_G)

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

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

      USE param
      USE param1, only: one

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

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Septadiagonal matrix A_W_g
      DOUBLE PRECISION, INTENT(INOUT) :: A_W_g(dimension_3, -3:3, 0:dimension_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_wcell_gvterms(u, v, ww)

! shear indicator:
      incr=0
      CALL calc_xsi (discretize(5), w_g, 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_t(ijk)) THEN

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

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

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


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


! Top face (i, j, k+1)
            a_w_g(ijk,top,0) = d_ft - xsi_t(ijk)*flux_t
            a_w_g(ijkp,bottom,0) = d_ft + (one - xsi_t(ijk))*flux_t
! Bottom face (i, j, k)
            IF (.NOT.flow_at_t(ijkm)) THEN
              a_w_g(ijk,bottom,0) = d_fb + (one - xsi_t(ijkm))*flux_b
            ENDIF

         ENDIF   ! end if flow_at_t
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE store_a_w_g1