conv_dif_w_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_DIF_W_s                                            C
!  Purpose: Determine convection diffusion terms for W_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_w_s                                             C
!                                                                      C
!  Author: M. Syamlal                                 Date: 24-DEC-96  C
!  Reviewer:                                          Date:            C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^CC
      SUBROUTINE conv_dif_w_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_z_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_z_eq(m)) THEN

             IF (def_cor) THEN
! USE DEFERRED CORRECTION TO SOLVE W_S
                CALL store_a_w_s0 (a_m(1,-3,m), m)
                IF (discretize(5) > 1)CALL store_a_w_sdc (m, b_m)

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

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

      RETURN
      END SUBROUTINE conv_dif_w_s

!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_svterms(U, V, WW, M)

! 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_s, v_s, w_s

      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
!---------------------------------------------------------------------//
! 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, 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_s(ijk,m) + &
                      theta_w_te(ijk) * u_s(ijkp,m))
            CALL get_interpolation_terms_s(ijk, m, '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_s(ijk,m) + &
                      theta_w_tn(ijk) * v_s(ijkp,m))
            CALL get_interpolation_terms_s(ijk, m, '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_s(ijk,m) + &
                      theta_wt(ijk) * w_s(ijkp,m))
            CALL get_interpolation_terms_s(ijk, m, 'W_MOMENTUM', &
               alpha_wt_c(ijk), aw, hw, velw)
            ww(ijk) = ww(ijk) * aw

         ELSE
            u(ijk) = avg_z(u_s(ijk,m),u_s(ijkp,m),k)
            v(ijk) = avg_z(v_s(ijk,m),v_s(ijkp,m),k)
            ww(ijk) = avg_z_t(w_s(ijk,m),w_s(ijkp,m))
         ENDIF   ! end if/else cut_w_treatment_at
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE get_wcell_svterms


!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_scflux_terms(FLUX_E, FLUX_W, FLUX_N, &
         flux_s, flux_t, flux_b, ijk, m)

! 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_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 :: 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_se(ijk,m) + &
                 theta_w_te(ijk) * flux_se(ijkp,m))
         CALL get_interpolation_terms_s(ijk, m, '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_se(imjk,m) + &
                   theta_w_te(imjk) * flux_se(imjkp,m))
         CALL get_interpolation_terms_s(ijk, m, '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_sn(ijk,m) + &
                 theta_w_tn(ijk) * flux_sn(ijkp,m))
         CALL get_interpolation_terms_s(ijk, m, '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_sn(ijmk,m) + &
                   theta_w_tn(ijmk) * flux_sn(ijmkp,m))
         CALL get_interpolation_terms_s(ijk, m, '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_st(ijk,m) + &
                   theta_wt(ijk) * flux_st(ijkp,m))
         CALL get_interpolation_terms_s(ijk, m, '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_st(ijkm,m) + &
                   theta_wt(ijkm) * flux_st(ijk,m))
         CALL get_interpolation_terms_s(ijk, m, 'W_MOMENTUM',&
            alpha_wt_c(ijkm), aw, hw, velw)
         flux_b = flux_b * aw
      ELSE
         flux_e = half * (flux_se(ijk,m)  + flux_se(ijkp,m))
         flux_w = half * (flux_se(imjk,m) + flux_se(imjkp,m))
         flux_n = half * (flux_sn(ijk,m)  + flux_sn(ijkp,m))
         flux_s = half * (flux_sn(ijmk,m) + flux_sn(ijmkp,m))
         flux_t = half * (flux_st(ijk,m)  + flux_st(ijkp,m))
         flux_b = half * (flux_st(ijkm,m) + flux_st(ijk,m))
      ENDIF   ! end if/else cut_w_treatment_at

      RETURN
      END SUBROUTINE get_wcell_scflux_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_sdiff_terms(D_FE, D_FW, D_FN, D_FS, &
         d_ft, d_fb, ijk, m)

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

      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 fun_avg, only: avg_x_h, avg_y_h, avg_z_h

      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_s, only: mu_s
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! 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
!---------------------------------------------------------------------//

      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(mu_s(ijkc,m),mu_s(ijke,m),i),&
                     avg_x_h(mu_s(ijkt,m),mu_s(ijkte,m),i),k)*&
               c_ae*ayz_w(ijk)
! West face (i-1/2, j, k+1/2)
      d_fw = avg_z_h(avg_x_h(mu_s(ijkw,m),mu_s(ijkc,m),im),&
                     avg_x_h(mu_s(ijkwt,m),mu_s(ijkt,m),im),k)*&
                c_aw*ayz_w(imjk)

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


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

      RETURN
      END SUBROUTINE get_wcell_sdiff_terms


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_W_s0                                            C
!  Purpose: Determine convection diffusion terms for W_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_w_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_w_s0(A_W_S, M)

