solve_vel_star.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_VEL_STAR                                          C
!  Purpose: Solve starred velocity components                          C
!                                                                      C
!  Author: M. Syamlal                                 Date: 25-APR-96  C
!                                                                      C
!  Revision Number: 2                                                  C
!  Purpose: allow multiparticle in D_E, D_N and D_T claculations       C
!           and account for the averaged Solid-Solid drag              C
!  Author: S. Dartevelle, LANL                        Date: 28-FEb-04  C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE solve_vel_star(IER)

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use discretelement, only: des_continuum_hybrid, discrete_element, des_continuum_coupled
      use fldvar, only: u_g, v_g, w_g, u_s, v_s, w_s
      use geometry, only: do_k, ijkmax2
      use leqsol, only: leq_sweep, leq_tol, leq_pc, leq_it, leq_method
      use param, only: dimension_3, dimension_m
      use param1, only: zero, dimension_lm
      use pgcor, only: d_e, d_t, d_n
      use ps, only: point_source
      use pscor, only: e_e, e_t, e_n, mcp
      use physprop,only: mmax
      use residual, only: resid_u, resid_v, resid_w, num_resid, den_resid, resid, max_resid, ijk_resid
      use run, only: kt_type_enum, ghd_2007
      use run, only: momentum_x_eq, momentum_y_eq, momentum_z_eq
      use ur_facs, only: ur_fac
      ! Use ambm
      ! Use tmp_array1,  VxF_gs => Arraym1
      ! Use tmp_array,  VxF_ss => ArrayLM
      USE qmom_kinetic_equation
      use usr_src, only: call_usr_source, calc_usr_source
      use usr_src, only: gas_u_mom, gas_v_mom, gas_w_mom
      use usr_src, only: solids_u_mom, solids_v_mom, solids_w_mom
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Error index
      INTEGER, INTENT(INOUT) :: IER

! Local variables
!---------------------------------------------------------------------//
! fluid cell index
      INTEGER :: IJK
! temporary velocity arrays
      DOUBLE PRECISION, DIMENSION(:), allocatable :: U_gtmp,  V_gtmp, W_gtmp
      DOUBLE PRECISION, DIMENSION(:,:), allocatable :: U_stmp, V_stmp, W_stmp
! linear equation solver method and iterations
      INTEGER :: LEQM, LEQI

      LOGICAL :: DO_SOLIDS

! temporary use of global arrays:
! arraym1 (locally vxf_gs)
! the volume x average gas-solids drag at momentum cell centers
!      DOUBLE PRECISION :: VXF_GS(DIMENSION_3, DIMENSION_M)
! arraylm (locally vxf_ss)
! the volume x average solids-solids drag at momentum cell centers
!      DOUBLE PRECISION :: VXF_SS(DIMENSION_3, DIMENSION_LM)
! Septadiagonal matrix A_m, vector b_m
!      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
!      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
!---------------------------------------------------------------------//

      allocate(u_gtmp(dimension_3))
      allocate(v_gtmp(dimension_3))
      allocate(w_gtmp(dimension_3))
      allocate(u_stmp(dimension_3,dimension_m))
      allocate(v_stmp(dimension_3,dimension_m))
      allocate(w_stmp(dimension_3,dimension_m))

      ! call lock_ambm        ! locks arrys a_m and b_m
      ! call lock_tmp_array1  ! locks arraym1 (locally vxf_gs)
      ! call lock_tmp_array2  ! locks arraylm (locally vxf_ss)

      CALL init(u_gtmp, v_gtmp, w_gtmp, u_stmp, v_stmp, w_stmp)

      do_solids = .NOT.(discrete_element .OR. qmomk) .OR. &
         des_continuum_hybrid

!!$omp parallel sections default(shared)
      CALL u_m_star(u_gtmp,u_stmp)
!!$omp section
      CALL v_m_star(v_gtmp,v_stmp)
!!$omp section
      CALL w_m_star(w_gtmp,w_stmp)
!!$omp end parallel sections

      CALL save(u_gtmp, v_gtmp, w_gtmp, u_stmp, v_stmp, w_stmp)

