File: RELATIVE:/../../../mfix.git/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     !  Reviewer:                                          Date:            C
8     !                                                                      C
9     !  Revision Number: 2                                                  C
10     !  Purpose: allow multiparticle in D_E, D_N and D_T claculations       C
11     !           and account for the averaged Solid-Solid drag              C
12     !  Author: S. Dartevelle, LANL                        Date: 28-FEb-04  C
13     !  Reviewer:                                          Date: dd-mmm-yy  C
14     !                                                                      C
15     !  Revision Number: 3                                                  C
16     !  Purpose: To flag the solids calculations: Continuum or DES          C
17     !           and to change the gas arrays incorporating the drag        C
18     !           when doing DES                                             C
19     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
20     !                                                                      C
21     !  Revision Number: 4                                                  C
22     !  Purpose: Incorporation of QMOM for the solution of the particle     C
23     !           kinetic equation                                           C
24     !  Author: Alberto Passalacqua - Fox Research Group   Date: 02-Dec-09  C
25     !                                                                      C
26     !  Literature/Document References:                                     C
27     !                                                                      C
28     !  Variables referenced:                                               C
29     !  Variables modified:                                                 C
30     !  Local variables:                                                    C
31     !                                                                      C
32     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
33     
34           SUBROUTINE SOLVE_VEL_STAR(IER)
35     
36     !-----------------------------------------------
37     ! Modules
38     !-----------------------------------------------
39           USE param
40           USE param1
41           USE toleranc
42           USE run
43           USE physprop
44           USE geometry
45           USE fldvar
46           USE ghdtheory
47           USE output
48           USE indices
49           USE drag
50           USE residual
51           USE ur_facs
52           USE pgcor
53           USE pscor
54           USE leqsol
55           Use ambm
56           Use tmp_array1,  VxF_gs => Arraym1
57           Use tmp_array,  VxF_ss => ArrayLM
58           USE compar
59           USE discretelement
60           USE qmom_kinetic_equation
61           use ps
62     
63           IMPLICIT NONE
64     !-----------------------------------------------
65     ! Dummy arguments
66     !-----------------------------------------------
67     ! Error index
68           INTEGER, INTENT(INOUT) :: IER
69     !-----------------------------------------------
70     ! Local variables
71     !-----------------------------------------------
72     ! solids phase index
73           INTEGER :: M
74     ! fluid cell index
75           INTEGER :: IJK
76     ! temporary velocity arrays
77           DOUBLE PRECISION, DIMENSION(:), allocatable :: U_gtmp,  V_gtmp, W_gtmp
78           DOUBLE PRECISION, DIMENSION(:,:), allocatable :: U_stmp, V_stmp, W_stmp
79     ! linear equation solver method and iterations
80           INTEGER :: LEQM, LEQI
81     
82           LOGICAL :: DO_SOLIDS
83     
84     ! temporary use of global arrays:
85     ! arraym1 (locally vxf_gs)
86     ! the volume x average gas-solids drag at momentum cell centers
87     !      DOUBLE PRECISION :: VXF_GS(DIMENSION_3, DIMENSION_M)
88     ! arraylm (locally vxf_ss)
89     ! the volume x average solids-solids drag at momentum cell centers
90     !      DOUBLE PRECISION :: VXF_SS(DIMENSION_3, DIMENSION_LM)
91     ! Septadiagonal matrix A_m, vector b_m
92     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
93     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
94     
95     !-----------------------------------------------
96     
97           allocate(U_gtmp(DIMENSION_3))
98           allocate(V_gtmp(DIMENSION_3))
99           allocate(W_gtmp(DIMENSION_3))
100           allocate(U_stmp(DIMENSION_3,DIMENSION_M))
101           allocate(V_stmp(DIMENSION_3,DIMENSION_M))
102           allocate(W_stmp(DIMENSION_3,DIMENSION_M))
103     
104           call lock_ambm        ! locks arrys a_m and b_m
105           call lock_tmp_array1  ! locks arraym1 (locally vxf_gs)
106           call lock_tmp_array2  ! locks arraylm (locally vxf_ss)
107     
108     ! Store the velocities so that the order of solving the momentum
109     ! equations does not matter
110           DO IJK = ijkstart3, ijkend3
111              U_gtmp(IJK) = U_g(IJK)
112              V_gtmp(IJK) = V_g(IJK)
113              W_gtmp(IJK) = W_g(IJK)
114           ENDDO
115           DO M = 1, MMAX
116              IF(TRIM(KT_TYPE) /= 'GHD' .OR. &
117                 (TRIM(KT_TYPE) == 'GHD' .AND. M==MMAX)) THEN
118                 DO IJK = ijkstart3, ijkend3
119                    U_stmp(IJK, M) = U_s(IJK, M)
120                    V_stmp(IJK, M) = V_s(IJK, M)
121                    W_stmp(IJK, M) = W_s(IJK, M)
122                 ENDDO
123              ENDIF
124           ENDDO
125     
126           DO_SOLIDS = .NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
127              DES_CONTINUUM_HYBRID
128     
129     
130     ! Calculate U_m_star and residuals
131     ! ---------------------------------------------------------------->>>
132           DO M = 0, MMAX
133              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
134              IF (M >= 1) VXF_GS(:,M) = ZERO
135           ENDDO
136     
137     ! calculate the convection-diffusion terms for the gas and solids phase
138     ! u-momentum equations
139           CALL CONV_DIF_U_G (A_M, B_M)
140           IF(DO_SOLIDS) CALL CONV_DIF_U_S (A_M, B_M)
141     
142     ! calculate the source terms for the gas and solids phase u-momentum
143     ! equations
144           CALL SOURCE_U_G (A_M, B_M)
145           IF(POINT_SOURCE) CALL POINT_SOURCE_U_G (A_M, B_M)
146           IF(DO_SOLIDS) THEN
147              CALL SOURCE_U_S (A_M, B_M)
148              IF(POINT_SOURCE) CALL POINT_SOURCE_U_S (A_M, B_M)
149           ENDIF
150     
151     ! evaluate local variable vxf_gs and vxf_ss.  both terms are sent to the
152     ! subroutine calc_d (pressure correction equation coefficients).  the
153     ! former is also used in the subroutine partial_elim_u while the latter
154     ! is effectively re-evaluated within said subroutine
155           CALL VF_GS_X (VXF_GS)
156           IF(DO_SOLIDS .AND. (TRIM(KT_TYPE) /= 'GHD')) THEN
157              IF (MMAX > 0) CALL VF_SS_X (VXF_SS)
158           ENDIF
159     
160     ! calculate coefficients for the pressure correction equation
161           IF(TRIM(KT_TYPE) == 'GHD') THEN
162              CALL CALC_D_GHD_E (A_M, VXF_GS, D_E)
163           ELSE
164              CALL CALC_D_E (A_M, VXF_GS, VXF_SS, D_E, IER)
165           ENDIF
166     
167           IF(DO_SOLIDS) THEN
168     ! calculate coefficients for a solids volume correction equation
169              IF (MMAX > 0) CALL CALC_E_E (A_M, MCP, E_E)
170     
171     ! calculate modifications to the A matrix center coefficient and B
172     ! source vector for partial elimination
173              IF(TRIM(KT_TYPE) == 'GHD') THEN
174                 IF (MMAX > 0) CALL PARTIAL_ELIM_GHD_U (U_G,U_S,VXF_GS,A_M,B_M)
175              ELSE
176                 IF (MMAX > 0) CALL PARTIAL_ELIM_U (U_G,U_S,VXF_GS,A_M,B_M)
177              ENDIF
178           ENDIF
179     
180     ! handle special case where center coefficient is zero
181           CALL ADJUST_A_U_G (A_M, B_M)
182           IF(DO_SOLIDS) CALL ADJUST_A_U_S (A_M, B_M)
183     
184     ! calculate modifications to the A matrix center coefficient and B
185     ! source vector for treating DEM drag terms
186           IF(DES_CONTINUUM_COUPLED) THEN
187              CALL GAS_DRAG_U(A_M, B_M, IER)
188              IF (DES_CONTINUUM_HYBRID) &
189                 CALL SOLID_DRAG_U(A_M, B_M)
190           ENDIF
191     
192           IF(QMOMK .AND. QMOMK_COUPLED) THEN
193              CALL QMOMK_GAS_DRAG(A_M, B_M, IER, 1, 0, 0)
194           ENDIF
195     
196           IF (MOMENTUM_X_EQ(0)) THEN
197              CALL CALC_RESID_U (U_G, V_G, W_G, A_M, B_M, 0, &
198                 NUM_RESID(RESID_U,0), DEN_RESID(RESID_U,0), &
199                 RESID(RESID_U,0), MAX_RESID(RESID_U,0), &
200                 IJK_RESID(RESID_U,0))
201              CALL UNDER_RELAX_U (U_G, A_M, B_M, 0, UR_FAC(3))
202     !         call check_ab_m(a_m, b_m, 0, .false., ier)
203     !         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
204     !         write(*,*) &
205     !            resid(resid_u, 0), max_resid(resid_u, 0), &
206     !            ijk_resid(resid_u, 0)
207           ENDIF
208     
209           IF(DO_SOLIDS) THEN
210              DO M = 1, MMAX
211                 IF(TRIM(KT_TYPE) /= 'GHD' .OR. &
212                   (TRIM(KT_TYPE) == 'GHD' .AND. M==MMAX)) THEN
213                    IF (MOMENTUM_X_EQ(M)) THEN
214                       CALL CALC_RESID_U (U_S(1,M), V_S(1,M), W_S(1,M), A_M,&
215                          B_M, M, NUM_RESID(RESID_U,M), &
216                          DEN_RESID(RESID_U,M), RESID(RESID_U,M), &
217                          MAX_RESID(RESID_U,M), IJK_RESID(RESID_U,M))
218                       CALL UNDER_RELAX_U (U_S(1,M), A_M, B_M, M, &
219                          UR_FAC(3))
220     !                  call check_ab_m(a_m, b_m, m, .false., ier)
221     !                  write(*,*) &
222     !                     resid(resid_u, m), max_resid(resid_u, m), &
223     !                     ijk_resid(resid_u, m)
224     !                  call write_ab_m(a_m, b_m, ijkmax2, m, ier)
225                    ENDIF   ! end if (momentum_x_eq(m))
226                 ENDIF ! end if check for GHD Theory
227              ENDDO   ! end do (m=1,mmax)
228           ENDIF
229     
230           IF (MOMENTUM_X_EQ(0)) THEN
231     !         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1,&
232     !            DO_K,ier)
233              CALL ADJUST_LEQ (RESID(RESID_U,0), LEQ_IT(3), LEQ_METHOD(3), &
234                 LEQI, LEQM)
235              CALL SOLVE_LIN_EQ ('U_g', 3, U_Gtmp, A_M, B_M, 0, LEQI, LEQM, &
236                 LEQ_SWEEP(3), LEQ_TOL(3),  LEQ_PC(3), IER)
237     !         call out_array(u_g, 'u_g')
238           ENDIF
239     
240           IF(DO_SOLIDS) THEN
241              DO M = 1, MMAX
242                 IF(TRIM(KT_TYPE) /= 'GHD' .OR. &
243                   (TRIM(KT_TYPE) == 'GHD' .AND. M==MMAX)) THEN
244                    IF (MOMENTUM_X_EQ(M)) THEN
245     !                  call test_lin_eq(ijkmax2, ijmax2, imax2, &
246     !                     a_m(1, -3, M), 1, DO_K,ier)
247                       CALL ADJUST_LEQ (RESID(RESID_U,M), LEQ_IT(3),&
248                          LEQ_METHOD(3), LEQI, LEQM)
249                       CALL SOLVE_LIN_EQ ('U_s', 3, U_Stmp(1,M), A_M, &
250                          B_M, M, LEQI, LEQM, LEQ_SWEEP(3), LEQ_TOL(3),&
251                          LEQ_PC(3), IER)
252     !                  call out_array(u_s(1,m), 'u_s')
253                    ENDIF   ! end if (momentum_x_eq(m))
254                 ENDIF ! end if check for GHD Theory
255              ENDDO   ! end do (m=1,mmax)
256           ENDIF
257     ! End U_m_star and residuals
258     ! ----------------------------------------------------------------<<<
259     
260     ! Calculate V_m_star and residuals
261     ! ---------------------------------------------------------------->>>
262           DO M = 0, MMAX
263              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
264              IF (M >= 1) VXF_GS(:,M) = ZERO
265           ENDDO
266     
267           CALL CONV_DIF_V_G (A_M, B_M, IER)
268     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
269           IF(DO_SOLIDS) CALL CONV_DIF_V_S (A_M, B_M, IER)
270     
271           CALL SOURCE_V_G (A_M, B_M)
272           IF(POINT_SOURCE) CALL POINT_SOURCE_V_G (A_M, B_M)
273     
274     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
275           IF(DO_SOLIDS) THEN
276              CALL SOURCE_V_S (A_M, B_M)
277              IF(POINT_SOURCE) CALL POINT_SOURCE_V_S (A_M, B_M)
278           END IF
279     
280           CALL VF_GS_Y (VXF_GS)
281           IF(DO_SOLIDS .AND. (TRIM(KT_TYPE) /= 'GHD')) THEN
282              IF (MMAX > 0) CALL VF_SS_Y (VXF_SS)
283           ENDIF
284     
285     ! calculate coefficients for the pressure correction equation
286           IF(TRIM(KT_TYPE) == 'GHD') THEN
287              CALL CALC_D_GHD_N (A_M, VXF_GS, D_N)
288           ELSE
289              CALL CALC_D_N (A_M, VXF_GS, VXF_SS, D_N, IER)
290           ENDIF
291     
292           IF(DO_SOLIDS) THEN
293     ! calculate coefficients for a solids pressure correction equation
294              IF (MMAX > 0) CALL CALC_E_N (A_M, MCP, E_N)
295     
296     ! calculate modifications to the A matrix center coefficient and B
297     ! source vector for partial elimination
298              IF(TRIM(KT_TYPE) == 'GHD') THEN
299                IF (MMAX > 0) CALL PARTIAL_ELIM_GHD_V (V_G,V_S,VXF_GS,A_M,B_M)
300              ELSE
301                IF (MMAX > 0) CALL PARTIAL_ELIM_V (V_G,V_S,VXF_GS,A_M,B_M)
302              ENDIF
303           ENDIF
304     
305     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
306           CALL ADJUST_A_V_G (A_M, B_M)
307     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
308           IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
309              DES_CONTINUUM_HYBRID) THEN
310              CALL ADJUST_A_V_S (A_M, B_M)
311           ENDIF
312     
313           IF(DES_CONTINUUM_COUPLED) THEN
314              CALL GAS_DRAG_V(A_M, B_M, IER)
315              IF (DES_CONTINUUM_HYBRID) &
316                 CALL SOLID_DRAG_V(A_M, B_M)
317           ENDIF
318     
319           IF(QMOMK .AND. QMOMK_COUPLED) THEN
320              CALL QMOMK_GAS_DRAG(A_M, B_M, IER, 0, 1, 0)
321           ENDIF
322     
323     
324           IF (MOMENTUM_Y_EQ(0)) THEN
325              CALL CALC_RESID_V (U_G, V_G, W_G, A_M, B_M, 0, &
326                 NUM_RESID(RESID_V,0), DEN_RESID(RESID_V,0), &
327                 RESID(RESID_V,0), MAX_RESID(RESID_V,0), &
328                 IJK_RESID(RESID_V,0))
329              CALL UNDER_RELAX_V (V_G, A_M, B_M, 0, UR_FAC(4))
330     !         call check_ab_m(a_m, b_m, 0, .false., ier)
331     !         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
332     !         write(*,*) &
333     !            resid(resid_v, 0), max_resid(resid_v, 0), &
334     !            ijk_resid(resid_v, 0)
335           ENDIF
336     
337           IF(DO_SOLIDS) THEN
338              DO M = 1, MMAX
339                 IF(TRIM(KT_TYPE) /= 'GHD' .OR. &
340                   (TRIM(KT_TYPE) == 'GHD' .AND. M==MMAX)) THEN
341                    IF (MOMENTUM_Y_EQ(M)) THEN
342                       CALL CALC_RESID_V (U_S(1,M), V_S(1,M), W_S(1,M), A_M,&
343                          B_M, M, NUM_RESID(RESID_V,M), &
344                          DEN_RESID(RESID_V,M),RESID(RESID_V,M), &
345                          MAX_RESID(RESID_V,M), IJK_RESID(RESID_V,M))
346                       CALL UNDER_RELAX_V (V_S(1,M),A_M,B_M,M,UR_FAC(4))
347     !                  call check_ab_m(a_m, b_m, m, .false., ier)
348     !                  write(*,*) &
349     !                     resid(resid_v, m), max_resid(resid_v, m),
350     !                     ijk_resid(resid_v, m)
351     !                  call write_ab_m(a_m, b_m, ijkmax2, m, ier)
352                    ENDIF   ! end if (momentum_y_eq(m))
353                 ENDIF ! end if check for GHD Theory
354              ENDDO   ! end do (m=1,mmax)
355           ENDIF
356     
357           IF (MOMENTUM_Y_EQ(0)) THEN
358     !         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), &
359     !            1, DO_K, ier)
360              CALL ADJUST_LEQ (RESID(RESID_V,0), LEQ_IT(4), LEQ_METHOD(4), &
361                 LEQI, LEQM)
362              CALL SOLVE_LIN_EQ ('V_g', 4, V_Gtmp, A_M, B_M, 0, LEQI, LEQM, &
363                 LEQ_SWEEP(4), LEQ_TOL(4),  LEQ_PC(4), IER)
364     !         call out_array(v_g, 'v_g')
365           ENDIF
366     
367           IF(DO_SOLIDS) THEN
368              DO M = 1, MMAX
369                 IF(TRIM(KT_TYPE) /= 'GHD' .OR. &
370                   (TRIM(KT_TYPE) == 'GHD' .AND. M==MMAX)) THEN
371                    IF (MOMENTUM_Y_EQ(M)) THEN
372     !                  call test_lin_eq(ijkmax2, ijmax2, imax2, &
373     !                     a_m(1, -3, M), 1, DO_K, ier)
374                       CALL ADJUST_LEQ (RESID(RESID_V,M), LEQ_IT(4), &
375                          LEQ_METHOD(4), LEQI, LEQM)
376                       CALL SOLVE_LIN_EQ ('V_s', 4, V_Stmp(1,M), A_M, &
377                          B_M, M, LEQI, LEQM, LEQ_SWEEP(4), LEQ_TOL(4), &
378                          LEQ_PC(4), IER)
379     !                  call out_array(v_s(1,m), 'v_s')
380                    ENDIF   ! end if (momentum_y_eq(m))
381                 ENDIF ! end if check for GHD Theory
382              ENDDO   ! end do (m=1,mmax)
383           ENDIF
384     ! End V_m_star and residuals
385     ! ----------------------------------------------------------------<<<
386     
387     
388     ! Calculate W_m_star and residuals
389     ! ---------------------------------------------------------------->>>
390           IF (DO_K)THEN
391              DO M = 0, MMAX
392                 CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
393                 IF (M >= 1) VXF_GS(:,M) = ZERO
394              ENDDO
395     
396              CALL CONV_DIF_W_G (A_M, B_M)
397              IF(DO_SOLIDS) CALL CONV_DIF_W_S (A_M, B_M)
398     
399              CALL SOURCE_W_G (A_M, B_M)
400              IF(POINT_SOURCE) CALL POINT_SOURCE_W_G (A_M, B_M)
401     !         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
402              IF(DO_SOLIDS) THEN
403                 CALL SOURCE_W_S (A_M, B_M)
404                 IF(POINT_SOURCE) CALL POINT_SOURCE_W_S (A_M, B_M)
405              ENDIF
406     !        call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
407     
408              CALL VF_GS_Z (VXF_GS)
409              IF(DO_SOLIDS .AND. (TRIM(KT_TYPE) /= 'GHD')) THEN
410                 IF (MMAX > 0) CALL VF_SS_Z (VXF_SS)
411              ENDIF
412     
413     ! calculate coefficients for the pressure correction equation
414              IF(TRIM(KT_TYPE) == 'GHD') THEN
415                 CALL CALC_D_GHD_T (A_M, VXF_GS, D_T)
416              ELSE
417                 CALL CALC_D_T (A_M, VXF_GS, VXF_SS, D_T, IER)
418              ENDIF
419     
420              IF(DO_SOLIDS) THEN
421     ! calculate coefficients for a solids pressure correction equation
422                 IF (MMAX > 0) CALL CALC_E_T (A_M, MCP, E_T)
423     
424     ! calculate modifications to the A matrix center coefficient and B
425     ! source vector for partial elimination
426                 IF(TRIM(KT_TYPE) == 'GHD') THEN
427                    IF (MMAX > 0) CALL PARTIAL_ELIM_GHD_W (W_G, W_S, VXF_GS, A_M, B_M)
428                 ELSE
429                    IF (MMAX > 0) CALL PARTIAL_ELIM_W (W_G, W_S, VXF_GS, A_M, B_M)
430                 ENDIF
431              ENDIF
432     
433              CALL ADJUST_A_W_G (A_M, B_M)
434              IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
435                 DES_CONTINUUM_HYBRID) THEN
436                 CALL ADJUST_A_W_S (A_M, B_M)
437              ENDIF
438     
439              IF(DES_CONTINUUM_COUPLED) THEN
440                 CALL GAS_DRAG_W(A_M, B_M, IER)
441                 IF (DISCRETE_ELEMENT .AND. DES_CONTINUUM_HYBRID) &
442                    CALL SOLID_DRAG_W(A_M, B_M)
443              ENDIF
444     
445              IF(QMOMK .AND. QMOMK_COUPLED) THEN
446                 CALL QMOMK_GAS_DRAG(A_M, B_M, IER, 0, 0, 1)
447              ENDIF
448     
449              IF (MOMENTUM_Z_EQ(0)) THEN
450                 CALL CALC_RESID_W (U_G, V_G, W_G, A_M, B_M, 0, &
451                    NUM_RESID(RESID_W,0), DEN_RESID(RESID_W,0), &
452                    RESID(RESID_W,0), MAX_RESID(RESID_W,0), &
453                    IJK_RESID(RESID_W,0))
454                 CALL UNDER_RELAX_W (W_G, A_M, B_M, 0, UR_FAC(5))
455     !            call check_ab_m(a_m, b_m, 0, .false., ier)
456     !            write(*,*) &
457     !               resid(resid_w, 0), max_resid(resid_w, 0), &
458     !               ijk_resid(resid_w, 0)
459              ENDIF
460     
461              IF(DO_SOLIDS) THEN
462                 DO M = 1, MMAX
463                    IF(TRIM(KT_TYPE) /= 'GHD' .OR. &
464                      (TRIM(KT_TYPE) == 'GHD' .AND. M==MMAX)) THEN
465                       IF (MOMENTUM_Z_EQ(M)) THEN
466                          CALL CALC_RESID_W (U_S(1,M), V_S(1,M), W_S(1,M),&
467                             A_M, B_M, M, NUM_RESID(RESID_W,M), &
468                             DEN_RESID(RESID_W,M), RESID(RESID_W,M), &
469                             MAX_RESID(RESID_W,M), IJK_RESID(RESID_W,M))
470                          CALL UNDER_RELAX_W (W_S(1,M), A_M, B_M, M, &
471                             UR_FAC(5))
472     !                     call check_ab_m(a_m, b_m, m, .false., ier)
473     !                     write(*,*) &
474     !                        resid(resid_w, m), max_resid(resid_w, m), &
475     !                        ijk_resid(resid_w, m)
476     !                     call write_ab_m(a_m, b_m, ijkmax2, m, ier)
477                       ENDIF   ! end if (momentum_z_eq(m))
478                    ENDIF ! end if check for GHD Theory
479                 ENDDO   ! end do (m=1,mmax)
480              ENDIF
481     
482              IF (MOMENTUM_Z_EQ(0)) THEN
483     !            call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), &
484     !               1, DO_K, ier)
485                 CALL ADJUST_LEQ (RESID(RESID_W,0), LEQ_IT(5), &
486                    LEQ_METHOD(5), LEQI, LEQM)
487                 CALL SOLVE_LIN_EQ ('W_g', 5, W_Gtmp, A_M, B_M, 0, LEQI,&
488                    LEQM, LEQ_SWEEP(5), LEQ_TOL(5), LEQ_PC(5), IER)
489     !            call out_array(w_g, 'w_g')
490              ENDIF
491     
492              IF(DO_SOLIDS) THEN
493                 DO M = 1, MMAX
494                    IF(TRIM(KT_TYPE) /= 'GHD' .OR. &
495                      (TRIM(KT_TYPE) == 'GHD' .AND. M==MMAX)) THEN
496                       IF (MOMENTUM_Z_EQ(M)) THEN
497     !                     call test_lin_eq(ijkmax2, ijmax2, imax2, &
498     !                        a_m(1, -3, M), 1, DO_K, ier)
499                          CALL ADJUST_LEQ (RESID(RESID_W,M), LEQ_IT(5), &
500                             LEQ_METHOD(5), LEQI, LEQM)
501                          CALL SOLVE_LIN_EQ ('W_s', 5, W_Stmp(1,M), &
502                             A_M, B_M, M, LEQI, LEQM, LEQ_SWEEP(5), &
503                             LEQ_TOL(5), LEQ_PC(5), IER)
504     !                     call out_array(w_s(1,m), 'w_s')
505                       ENDIF   ! end if (momentum_z_eq(m))
506                    ENDIF ! end if check for GHD Theory
507                 ENDDO   ! end do (m=1,mmax)
508              ENDIF
509           ENDIF   ! end if (do_k)
510     ! End W_m_star and residuals
511     ! ----------------------------------------------------------------<<<
512     
513     ! Now update all velocity components
514           DO IJK = ijkstart3, ijkend3
515              U_g(IJK) = U_gtmp(IJK)
516              V_g(IJK) = V_gtmp(IJK)
517              W_g(IJK) = W_gtmp(IJK)
518           ENDDO
519           DO M = 1, MMAX
520              IF(TRIM(KT_TYPE) /= 'GHD' .OR. &
521                (TRIM(KT_TYPE) == 'GHD' .AND. M==MMAX)) THEN
522                 DO IJK = ijkstart3, ijkend3
523                    U_s(IJK, M) = U_stmp(IJK, M)
524                    V_s(IJK, M) = V_stmp(IJK, M)
525                    W_s(IJK, M) = W_stmp(IJK, M)
526                 ENDDO
527              ENDIF
528           ENDDO
529     
530     ! modification for GHD theory to compute species velocity: Ui = Joi/(mi ni) + U.
531           IF(TRIM(KT_TYPE) == 'GHD') THEN
532              CALL calc_external_forces()
533              CALL GHDMassFlux() ! to compute solid species mass flux
534              CALL UpdateSpeciesVelocities() ! located at end of ghdMassFlux.f file
535           ENDIF
536     
537           call unlock_ambm
538           call unlock_tmp_array1
539           call unlock_tmp_array2
540     
541           deallocate(U_gtmp)
542           deallocate(V_gtmp)
543           deallocate(W_gtmp)
544           deallocate(U_stmp)
545           deallocate(V_stmp)
546           deallocate(W_stmp)
547     
548           RETURN
549           END SUBROUTINE SOLVE_VEL_STAR
550