! 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
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER, INTENT(IN) :: M
! Septadiagonal matrix A_W_s
      DOUBLE PRECISION, INTENT(INOUT) :: A_W_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, m, a_w_s)

      DO ijk = ijkstart3, ijkend3

         IF (flow_at_t(ijk)) THEN

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

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

            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_s(ijk,east,m) = d_fe
               a_w_s(ipjk,west,m) = d_fe + flux_e
            ELSE
               a_w_s(ijk,east,m) = d_fe - flux_e
               a_w_s(ipjk,west,m) = 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_s(ijk,west,m) = d_fw + flux_w
               ELSE
                  a_w_s(ijk,west,m) = d_fw
               ENDIF
            ENDIF


! North face (i, j+1/2, k+1/2)
            IF (flux_n >= zero) THEN
               a_w_s(ijk,north,m) = d_fn
               a_w_s(ijpk,south,m) = d_fn + flux_n
            ELSE
               a_w_s(ijk,north,m) = d_fn - flux_n
               a_w_s(ijpk,south,m) = 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_s(ijk,south,m) = d_fs + flux_s
               ELSE
                  a_w_s(ijk,south,m) = d_fs
               ENDIF
            ENDIF


! Top face (i, j, k+1)
            IF (flux_t >= zero) THEN
               a_w_s(ijk,top,m) = d_ft
               a_w_s(ijkp,bottom,m) = d_ft + flux_t
            ELSE
               a_w_s(ijk,top,m) = d_ft - flux_t
               a_w_s(ijkp,bottom,m) = d_ft
            ENDIF
! Bottom face (i, j, k)
            IF (.NOT.flow_at_t(ijkm)) THEN
               IF (flux_b >= zero) THEN
                  a_w_s(ijk,bottom,m) = d_fb + flux_b
               ELSE
                  a_w_s(ijk,bottom,m) = d_fb
               ENDIF
            ENDIF

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

      RETURN
      END SUBROUTINE store_a_w_s0


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_W_sdc                                           C
!  Purpose: To use deferred correction method to solve the w-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_w_sdc(M, B_M)

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

      USE discretization, only: fpfoi_of

      USE fldvar, only: w_s

      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_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 :: 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_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) = w_s(ijk,m)
         ENDDO
         CALL send_recv3(tmp4)
      ENDIF

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

! Calculate convection fluxes through each of the faces
            CALL get_wcell_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, k+1/2)
            IF(u(ijk) >= zero)THEN
               mom_lo = w_s(ijk,m)
               IF (fpfoi) mom_ho = fpfoi_of(w_s(ipjk,m), w_s(ijk,m), &
                                   w_s(imjk,m), w_s(im_of(imjk),m))
            ELSE
               mom_lo = w_s(ipjk,m)
               IF (fpfoi) mom_ho = fpfoi_of(w_s(ijk,m), w_s(ipjk,m), &
                                   w_s(ip_of(ipjk),m), tmp4(ippp4))
            ENDIF
            IF (.NOT. fpfoi) mom_ho = xsi_e(ijk)*w_s(ipjk,m)+ &
                                      (1.0-xsi_e(ijk))*w_s(ijk,m)
            east_dc = flux_e*(mom_lo-mom_ho)

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


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

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


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

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


! 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_t
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE store_a_w_sdc

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STORE_A_W_s1                                            C
!  Purpose: Determine convection diffusion terms for W_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_w_s                                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 19-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_w_s1(A_W_S, M)

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

      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
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! Septadiagonal matrix A_W_s
      DOUBLE PRECISION, INTENT(INOUT) :: A_W_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_wcell_svterms(u, v, ww, m)


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

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

            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_s(ijk,east,m) = d_fe - xsi_e(ijk)*flux_e
            a_w_s(ipjk,west,m) = 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_s(ijk,west,m) = d_fw + (one - xsi_e(imjk))*flux_w
            ENDIF


! North face (i, j+1/2, k+1/2)
            a_w_s(ijk,north,m) = d_fn - xsi_n(ijk)*flux_n
            a_w_s(ijpk,south,m) = 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_s(ijk,south,m) = d_fs + (one - xsi_n(ijmk))*flux_s
            ENDIF


! Top face (i, j, k+1)
            a_w_s(ijk,top,m) = d_ft - xsi_t(ijk)*flux_t
            a_w_s(ijkp,bottom,m) = d_ft + (one - xsi_t(ijk))*flux_t
! Bottom face (i, j, k)
            IF (.NOT.flow_at_t(ijkm)) THEN
              a_w_s(ijk,bottom,m) = 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_s1