! modification for GHD theory to compute species velocity: Ui = Joi/(mi ni) + U.
      IF(kt_type_enum == ghd_2007) THEN
         CALL calc_external_forces()
         CALL ghdmassflux() ! to compute solid species mass flux
         CALL updatespeciesvelocities() ! located at end of ghdMassFlux.f file
      ENDIF

      ! call unlock_ambm
      ! call unlock_tmp_array1
      ! call unlock_tmp_array2

      deallocate(u_gtmp)
      deallocate(v_gtmp)
      deallocate(w_gtmp)
      deallocate(u_stmp)
      deallocate(v_stmp)
      deallocate(w_stmp)

      RETURN

    CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE init(U_g_tmp, V_g_tmp, W_g_tmp, U_s_tmp, V_s_tmp, W_s_tmp)
        IMPLICIT NONE
!---------------------------------------------------------------------//
        DOUBLE PRECISION, DIMENSION(:), intent(out) :: U_g_tmp,  V_g_tmp, W_g_tmp
        DOUBLE PRECISION, DIMENSION(:,:), intent(out) :: U_s_tmp, V_s_tmp, W_s_tmp
!---------------------------------------------------------------------//
        ! solids phase index
        INTEGER :: M
!---------------------------------------------------------------------//

        ! Store the velocities so that the order of solving the momentum
        ! equations does not matter
        DO ijk = ijkstart3, ijkend3
           u_g_tmp(ijk) = u_g(ijk)
           v_g_tmp(ijk) = v_g(ijk)
           w_g_tmp(ijk) = w_g(ijk)
        ENDDO
        DO m = 1, mmax
           IF(kt_type_enum /= ghd_2007 .OR. &
                (kt_type_enum == ghd_2007 .AND. m==mmax)) THEN
              DO ijk = ijkstart3, ijkend3
                 u_s_tmp(ijk, m) = u_s(ijk, m)
                 v_s_tmp(ijk, m) = v_s(ijk, m)
                 w_s_tmp(ijk, m) = w_s(ijk, m)
              ENDDO
           ENDIF
        ENDDO

      END SUBROUTINE init

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE save(U_g_tmp, V_g_tmp, W_g_tmp, U_s_tmp, V_s_tmp, W_s_tmp)
        IMPLICIT NONE
!---------------------------------------------------------------------//
        DOUBLE PRECISION, DIMENSION(:), intent(in) :: U_g_tmp,  V_g_tmp, W_g_tmp
        DOUBLE PRECISION, DIMENSION(:,:), intent(in) :: U_s_tmp, V_s_tmp, W_s_tmp
!---------------------------------------------------------------------//
        ! solids phase index
        INTEGER :: M
!---------------------------------------------------------------------//

        ! Now update all velocity components
        DO ijk = ijkstart3, ijkend3
           u_g(ijk) = u_g_tmp(ijk)
           v_g(ijk) = v_g_tmp(ijk)
           w_g(ijk) = w_g_tmp(ijk)
        ENDDO
        DO m = 1, mmax
           IF(kt_type_enum /= ghd_2007 .OR. &
                (kt_type_enum == ghd_2007 .AND. m==mmax)) THEN
              DO ijk = ijkstart3, ijkend3
                 u_s(ijk, m) = u_s_tmp(ijk, m)
                 v_s(ijk, m) = v_s_tmp(ijk, m)
                 w_s(ijk, m) = w_s_tmp(ijk, m)
              ENDDO
           ENDIF
        ENDDO

      END SUBROUTINE save

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE u_m_star(U_g_tmp, U_s_tmp)
        IMPLICIT NONE
!---------------------------------------------------------------------//
        DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: U_g_tmp
        DOUBLE PRECISION, DIMENSION(:, :), INTENT(OUT) :: U_s_tmp
!---------------------------------------------------------------------//
        ! solids phase index
        INTEGER :: M
        DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: VXF_GS, VXF_SS, B_M
        DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: A_M
