File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/ghdmassflux.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: GHDMASSFLUX                                             C
4     !  Purpose: Calculate the species mass flux 3-components of Joi at cellC
5     !           faces to compute species velocities and source terms in T  C
6     !           equation.                                                  C
7     !                                                                      C
8     !  Author: S. Benyahia                              Date: 13-MAR-09    C
9     !  Reviewer:                                          Date:            C
10     !                                                                      C
11     !  Literature/Document References:                                     C
12     !     C. Hrenya handnotes and Garzo, Hrenya, Dufty papers (PRE, 2007)  C
13     !                                                                      C
14     !  Variables modified:   JoiX  JoiY   JoiZ                             C
15     !                                                                      C
16     !  Local variables: all terms in mass flux                             C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     
20           SUBROUTINE GHDMASSFLUX ()
21     
22     !-----------------------------------------------
23     ! Modules
24     !-----------------------------------------------
25           USE param
26           USE param1
27           USE geometry
28           USE compar
29           USE fldvar
30           USE indices
31           USE visc_s
32           USE ghdtheory
33           USE physprop
34           USE run
35           USE constant
36           USE toleranc
37           USE drag
38           USE fun_avg
39           USE functions
40           IMPLICIT NONE
41     !-----------------------------------------------
42     ! Local variables
43     !-----------------------------------------------
44     ! Index
45           INTEGER :: IJK, I, J, K
46     ! Index
47           INTEGER :: IJKE, IJKN, IJKT
48     ! Solids phase
49           INTEGER :: M, L
50     ! number densities to compute del(Nj)
51           DOUBLE PRECISION NjC, NjE, NjN, NjT
52     
53     ! mass of species
54           DOUBLE PRECISION :: Mi, Mj
55     ! number density of species
56           DOUBLE PRECISION :: Ni
57     ! mixture density and temperature at cell faces
58           DOUBLE PRECISION :: ropsE, ropsN, ropsT, ThetaE, ThetaN, ThetaT
59           DOUBLE PRECISION :: EPSA
60     ! transport coefficient at cell faces
61           DOUBLE PRECISION :: DiTE, DiTN, DiTT
62           DOUBLE PRECISION :: DijE, DijN, DijT
63           DOUBLE PRECISION :: DijFE, DijFN, DijFT
64     ! Terms in the calculation of Joi-X,Y,Z
65           DOUBLE PRECISION :: ordinDiffTermX, ordinDiffTermY, ordinDiffTermZ
66           DOUBLE PRECISION :: massMobilityTermX, massMobilityTermY, &
67                               massMobilityTermZ
68           DOUBLE PRECISION :: massMobilityTermXvelupdate, &
69                               massMobilityTermYvelupdate, &
70                               massMobilityTermZvelupdate
71           DOUBLE PRECISION :: thermalDiffTermX, thermalDiffTermY, &
72                               thermalDiffTermZ
73           DOUBLE PRECISION :: ropsme, ropsmn, ropsmt
74     
75           DOUBLE PRECISION :: addtermx, &
76                               addtermy, addtermz
77           DOUBLE PRECISION :: massMobilityTermNoDragX, &
78                               massMobilityTermNoDragY, &
79                               massMobilityTermNoDragZ
80           DOUBLE PRECISION :: DiTE_H, DiTE_A, DiTN_H, DiTN_A, DiTT_H, DiTT_A
81           DOUBLE PRECISION :: DijE_H, DijE_A, DijN_H, DijN_A, DijT_H, DijT_A
82           DOUBLE PRECISION :: DijFE_H, DijFE_A, DijFN_H, DijFN_A, DijFT_H, DijFT_A
83     !-----------------------------------------------
84     !     Function subroutines
85     !-----------------------------------------------
86     
87           DO M = 1, SMAX
88             DO 200 IJK = ijkstart3, ijkend3
89               I = I_OF(IJK)
90               J = J_OF(IJK)
91               K = K_OF(IJK)
92     
93               IF ( FLUID_AT(IJK) ) THEN
94                 Mi = (PI/6.d0)*D_P(IJK,M)**3 * RO_S(IJK,M)
95                 Ni = ROP_s(IJK,M) / Mi
96     
97                 IJKE = EAST_OF(IJK)
98                 IJKN = NORTH_OF(IJK)
99                 IJKT = TOP_OF(IJK)
100     
101     ! mixture density and temperature evaluated at cell faces
102                 ropsE = AVG_X(ROP_S(IJK,MMAX),ROP_S(IJKE,MMAX),I)
103                 ropsN = AVG_Y(ROP_S(IJK,MMAX),ROP_S(IJKN,MMAX),J)
104                 ropsT = AVG_Z(ROP_S(IJK,MMAX),ROP_S(IJKT,MMAX),K)
105     
106                 ThetaE = AVG_X(THETA_M(IJK,MMAX),THETA_M(IJKE,MMAX),I)
107                 ThetaN = AVG_Y(THETA_M(IJK,MMAX),THETA_M(IJKN,MMAX),J)
108                 ThetaT = AVG_Z(THETA_M(IJK,MMAX),THETA_M(IJKT,MMAX),K)
109     
110     ! Thermal diffusion evaluated at cell faces (all used transport coef.
111     ! will be evaluated this way)
112                 DiTE_H = AVG_X_S(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
113                    DiT(IJKE,M)*ROP_S(IJKE,MMAX)/Theta_m(IJKE,MMAX),I)
114                 DiTE_A = AVG_X(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
115                     DiT(IJKE,M)*ROP_S(IJKE,MMAX)/Theta_m(IJKE,MMAX),I)
116     
117                 DiTN_H = AVG_Y_S(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
118                     DiT(IJKN,M)*ROP_S(IJKN,MMAX)/Theta_m(IJKN,MMAX),J)
119                 DiTN_A = AVG_Y(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
120                      DiT(IJKN,M)*ROP_S(IJKN,MMAX)/Theta_m(IJKN,MMAX),J)
121     
122                 DiTT_H = AVG_Z_S(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
123                     DiT(IJKT,M)*ROP_S(IJKT,MMAX)/Theta_m(IJKT,MMAX),K)
124                 DiTT_A = AVG_Z(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
125                     DiT(IJKT,M)*ROP_S(IJKT,MMAX)/Theta_m(IJKT,MMAX),K)
126     
127                 IF(M .eq. 1)THEN
128                   IF(MIN(ABS(DiTE_H),ABS(DiTE_A)) .eq. ABS(DiTE_H))THEN
129                     DiTE = DiTE_H
130                     DiT_HarmE(IJK) = .TRUE.
131                   ELSE
132                     DiTE = DiTE_A
133                     DiT_HarmE(IJK) = .FALSE.
134                   ENDIF
135     
136                   IF(MIN(ABS(DiTN_H),ABS(DiTN_A)) .eq. ABS(DiTN_H))THEN
137                     DiTN = DiTN_H
138                     DiT_HarmN(IJK) = .TRUE.
139                   ELSE
140                     DiTN = DiTN_A
141                     DiT_HarmN(IJK) = .FALSE.
142                   ENDIF
143     
144                   IF(MIN(ABS(DiTT_H),ABS(DiTT_A)) .eq. ABS(DiTT_H))THEN
145                     DiTT = DiTT_H
146                     DiT_HarmT(IJK) = .TRUE.
147                   ELSE
148                     DiTT = DiTT_A
149                     DiT_HarmT(IJK) = .FALSE.
150                   ENDIF
151                 ELSE
152                   IF(DiT_HarmE(IJK))THEN
153                     DiTE = DiTE_H
154                   ELSE
155                     DiTE = DiTE_A
156                   ENDIF
157     
158                   IF(DiT_HarmN(IJK))THEN
159                     DiTN = DiTN_H
160                   ELSE
161                     DiTN = DiTN_A
162                   ENDIF
163     
164                   IF(DiT_HarmT(IJK))THEN
165                     DiTT = DiTT_H
166                   ELSE
167                     DiTT = DiTT_A
168                   ENDIF
169                 ENDIF   ! end if/else M=1
170     
171     ! initializing variables for summation over L
172                 ordinDiffTermX = ZERO
173                 ordinDiffTermY = ZERO
174                 ordinDiffTermZ = ZERO
175                 massMobilityTermX = ZERO
176                 massMobilityTermY = ZERO
177                 massMobilityTermZ = ZERO
178                 massMobilityTermXvelUpdate = ZERO
179                 massMobilityTermYvelUpdate = ZERO
180                 massMobilityTermZvelUpdate = ZERO
181                 addtermx = ZERO
182                 addtermy = ZERO
183                 addtermz = ZERO
184                 massMobilityTermNoDragX = ZERO
185                 massMobilityTermNoDragY = ZERO
186                 massMobilityTermNoDragZ = ZERO
187     
188                 DO L = 1, SMAX
189                   Mj  = (PI/6.d0)*D_P(IJK,L)**3 * RO_S(IJK,L)
190     
191                   NjC = ROP_s(IJK,L) / Mj
192                   NjE = ROP_S(IJKE,L) / Mj
193                   NjN = ROP_S(IJKN,L) / Mj
194                   NjT = ROP_S(IJKT,L) / Mj
195     
196                   IF((ROP_S(IJK,MMAX)/RO_S(IJK,M) > DIL_EP_S) .and. &
197                     (ROP_S(IJKE,MMAX)/RO_S(IJK,M) > DIL_EP_S))THEN
198                     DijE_H = AVG_X_S(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
199                        Dij(IJKE,M,L)*Mi*Mj/ROP_S(IJKE,MMAX),I)
200                     DijE_A = AVG_X(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
201                        Dij(IJKE,M,L)*Mi*Mj/ROP_S(IJKE,MMAX),I)
202     
203                     IF(M .eq. 1)THEN
204                       IF(MIN(ABS(DijE_H),ABS(DijE_A)) .eq. ABS(DijE_H))THEN
205                         DijE = DijE_H
206                         Dij_HarmE(IJK,L) = .TRUE.
207                       ELSE
208                         DijE = DijE_A
209                         Dij_HarmE(IJK,L) = .FALSE.
210                       ENDIF
211                     ELSE
212                       IF(Dij_HarmE(IJK,L))THEN
213                         DijE = DijE_H
214                       ELSE
215                         DijE = DijE_A
216                       ENDIF
217                     ENDIF
218                   ELSE
219                     DijE = ZERO
220                   ENDIF
221     
222                   IF((ROP_S(IJK,MMAX)/RO_S(IJK,M) > DIL_EP_S) .and. &
223                      (ROP_S(IJKN,MMAX)/RO_S(IJK,M) > DIL_EP_S))THEN
224                     DijN_H = AVG_Y_S(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
225                       Dij(IJKN,M,L)*Mi*Mj/ROP_S(IJKN,MMAX),J)
226                     DijN_A = AVG_Y(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
227                       Dij(IJKN,M,L)*Mi*Mj/ROP_S(IJKN,MMAX),J)
228                     IF(M .eq. 1)THEN
229                       IF(MIN(ABS(DijN_H),ABS(DijN_A)) .eq. ABS(DijN_H))THEN
230                         DijN = DijN_H
231                         Dij_HarmN(IJK,L) = .TRUE.
232                       ELSE
233                         DijN = DijN_A
234                         Dij_HarmN(IJK,L) = .FALSE.
235                       ENDIF
236                     ELSE
237                       IF(Dij_HarmN(IJK,L))THEN
238                         DijN = DijN_H
239                       ELSE
240                         DijN = DijN_A
241                       ENDIF
242                     ENDIF
243                   ELSE
244                     DijN = ZERO
245                   ENDIF
246     
247                   IF((ROP_S(IJK,MMAX)/RO_S(IJK,M) > DIL_EP_S) .and. &
248                      (ROP_S(IJKT,MMAX)/RO_S(IJK,M) > DIL_EP_S))THEN
249                     DijT_H = AVG_Z_S(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
250                       Dij(IJKT,M,L)*Mi*Mj/ROP_S(IJKT,MMAX),K)
251                     DijT_A = AVG_Z(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
252                       Dij(IJKT,M,L)*Mi*Mj/ROP_S(IJKT,MMAX),K)
253                     IF(M .eq. 1)THEN
254                       IF(MIN(ABS(DijT_H),ABS(DijT_A)) .eq. ABS(DijT_H))THEN
255                         DijT = DijT_H
256                         Dij_HarmT(IJK,L) = .TRUE.
257                       ELSE
258                         DijT = DijT_A
259                         Dij_HarmT(IJK,L) = .FALSE.
260                       ENDIF
261                     ELSE
262                       IF(Dij_HarmT(IJK,L))THEN
263                         DijT = DijT_H
264                       ELSE
265                         DijT = DijT_A
266                       ENDIF
267                     ENDIF
268                   ELSE
269                     DijT = ZERO
270                   ENDIF
271     
272     
273                   DijFE_H = AVG_X_S(DijF(IJK,M,L),DijF(IJKE,M,L),I)
274                   DijFE_A = AVG_X(DijF(IJK,M,L),DijF(IJKE,M,L),I)
275                   DijFN_H = AVG_Y_S(DijF(IJK,M,L),DijF(IJKN,M,L),J)
276                   DijFN_A = AVG_Y(DijF(IJK,M,L),DijF(IJKN,M,L),J)
277                   DijFT_H = AVG_Z_S(DijF(IJK,M,L),DijF(IJKT,M,L),K)
278                   DijFT_A = AVG_Z(DijF(IJK,M,L),DijF(IJKT,M,L),K)
279     
280                   IF(M .eq. 1)THEN
281                     IF(MIN(ABS(DijFE_H),ABS(DijFE_A)) .eq. ABS(DijFE_H))THEN
282                       DijFE = DijFE_H
283                       DijF_HarmE(IJK,L) = .TRUE.
284                     ELSE
285                       DijFE = DijFE_A
286                       DijF_HarmE(IJK,L) = .FALSE.
287                     ENDIF
288     
289                     IF(MIN(ABS(DijFN_H),ABS(DijFN_A)) .eq. ABS(DijFN_H))THEN
290                       DijFN = DijFN_H
291                       DijF_HarmN(IJK,L) = .TRUE.
292                     ELSE
293                       DijFN = DijFN_A
294                       DijF_HarmN(IJK,L) = .FALSE.
295                     ENDIF
296     
297                     IF(MIN(ABS(DijFT_H),ABS(DijFT_A)) .eq. ABS(DijFT_H))THEN
298                       DijFT = DijFT_H
299                       DijF_HarmT(IJK,L) = .TRUE.
300                     ELSE
301                       DijFT = DijFT_A
302                       DijF_HarmT(IJK,L) = .FALSE.
303                     ENDIF
304                   ELSE
305                     IF(DijF_HarmE(IJK,L))THEN
306                       DijFE = DijFE_H
307                     ELSE
308                       DijFE = DijFE_A
309                     ENDIF
310     
311                     IF(DijF_HarmN(IJK,L))THEN
312                       DijFN = DijFN_H
313                     ELSE
314                       DijFN = DijFN_A
315                     ENDIF
316     
317                     IF(DijF_HarmT(IJK,L))THEN
318                       DijFT = DijFT_H
319                     ELSE
320                       DijFT = DijFT_A
321                     ENDIF
322                   ENDIF   ! end if/else (M=1)
323     
324                   ordinDiffTermX = ordinDiffTermX + DijE * (NjE - NjC) * oDX_E(I)
325                   ordinDiffTermY = ordinDiffTermY + DijN * (NjN - NjC) * oDY_N(J)
326                   ordinDiffTermZ = ordinDiffTermZ + DijT * (NjT - NjC) * (oX_E(I)*oDZ_T(K))
327     
328                   massMobilityTermX = massMobilityTermX + DijFE * FiX(IJK,L)
329                   massMobilityTermY = massMobilityTermY + DijFN * FiY(IJK,L)
330                   massMobilityTermZ = massMobilityTermZ + DijFT * FiZ(IJK,L)
331     
332                   massMobilityTermXvelUpdate = massMobilityTermXvelUpdate + DijFE * FiXvel(IJK,L)
333                   massMobilityTermYvelUpdate = massMobilityTermYvelUpdate + DijFN * FiYvel(IJK,L)
334                   massMobilityTermZvelUpdate = massMobilityTermZvelUpdate + DijFT * FiZvel(IJK,L)
335     
336                   massMobilityTermNoDragX = massMobilityTermNoDragX + DijFE * FiMinusDragX(IJK,L)
337                   massMobilityTermNoDragY = massMobilityTermNoDragY + DijFN * FiMinusDragY(IJK,L)
338                   massMobilityTermNoDragZ = massMobilityTermNoDragZ + DijFT * FiMinusDragZ(IJK,L)
339     
340                 ENDDO   ! end do (l=1,smax)
341     
342     
343                 thermalDiffTermX = DiTE * ( THETA_M(IJKE,MMAX) - THETA_M(IJK,MMAX) )  * oDX_E(I)
344                 thermalDiffTermY = DiTN * ( THETA_M(IJKN,MMAX) - THETA_M(IJK,MMAX) )  * oDY_N(J)
345                 thermalDiffTermZ = DiTT * ( THETA_M(IJKT,MMAX) - THETA_M(IJK,MMAX) )  * (oX_E(I)*oDZ_T(K))
346     
347                 JoiX(IJK,M) = -(ordinDiffTermX + thermalDiffTermX + massMobilityTermX)
348                 JoiY(IJK,M) = -(ordinDiffTermY + thermalDiffTermY + massMobilityTermY)
349                 JoiZ(IJK,M) = -(ordinDiffTermZ + thermalDiffTermZ + massMobilityTermZ)
350     
351                 JoiMinusDragX(IJK,M) = (ordinDiffTermX + thermalDiffTermX + massMobilityTermNoDragX)
352                 JoiMinusDragY(IJK,M) = (ordinDiffTermY + thermalDiffTermY + massMobilityTermNoDragY)
353                 JoiMinusDragZ(IJK,M) = (ordinDiffTermZ + thermalDiffTermZ + massMobilityTermNoDragZ)
354     
355                 ropsME=AVG_X(ROP_S(IJK,M),ROP_S(IJKE,M),I)
356                 ropsMN=AVG_Y(ROP_S(IJK,M),ROP_S(IJKN,M),J)
357                 ropsMT=AVG_Z(ROP_S(IJK,M),ROP_S(IJKT,M),K)
358     
359                 DELTAU(IJK,M) = -(ordinDiffTermX+thermalDiffTermX+massMobilityTermXvelupdate)
360                 DELTAV(IJK,M) = -(ordinDiffTermY+thermalDiffTermY+massMobilityTermYvelupdate)
361                 DELTAW(IJK,M) = -(ordinDiffTermZ+thermalDiffTermz+massMobilityTermZvelUpdate)
362     
363     ! set fluxes to zero in case very dilute conditions
364                 EPSA = AVG_X(ROP_S(IJK,M),ROP_S(IJKE,M),I) / RO_S(IJK,M)
365                 IF(EPSA <= ZERO_EP_S) JoiX(IJK,M) = ZERO
366     
367                 EPSA = AVG_Y(ROP_S(IJK,M),ROP_S(IJKN,M),J) / RO_S(IJK,M)
368                 IF(EPSA <= ZERO_EP_S) JoiY(IJK,M) = ZERO
369     
370                 EPSA = AVG_Z(ROP_S(IJK,M),ROP_S(IJKT,M),K) / RO_S(IJK,M)
371                 IF(EPSA <= ZERO_EP_S) JoiZ(IJK,M) = ZERO
372     
373     ! set flux components to zero in case of walls
374                 IF (IP_AT_E(IJK))  JoiX(IJK,M) = ZERO
375                 IF (IP_AT_N(IJK))  JoiY(IJK,M) = ZERO
376                 IF (IP_AT_T(IJK))  JoiZ(IJK,M) = ZERO
377     
378               ELSE
379                 JoiX(IJK,M) = ZERO
380                 JoiY(IJK,M) = ZERO
381                 JoiZ(IJK,M) = ZERO
382               ENDIF   ! end if/else (fluid_at(ijk))
383     
384      200  CONTINUE   ! outer IJK loop
385           ENDDO   ! for m = 1,smax
386     
387           RETURN
388           END SUBROUTINE GHDMASSFLUX
389     
390     
391     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
392     !                                                                      C
393     !  Subroutine: UpdateSpeciesVelocities                                 C
394     !  Purpose: Update solids velocities at celll faces based on the       C
395     !           formula Ui = U + Joi/(mi ni); also calculate averaged      C
396     !           velocities for dilute conditions.                          C
397     !                                                                      C
398     !  Author: S. Benyahia                              Date: 6-MAR-09     C
399     !  Reviewer:                                          Date:            C
400     !                                                                      C
401     !  Literature/Document References:                                     C
402     !     C. Hrenya handnotes and Garzo, Hrenya, Dufty papers (PRE, 2007)  C
403     !                                                                      C
404     !  Variables modified: solid species velocity components: Us, Vs, Ws   C
405     !                                                                      C
406     !                                                                      C
407     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
408     
409           SUBROUTINE UpdateSpeciesVelocities ()
410     
411     !-----------------------------------------------
412     ! Modules
413     !-----------------------------------------------
414           USE param
415           USE param1
416           USE geometry
417           USE compar
418           USE fldvar
419           USE indices
420           USE is
421           USE drag
422           USE visc_s
423           USE ghdtheory
424           USE physprop
425           USE run
426           USE constant
427           USE toleranc
428           USE fun_avg
429           USE functions
430           IMPLICIT NONE
431     !-----------------------------------------------
432     ! Local variables
433     !-----------------------------------------------
434     ! Index
435           INTEGER :: IJK, I, J, K
436     ! Index
437           INTEGER :: IJKE, IJKN, IJKT, IJKW, IJKS, IJKB, IMJK, IJMK, &
438                      IJKM
439     ! Solids phase
440           INTEGER :: S
441     ! species density at cell faces
442           integer :: kk, maxFluxS
443           double precision :: epgN, rogN, mugN, Vg
444           double precision :: Ur(smax), vrelSq(smax), vel, velup(smax)
445           double precision :: maxFlux
446           double precision :: rosN(smax), dp(smax)
447           double precision :: DijN(smax,smax), JoiM(smax), &
448                               DijN_H(smax,smax), DijN_A(smax,smax)
449           double precision :: beta_cell(smax), beta_ij_cell(smax,smax)
450     
451           integer :: ntrial
452           double precision tolx, tolf
453     !-----------------------------------------------
454     !     Function subroutines
455     !-----------------------------------------------
456     
457     ! for Newton method
458           ntrial = 100
459           tolx = 1d-14
460           tolf = 1d-14
461     
462           DO 200 IJK = ijkstart3, ijkend3
463             I = I_OF(IJK)
464             J = J_OF(IJK)
465             K = K_OF(IJK)
466     
467             IF ( FLUID_AT(IJK) ) THEN
468     
469               IJKE = EAST_OF(IJK)
470               IJKW = WEST_OF(IJK)
471               IJKN = NORTH_OF(IJK)
472               IJKS = SOUTH_OF(IJK)
473               IJKT = TOP_OF(IJK)
474               IJKB = BOTTOM_OF(IJK)
475               IMJK = IM_OF(IJK)
476               IJMK = JM_OF(IJK)
477               IJKM = KM_OF(IJK)
478     
479     ! First Compute Us
480     ! ---------------------------------------------------------------->>>
481               IF(.NOT.IP_AT_E(IJK) .OR. .NOT.SIP_AT_E(IJK)) THEN
482                 IF(RO_g0 == ZERO) THEN ! special case of a granular flow
483                   do s = 1, smax
484                     rosN(s) = AVG_X(ROP_S(IJK,s),ROP_S(IJKE,s),I)/ RO_S(IJK,s) ! this is ep_s
485                     if(rosN(s) > zero_ep_s) then
486                       u_s(ijk,s) = u_s(ijk,mmax) + JoiX(IJK,s)/(rosN(s)*ro_s(IJK,s))
487                     else
488                       u_s(ijk,s) = u_s(ijk,mmax) ! case of zero flux
489                     endif
490                   enddo
491                 ELSE
492                   Vg   = U_g(ijk) - u_s(ijk,mmax)
493                   epgN = AVG_X(EP_g(IJK),EP_g(IJKE),I)
494                   rogN = AVG_X(ROP_g(IJK),ROP_g(IJKE),I)
495                   mugN = AVG_X(MU_g(IJK),MU_g(IJKE),I)
496     
497                   do s = 1, smax
498                     Ur(s) = u_g(ijk)-u_s(ijk,s)
499     ! vrel must not include Ur, which is being solved for iteratively and must be updated.
500                     vrelSq(s) = (v_g(ijk)-v_s(ijk,s))**2 + (w_g(ijk)-w_s(ijk,s))**2
501                     rosN(s) = AVG_X(ROP_S(IJK,s),ROP_S(IJKE,s),I)
502                     velup(s) = 0.d0
503                     beta_cell(s) = beta_cell_X(IJK,s)
504                     dp(s)   = D_P(IJK,s)
505     
506                     IF(DRAG_TYPE_ENUM .eq. HYS)THEN
507                       JoiM(s) = DELTAU(IJK,s)
508                     ELSE
509                       JoiM(s) = JoiMinusDragX(ijk,s)
510                     ENDIF
511     
512                     do kk = 1, smax
513                       DijN_H(s,kk) = AVG_X_S(DijF(IJK,s,kk),DijF(IJKE,s,kk),I)
514                       DijN_A(s,kk) = AVG_X(DijF(IJK,s,kk),DijF(IJKE,s,kk),I)
515                       if(DijF_HarmE(IJK,kk))THEN
516                         DijN(s,kk) = DijN_H(s,kk)
517                       ELSE
518                         DijN(s,kk) = DijN_A(s,kk)
519                       ENDIF
520                       if(s .eq. kk)then
521                         beta_ij_cell(s,kk)=0.d0
522                       else
523                         beta_ij_cell(s,kk)=beta_ij_cell_X(IJK,s,kk)
524                       endif
525                     enddo
526                   enddo
527     
528                   IF(DRAG_TYPE_ENUM .eq. HYS)THEN
529                     vel=U_S(IJK,MMAX)
530                     CALL VELOCITY_UPDATE(velup, smax, rosN, DijN, &
531                       JoiM, beta_cell, beta_ij_cell,vel)
532                   ELSE
533                     CALL UrNEWT(ntrial, Ur, smax, ijk, tolx, tolf, &
534                       epgN, rogN, mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
535                   ENDIF
536     
537     ! species velocity and flux update.
538                   do s = 1, smax
539                     IF(DRAG_TYPE_ENUM .eq. HYS)THEN
540                       U_S(IJK,s)=velup(s)
541                       JoiX(IJK,s) = rosN(s) * (u_s(ijk,s)-u_s(ijk,mmax))
542                     ELSE
543                       u_s(ijk,s) =  u_g(ijk) - Ur(s)
544                       JoiX(IJK,s) = rosN(s) * (u_s(ijk,s)-u_s(ijk,mmax))
545                     ENDIF
546                   enddo
547                 ENDIF ! for a granular case (no gas and no drag)
548               ENDIF   ! end if (.not.ip_at_e or .not.sip_at_e)
549     
550               if(smax==2) then ! only for binary, how to implement for smax > 2?
551                 JoiX(IJK,2)=-JoiX(IJK,1)
552               elseif(smax > 2) then
553                 maxFlux = JoiX(IJK,1)
554                 maxFluxS = 1
555                 do s = 2, smax  ! finding species with maximum flux in a cell
556                   if( abs(JoiX(IJK,s)) > abs(maxFlux) ) then
557                     maxFlux = JoiX(IJK,s)
558                     maxFluxS = s
559                   endif
560                 enddo
561                 JoiX(IJK,maxFluxS) = 0d0 ! reset max. flux to zero
562                 do s = 1, smax ! re-calc species with max. flux to satisfy SUM(fluxes) = 0
563                   if(s /= maxFluxS) JoiX(IJK,maxFluxS) = JoiX(IJK,maxFluxS) - JoiX(IJK,s)
564                 enddo
565               endif
566     ! End Compute Us
567     ! ----------------------------------------------------------------<<<
568     
569     
570     
571     ! Now Compute Vs
572     ! ---------------------------------------------------------------->>>
573               IF (.NOT.IP_AT_N(IJK) .OR. .NOT.SIP_AT_N(IJK)) THEN
574                 IF(RO_g0 == ZERO) THEN ! special case of a granular flow
575                   do s = 1, smax
576                     rosN(s) = AVG_Y(ROP_S(IJK,s),ROP_S(IJKN,s),J)/ RO_S(IJK,s) ! this is ep_s
577                     if(rosN(s) > zero_ep_s) then
578                       v_s(ijk,s) = v_s(ijk,mmax) + JoiY(IJK,s)/(rosN(s)*ro_s(IJK,s))
579                     else
580                        v_s(ijk,s) = v_s(ijk,mmax) ! case of zero flux
581                     endif
582                   enddo
583                 ELSE
584                   Vg   = V_g(ijk) - v_s(ijk,mmax)
585                   epgN = AVG_Y(EP_g(IJK),EP_g(IJKN),J)
586                   rogN = AVG_Y(ROP_g(IJK),ROP_g(IJKN),J)
587                   mugN = AVG_Y(MU_g(IJK),MU_g(IJKN),J)
588     
589                   do s = 1, smax
590                     Ur(s) = v_g(ijk)-v_s(ijk,s)
591     ! vrel must not include Ur, which is being solved for iteratively and must be updated.
592                     vrelSq(s) = (u_g(ijk)-u_s(ijk,s))**2 + (w_g(ijk)-w_s(ijk,s))**2
593                     rosN(s) = AVG_Y(ROP_S(IJK,s),ROP_S(IJKN,s),J)
594                     velup(s) = 0.d0
595                     beta_cell(s) = beta_cell_Y(IJK,s)
596                     dp(s)   = D_P(IJK,s)
597     
598                     IF(DRAG_TYPE_ENUM .eq. HYS)THEN
599                       JoiM(s) = DELTAV(IJK,s)
600                     ELSE
601                       JoiM(s) = JoiMinusDragY(ijk,s)
602                     ENDIF
603     
604                     do kk = 1, smax
605     
606                       DijN_H(s,kk) = AVG_Y_S(DijF(IJK,s,kk),DijF(IJKN,s,kk),J)
607                       DijN_A(s,kk) = AVG_Y(DijF(IJK,s,kk),DijF(IJKN,s,kk),J)
608     
609                       if(DijF_HarmN(IJK,kk))THEN
610                         DijN(s,kk) = DijN_H(s,kk)
611                       ELSE
612                         DijN(s,kk) = DijN_A(s,kk)
613                       ENDIF
614     
615                       if(s .eq. kk)then
616                          beta_ij_cell(s,kk)=0.d0
617                       else
618                         beta_ij_cell(s,kk)=beta_ij_cell_Y(IJK,s,kk)
619                       endif
620                     enddo
621                   enddo
622     
623                   IF(DRAG_TYPE_ENUM .eq. HYS)THEN
624                     vel=V_S(IJK,MMAX)
625                     CALL VELOCITY_UPDATE(velup, smax, rosN, DijN, JoiM, &
626                       beta_cell, beta_ij_cell,vel)
627                   ELSE
628                     CALL UrNEWT(ntrial, Ur, smax, ijk, tolx, tolf, &
629                       epgN, rogN, mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
630                   ENDIF
631     
632     ! species velocity and flux update
633                   do s = 1, smax
634                     IF(DRAG_TYPE_ENUM .eq. HYS)THEN
635                       V_S(IJK,s)=velup(s)
636                       JoiY(IJK,s) = rosN(s) * (v_s(ijk,s)-v_s(ijk,mmax))
637                     ELSE
638                       v_s(ijk,s) =  v_g(ijk) - Ur(s)
639                       JoiY(IJK,s) = rosN(s) * (v_s(ijk,s)-v_s(ijk,mmax))
640                     ENDIF
641                   enddo
642                 ENDIF ! for a granular case (no gas and no drag)
643               ENDIF   ! end if (.not.ip_at_n or .not.sip_at_n)
644     
645               if(smax==2) then ! only for binary, how to implement for smax > 2?
646                 JoiY(IJK,2)=-JoiY(IJK,1)
647               elseif(smax > 2) then
648                 maxFlux = JoiY(IJK,1)
649                 maxFluxS = 1
650                 do s = 2, smax  ! finding species with maximum flux in a cell
651                   if( abs(JoiY(IJK,s)) > abs(maxFlux) ) then
652                     maxFlux = JoiY(IJK,s)
653                     maxFluxS = s
654                   endif
655                 enddo
656                 JoiY(IJK,maxFluxS) = 0d0 ! reset max. flux to zero
657                 do s = 1, smax ! re-calc species with max. flux to satisfy SUM(fluxes) = 0
658                   if(s /= maxFluxS) JoiY(IJK,maxFluxS) = JoiY(IJK,maxFluxS) - JoiY(IJK,s)
659                 enddo
660               endif
661     ! End Compute Vs
662     ! ----------------------------------------------------------------<<<
663     
664     
665     ! Finaly Compute Ws
666     ! ---------------------------------------------------------------->>>
667               IF(.NOT.NO_K .AND. (.NOT.IP_AT_T(IJK) .OR. .NOT.SIP_AT_T(IJK))) THEN
668                 IF(RO_g0 == ZERO) THEN ! special case of a granular flow
669                   do s = 1, smax
670                     rosN(s) = AVG_Z(ROP_S(IJK,s),ROP_S(IJKT,s),K)/ RO_S(IJK,s) ! this is ep_s
671                     if(rosN(s) > zero_ep_s) then
672                       w_s(ijk,s) = w_s(ijk,mmax) + JoiZ(IJK,s)/(rosN(s)*ro_s(IJK,s))
673                     else
674                       w_s(ijk,s) = w_s(ijk,mmax) ! case of zero flux
675                     endif
676                   enddo
677                 ELSE
678                   Vg   = W_g(ijk) - W_s(ijk,mmax)
679                   epgN = AVG_Z(EP_g(IJK),EP_g(IJKT),K)
680                   rogN = AVG_Z(ROP_g(IJK),ROP_g(IJKT),K)
681                   mugN = AVG_Z(MU_g(IJK),MU_g(IJKT),K)
682     
683                   do s = 1, smax
684                     Ur(s) = w_g(ijk)-w_s(ijk,s)
685     ! vrel must not include Ur, which is being solved for iteratively and must be updated.
686                     vrelSq(s) = (u_g(ijk)-u_s(ijk,s))**2 + (v_g(ijk)-v_s(ijk,s))**2
687                     rosN(s) = AVG_Z(ROP_S(IJK,s),ROP_S(IJKT,s),K)
688                     velup(s) = 0.d0
689                     beta_cell(s) = beta_cell_Z(IJK,s)
690                     dp(s)   = D_P(IJK,s)
691     
692                     IF(DRAG_TYPE_ENUM .eq. HYS)THEN
693                       JoiM(s) = DELTAW(IJK,s)
694                     ELSE
695                       JoiM(s) = JoiMinusDragZ(ijk,s)
696                     ENDIF
697     
698                     do kk = 1, smax
699                       DijN_H(s,kk) = AVG_Z_S(DijF(IJK,s,kk),DijF(IJKT,s,kk),K)
700                       DijN_A(s,kk) = AVG_Z(DijF(IJK,s,kk),DijF(IJKT,s,kk),K)
701                       if(DijF_HarmT(IJK,kk))THEN
702                         DijN(s,kk) = DijN_H(s,kk)
703                       ELSE
704                         DijN(s,kk) = DijN_A(s,kk)
705                       ENDIF
706                       if(s .eq. kk)then
707                          beta_ij_cell(s,kk)=0.d0
708                       else
709                         beta_ij_cell(s,kk)=beta_ij_cell_Z(IJK,s,kk)
710                       endif
711                     enddo
712                   enddo
713     
714                   IF(DRAG_TYPE_ENUM .eq. HYS)THEN
715                     vel=W_S(IJK,MMAX)
716                     CALL VELOCITY_UPDATE(velup, smax, rosN, DijN, &
717                       JoiM, beta_cell, beta_ij_cell,vel)
718                   ELSE
719                     CALL UrNEWT(ntrial, Ur, smax, ijk, tolx, tolf, &
720                       epgN, rogN, mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
721                   ENDIF
722     
723     ! species velocity and flux update
724                   do s = 1, smax
725                     IF(DRAG_TYPE_ENUM .eq. HYS)THEN
726                       W_S(IJK,s)=velup(s)
727                       JoiZ(IJK,s) = rosN(s) * (w_s(ijk,s)-w_s(ijk,mmax))
728                     ELSE
729                       w_s(ijk,s) =  w_g(ijk) - Ur(s)
730                       JoiZ(IJK,s) = rosN(s) * (w_s(ijk,s)-w_s(ijk,mmax))
731                     ENDIF
732                   enddo
733                 ENDIF ! for a granular case (no gas and no drag)
734               ENDIF   ! end if (.not.ip_at_t or .not.sip_at_t)
735     
736               if(smax==2) then ! only for binary, how to implement for smax > 2?
737                 JoiZ(IJK,2)=-JoiZ(IJK,1)
738               elseif(smax > 2) then
739                 maxFlux = JoiZ(IJK,1)
740                 maxFluxS = 1
741                 do s = 2, smax  ! finding species with maximum flux in a cell
742                   if( abs(JoiZ(IJK,s)) > abs(maxFlux) ) then
743                     maxFlux = JoiZ(IJK,s)
744                     maxFluxS = s
745                   endif
746                 enddo
747                 JoiZ(IJK,maxFluxS) = 0d0 ! reset max. flux to zero
748                 do s = 1, smax ! re-calc species with max. flux to satisfy SUM(fluxes) = 0
749                   if(s /= maxFluxS) JoiZ(IJK,maxFluxS) = JoiZ(IJK,maxFluxS) - JoiZ(IJK,s)
750                 enddo
751               endif
752     ! End Compute Ws
753     ! ----------------------------------------------------------------<<<
754     
755     ! if .not. fluid_at(ijk), do nothing (keep old values of velocities).
756     
757               ENDIF     ! Fluid_at
758      200  CONTINUE     ! outer IJK loop
759     
760           RETURN
761           END SUBROUTINE UpdateSpeciesVelocities
762     
763     
764     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
765     !                                                                      C
766     !  Subroutine: UrNEWT                                                  C
767     !                                                                      C
768     !                                                                      C
769     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
770     
771           subroutine UrNEWT(ntrial, x, s, ijk, tolx, tolf, epgN, rogN, &
772                             mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
773     
774           Implicit NONE
775     
776     !-----------------------------------------------
777     ! Dummy arguments
778     !-----------------------------------------------
779           integer, intent(in) :: s
780           integer, intent(in) :: ijk
781           integer, intent(in) :: ntrial
782           double precision, intent(in) :: tolx, tolf
783     
784           DOUBLE PRECISION, intent(inout) :: X(s)
785           double precision, intent(in) :: epgN, rogN, mugN, Vg
786           double precision, intent(in) :: vrelSq(s), rosN(s), dp(s)
787           double precision, intent(in) :: DijN(s,s), JoiM(s)
788     !-----------------------------------------------
789     ! Local parameters
790     !-----------------------------------------------
791           INTEGER :: NP
792           PARAMETER (NP=15)  ! solves up to NP variable/equations;
793     !-----------------------------------------------
794     ! Local variables
795     !-----------------------------------------------
796           INTEGER :: I, K, INDX(s)
797           DOUBLE PRECISION :: D, ERRF, ERRX, FJAC(s,s), FVEC(s), P(s)
798     !-----------------------------------------------
799     
800           DO K = 1, NTRIAL
801             CALL Ur_JACOBI_EVAL(X, s, ijk, FVEC, FJAC, epgN, rogN, mugN, Vg, &
802                                  vrelSq, rosN, dp, DijN, JoiM)
803             ERRF = 0d0
804             DO I = 1, s
805                  ERRF = ERRF + DABS(FVEC(I))
806             ENDDO
807     
808             IF(ERRF <= TOLF) RETURN
809             DO I = 1, s
810                 P(I) = -FVEC(I) ! RHS of linear equations.
811                 IF(rosN(I) .eq. 0.d0)THEN
812                   P(I) = 0.d0
813                 ENDIF
814             ENDDO
815     
816             CALL LUDCMP(fjac, s, NP, indx, d, 'UrNewt')  ! solve system of s linear equations using
817             CALL LUBKSB(fjac, s, NP, indx, p)  ! LU decomposition
818     
819             ERRX = 0d0
820             DO I = 1, s
821                 ERRX = ERRX + DABS(P(I))
822                 X(I) = X(I) + P(I)
823             ENDDO
824             IF(ERRX <= TOLX) RETURN
825           ENDDO  ! for K trials
826     
827           RETURN
828           END SUBROUTINE UrNEWT
829     
830     
831     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
832     !                                                                      C
833     !  Subroutine: Ur_JACOBI_EVAL                                          C
834     !                                                                      C
835     !                                                                      C
836     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
837     
838           SUBROUTINE Ur_JACOBI_EVAL(X, s, ijk, FVEC, FJAC, epgN, rogN, &
839                                     mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
840           Implicit NONE
841     
842     !-----------------------------------------------
843     ! Dummy arguments
844     !-----------------------------------------------
845           integer, intent(in) :: s
846           integer, intent(in) :: ijk
847           double precision, intent(in) :: X(s)
848     ! vector function and matrix jacobian
849           DOUBLE PRECISION, INTENT(OUT) :: FVEC(s), FJAC(s,s)
850           DOUBLE PRECISION, INTENT(IN) :: epgN, rogN, mugN, Vg
851           DOUBLE PRECISION, INTENT(IN) :: vrelSq(s), rosN(s), dp(s)
852           DOUBLE PRECISION, INTENT(IN) :: DijN(s,s), JoiM(s)
853     !-----------------------------------------------
854     ! Local parameters
855     !-----------------------------------------------
856           double precision :: pi
857           parameter (pi=3.14159265458979323846d0)
858           double precision :: one
859           parameter (one=1.d0)
860           double precision :: zero
861           parameter (zero=0.d0)
862     !-----------------------------------------------
863     ! Local variables
864     !-----------------------------------------------
865     ! Indices
866           INTEGER :: I, J
867     ! various quantities
868           DOUBLE PRECISION :: RE_G, C_d, DgA, Vi, vrel
869           DOUBLE PRECISION :: FgsOni(s), dFgsdVi(s), sum(s)
870     !-----------------------------------------------
871     
872           DO i = 1, s
873             vrel = dsqrt(vrelSq(i) + x(i)**2)
874             Vi = pi * dp(i)**3 / 6d0
875             RE_G = dp(i)*vrel*rogN/mugN
876             IF(Re_G <= 1000D0 .and. Re_G> zero)THEN
877                C_d = (24.D0/Re_G) * (ONE + 0.15D0 * Re_G**0.687D0)
878             ELSE
879                C_d = 0.44D0
880             ENDIF
881             DgA = 0.75D0 * C_d * vrel * rogN / (epgN**2.65D0 * dp(i))
882             FgsOni(i) = DgA * Vi  ! this is the wen-Yu drag
883             IF(vrel == ZERO)THEN
884                dFgsdVi(i) = ZERO
885             ELSEIF(Re_G <= 1000D0)THEN
886                dFgsdVi(i) = 1.8549d0*mugN*Re_G**0.687d0*Vi*dabs(x(i))/ &
887                                 (dp(i)**2*epgN**2.65d0*vrel**2)
888             ELSE
889                dFgsdVi(i) = 0.33d0*rogN*Vi*dabs(x(i))/ &
890                                  (dp(i)*epgN**2.65d0*vrel)
891             ENDIF
892           ENDDO
893     
894     ! Start computing values of the function FVEC(i)
895           DO i = 1, s
896              sum(i) = zero
897              DO j = 1, s
898                 sum(i) = sum(i) + DijN(i,j) * FgsOni(j) * x(j)
899              ENDDO
900           ENDDO
901           DO i = 1, s
902              FVEC(i) = sum(i) - rosN(i)*x(i) + rosN(i)*Vg + JoiM(i)
903           ENDDO
904     
905     !
906     ! Evaluation of functions FVEC(i) is done above, and their Jacobian is computed below.
907     !
908           DO i = 1, s
909              DO j = 1, s
910                 if(j == i) THEN
911                    FJAC(i,j) = DijN(i,i) * (FgsOni(i) + dFgsdVi(i)*x(i)) - rosN(i)
912                    IF(rosN(i) .eq. 0.d0)THEN
913                      FJAC(i,j) = 1.d0
914                    ENDIF
915                 else
916                    FJAC(i,j) = DijN(i,j) * (FgsOni(j) + dFgsdVi(j)*x(j))
917                 endif
918              ENDDO ! for j
919           ENDDO ! for i
920     !
921     ! End of computing jacobian (J_ij).
922     !
923           RETURN
924           END SUBROUTINE Ur_JACOBI_EVAL
925     
926     
927     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
928     !                                                                      C
929     !  Subroutine: Velocity_update                                         C
930     !                                                                      C
931     !                                                                      C
932     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
933     
934           SUBROUTINE VELOCITY_UPDATE(FVEC, s, rosi, Diji, Joii, &
935                                      beta_celli, beta_ij_celli, &
936                                      velocity)
937           USE toleranc
938           USE fldvar
939           Implicit NONE
940     
941     !-----------------------------------------------
942     ! Dummy arguments
943     !-----------------------------------------------
944           integer, intent (in) :: s
945     ! various quantities
946           DOUBLE PRECISION, INTENT(IN) :: rosi(s), Diji(s,s), &
947                                           Joii(s), beta_celli(s), &
948                                           beta_ij_celli(s,s), velocity
949     ! vector function
950           DOUBLE PRECISION, INTENT(OUT) :: FVEC(s)
951     !-----------------------------------------------
952     ! Local parameters
953     !-----------------------------------------------
954           integer :: NP
955           PARAMETER (NP=15)  ! solves up to NP variable/equations;
956     !-----------------------------------------------
957     ! Local variables
958     !-----------------------------------------------
959     ! Indices
960           INTEGER :: i, k, l
961           INTEGER :: indx(s)
962           DOUBLE PRECISION :: D
963     ! matrix jacobian
964           DOUBLE PRECISION :: FJAC(s,s)
965     !-----------------------------------------------
966     
967           DO i=1,s
968              FVEC(i)= 0.d0
969              DO k=1,s
970                FJAC(i,k) = 0.d0
971              ENDDO
972           ENDDO
973     
974           DO i=1,s
975     
976             IF(rosi(i) .ne. 0.d0)THEN
977                FVEC(i) = rosi(i)*velocity+Joii(i)
978             ELSE
979                FVEC(i) = 0.d0
980             ENDIF
981     
982             DO l=1,s
983     
984               if(i .eq. l)then
985                 FJAC(i,l)=FJAC(i,l)+rosi(i)
986                 IF(rosi(i) .eq. 0.d0)THEN   ! This is done to avoid a singular matrix if ros(i) is ZERO
987                   FJAC(i,l) = 1.d0
988                 ENDIF
989               endif
990     
991               FJAC(i,l)=FJAC(i,l)-Diji(i,l)*beta_celli(l)
992     
993               DO k=1,s
994     
995                if(l .ne. k)then
996                   FJAC(i,k)=FJAC(i,k)+Diji(i,l)*beta_ij_celli(l,k)
997                endif
998     
999               enddo
1000     
1001             enddo
1002           enddo
1003     
1004           CALL LUDCMP(fjac, s, NP, indx, d, 'velocity_update')  ! solve system of s linear equations using
1005           CALL LUBKSB(fjac, s, NP, indx, FVEC)  ! LU decomposition
1006     
1007           RETURN
1008     
1009           END SUBROUTINE VELOCITY_UPDATE
1010     
1011