File: N:\mfix\model\solve_vel_star.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_VEL_STAR                                          C
4     !  Purpose: Solve starred velocity components                          C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 25-APR-96  C
7     !                                                                      C
8     !  Revision Number: 2                                                  C
9     !  Purpose: allow multiparticle in D_E, D_N and D_T claculations       C
10     !           and account for the averaged Solid-Solid drag              C
11     !  Author: S. Dartevelle, LANL                        Date: 28-FEb-04  C
12     !                                                                      C
13     !                                                                      C
14     !  Literature/Document References:                                     C
15     !                                                                      C
16     !                                                                      C
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
18           SUBROUTINE SOLVE_VEL_STAR(IER)
19     
20     ! Modules
21     !---------------------------------------------------------------------//
22           use compar, only: ijkstart3, ijkend3
23           use discretelement, only: des_continuum_hybrid, discrete_element, des_continuum_coupled
24           use fldvar, only: u_g, v_g, w_g, u_s, v_s, w_s
25           use geometry, only: do_k, ijkmax2
26           use leqsol, only: leq_sweep, leq_tol, leq_pc, leq_it, leq_method
27           use param, only: dimension_3, dimension_m
28           use param1, only: zero, dimension_lm
29           use pgcor, only: d_e, d_t, d_n
30           use ps, only: point_source
31           use pscor, only: e_e, e_t, e_n, mcp
32           use physprop,only: mmax
33           use residual, only: resid_u, resid_v, resid_w, num_resid, den_resid, resid, max_resid, ijk_resid
34           use run, only: kt_type_enum, GHD_2007
35           use run, only: momentum_x_eq, momentum_y_eq, momentum_z_eq
36           use ur_facs, only: ur_fac
37           ! Use ambm
38           ! Use tmp_array1,  VxF_gs => Arraym1
39           ! Use tmp_array,  VxF_ss => ArrayLM
40           USE qmom_kinetic_equation
41           use usr_src, only: call_usr_source, calc_usr_source
42           use usr_src, only: gas_u_mom, gas_v_mom, gas_w_mom
43           use usr_src, only: solids_u_mom, solids_v_mom, solids_w_mom
44           IMPLICIT NONE
45     
46     ! Dummy arguments
47     !---------------------------------------------------------------------//
48     ! Error index
49           INTEGER, INTENT(INOUT) :: IER
50     
51     ! Local variables
52     !---------------------------------------------------------------------//
53     ! fluid cell index
54           INTEGER :: IJK
55     ! temporary velocity arrays
56           DOUBLE PRECISION, DIMENSION(:), allocatable :: U_gtmp,  V_gtmp, W_gtmp
57           DOUBLE PRECISION, DIMENSION(:,:), allocatable :: U_stmp, V_stmp, W_stmp
58     ! linear equation solver method and iterations
59           INTEGER :: LEQM, LEQI
60     
61           LOGICAL :: DO_SOLIDS
62     
63     ! temporary use of global arrays:
64     ! arraym1 (locally vxf_gs)
65     ! the volume x average gas-solids drag at momentum cell centers
66     !      DOUBLE PRECISION :: VXF_GS(DIMENSION_3, DIMENSION_M)
67     ! arraylm (locally vxf_ss)
68     ! the volume x average solids-solids drag at momentum cell centers
69     !      DOUBLE PRECISION :: VXF_SS(DIMENSION_3, DIMENSION_LM)
70     ! Septadiagonal matrix A_m, vector b_m
71     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
72     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
73     !---------------------------------------------------------------------//
74     
75           allocate(U_gtmp(DIMENSION_3))
76           allocate(V_gtmp(DIMENSION_3))
77           allocate(W_gtmp(DIMENSION_3))
78           allocate(U_stmp(DIMENSION_3,DIMENSION_M))
79           allocate(V_stmp(DIMENSION_3,DIMENSION_M))
80           allocate(W_stmp(DIMENSION_3,DIMENSION_M))
81     
82           ! call lock_ambm        ! locks arrys a_m and b_m
83           ! call lock_tmp_array1  ! locks arraym1 (locally vxf_gs)
84           ! call lock_tmp_array2  ! locks arraylm (locally vxf_ss)
85     
86           CALL init(U_gtmp, V_gtmp, W_gtmp, U_stmp, V_stmp, W_stmp)
87     
88           DO_SOLIDS = .NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
89              DES_CONTINUUM_HYBRID
90     
91     !!$omp parallel sections default(shared)
92           CALL U_m_star(U_gtmp,U_stmp)
93     !!$omp section
94           CALL V_m_star(V_gtmp,V_stmp)
95     !!$omp section
96           CALL W_m_star(W_gtmp,W_stmp)
97     !!$omp end parallel sections
98     
99           CALL save(U_gtmp, V_gtmp, W_gtmp, U_stmp, V_stmp, W_stmp)
100     
101     ! modification for GHD theory to compute species velocity: Ui = Joi/(mi ni) + U.
102           IF(KT_TYPE_ENUM == GHD_2007) THEN
103              CALL calc_external_forces()
104              CALL GHDMassFlux() ! to compute solid species mass flux
105              CALL UpdateSpeciesVelocities() ! located at end of ghdMassFlux.f file
106           ENDIF
107     
108           ! call unlock_ambm
109           ! call unlock_tmp_array1
110           ! call unlock_tmp_array2
111     
112           deallocate(U_gtmp)
113           deallocate(V_gtmp)
114           deallocate(W_gtmp)
115           deallocate(U_stmp)
116           deallocate(V_stmp)
117           deallocate(W_stmp)
118     
119           RETURN
120     
121         CONTAINS
122     
123     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
124     !                                                                      C
125     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
126           SUBROUTINE init(U_g_tmp, V_g_tmp, W_g_tmp, U_s_tmp, V_s_tmp, W_s_tmp)
127             IMPLICIT NONE
128     !---------------------------------------------------------------------//
129             DOUBLE PRECISION, DIMENSION(:), intent(out) :: U_g_tmp,  V_g_tmp, W_g_tmp
130             DOUBLE PRECISION, DIMENSION(:,:), intent(out) :: U_s_tmp, V_s_tmp, W_s_tmp
131     !---------------------------------------------------------------------//
132             ! solids phase index
133             INTEGER :: M
134     !---------------------------------------------------------------------//
135     
136             ! Store the velocities so that the order of solving the momentum
137             ! equations does not matter
138             DO IJK = ijkstart3, ijkend3
139                U_g_tmp(IJK) = U_g(IJK)
140                V_g_tmp(IJK) = V_g(IJK)
141                W_g_tmp(IJK) = W_g(IJK)
142             ENDDO
143             DO M = 1, MMAX
144                IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
145                     (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
146                   DO IJK = ijkstart3, ijkend3
147                      U_s_tmp(IJK, M) = U_s(IJK, M)
148                      V_s_tmp(IJK, M) = V_s(IJK, M)
149                      W_s_tmp(IJK, M) = W_s(IJK, M)
150                   ENDDO
151                ENDIF
152             ENDDO
153     
154           END SUBROUTINE init
155     
156     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
157     !                                                                      C
158     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
159           SUBROUTINE save(U_g_tmp, V_g_tmp, W_g_tmp, U_s_tmp, V_s_tmp, W_s_tmp)
160             IMPLICIT NONE
161     !---------------------------------------------------------------------//
162             DOUBLE PRECISION, DIMENSION(:), intent(in) :: U_g_tmp,  V_g_tmp, W_g_tmp
163             DOUBLE PRECISION, DIMENSION(:,:), intent(in) :: U_s_tmp, V_s_tmp, W_s_tmp
164     !---------------------------------------------------------------------//
165             ! solids phase index
166             INTEGER :: M
167     !---------------------------------------------------------------------//
168     
169             ! Now update all velocity components
170             DO IJK = ijkstart3, ijkend3
171                U_g(IJK) = U_g_tmp(IJK)
172                V_g(IJK) = V_g_tmp(IJK)
173                W_g(IJK) = W_g_tmp(IJK)
174             ENDDO
175             DO M = 1, MMAX
176                IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
177                     (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
178                   DO IJK = ijkstart3, ijkend3
179                      U_s(IJK, M) = U_s_tmp(IJK, M)
180                      V_s(IJK, M) = V_s_tmp(IJK, M)
181                      W_s(IJK, M) = W_s_tmp(IJK, M)
182                   ENDDO
183                ENDIF
184             ENDDO
185     
186           END SUBROUTINE save
187     
188     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
189     !                                                                      C
190     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
191           SUBROUTINE U_m_star(U_g_tmp, U_s_tmp)
192             IMPLICIT NONE
193     !---------------------------------------------------------------------//
194             DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: U_g_tmp
195             DOUBLE PRECISION, DIMENSION(:, :), INTENT(OUT) :: U_s_tmp
196     !---------------------------------------------------------------------//
197             ! solids phase index
198             INTEGER :: M
199             DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: VXF_GS, VXF_SS, B_M
200             DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: A_M
201     !---------------------------------------------------------------------//
202     
203             allocate(vxf_gs(DIMENSION_3,DIMENSION_M))
204             allocate(vxf_ss(DIMENSION_3,DIMENSION_LM))
205             ALLOCATE(A_M(DIMENSION_3, -3:3, 0:DIMENSION_M))
206             ALLOCATE(B_M(DIMENSION_3, 0:DIMENSION_M))
207     
208     ! Calculate U_m_star and residuals
209     ! ---------------------------------------------------------------->>>
210           DO M = 0, MMAX
211              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
212              IF (M >= 1) VXF_GS(:,M) = ZERO
213           ENDDO
214     
215     ! calculate the convection-diffusion terms for the gas and solids phase
216     ! u-momentum equations
217           CALL CONV_DIF_U_G (A_M, B_M)
218           IF(DO_SOLIDS) CALL CONV_DIF_U_S (A_M, B_M)
219     
220     ! calculate the source terms for the gas and solids phase u-momentum
221     ! equations
222           CALL SOURCE_U_G (A_M, B_M)
223           IF(POINT_SOURCE) CALL POINT_SOURCE_U_G (A_M, B_M)
224           IF(CALL_USR_SOURCE(3)) CALL CALC_USR_SOURCE(GAS_U_MOM, A_M, B_M)
225           IF(DO_SOLIDS) THEN
226              CALL SOURCE_U_S (A_M, B_M)
227              IF(POINT_SOURCE) CALL POINT_SOURCE_U_S (A_M, B_M)
228              IF(CALL_USR_SOURCE(3)) CALL CALC_USR_SOURCE(SOLIDS_U_MOM, A_M, B_M)
229           ENDIF
230     
231     ! evaluate local variable vxf_gs and vxf_ss.  both terms are sent to the
232     ! subroutine calc_d (pressure correction equation coefficients).  the
233     ! former is also used in the subroutine partial_elim_u while the latter
234     ! is effectively re-evaluated within said subroutine
235           CALL VF_GS_X (VXF_GS)
236           IF(DO_SOLIDS .AND. (KT_TYPE_ENUM /= GHD_2007)) THEN
237              IF (MMAX > 0) CALL VF_SS_X (VXF_SS)
238           ENDIF
239     
240     ! calculate coefficients for the pressure correction equation
241           IF(KT_TYPE_ENUM == GHD_2007) THEN
242              CALL CALC_D_GHD_E (A_M, VXF_GS, D_E)
243           ELSE
244              CALL CALC_D_E (A_M, VXF_GS, VXF_SS, D_E, IER)
245           ENDIF
246     
247           IF(DO_SOLIDS) THEN
248     ! calculate coefficients for a solids volume correction equation
249              IF (MMAX > 0) CALL CALC_E_E (A_M, MCP, E_E)
250     
251     ! calculate modifications to the A matrix center coefficient and B
252     ! source vector for partial elimination
253              IF(KT_TYPE_ENUM == GHD_2007) THEN
254                 IF (MMAX > 0) CALL PARTIAL_ELIM_GHD_U (U_G,U_S,VXF_GS,A_M,B_M)
255              ELSE
256                 IF (MMAX > 0) CALL PARTIAL_ELIM_U (U_G,U_S,VXF_GS,A_M,B_M)
257              ENDIF
258           ENDIF
259     
260     ! handle special case where center coefficient is zero
261           CALL ADJUST_A_U_G (A_M, B_M)
262           IF(DO_SOLIDS) CALL ADJUST_A_U_S (A_M, B_M)
263     
264     ! calculate modifications to the A matrix center coefficient and B
265     ! source vector for treating DEM drag terms
266           IF(DES_CONTINUUM_COUPLED) THEN
267              CALL GAS_DRAG_U(A_M, B_M, IER)
268              IF (DES_CONTINUUM_HYBRID) &
269                 CALL SOLID_DRAG_U(A_M, B_M)
270           ENDIF
271     
272           IF(QMOMK .AND. QMOMK_COUPLED) THEN
273              CALL QMOMK_GAS_DRAG(A_M, B_M, IER, 1, 0, 0)
274           ENDIF
275     
276           IF (MOMENTUM_X_EQ(0)) THEN
277              CALL CALC_RESID_U (U_G, V_G, W_G, A_M, B_M, 0, &
278                 NUM_RESID(RESID_U,0), DEN_RESID(RESID_U,0), &
279                 RESID(RESID_U,0), MAX_RESID(RESID_U,0), &
280                 IJK_RESID(RESID_U,0))
281              CALL UNDER_RELAX_U (U_G, A_M, B_M, 0, UR_FAC(3))
282     !         call check_ab_m(a_m, b_m, 0, .false., ier)
283     !         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
284     !         write(*,*) &
285     !            resid(resid_u, 0), max_resid(resid_u, 0), &
286     !            ijk_resid(resid_u, 0)
287           ENDIF
288     
289           IF(DO_SOLIDS) THEN
290              DO M = 1, MMAX
291                 IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
292                   (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
293                    IF (MOMENTUM_X_EQ(M)) THEN
294                       CALL CALC_RESID_U (U_S(1,M), V_S(1,M), W_S(1,M), A_M,&
295                          B_M, M, NUM_RESID(RESID_U,M), &
296                          DEN_RESID(RESID_U,M), RESID(RESID_U,M), &
297                          MAX_RESID(RESID_U,M), IJK_RESID(RESID_U,M))
298                       CALL UNDER_RELAX_U (U_S(1,M), A_M, B_M, M, &
299                          UR_FAC(3))
300     !                  call check_ab_m(a_m, b_m, m, .false., ier)
301     !                  write(*,*) &
302     !                     resid(resid_u, m), max_resid(resid_u, m), &
303     !                     ijk_resid(resid_u, m)
304     !                  call write_ab_m(a_m, b_m, ijkmax2, m, ier)
305                    ENDIF   ! end if (momentum_x_eq(m))
306                 ENDIF ! end if check for GHD Theory
307              ENDDO   ! end do (m=1,mmax)
308           ENDIF
309     
310           IF (MOMENTUM_X_EQ(0)) THEN
311     !         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1,&
312     !            DO_K,ier)
313              CALL ADJUST_LEQ (RESID(RESID_U,0), LEQ_IT(3), LEQ_METHOD(3), &
314                 LEQI, LEQM)
315              CALL SOLVE_LIN_EQ ('U_g', 3, U_G_tmp, A_M, B_M, 0, LEQI, LEQM, &
316                 LEQ_SWEEP(3), LEQ_TOL(3),  LEQ_PC(3), IER)
317     !         call out_array(u_g, 'u_g')
318           ENDIF
319     
320           IF(DO_SOLIDS) THEN
321              DO M = 1, MMAX
322                 IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
323                   (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
324                    IF (MOMENTUM_X_EQ(M)) THEN
325     !                  call test_lin_eq(ijkmax2, ijmax2, imax2, &
326     !                     a_m(1, -3, M), 1, DO_K,ier)
327                       CALL ADJUST_LEQ (RESID(RESID_U,M), LEQ_IT(3),&
328                          LEQ_METHOD(3), LEQI, LEQM)
329                       CALL SOLVE_LIN_EQ ('U_s', 3, U_S_tmp(1,M), A_M, &
330                          B_M, M, LEQI, LEQM, LEQ_SWEEP(3), LEQ_TOL(3),&
331                          LEQ_PC(3), IER)
332     !                  call out_array(u_s(1,m), 'u_s')
333                    ENDIF   ! end if (momentum_x_eq(m))
334                 ENDIF ! end if check for GHD Theory
335              ENDDO   ! end do (m=1,mmax)
336           ENDIF
337     ! End U_m_star and residuals
338     ! ----------------------------------------------------------------<<<
339     
340           deallocate(vxf_gs)
341           deallocate(vxf_ss)
342           deallocate(a_m)
343           deallocate(b_m)
344     
345         END SUBROUTINE U_m_star
346     
347     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
348     !                                                                      C
349     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
350           SUBROUTINE V_m_star(V_g_tmp, V_s_tmp)
351             IMPLICIT NONE
352     !---------------------------------------------------------------------//
353             DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: V_g_tmp
354             DOUBLE PRECISION, DIMENSION(:, :), INTENT(OUT) :: V_s_tmp
355     !---------------------------------------------------------------------//
356             ! solids phase index
357             INTEGER :: M
358             DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: VXF_GS, VXF_SS, B_M
359             DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: A_M
360     !---------------------------------------------------------------------//
361     
362             allocate(vxf_gs(DIMENSION_3,DIMENSION_M))
363             allocate(vxf_ss(DIMENSION_3,DIMENSION_LM))
364             ALLOCATE(A_M(DIMENSION_3, -3:3, 0:DIMENSION_M))
365             ALLOCATE(B_M(DIMENSION_3, 0:DIMENSION_M))
366     
367     ! Calculate V_m_star and residuals
368     ! ---------------------------------------------------------------->>>
369           DO M = 0, MMAX
370              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
371              IF (M >= 1) VXF_GS(:,M) = ZERO
372           ENDDO
373     
374     ! convection-diffusion terms
375           CALL CONV_DIF_V_G (A_M, B_M, IER)
376     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
377           IF(DO_SOLIDS) CALL CONV_DIF_V_S (A_M, B_M, IER)
378     
379     ! source terms
380           CALL SOURCE_V_G (A_M, B_M)
381           IF(POINT_SOURCE) CALL POINT_SOURCE_V_G (A_M, B_M)
382           IF(CALL_USR_SOURCE(4)) CALL CALC_USR_SOURCE(GAS_V_MOM, A_M, B_M)
383     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
384           IF(DO_SOLIDS) THEN
385              CALL SOURCE_V_S (A_M, B_M)
386              IF(POINT_SOURCE) CALL POINT_SOURCE_V_S (A_M, B_M)
387              IF(CALL_USR_SOURCE(4)) CALL CALC_USR_SOURCE(SOLIDS_V_MOM, A_M, B_M)
388           ENDIF
389     
390     ! evaluate local vxf_gs and vxf_ss
391           CALL VF_GS_Y (VXF_GS)
392           IF(DO_SOLIDS .AND. (KT_TYPE_ENUM /= GHD_2007)) THEN
393              IF (MMAX > 0) CALL VF_SS_Y (VXF_SS)
394           ENDIF
395     
396     ! calculate coefficients for the pressure correction equation
397           IF(KT_TYPE_ENUM == GHD_2007) THEN
398              CALL CALC_D_GHD_N (A_M, VXF_GS, D_N)
399           ELSE
400              CALL CALC_D_N (A_M, VXF_GS, VXF_SS, D_N, IER)
401           ENDIF
402     
403           IF(DO_SOLIDS) THEN
404     ! calculate coefficients for a solids pressure correction equation
405              IF (MMAX > 0) CALL CALC_E_N (A_M, MCP, E_N)
406     
407     ! calculate modifications to the A matrix center coefficient and B
408     ! source vector for partial elimination
409              IF(KT_TYPE_ENUM == GHD_2007) THEN
410                IF (MMAX > 0) CALL PARTIAL_ELIM_GHD_V (V_G,V_S,VXF_GS,A_M,B_M)
411              ELSE
412                IF (MMAX > 0) CALL PARTIAL_ELIM_V (V_G,V_S,VXF_GS,A_M,B_M)
413              ENDIF
414           ENDIF
415     
416     ! handle special case where center coefficient is zero
417           CALL ADJUST_A_V_G (A_M, B_M)
418     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
419           IF(DO_SOLIDS) CALL ADJUST_A_V_S (A_M, B_M)
420     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
421     
422     ! modification to matrix equation for DEM drag terms
423           IF(DES_CONTINUUM_COUPLED) THEN
424              CALL GAS_DRAG_V(A_M, B_M, IER)
425              IF (DES_CONTINUUM_HYBRID) &
426                 CALL SOLID_DRAG_V(A_M, B_M)
427           ENDIF
428     
429           IF(QMOMK .AND. QMOMK_COUPLED) THEN
430              CALL QMOMK_GAS_DRAG(A_M, B_M, IER, 0, 1, 0)
431           ENDIF
432     
433     
434           IF (MOMENTUM_Y_EQ(0)) THEN
435              CALL CALC_RESID_V (U_G, V_G, W_G, A_M, B_M, 0, &
436                 NUM_RESID(RESID_V,0), DEN_RESID(RESID_V,0), &
437                 RESID(RESID_V,0), MAX_RESID(RESID_V,0), &
438                 IJK_RESID(RESID_V,0))
439              CALL UNDER_RELAX_V (V_G, A_M, B_M, 0, UR_FAC(4))
440     !         call check_ab_m(a_m, b_m, 0, .false., ier)
441     !         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
442     !         write(*,*) &
443     !            resid(resid_v, 0), max_resid(resid_v, 0), &
444     !            ijk_resid(resid_v, 0)
445           ENDIF
446     
447           IF(DO_SOLIDS) THEN
448              DO M = 1, MMAX
449                 IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
450                   (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
451                    IF (MOMENTUM_Y_EQ(M)) THEN
452                       CALL CALC_RESID_V (U_S(1,M), V_S(1,M), W_S(1,M), A_M,&
453                          B_M, M, NUM_RESID(RESID_V,M), &
454                          DEN_RESID(RESID_V,M),RESID(RESID_V,M), &
455                          MAX_RESID(RESID_V,M), IJK_RESID(RESID_V,M))
456                       CALL UNDER_RELAX_V (V_S(1,M),A_M,B_M,M,UR_FAC(4))
457     !                  call check_ab_m(a_m, b_m, m, .false., ier)
458     !                  write(*,*) &
459     !                     resid(resid_v, m), max_resid(resid_v, m),
460     !                     ijk_resid(resid_v, m)
461     !                  call write_ab_m(a_m, b_m, ijkmax2, m, ier)
462                    ENDIF   ! end if (momentum_y_eq(m))
463                 ENDIF ! end if check for GHD Theory
464              ENDDO   ! end do (m=1,mmax)
465           ENDIF
466     
467           IF (MOMENTUM_Y_EQ(0)) THEN
468     !         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), &
469     !            1, DO_K, ier)
470              CALL ADJUST_LEQ (RESID(RESID_V,0), LEQ_IT(4), LEQ_METHOD(4), &
471                 LEQI, LEQM)
472              CALL SOLVE_LIN_EQ ('V_g', 4, V_G_tmp, A_M, B_M, 0, LEQI, LEQM, &
473                 LEQ_SWEEP(4), LEQ_TOL(4),  LEQ_PC(4), IER)
474     !         call out_array(v_g, 'v_g')
475           ENDIF
476     
477           IF(DO_SOLIDS) THEN
478              DO M = 1, MMAX
479                 IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
480                   (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
481                    IF (MOMENTUM_Y_EQ(M)) THEN
482     !                  call test_lin_eq(ijkmax2, ijmax2, imax2, &
483     !                     a_m(1, -3, M), 1, DO_K, ier)
484                       CALL ADJUST_LEQ (RESID(RESID_V,M), LEQ_IT(4), &
485                          LEQ_METHOD(4), LEQI, LEQM)
486                       CALL SOLVE_LIN_EQ ('V_s', 4, V_S_tmp(1,M), A_M, &
487                          B_M, M, LEQI, LEQM, LEQ_SWEEP(4), LEQ_TOL(4), &
488                          LEQ_PC(4), IER)
489     !                  call out_array(v_s(1,m), 'v_s')
490                    ENDIF   ! end if (momentum_y_eq(m))
491                 ENDIF ! end if check for GHD Theory
492              ENDDO   ! end do (m=1,mmax)
493           ENDIF
494     ! End V_m_star and residuals
495     ! ----------------------------------------------------------------<<<
496     
497           deallocate(vxf_gs)
498           deallocate(vxf_ss)
499           deallocate(a_m)
500           deallocate(b_m)
501     
502         END SUBROUTINE V_m_star
503     
504     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
505     !                                                                      C
506     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
507           SUBROUTINE W_m_star(W_g_tmp, W_s_tmp)
508             IMPLICIT NONE
509     !---------------------------------------------------------------------//
510             DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: W_g_tmp
511             DOUBLE PRECISION, DIMENSION(:, :), INTENT(OUT) :: W_s_tmp
512     !---------------------------------------------------------------------//
513             ! solids phase index
514             INTEGER :: M
515             DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: VXF_GS, VXF_SS, B_M
516             DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: A_M
517     !---------------------------------------------------------------------//
518     
519             ALLOCATE(VXF_GS(DIMENSION_3,DIMENSION_M))
520             ALLOCATE(VXF_SS(DIMENSION_3,DIMENSION_LM))
521             ALLOCATE(A_M(DIMENSION_3, -3:3, 0:DIMENSION_M))
522             ALLOCATE(B_M(DIMENSION_3, 0:DIMENSION_M))
523     
524     ! Calculate W_m_star and residuals
525     ! ---------------------------------------------------------------->>>
526           IF (DO_K)THEN
527              DO M = 0, MMAX
528                 CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
529                 IF (M >= 1) VXF_GS(:,M) = ZERO
530              ENDDO
531     
532     ! convection diffusion
533              CALL CONV_DIF_W_G (A_M, B_M)
534              IF(DO_SOLIDS) CALL CONV_DIF_W_S (A_M, B_M)
535     
536     ! source terms
537              CALL SOURCE_W_G (A_M, B_M)
538              IF(POINT_SOURCE) CALL POINT_SOURCE_W_G (A_M, B_M)
539              IF(CALL_USR_SOURCE(5)) CALL CALC_USR_SOURCE(GAS_W_MOM, A_M, B_M)
540     !         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
541              IF(DO_SOLIDS) THEN
542                 CALL SOURCE_W_S (A_M, B_M)
543                 IF(POINT_SOURCE) CALL POINT_SOURCE_W_S (A_M, B_M)
544                 IF(CALL_USR_SOURCE(5)) CALL CALC_USR_SOURCE(SOLIDS_W_MOM, A_M, B_M)
545              ENDIF
546     !        call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
547     
548     ! evaluate local variable vxf_gs and vxf_ss
549              CALL VF_GS_Z (VXF_GS)
550              IF(DO_SOLIDS .AND. (KT_TYPE_ENUM /= GHD_2007)) THEN
551                 IF (MMAX > 0) CALL VF_SS_Z (VXF_SS)
552              ENDIF
553     
554     ! calculate coefficients for the pressure correction equation
555              IF(KT_TYPE_ENUM == GHD_2007) THEN
556                 CALL CALC_D_GHD_T (A_M, VXF_GS, D_T)
557              ELSE
558                 CALL CALC_D_T (A_M, VXF_GS, VXF_SS, D_T, IER)
559              ENDIF
560     
561              IF(DO_SOLIDS) THEN
562     ! calculate coefficients for a solids pressure correction equation
563                 IF (MMAX > 0) CALL CALC_E_T (A_M, MCP, E_T)
564     
565     ! calculate modifications to the A matrix center coefficient and B
566     ! source vector for partial elimination
567                 IF(KT_TYPE_ENUM == GHD_2007) THEN
568                    IF (MMAX > 0) CALL PARTIAL_ELIM_GHD_W (W_G, W_S, VXF_GS, A_M, B_M)
569                 ELSE
570                    IF (MMAX > 0) CALL PARTIAL_ELIM_W (W_G, W_S, VXF_GS, A_M, B_M)
571                 ENDIF
572              ENDIF
573     
574     ! handle special case where center coefficient is zero
575              CALL ADJUST_A_W_G (A_M, B_M)
576              IF(DO_SOLIDS) THEN
577                 CALL ADJUST_A_W_S (A_M, B_M)
578              ENDIF
579     
580     ! modifications to matrix equation for DEM
581              IF(DES_CONTINUUM_COUPLED) THEN
582                 CALL GAS_DRAG_W(A_M, B_M, IER)
583                 IF (DISCRETE_ELEMENT .AND. DES_CONTINUUM_HYBRID) &
584                    CALL SOLID_DRAG_W(A_M, B_M)
585              ENDIF
586     
587              IF(QMOMK .AND. QMOMK_COUPLED) THEN
588                 CALL QMOMK_GAS_DRAG(A_M, B_M, IER, 0, 0, 1)
589              ENDIF
590     
591              IF (MOMENTUM_Z_EQ(0)) THEN
592                 CALL CALC_RESID_W (U_G, V_G, W_G, A_M, B_M, 0, &
593                    NUM_RESID(RESID_W,0), DEN_RESID(RESID_W,0), &
594                    RESID(RESID_W,0), MAX_RESID(RESID_W,0), &
595                    IJK_RESID(RESID_W,0))
596                 CALL UNDER_RELAX_W (W_G, A_M, B_M, 0, UR_FAC(5))
597     !            call check_ab_m(a_m, b_m, 0, .false., ier)
598     !            write(*,*) &
599     !               resid(resid_w, 0), max_resid(resid_w, 0), &
600     !               ijk_resid(resid_w, 0)
601              ENDIF
602     
603              IF(DO_SOLIDS) THEN
604                 DO M = 1, MMAX
605                    IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
606                      (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
607                       IF (MOMENTUM_Z_EQ(M)) THEN
608                          CALL CALC_RESID_W (U_S(1,M), V_S(1,M), W_S(1,M),&
609                             A_M, B_M, M, NUM_RESID(RESID_W,M), &
610                             DEN_RESID(RESID_W,M), RESID(RESID_W,M), &
611                             MAX_RESID(RESID_W,M), IJK_RESID(RESID_W,M))
612                          CALL UNDER_RELAX_W (W_S(1,M), A_M, B_M, M, &
613                             UR_FAC(5))
614     !                     call check_ab_m(a_m, b_m, m, .false., ier)
615     !                     write(*,*) &
616     !                        resid(resid_w, m), max_resid(resid_w, m), &
617     !                        ijk_resid(resid_w, m)
618     !                     call write_ab_m(a_m, b_m, ijkmax2, m, ier)
619                       ENDIF   ! end if (momentum_z_eq(m))
620                    ENDIF ! end if check for GHD Theory
621                 ENDDO   ! end do (m=1,mmax)
622              ENDIF
623     
624              IF (MOMENTUM_Z_EQ(0)) THEN
625     !            call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), &
626     !               1, DO_K, ier)
627                 CALL ADJUST_LEQ (RESID(RESID_W,0), LEQ_IT(5), &
628                    LEQ_METHOD(5), LEQI, LEQM)
629                 CALL SOLVE_LIN_EQ ('W_g', 5, W_G_tmp, A_M, B_M, 0, LEQI,&
630                    LEQM, LEQ_SWEEP(5), LEQ_TOL(5), LEQ_PC(5), IER)
631     !            call out_array(w_g, 'w_g')
632              ENDIF
633     
634              IF(DO_SOLIDS) THEN
635                 DO M = 1, MMAX
636                    IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
637                      (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
638                       IF (MOMENTUM_Z_EQ(M)) THEN
639     !                     call test_lin_eq(ijkmax2, ijmax2, imax2, &
640     !                        a_m(1, -3, M), 1, DO_K, ier)
641                          CALL ADJUST_LEQ (RESID(RESID_W,M), LEQ_IT(5), &
642                             LEQ_METHOD(5), LEQI, LEQM)
643                          CALL SOLVE_LIN_EQ ('W_s', 5, W_S_tmp(1,M), &
644                             A_M, B_M, M, LEQI, LEQM, LEQ_SWEEP(5), &
645                             LEQ_TOL(5), LEQ_PC(5), IER)
646     !                     call out_array(w_s(1,m), 'w_s')
647                       ENDIF   ! end if (momentum_z_eq(m))
648                    ENDIF ! end if check for GHD Theory
649                 ENDDO   ! end do (m=1,mmax)
650              ENDIF
651           ENDIF   ! end if (do_k)
652     ! End W_m_star and residuals
653     ! ----------------------------------------------------------------<<<
654     
655           deallocate(vxf_gs)
656           deallocate(vxf_ss)
657           deallocate(a_m)
658           deallocate(b_m)
659     
660         END SUBROUTINE W_M_STAR
661     
662       END SUBROUTINE SOLVE_VEL_STAR
663