!---------------------------------------------------------------------//

        allocate(vxf_gs(dimension_3,dimension_m))
        allocate(vxf_ss(dimension_3,dimension_lm))
        ALLOCATE(a_m(dimension_3, -3:3, 0:dimension_m))
        ALLOCATE(b_m(dimension_3, 0:dimension_m))

! Calculate U_m_star and residuals
! ---------------------------------------------------------------->>>
      DO m = 0, mmax
         CALL init_ab_m (a_m, b_m, ijkmax2, m)
         IF (m >= 1) vxf_gs(:,m) = zero
      ENDDO

! calculate the convection-diffusion terms for the gas and solids phase
! u-momentum equations
      CALL conv_dif_u_g (a_m, b_m)
      IF(do_solids) CALL conv_dif_u_s (a_m, b_m)

! calculate the source terms for the gas and solids phase u-momentum
! equations
      CALL source_u_g (a_m, b_m)
      IF(point_source) CALL point_source_u_g (a_m, b_m)
      IF(call_usr_source(3)) CALL calc_usr_source(gas_u_mom, a_m, b_m)
      IF(do_solids) THEN
         CALL source_u_s (a_m, b_m)
         IF(point_source) CALL point_source_u_s (a_m, b_m)
         IF(call_usr_source(3)) CALL calc_usr_source(solids_u_mom, a_m, b_m)
      ENDIF

! evaluate local variable vxf_gs and vxf_ss.  both terms are sent to the
! subroutine calc_d (pressure correction equation coefficients).  the
! former is also used in the subroutine partial_elim_u while the latter
! is effectively re-evaluated within said subroutine
      CALL vf_gs_x (vxf_gs)
      IF(do_solids .AND. (kt_type_enum /= ghd_2007)) THEN
         IF (mmax > 0) CALL vf_ss_x (vxf_ss)
      ENDIF

! calculate coefficients for the pressure correction equation
      IF(kt_type_enum == ghd_2007) THEN
         CALL calc_d_ghd_e (a_m, vxf_gs, d_e)
      ELSE
         CALL calc_d_e (a_m, vxf_gs, vxf_ss, d_e, ier)
      ENDIF

      IF(do_solids) THEN
! calculate coefficients for a solids volume correction equation
         IF (mmax > 0) CALL calc_e_e (a_m, mcp, e_e)

! calculate modifications to the A matrix center coefficient and B
! source vector for partial elimination
         IF(kt_type_enum == ghd_2007) THEN
            IF (mmax > 0) CALL partial_elim_ghd_u (u_g,u_s,vxf_gs,a_m,b_m)
         ELSE
            IF (mmax > 0) CALL partial_elim_u (u_g,u_s,vxf_gs,a_m,b_m)
         ENDIF
      ENDIF

! handle special case where center coefficient is zero
      CALL adjust_a_u_g (a_m, b_m)
      IF(do_solids) CALL adjust_a_u_s (a_m, b_m)

! calculate modifications to the A matrix center coefficient and B
! source vector for treating DEM drag terms
      IF(des_continuum_coupled) THEN
         CALL gas_drag_u(a_m, b_m, ier)
         IF (des_continuum_hybrid) &
            CALL solid_drag_u(a_m, b_m)
      ENDIF

      IF(qmomk .AND. qmomk_coupled) THEN
         CALL qmomk_gas_drag(a_m, b_m, ier, 1, 0, 0)
      ENDIF

      IF (momentum_x_eq(0)) THEN
         CALL calc_resid_u (u_g, v_g, w_g, a_m, b_m, 0, &
            num_resid(resid_u,0), den_resid(resid_u,0), &
            resid(resid_u,0), max_resid(resid_u,0), &
            ijk_resid(resid_u,0))
         CALL under_relax_u (u_g, a_m, b_m, 0, ur_fac(3))
!         call check_ab_m(a_m, b_m, 0, .false., ier)
!         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
!         write(*,*) &
!            resid(resid_u, 0), max_resid(resid_u, 0), &
!            ijk_resid(resid_u, 0)
      ENDIF

      IF(do_solids) THEN
         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
                  CALL calc_resid_u (u_s(1,m), v_s(1,m), w_s(1,m), a_m,&
                     b_m, m, num_resid(resid_u,m), &
                     den_resid(resid_u,m), resid(resid_u,m), &
                     max_resid(resid_u,m), ijk_resid(resid_u,m))
                  CALL under_relax_u (u_s(1,m), a_m, b_m, m, &
                     ur_fac(3))
!                  call check_ab_m(a_m, b_m, m, .false., ier)
!                  write(*,*) &
!                     resid(resid_u, m), max_resid(resid_u, m), &
!                     ijk_resid(resid_u, m)
!                  call write_ab_m(a_m, b_m, ijkmax2, m, ier)
               ENDIF   ! end if (momentum_x_eq(m))
            ENDIF ! end if check for GHD Theory
         ENDDO   ! end do (m=1,mmax)
      ENDIF

      IF (momentum_x_eq(0)) THEN
!         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1,&
!            DO_K,ier)
         CALL adjust_leq (resid(resid_u,0), leq_it(3), leq_method(3), &
            leqi, leqm)
         CALL solve_lin_eq ('U_g', 3, u_g_tmp, a_m, b_m, 0, leqi, leqm, &
            leq_sweep(3), leq_tol(3),  leq_pc(3), ier)
!         call out_array(u_g, 'u_g')
      ENDIF

      IF(do_solids) THEN
         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
!                  call test_lin_eq(ijkmax2, ijmax2, imax2, &
!                     a_m(1, -3, M), 1, DO_K,ier)
                  CALL adjust_leq (resid(resid_u,m), leq_it(3),&
                     leq_method(3), leqi, leqm)
                  CALL solve_lin_eq ('U_s', 3, u_s_tmp(1,m), a_m, &
                     b_m, m, leqi, leqm, leq_sweep(3), leq_tol(3),&
                     leq_pc(3), ier)
!                  call out_array(u_s(1,m), 'u_s')
               ENDIF   ! end if (momentum_x_eq(m))
            ENDIF ! end if check for GHD Theory
         ENDDO   ! end do (m=1,mmax)
      ENDIF
! End U_m_star and residuals
! ----------------------------------------------------------------<<<

      deallocate(vxf_gs)
      deallocate(vxf_ss)
      deallocate(a_m)
      deallocate(b_m)

    END SUBROUTINE u_m_star

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE v_m_star(V_g_tmp, V_s_tmp)
        IMPLICIT NONE
!---------------------------------------------------------------------//
        DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: V_g_tmp
        DOUBLE PRECISION, DIMENSION(:, :), INTENT(OUT) :: V_s_tmp
!---------------------------------------------------------------------//
        ! solids phase index
        INTEGER :: M
        DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: VXF_GS, VXF_SS, B_M
        DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: A_M
!---------------------------------------------------------------------//

        allocate(vxf_gs(dimension_3,dimension_m))
        allocate(vxf_ss(dimension_3,dimension_lm))
        ALLOCATE(a_m(dimension_3, -3:3, 0:dimension_m))
        ALLOCATE(b_m(dimension_3, 0:dimension_m))

! Calculate V_m_star and residuals
! ---------------------------------------------------------------->>>
      DO m = 0, mmax
         CALL init_ab_m (a_m, b_m, ijkmax2, m)
         IF (m >= 1) vxf_gs(:,m) = zero
      ENDDO

! convection-diffusion terms
      CALL conv_dif_v_g (a_m, b_m, ier)
!      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
      IF(do_solids) CALL conv_dif_v_s (a_m, b_m, ier)

! source terms
      CALL source_v_g (a_m, b_m)
      IF(point_source) CALL point_source_v_g (a_m, b_m)
      IF(call_usr_source(4)) CALL calc_usr_source(gas_v_mom, a_m, b_m)
!      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
      IF(do_solids) THEN
         CALL source_v_s (a_m, b_m)
         IF(point_source) CALL point_source_v_s (a_m, b_m)
         IF(call_usr_source(4)) CALL calc_usr_source(solids_v_mom, a_m, b_m)
      ENDIF

! evaluate local vxf_gs and vxf_ss
      CALL vf_gs_y (vxf_gs)
      IF(do_solids .AND. (kt_type_enum /= ghd_2007)) THEN
         IF (mmax > 0) CALL vf_ss_y (vxf_ss)
      ENDIF

! calculate coefficients for the pressure correction equation
      IF(kt_type_enum == ghd_2007) THEN
         CALL calc_d_ghd_n (a_m, vxf_gs, d_n)
      ELSE
         CALL calc_d_n (a_m, vxf_gs, vxf_ss, d_n, ier)
      ENDIF

      IF(do_solids) THEN
! calculate coefficients for a solids pressure correction equation
         IF (mmax > 0) CALL calc_e_n (a_m, mcp, e_n)

! calculate modifications to the A matrix center coefficient and B
! source vector for partial elimination
         IF(kt_type_enum == ghd_2007) THEN
           IF (mmax > 0) CALL partial_elim_ghd_v (v_g,v_s,vxf_gs,a_m,b_m)
         ELSE
           IF (mmax > 0) CALL partial_elim_v (v_g,v_s,vxf_gs,a_m,b_m)
         ENDIF
      ENDIF

! handle special case where center coefficient is zero
      CALL adjust_a_v_g (a_m, b_m)
!      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
      IF(do_solids) CALL adjust_a_v_s (a_m, b_m)
!      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)

! modification to matrix equation for DEM drag terms
      IF(des_continuum_coupled) THEN
         CALL gas_drag_v(a_m, b_m, ier)
         IF (des_continuum_hybrid) &
            CALL solid_drag_v(a_m, b_m)
      ENDIF

      IF(qmomk .AND. qmomk_coupled) THEN
         CALL qmomk_gas_drag(a_m, b_m, ier, 0, 1, 0)
      ENDIF


      IF (momentum_y_eq(0)) THEN
         CALL calc_resid_v (u_g, v_g, w_g, a_m, b_m, 0, &
            num_resid(resid_v,0), den_resid(resid_v,0), &
            resid(resid_v,0), max_resid(resid_v,0), &
            ijk_resid(resid_v,0))
         CALL under_relax_v (v_g, a_m, b_m, 0, ur_fac(4))
!         call check_ab_m(a_m, b_m, 0, .false., ier)
!         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
!         write(*,*) &
!            resid(resid_v, 0), max_resid(resid_v, 0), &
!            ijk_resid(resid_v, 0)
      ENDIF

      IF(do_solids) THEN
         DO m = 1, mmax
            IF(kt_type_enum /= ghd_2007 .OR. &
              (kt_type_enum == ghd_2007 .AND. m==mmax)) THEN
               IF (momentum_y_eq(m)) THEN
                  CALL calc_resid_v (u_s(1,m), v_s(1,m), w_s(1,m), a_m,&
                     b_m, m, num_resid(resid_v,m), &
                     den_resid(resid_v,m),resid(resid_v,m), &
                     max_resid(resid_v,m), ijk_resid(resid_v,m))
                  CALL under_relax_v (v_s(1,m),a_m,b_m,m,ur_fac(4))
!                  call check_ab_m(a_m, b_m, m, .false., ier)
!                  write(*,*) &
!                     resid(resid_v, m), max_resid(resid_v, m),
!                     ijk_resid(resid_v, m)
!                  call write_ab_m(a_m, b_m, ijkmax2, m, ier)
               ENDIF   ! end if (momentum_y_eq(m))
            ENDIF ! end if check for GHD Theory
         ENDDO   ! end do (m=1,mmax)
      ENDIF

      IF (momentum_y_eq(0)) THEN
!         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), &
!            1, DO_K, ier)
         CALL adjust_leq (resid(resid_v,0), leq_it(4), leq_method(4), &
            leqi, leqm)
         CALL solve_lin_eq ('V_g', 4, v_g_tmp, a_m, b_m, 0, leqi, leqm, &
            leq_sweep(4), leq_tol(4),  leq_pc(4), ier)
!         call out_array(v_g, 'v_g')
      ENDIF

      IF(do_solids) THEN
         DO m = 1, mmax
            IF(kt_type_enum /= ghd_2007 .OR. &
              (kt_type_enum == ghd_2007 .AND. m==mmax)) THEN
               IF (momentum_y_eq(m)) THEN
!                  call test_lin_eq(ijkmax2, ijmax2, imax2, &
!                     a_m(1, -3, M), 1, DO_K, ier)
                  CALL adjust_leq (resid(resid_v,m), leq_it(4), &
                     leq_method(4), leqi, leqm)
                  CALL solve_lin_eq ('V_s', 4, v_s_tmp(1,m), a_m, &
                     b_m, m, leqi, leqm, leq_sweep(4), leq_tol(4), &
                     leq_pc(4), ier)
!                  call out_array(v_s(1,m), 'v_s')
               ENDIF   ! end if (momentum_y_eq(m))
            ENDIF ! end if check for GHD Theory
         ENDDO   ! end do (m=1,mmax)
      ENDIF
! End V_m_star and residuals
! ----------------------------------------------------------------<<<

      deallocate(vxf_gs)
      deallocate(vxf_ss)
      deallocate(a_m)
      deallocate(b_m)

    END SUBROUTINE v_m_star

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE w_m_star(W_g_tmp, W_s_tmp)
        IMPLICIT NONE
!---------------------------------------------------------------------//
        DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: W_g_tmp
        DOUBLE PRECISION, DIMENSION(:, :), INTENT(OUT) :: W_s_tmp
!---------------------------------------------------------------------//
        ! solids phase index
        INTEGER :: M
        DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: VXF_GS, VXF_SS, B_M
        DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: A_M
!---------------------------------------------------------------------//

        ALLOCATE(vxf_gs(dimension_3,dimension_m))
        ALLOCATE(vxf_ss(dimension_3,dimension_lm))
        ALLOCATE(a_m(dimension_3, -3:3, 0:dimension_m))
        ALLOCATE(b_m(dimension_3, 0:dimension_m))

! Calculate W_m_star and residuals
! ---------------------------------------------------------------->>>
      IF (do_k)THEN
         DO m = 0, mmax
            CALL init_ab_m (a_m, b_m, ijkmax2, m)
            IF (m >= 1) vxf_gs(:,m) = zero
         ENDDO

! convection diffusion
         CALL conv_dif_w_g (a_m, b_m)
         IF(do_solids) CALL conv_dif_w_s (a_m, b_m)

! source terms
         CALL source_w_g (a_m, b_m)
         IF(point_source) CALL point_source_w_g (a_m, b_m)
         IF(call_usr_source(5)) CALL calc_usr_source(gas_w_mom, a_m, b_m)
!         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
         IF(do_solids) THEN
            CALL source_w_s (a_m, b_m)
            IF(point_source) CALL point_source_w_s (a_m, b_m)
            IF(call_usr_source(5)) CALL calc_usr_source(solids_w_mom, a_m, b_m)
         ENDIF
!        call write_ab_m(a_m, b_m, ijkmax2, 0, ier)

! evaluate local variable vxf_gs and vxf_ss
         CALL vf_gs_z (vxf_gs)
         IF(do_solids .AND. (kt_type_enum /= ghd_2007)) THEN
            IF (mmax > 0) CALL vf_ss_z (vxf_ss)
         ENDIF

! calculate coefficients for the pressure correction equation
         IF(kt_type_enum == ghd_2007) THEN
            CALL calc_d_ghd_t (a_m, vxf_gs, d_t)
         ELSE
            CALL calc_d_t (a_m, vxf_gs, vxf_ss, d_t, ier)
         ENDIF

         IF(do_solids) THEN
! calculate coefficients for a solids pressure correction equation
            IF (mmax > 0) CALL calc_e_t (a_m, mcp, e_t)

! calculate modifications to the A matrix center coefficient and B
! source vector for partial elimination
            IF(kt_type_enum == ghd_2007) THEN
               IF (mmax > 0) CALL partial_elim_ghd_w (w_g, w_s, vxf_gs, a_m, b_m)
            ELSE
               IF (mmax > 0) CALL partial_elim_w (w_g, w_s, vxf_gs, a_m, b_m)
            ENDIF
         ENDIF

! handle special case where center coefficient is zero
         CALL adjust_a_w_g (a_m, b_m)
         IF(do_solids) THEN
            CALL adjust_a_w_s (a_m, b_m)
         ENDIF

! modifications to matrix equation for DEM
         IF(des_continuum_coupled) THEN
            CALL gas_drag_w(a_m, b_m, ier)
            IF (discrete_element .AND. des_continuum_hybrid) &
               CALL solid_drag_w(a_m, b_m)
         ENDIF

         IF(qmomk .AND. qmomk_coupled) THEN
            CALL qmomk_gas_drag(a_m, b_m, ier, 0, 0, 1)
         ENDIF

         IF (momentum_z_eq(0)) THEN
            CALL calc_resid_w (u_g, v_g, w_g, a_m, b_m, 0, &
               num_resid(resid_w,0), den_resid(resid_w,0), &
               resid(resid_w,0), max_resid(resid_w,0), &
               ijk_resid(resid_w,0))
            CALL under_relax_w (w_g, a_m, b_m, 0, ur_fac(5))
!            call check_ab_m(a_m, b_m, 0, .false., ier)
!            write(*,*) &
!               resid(resid_w, 0), max_resid(resid_w, 0), &
!               ijk_resid(resid_w, 0)
         ENDIF

         IF(do_solids) THEN
            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
                     CALL calc_resid_w (u_s(1,m), v_s(1,m), w_s(1,m),&
                        a_m, b_m, m, num_resid(resid_w,m), &
                        den_resid(resid_w,m), resid(resid_w,m), &
                        max_resid(resid_w,m), ijk_resid(resid_w,m))
                     CALL under_relax_w (w_s(1,m), a_m, b_m, m, &
                        ur_fac(5))
!                     call check_ab_m(a_m, b_m, m, .false., ier)
!                     write(*,*) &
!                        resid(resid_w, m), max_resid(resid_w, m), &
!                        ijk_resid(resid_w, m)
!                     call write_ab_m(a_m, b_m, ijkmax2, m, ier)
                  ENDIF   ! end if (momentum_z_eq(m))
               ENDIF ! end if check for GHD Theory
            ENDDO   ! end do (m=1,mmax)
         ENDIF

         IF (momentum_z_eq(0)) THEN
!            call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), &
!               1, DO_K, ier)
            CALL adjust_leq (resid(resid_w,0), leq_it(5), &
               leq_method(5), leqi, leqm)
            CALL solve_lin_eq ('W_g', 5, w_g_tmp, a_m, b_m, 0, leqi,&
               leqm, leq_sweep(5), leq_tol(5), leq_pc(5), ier)
!            call out_array(w_g, 'w_g')
         ENDIF

         IF(do_solids) THEN
            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
!                     call test_lin_eq(ijkmax2, ijmax2, imax2, &
!                        a_m(1, -3, M), 1, DO_K, ier)
                     CALL adjust_leq (resid(resid_w,m), leq_it(5), &
                        leq_method(5), leqi, leqm)
                     CALL solve_lin_eq ('W_s', 5, w_s_tmp(1,m), &
                        a_m, b_m, m, leqi, leqm, leq_sweep(5), &
                        leq_tol(5), leq_pc(5), ier)
!                     call out_array(w_s(1,m), 'w_s')
                  ENDIF   ! end if (momentum_z_eq(m))
               ENDIF ! end if check for GHD Theory
            ENDDO   ! end do (m=1,mmax)
         ENDIF
      ENDIF   ! end if (do_k)
! End W_m_star and residuals
! ----------------------------------------------------------------<<<

      deallocate(vxf_gs)
      deallocate(vxf_ss)
      deallocate(a_m)
      deallocate(b_m)

    END SUBROUTINE w_m_star

  END SUBROUTINE solve_vel_star