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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_RESID_C                                            C
4     !  Purpose: Calculate residuals for continuity equations               C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
7     !  Reviewer:                                          Date:            C
8     !                                                                      C
9     !  Literature/Document References:                                     C
10     !  Variables referenced:                                               C
11     !  Variables modified:                                                 C
12     !  Local variables:                                                    C
13     !                                                                      C
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15     
16           SUBROUTINE CALC_RESID_C(VAR, A_M, B_M, M, NUM, DEN, &
17              RESID, MAX_RESID, IJK_RESID)
18     
19     !-----------------------------------------------
20     ! Modules
21     !-----------------------------------------------
22           USE param
23           USE param1
24           USE matrix
25           USE parallel
26           USE geometry
27           USE indices
28           USE compar
29           USE mpi_utility
30           USE run
31           USE functions
32           IMPLICIT NONE
33     !-----------------------------------------------
34     ! Dummy arguments
35     !-----------------------------------------------
36     ! Variable being evaluated
37           DOUBLE PRECISION, INTENT(IN) :: Var(DIMENSION_3)
38     ! Septadiagonal matrix A_m
39           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
40     ! Vector b_m
41           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3, 0:DIMENSION_M)
42     ! Phase index
43           INTEGER, INTENT(IN) :: M
44     ! Numerator and denominator
45           DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
46     ! Average value of Residual
47           DOUBLE PRECISION, INTENT(OUT) :: RESID
48     ! Maximum value of Residual
49           DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
50     ! IJK of Maximum value of Residual
51           INTEGER, INTENT(OUT) :: IJK_RESID
52     !-----------------------------------------------
53     ! Local variables
54     !-----------------------------------------------
55     ! Indices
56           INTEGER :: IJK, IJKW, IJKS, IJKB, IJKE, IJKN, IJKT
57     ! Numerators and denominators
58           DOUBLE PRECISION :: NUM1, DEN1
59     ! Number of fluid cells
60           INTEGER :: NCELLS
61     ! New local variables for DMP version
62           DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
63           DOUBLE PRECISION :: MAX_RESID_GL(0:numPEs-1), MAX_RESID_L(0:numPEs-1)
64           INTEGER :: IJK_RESID_GL(0:numPEs-1), IJK_RESID_L(0:numPEs-1)
65           INTEGER :: nproc
66     !-----------------------------------------------
67     
68     ! initializing values
69           NUM = ZERO
70           DEN = ZERO
71           MAX_RESID = -ONE
72           NCELLS = 0
73     
74     !!$omp parallel do private( IJK )
75           DO IJK = ijkstart3, ijkend3
76               RESID_IJK(IJK) = ZERO
77           ENDDO
78     
79     !!$omp  parallel do private( IJK, IJKW, IJKS, IJKB, IJKE, IJKN, IJKT,  &
80     !!$omp&  NUM1, DEN1) &
81     !!$omp&  REDUCTION(+:NUM,DEN,NCELLS)
82           DO IJK = ijkstart3, ijkend3
83              IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
84              IF (FLUID_AT(IJK)) THEN
85                 IJKW = WEST_OF(IJK)
86                 IJKS = SOUTH_OF(IJK)
87                 IJKE = EAST_OF(IJK)
88                 IJKN = NORTH_OF(IJK)
89     
90     ! evaluating the residual at cell ijk:
91     !   RESp = B-sum(Anb*VARnb)-Ap*VARp
92     !   (where nb = neighbor cells and p = center/0 cell)
93                 NUM1 = B_M(IJK,M) - (A_M(IJK,0,M)*VAR(IJK)+A_M(IJK,E,M)*VAR(IJKE)+&
94                    A_M(IJK,W,M)*VAR(IJKW)+A_M(IJK,N,M)*VAR(IJKN)+A_M(IJK,S,M)*VAR(&
95                    IJKS))
96     
97                 IF (DO_K) THEN
98                    IJKB = BOTTOM_OF(IJK)
99                    IJKT = TOP_OF(IJK)
100                    NUM1 = NUM1 - (A_M(IJK,T,M)*VAR(IJKT)+A_M(IJK,B,M)*VAR(IJKB))
101                 ENDIF
102     
103                 NUM1 = ABS(NUM1)
104                 DEN1 = ABS(A_M(IJK,0,M)*VAR(IJK))
105     ! storing value of residual at each ijk location
106                 RESID_IJK(IJK) = NUM1
107     
108     ! adding to terms that are accumulated
109                 NCELLS = NCELLS + 1
110                 NUM = NUM + NUM1
111                 DEN = DEN + DEN1
112              ENDIF
113           ENDDO
114     
115           IF(.not.debug_resid) RETURN
116     
117     ! Collecting all the information among all the procesors -
118     ! determining the global sum
119           call global_all_sum(NUM)
120           call global_all_sum(DEN)
121           call global_all_sum(NCELLS)
122     
123           IJK_RESID = 1
124           MAX_RESID = RESID_IJK( IJK_RESID )
125           DO IJK = ijkstart3, ijkend3
126              IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
127              IF (RESID_IJK(IJK) > MAX_RESID) THEN
128                 IJK_RESID = IJK
129                 MAX_RESID = RESID_IJK( IJK_RESID )
130              ENDIF
131           ENDDO
132     
133     ! Determining the max residual
134           do nproc=0,NumPEs-1
135              if(nproc.eq.myPE) then
136                 MAX_RESID_L(nproc) = MAX_RESID
137                 IJK_RESID_L(nproc) = FUNIJK_GL(I_OF(IJK_RESID), J_OF(IJK_RESID), K_OF(IJK_RESID))
138              else
139                 MAX_RESID_L(nproc) = 0.0
140                 IJK_RESID_L(nproc) = 0
141              endif
142           enddo
143     
144     ! Determining the maximum among all the procesors
145           call global_all_max(MAX_RESID)
146     
147     ! Collecting all the information among all the procesors
148           call global_all_sum(MAX_RESID_L, MAX_RESID_GL)
149           call global_all_sum(IJK_RESID_L, IJK_RESID_GL)
150     
151     ! Calling to determine the global IJK location w.r.t. serial version
152           IJK_RESID = IJKMAX2
153           do nproc=0,NumPEs-1
154              if(MAX_RESID_GL(nproc).eq.MAX_RESID.and.IJK_RESID_GL(nproc).lt.IJK_RESID) then
155                 IJK_RESID = IJK_RESID_GL(nproc)
156              endif
157           enddo
158     
159     ! Normalizing the residual
160           IF (DEN > ZERO) THEN
161              RESID = NUM/DEN
162              MAX_RESID = NCELLS*MAX_RESID/DEN
163           ELSEIF (NUM == ZERO) THEN
164              RESID = ZERO
165              MAX_RESID = ZERO
166              IJK_RESID = 0
167           ELSE
168              RESID = UNDEFINED
169              MAX_RESID = UNDEFINED
170     !         WRITE (LINE, *) 'Warning: All center coefficients are zero.'
171     !         CALL WRITE_ERROR ('CALC_RESID_C', LINE, 1)
172           ENDIF
173     
174           RETURN
175           END SUBROUTINE CALC_RESID_C
176     
177     
178     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
179     !                                                                      C
180     !  Subroutine: CALC_RESID_S                                            C
181     !  Purpose: Calculate residuals for scalar equations                   C
182     !                                                                      C
183     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
184     !  Reviewer:                                          Date:            C
185     !                                                                      C
186     !  Literature/Document References:                                     C
187     !  Variables referenced:                                               C
188     !  Variables modified:                                                 C
189     !  Local variables:                                                    C
190     !                                                                      C
191     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
192     
193           SUBROUTINE CALC_RESID_S(VAR, A_M, B_M, M, NUM, DEN, &
194              RESID, MAX_RESID, IJK_RESID, TOL)
195     
196     !-----------------------------------------------
197     ! Modules
198     !-----------------------------------------------
199           USE param
200           USE param1
201           USE matrix
202           USE parallel
203           USE geometry
204           USE indices
205           USE compar
206           USE mpi_utility
207           USE run
208     
209           USE fldvar
210           USE physprop
211           USE toleranc
212           USE functions
213     
214           IMPLICIT NONE
215     !-----------------------------------------------
216     ! Dummy arguments
217     !-----------------------------------------------
218     ! Variable being evaluated
219           DOUBLE PRECISION, INTENT(IN) :: Var(DIMENSION_3)
220     ! Septadiagonal matrix A_m
221           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
222     ! Vector b_m
223           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3, 0:DIMENSION_M)
224     ! Phase index
225           INTEGER, INTENT(IN) :: M
226     ! Numerator and denominator
227           DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
228     ! Average value of Residual
229           DOUBLE PRECISION, INTENT(OUT) :: RESID
230     ! Maximum value of Residual
231           DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
232     ! IJK of Maximum value of Residual
233           INTEGER, INTENT(OUT) :: IJK_RESID
234     ! Ignore residual calculation for scalar values below this
235           DOUBLE PRECISION, INTENT(IN) :: TOL
236     !-----------------------------------------------
237     ! Local variables
238     !-----------------------------------------------
239     ! Indices
240           INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
241     ! Numerators and denominators
242           DOUBLE PRECISION :: NUM1, DEN1
243     ! Number of fluid cells
244           INTEGER :: NCELLS
245     ! New local variables for DMP version
246           DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
247           DOUBLE PRECISION :: MAX_RESID_GL(0:numPEs-1), MAX_RESID_L(0:numPEs-1)
248           INTEGER :: IJK_RESID_GL(0:numPEs-1), IJK_RESID_L(0:numPEs-1)
249           INTEGER :: nproc
250     !-----------------------------------------------
251     
252     ! initializing
253           NUM = ZERO
254           DEN = ZERO
255           MAX_RESID = -ONE
256           NCELLS = 0
257     
258     !!$omp parallel do private( IJK )
259           DO IJK = ijkstart3, ijkend3
260              RESID_IJK(IJK) = ZERO
261           ENDDO
262     
263     !!$omp    parallel do &
264     !!$omp&   private(   IJK, IMJK, IJMK, IPJK, IJPK, IJKM, IJKP, &
265     !!$omp&   NUM1, DEN1) &
266     !!$omp&   REDUCTION(+:NUM, DEN,NCELLS)
267     
268           DO IJK = ijkstart3, ijkend3
269              IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
270     
271              IF (FLUID_AT(IJK) .AND. ABS(VAR(IJK)) > TOL) THEN
272                 IMJK = IM_OF(IJK)
273                 IJMK = JM_OF(IJK)
274                 IPJK = IP_OF(IJK)
275                 IJPK = JP_OF(IJK)
276     
277                 if(M/=0) then
278                    if(EP_S(IJK,M) <= DIL_EP_s) CYCLE
279                 endif
280     
281     
282     ! evaluating the residual at cell ijk:
283     !   RESp = B-sum(Anb*VARnb)-Ap*VARp
284     !   (where nb = neighbor cells and p = center/0 cell)
285                 NUM1 = B_M(IJK,M) - (A_M(IJK,0,M)*VAR(IJK)+A_M(IJK,E,M)*VAR(IPJK)+&
286                    A_M(IJK,W,M)*VAR(IMJK)+A_M(IJK,N,M)*VAR(IJPK)+A_M(IJK,S,M)*VAR(&
287                    IJMK))
288                 IF (DO_K) THEN
289                    IJKM = KM_OF(IJK)
290                    IJKP = KP_OF(IJK)
291                    NUM1 = NUM1 - (A_M(IJK,T,M)*VAR(IJKP)+A_M(IJK,B,M)*VAR(IJKM))
292                 ENDIF
293     
294                 NUM1 = ABS(NUM1)
295                 DEN1 = ABS(A_M(IJK,0,M)*VAR(IJK))
296     ! storing value of residual at each ijk location
297                 RESID_IJK(IJK) = NUM1
298     
299     ! adding to terms that are accumulated
300                 NCELLS = NCELLS + 1
301                 NUM = NUM + NUM1
302                 DEN = DEN + DEN1
303              ENDIF
304           ENDDO
305     
306           IF(.not.debug_resid) RETURN
307     
308     ! Collecting all the information among all the procesors -
309     ! determining the global sum
310           call global_all_sum(NUM)
311           call global_all_sum(DEN)
312           call global_all_sum(NCELLS)
313     
314           IJK_RESID = 1
315           MAX_RESID = RESID_IJK( IJK_RESID )
316           DO IJK = ijkstart3, ijkend3
317           IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
318              IF (RESID_IJK(IJK) > MAX_RESID) THEN
319                    IJK_RESID = IJK
320                    MAX_RESID = RESID_IJK( IJK_RESID )
321              ENDIF
322           ENDDO
323     
324     ! Determining the max residual
325           do nproc=0,NumPEs-1
326              if(nproc.eq.myPE) then
327                 MAX_RESID_L(nproc) = MAX_RESID
328                 IJK_RESID_L(nproc) = FUNIJK_GL(I_OF(IJK_RESID), J_OF(IJK_RESID), K_OF(IJK_RESID))
329              else
330                 MAX_RESID_L(nproc) = 0.0
331                 IJK_RESID_L(nproc) = 0
332              endif
333           enddo
334     
335     ! Determining the maximum among all the procesors
336           call global_all_max(MAX_RESID)
337     
338     ! Collecting all the information among all the procesors
339           call global_all_sum(MAX_RESID_L, MAX_RESID_GL)
340           call global_all_sum(IJK_RESID_L, IJK_RESID_GL)
341     
342     ! Determining the global IJK location w.r.t. serial version
343           IJK_RESID = IJKMAX2
344           do nproc=0,NumPEs-1
345             if(MAX_RESID_GL(nproc).eq.MAX_RESID.and.IJK_RESID_GL(nproc).lt.IJK_RESID) then
346                IJK_RESID = IJK_RESID_GL(nproc)
347             endif
348           enddo
349     
350     ! Normalizing the residual
351           IF (DEN > ZERO) THEN
352              RESID = NUM/DEN
353              MAX_RESID = NCELLS*MAX_RESID/DEN
354           ELSEIF (NUM == ZERO) THEN
355              RESID = ZERO
356              MAX_RESID = ZERO
357              IJK_RESID = 0
358           ELSE
359              RESID = UNDEFINED
360              MAX_RESID = UNDEFINED
361     !         WRITE(LINE,*)'Message: All center coefficients are zero.'
362     !         CALL WRITE_ERROR('CALC_RESID_S', LINE, 1)
363           ENDIF
364     
365           RETURN
366           END SUBROUTINE CALC_RESID_S
367     
368     
369     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
370     !                                                                      C
371     !  Subroutine: CALC_RESID_pp                                           C
372     !  Purpose: Calculate residuals for pressure correction equation       C
373     !                                                                      C
374     !  Author: M. Syamlal                                 Date: 21-JUN-96  C
375     !  Reviewer:                                          Date:            C
376     !                                                                      C
377     !  Comments:                                                           C
378     !  for correction equations the convergence for the corrections must   C
379     !  go to zero, therefore the vector b must go to zero. this value      C
380     !  cannot be normalized as the other equations are since the           C
381     !  denominator here will vanish.  thus the residual is normalized      C
382     !  based on its value in the first iteration                           C
383     !                                                                      C
384     !                                                                      C
385     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
386     
387           SUBROUTINE CALC_RESID_PP(B_M, NORM, NUM, DEN, RESID, MAX_RESID, &
388              IJK_RESID)
389     
390     !-----------------------------------------------
391     ! Modules
392     !-----------------------------------------------
393           USE param
394           USE param1
395           USE matrix
396           USE parallel
397           USE geometry
398           USE indices
399           USE compar
400           USE mpi_utility
401           USE run
402           USE functions
403           IMPLICIT NONE
404     !-----------------------------------------------
405     ! Dummy arguments
406     !-----------------------------------------------
407     ! Vector b_m
408           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3, 0:DIMENSION_M)
409     ! Normalization factor
410           DOUBLE PRECISION, INTENT(IN) :: NORM
411     ! Numerator and denominator
412           DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
413     ! Average value of Residual
414           DOUBLE PRECISION, INTENT(OUT) :: RESID
415     ! Maximum value of Residual
416           DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
417     ! IJK of Maximum value of Residual
418           INTEGER, INTENT(OUT) :: IJK_RESID
419     !-----------------------------------------------
420     ! Local variables
421     !-----------------------------------------------
422     ! Indices
423           INTEGER :: IJK
424     ! Number of fluid cells
425           INTEGER :: NCELLS
426     ! Numerators and denominators
427           DOUBLE PRECISION :: NUM1, DEN1
428     ! New local variables for DMP version
429           DOUBLE PRECISION :: MAX_RESID_GL(0:numPEs-1), MAX_RESID_L(0:numPEs-1)
430           INTEGER :: IJK_RESID_GL(0:numPEs-1), IJK_RESID_L(0:numPEs-1)
431           INTEGER :: nproc
432     !-----------------------------------------------
433     
434     ! initializing values
435           NUM = ZERO
436           DEN = ZERO
437           MAX_RESID = -ONE
438           NCELLS = 0
439           DEN1 = ONE
440           IJK_RESID = 1
441     
442           DO IJK = ijkstart3, ijkend3
443              IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
444              IF (FLUID_AT(IJK)) THEN
445     
446     ! evaluating the residual at cell ijk:
447                 NUM1 = ABS(B_M(IJK,0))
448                 IF (NUM1 > MAX_RESID) THEN
449                    MAX_RESID = NUM1
450                    IJK_RESID = IJK
451                 ENDIF
452     
453     ! adding to terms that are accumulated
454                 NCELLS = NCELLS + 1
455                 NUM = NUM + NUM1
456                 DEN = DEN + DEN1
457              ENDIF
458           ENDDO
459     
460     
461           IF(.not.debug_resid) THEN
462     ! Collecting all the information among all the procesors
463              call global_all_sum(NUM)
464              call global_all_sum(DEN)
465     
466     ! Normalizing the residual
467              IF (DEN*NORM > ZERO) THEN
468     ! if norm=1 then this simply becomes an unscaled 'average' residual
469                 RESID = NUM/(DEN*NORM)
470              ELSEIF (NUM == ZERO) THEN
471                 RESID = ZERO
472              ELSE
473                 RESID = LARGE_NUMBER
474              ENDIF
475           ELSE   ! if(debug_resid) branch
476     
477     ! Collecting all the information among all the procesors -
478     ! determining the global sum
479              call global_all_sum(NUM)
480              call global_all_sum(DEN)
481              call global_all_sum(NCELLS)
482     
483     ! Determining the max residual
484              do nproc=0,NumPEs-1
485                 if(nproc.eq.myPE) then
486                    MAX_RESID_L(nproc) = MAX_RESID
487                    IJK_RESID_L(nproc) = FUNIJK_GL(I_OF(IJK_RESID), J_OF(IJK_RESID), K_OF(IJK_RESID))
488                 else
489                    MAX_RESID_L(nproc) = 0.0
490                    IJK_RESID_L(nproc) = 0
491                 endif
492              enddo
493     
494     ! Determining the maximum among all the procesors
495              call global_all_max(MAX_RESID)
496     
497     ! Collecting all the information among all the procesors
498              call global_all_sum(MAX_RESID_L, MAX_RESID_GL)
499              call global_all_sum(IJK_RESID_L, IJK_RESID_GL)
500     
501     ! Determining the global IJK location w.r.t. serial version
502              IJK_RESID = IJKMAX2
503              do nproc=0,NumPEs-1
504                 if(MAX_RESID_GL(nproc).eq.MAX_RESID.and.IJK_RESID_GL(nproc).lt.IJK_RESID) then
505                    IJK_RESID = IJK_RESID_GL(nproc)
506                 endif
507              enddo
508     
509     ! Normalizing the residual
510              IF (DEN*NORM > ZERO) THEN
511                 RESID = NUM/(DEN*NORM)
512                 MAX_RESID = NCELLS*MAX_RESID/(DEN*NORM)
513              ELSEIF (NUM == ZERO) THEN
514                 RESID = ZERO
515                 MAX_RESID = ZERO
516                 IJK_RESID = 0
517              ELSE
518                 RESID = LARGE_NUMBER
519                 MAX_RESID = LARGE_NUMBER
520     !            WRITE (LINE, *) 'Warning: All center coefficients are zero.'
521     !             CALL WRITE_ERROR ('CALC_RESID_pp', LINE, 1)
522              ENDIF
523     
524           ENDIF   ! end if/else debug_resid branch
525     
526           RETURN
527           END SUBROUTINE CALC_RESID_PP
528     
529     
530     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
531     !                                                                      C
532     !  Subroutine: CALC_RESID_mb                                           C
533     !  Purpose: Calculate overall mass balance error                       C
534     !                                                                      C
535     !                                                                      C
536     !  Author: M. Syamlal                                 Date: 9-DEC-02   C
537     !  Reviewer:                                          Date:            C
538     !                                                                      C
539     !                                                                      C
540     !                                                                      C
541     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
542     
543           SUBROUTINE CALC_RESID_MB(INIT, ErrorPercent)
544     
545     !-----------------------------------------------
546     ! Modules
547     !-----------------------------------------------
548           USE param
549           USE param1
550           USE fldvar
551           USE parallel
552           USE geometry
553           USE indices
554           USE run
555           USE bc
556           USE constant
557           USE physprop
558           USE compar
559           USE mpi_utility
560           USE residual
561           USE rxns
562           USE mflux
563           USE functions
564           IMPLICIT NONE
565     !-----------------------------------------------
566     ! Dummy arguments
567     !-----------------------------------------------
568     ! Flag to check whether this is an initialization call
569     ! 0 -> initialize old accumulation; 1 -> calc residual
570           INTEGER, INTENT(IN) :: init
571     ! Total mass balance error as a % of inflow
572           DOUBLE PRECISION, INTENT(OUT) :: ErrorPercent(0:MMAX)
573     !-----------------------------------------------
574     ! Local variables
575     !-----------------------------------------------
576     ! phase index
577           INTEGER :: M
578     ! boundary index
579           INTEGER :: L
580     ! locally define dt, so that this routine works when dt is not defined
581           DOUBLE PRECISION :: dt_local
582     ! error index
583           INTEGER :: IER
584           DOUBLE PRECISION :: flux_in, flux_out, fin, fout
585           DOUBLE PRECISION :: err, accum_new, denom
586     !-----------------------------------------------
587     ! Functions
588     !-----------------------------------------------
589           DOUBLE PRECISION Accumulation
590     !-----------------------------------------------
591     
592           if(dt == UNDEFINED)then
593             dt_local = ONE
594           else
595             dt_local = dt
596           endif
597     
598           IF(init == 0) THEN
599     ! Initialize this routine
600     ! ---------------------------------------------------------------->>>
601     
602     ! Accumulation
603             if(dt == UNDEFINED)then
604               Accum_resid_g = ZERO
605             else
606               Accum_resid_g = Accumulation(ROP_g)
607             endif
608             DO M=1, MMAX
609               if(dt == UNDEFINED)then
610                 Accum_resid_s(M) = ZERO
611               else
612                 Accum_resid_s(M) = Accumulation(ROP_s(1,M))
613               endif
614             ENDDO
615             RETURN
616     ! end initialization
617     ! ----------------------------------------------------------------<<<
618     
619           ELSE
620     
621     ! Calculate residual
622     ! ---------------------------------------------------------------->>>
623             if(dt == UNDEFINED)then
624               Accum_new = - Accumulation(SUM_R_g) * dt_local
625             else
626               Accum_new = Accumulation(ROP_g) - Accumulation(SUM_R_g) * dt_local
627             endif
628     
629             flux_out = zero
630             flux_in = zero
631             DO L = 1, DIMENSION_BC
632               IF (BC_DEFINED(L)) THEN
633     !            call Calc_mass_flux(BC_I_W(L), BC_I_E(L), BC_J_S(L), &
634     !            BC_J_N(L), BC_K_B(L), BC_K_T(L), BC_PLANE(L), U_g, V_g, W_g, &
635     !            ROP_g, fin, fout, IER)
636                 IF(.NOT.Added_Mass) THEN
637                   call Calc_mass_fluxHR(BC_I_W(L), BC_I_E(L), BC_J_S(L), &
638                      BC_J_N(L), BC_K_B(L), BC_K_T(L), BC_PLANE(L), &
639                      Flux_gE, Flux_gN, Flux_gT, fin, fout, IER)
640                 ELSE
641                   call Calc_mass_fluxHR(BC_I_W(L), BC_I_E(L), BC_J_S(L), &
642                      BC_J_N(L), BC_K_B(L), BC_K_T(L), BC_PLANE(L), Flux_gSE,&
643                      Flux_gSN, Flux_gST, fin, fout, IER)
644                 ENDIF
645                 flux_out = flux_out + fout  * dt_local
646                 flux_in = flux_in + fin * dt_local
647               ENDIF
648             END DO
649     
650             Err = (accum_new - Accum_resid_g) - (flux_in - flux_out)
651             denom = max(abs(accum_new), abs(Accum_resid_g), abs(flux_in), abs(flux_out))
652             IF (denom /= ZERO) THEN
653                ErrorPercent(0) = err*100./denom
654             ELSE
655                ErrorPercent(0) = err*100./SMALL_NUMBER
656             ENDIF
657     
658             DO M =1, MMAX
659               if(dt == UNDEFINED)then
660                 Accum_new =  - Accumulation(SUM_R_s(1,M)) * dt_local
661               else
662                 Accum_new = Accumulation(ROP_s(1,M)) - Accumulation(SUM_R_s(1,M)) * dt_local
663               endif
664     
665               flux_out = zero
666               flux_in = zero
667               DO L = 1, DIMENSION_BC
668                 IF (BC_DEFINED(L)) THEN
669     !              call Calc_mass_flux(BC_I_W(L), BC_I_E(L), BC_J_S(L), BC_J_N(L), &
670     !              BC_K_B(L), BC_K_T(L), BC_PLANE(L), U_s(1,M), V_s(1,M), W_s(1,M), &
671     !              ROP_s(1,M), fin, fout, IER)
672                   IF(.NOT.Added_Mass .OR. M /= M_AM) THEN
673                     call Calc_mass_fluxHR(BC_I_W(L), BC_I_E(L), BC_J_S(L), &
674                        BC_J_N(L), BC_K_B(L), BC_K_T(L), BC_PLANE(L), &
675                        Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M), fin, fout, IER)
676                   ELSE
677                     call Calc_mass_fluxHR(BC_I_W(L), BC_I_E(L), BC_J_S(L), &
678                        BC_J_N(L), BC_K_B(L), BC_K_T(L), BC_PLANE(L), &
679                        Flux_sSE, Flux_sSN, Flux_sST, fin, fout, IER)
680                   ENDIF
681                   flux_out = flux_out + fout  * dt_local
682                   flux_in = flux_in + fin * dt_local
683                 ENDIF
684               ENDDO
685     
686               Err = (accum_new - Accum_resid_s(M)) - (flux_in - flux_out)
687               denom = max(abs(accum_new), abs(Accum_resid_s(M)), abs(flux_in), abs(flux_out))
688               if(denom /= ZERO) THEN
689                 ErrorPercent(M) = err*100./denom
690               else
691                 ErrorPercent(M) = err*100./SMALL_NUMBER
692               endif
693             ENDDO  ! end do m=1,mmax
694     ! end calculate residual
695     ! ----------------------------------------------------------------<<<
696     
697           ENDIF   ! end if/else (init==0)
698     
699           RETURN
700           END SUBROUTINE CALC_RESID_MB
701     
702     
703     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
704     !                                                                      C
705     !  Subroutine: CALC_RESID_U                                            C
706     !  Purpose: Calculate residuals for u-momentum equations               C
707     !                                                                      C
708     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
709     !  Reviewer:                                          Date:            C
710     !                                                                      C
711     !                                                                      C
712     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
713     
714           SUBROUTINE CALC_RESID_U(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, &
715              RESID, MAX_RESID, IJK_RESID)
716     
717     !-----------------------------------------------
718     ! Modules
719     !-----------------------------------------------
720           USE param
721           USE param1
722           USE matrix
723           USE parallel
724           USE geometry
725           USE indices
726           USE compar
727           USE mpi_utility
728           USE run
729           USE fldvar
730           USE physprop
731           USE toleranc
732           USE fun_avg
733           USE functions
734     
735           IMPLICIT NONE
736     !-----------------------------------------------
737     ! Dummy arguments
738     !-----------------------------------------------
739     ! U velocity (x-dir)
740           DOUBLE PRECISION, INTENT(IN) :: U_m(DIMENSION_3)
741     ! V velocity (y-dir), used here for scaling
742           DOUBLE PRECISION, INTENT(IN) :: V_m(DIMENSION_3)
743     ! W velocity (z-dir), used here for scaling
744           DOUBLE PRECISION, INTENT(IN) :: W_m(DIMENSION_3)
745     ! Septadiagonal matrix A_m
746           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
747     ! Vector b_m
748           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3, 0:DIMENSION_M)
749     ! Phase index
750           INTEGER, INTENT(IN) :: M
751     ! Numerator and denominator
752           DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
753     ! Average value of Residual
754           DOUBLE PRECISION, INTENT(OUT) :: RESID
755     ! Maximum value of Residual
756           DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
757     ! IJK of Maximum value of Residual
758           INTEGER, INTENT(OUT) :: IJK_RESID
759     !-----------------------------------------------
760     !     Local variables
761     !-----------------------------------------------
762     ! Velocity magnitude
763           DOUBLE PRECISION :: VEL
764     ! Indices
765           INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
766     ! Numerators and denominators
767           DOUBLE PRECISION :: NUM1, DEN1
768     ! Number of fluid cells
769           INTEGER :: NCELLS
770     ! New local variables for DMP version
771           DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
772           DOUBLE PRECISION :: MAX_RESID_GL(0:numPEs-1), MAX_RESID_L(0:numPEs-1)
773           INTEGER :: IJK_RESID_GL(0:numPEs-1), IJK_RESID_L(0:numPEs-1)
774           INTEGER :: nproc
775     
776     ! Solids volume fraction at face
777           DOUBLE PRECISION :: EPSA
778     
779     !-----------------------------------------------
780     
781     ! initializing
782           NUM = ZERO
783           DEN = ZERO
784           MAX_RESID = -ONE
785           NCELLS = 0
786     
787     !$omp  parallel default(none) &
788     !$omp  private( IMJK, IPJK, IJMK, IJPK, IJKM, IJKP, &
789     !$omp           VEL,  NUM1, DEN1,EPSA) &
790     !$omp  shared (ijkstart3, ijkend3, resid_ijk, i_of, j_of, k_of,m,a_m,b_m,w_m,do_k,u_m,v_m,num,den,ncells)
791     !$omp do reduction(+:num, DEN, NCELLS )
792           DO IJK = ijkstart3, ijkend3
793             RESID_IJK(IJK) = ZERO
794     
795             IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
796     
797     ! Skip walls where some values are undefined.
798             IF(WALL_AT(IJK)) cycle
799     
800              if(m/=0) then
801                 EPSA = AVG_X(EP_S(IJK,M),EP_S(EAST_OF(IJK),M),I_OF(IJK))
802                 if(EPSA <= DIL_EP_s) CYCLE
803              endif
804     
805              IF (.NOT.IP_AT_E(IJK)) THEN
806                 IMJK = IM_OF(IJK)
807                 IJMK = JM_OF(IJK)
808                 IPJK = IP_OF(IJK)
809                 IJPK = JP_OF(IJK)
810     
811     ! evaluating the residual at cell ijk:
812     !   RESp = B-sum(Anb*VARnb)-Ap*VARp
813     !   (where nb = neighbor cells and p = center/0 cell)
814                 NUM1 = B_M(IJK,M) - (A_M(IJK,0,M)*U_M(IJK)+&
815                    A_M(IJK,E,M)*U_M(IPJK)+A_M(IJK,W,M)*U_M(IMJK)+&
816                    A_M(IJK,N,M)*U_M(IJPK)+A_M(IJK,S,M)*U_M(IJMK))
817                 IF (DO_K) THEN
818                    IJKM = KM_OF(IJK)
819                    IJKP = KP_OF(IJK)
820                    NUM1 = NUM1 - (A_M(IJK,T,M)*U_M(IJKP)+A_M(IJK,B,M)*U_M(IJKM))
821                 ENDIF
822     
823     ! Ignore momentum residual in stagnant regions.  Need an alternative
824     ! criteria for residual scaling for such cases.
825                 VEL = SQRT(U_M(IJK)**2+V_M(IJK)**2+W_M(IJK)**2)
826                 IF (VEL > SMALL_NUMBER) THEN
827                    NUM1 = ABS(NUM1)
828                    DEN1 = ABS(A_M(IJK,0,M)*VEL)
829     ! storing value of residual at each ijk location
830                    RESID_IJK(IJK) = NUM1
831     ! adding to terms that are accumulated
832                    NCELLS = NCELLS + 1
833                    NUM = NUM + NUM1
834                    DEN = DEN + DEN1
835                 ENDIF
836              ENDIF
837           ENDDO
838     !$omp end parallel
839     
840           IF(.not.debug_resid) RETURN
841     
842     ! Collecting all the information among all the procesors -
843     ! determining the global sum
844           call global_all_sum(NUM)
845           call global_all_sum(DEN)
846           call global_all_sum(NCELLS)
847     
848           IJK_RESID = 1
849           MAX_RESID = RESID_IJK( IJK_RESID )
850           DO IJK = ijkstart3, ijkend3
851              IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
852              IF (RESID_IJK( IJK ) > MAX_RESID) THEN
853                 IJK_RESID = IJK
854                 MAX_RESID = RESID_IJK( IJK_RESID )
855              ENDIF
856           ENDDO
857     
858     ! Determining the max residual
859           do nproc=0,NumPEs-1
860              if(nproc.eq.myPE) then
861                 MAX_RESID_L(nproc) = MAX_RESID
862                 IJK_RESID_L(nproc) = FUNIJK_GL(I_OF(IJK_RESID), J_OF(IJK_RESID), K_OF(IJK_RESID))
863              else
864                 MAX_RESID_L(nproc) = 0.0
865                 IJK_RESID_L(nproc) = 0
866              endif
867           enddo
868     
869     ! Determining the maximum among all the procesors
870           call global_all_max(MAX_RESID)
871     
872     ! Collecting all the information among all the procesors
873           call global_all_sum(MAX_RESID_L, MAX_RESID_GL)
874           call global_all_sum(IJK_RESID_L, IJK_RESID_GL)
875     
876     ! Determining the global IJK location w.r.t. serial version
877           IJK_RESID = IJKMAX2
878           do nproc=0,NumPEs-1
879             IF(MAX_RESID_GL(nproc).eq.MAX_RESID.and.&
880                IJK_RESID_GL(nproc).lt.IJK_RESID) THEN
881               IJK_RESID = IJK_RESID_GL(nproc)
882             ENDIF
883           ENDDO
884     
885     ! Normalizing the residual
886           IF (DEN > ZERO) THEN
887              RESID = NUM/DEN
888              MAX_RESID = NCELLS*MAX_RESID/DEN
889           ELSEIF (NUM == ZERO) THEN
890              RESID = ZERO
891              MAX_RESID = ZERO
892              IJK_RESID = 0
893           ELSE
894              RESID = LARGE_NUMBER
895              MAX_RESID = LARGE_NUMBER
896     !         WRITE (LINE, *) 'Warning: All center coefficients are zero.'
897     !         CALL WRITE_ERROR ('CALC_RESID_U', LINE, 1)
898           ENDIF
899     
900           RETURN
901           END SUBROUTINE CALC_RESID_U
902     
903     
904     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
905     !                                                                      C
906     !  Subroutine: CALC_RESID_V                                            C
907     !  Purpose: Calculate residuals for v-momentum equations               C
908     !                                                                      C
909     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
910     !  Reviewer:                                          Date:            C
911     !                                                                      C
912     !                                                                      C
913     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
914     
915           SUBROUTINE CALC_RESID_V(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, &
916              RESID, MAX_RESID, IJK_RESID)
917     
918     !-----------------------------------------------
919     ! Modules
920     !-----------------------------------------------
921           USE param
922           USE param1
923           USE matrix
924           USE parallel
925           USE geometry
926           USE indices
927           USE compar
928           USE mpi_utility
929           USE run
930           USE fldvar
931           USE physprop
932           USE toleranc
933           USE fun_avg
934           USE functions
935     
936           IMPLICIT NONE
937     !-----------------------------------------------
938     ! Dummy arguments
939     !-----------------------------------------------
940     ! U velocity (x-dir)
941           DOUBLE PRECISION, INTENT(IN) :: U_m(DIMENSION_3)
942     ! V velocity (y-dir), used here for scaling
943           DOUBLE PRECISION, INTENT(IN) :: V_m(DIMENSION_3)
944     ! W velocity (z-dir), used here for scaling
945           DOUBLE PRECISION, INTENT(IN) :: W_m(DIMENSION_3)
946     ! Septadiagonal matrix A_m
947           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
948     ! Vector b_m
949           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3, 0:DIMENSION_M)
950     ! Phase index
951           INTEGER, INTENT(IN) :: M
952     ! Numerator and denominator
953           DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
954     ! Average value of Residual
955           DOUBLE PRECISION, INTENT(OUT) :: RESID
956     ! Maximum value of Residual
957           DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
958     ! IJK of Maximum value of Residual
959           INTEGER, INTENT(OUT) :: IJK_RESID
960     !-----------------------------------------------
961     !     Local variables
962     !-----------------------------------------------
963     ! Velocity magnitude
964           DOUBLE PRECISION :: VEL
965     ! Indices
966           INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
967     ! Numerators and denominators
968           DOUBLE PRECISION :: NUM1, DEN1
969     ! Number of fluid cells
970           INTEGER :: NCELLS
971     ! New local variables for DMP version
972           DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
973           DOUBLE PRECISION :: MAX_RESID_GL(0:numPEs-1), MAX_RESID_L(0:numPEs-1)
974           INTEGER :: IJK_RESID_GL(0:numPEs-1), IJK_RESID_L(0:numPEs-1)
975           INTEGER :: nproc
976     ! Solids volume fraction at face
977           DOUBLE PRECISION :: EPSA
978     
979     !-----------------------------------------------
980     
981     ! initializing
982           NUM = ZERO
983           DEN = ZERO
984           MAX_RESID = -ONE
985           NCELLS = 0
986     
987     !$omp  parallel default(none) &
988     !$omp  private( IMJK, IPJK, IJMK, IJPK, IJKM, IJKP, &
989     !$omp           VEL,  NUM1, DEN1,EPSA) &
990     !$omp  shared (ijkstart3, ijkend3, resid_ijk, i_of, j_of, k_of,m,a_m,b_m,w_m,do_k,u_m,v_m,num,den,ncells)
991     !$omp do reduction(+:num, DEN, NCELLS )
992           DO IJK = ijkstart3, ijkend3
993             RESID_IJK(IJK) = ZERO
994     
995             IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
996     
997     ! Skip walls where some values are undefined.
998             IF(WALL_AT(IJK)) cycle
999     
1000     
1001              if(m/=0) then
1002                 EPSA = AVG_Y(EP_S(IJK,M),EP_S(NORTH_OF(IJK),M),J_OF(IJK))
1003                 if(EPSA <= DIL_EP_s) CYCLE
1004              endif
1005     
1006              IF (.NOT.IP_AT_N(IJK)) THEN
1007                 IMJK = IM_OF(IJK)
1008                 IJMK = JM_OF(IJK)
1009                 IPJK = IP_OF(IJK)
1010                 IJPK = JP_OF(IJK)
1011     
1012     ! evaluating the residual at cell ijk:
1013     !   RESp = B-sum(Anb*VARnb)-Ap*VARp
1014     !   (where nb = neighbor cells and p = center/0 cell)
1015                 NUM1 = B_M(IJK,M) - (A_M(IJK,0,M)*V_M(IJK)+&
1016                    A_M(IJK,E,M)*V_M(IPJK)+A_M(IJK,W,M)*V_M(IMJK)+&
1017                    A_M(IJK,N,M)*V_M(IJPK)+A_M(IJK,S,M)*V_M(IJMK))
1018                 IF (DO_K) THEN
1019                    IJKM = KM_OF(IJK)
1020                    IJKP = KP_OF(IJK)
1021                    NUM1 = NUM1 - (A_M(IJK,T,M)*V_M(IJKP)+A_M(IJK,B,M)*V_M(IJKM))
1022                 ENDIF
1023     
1024     ! Ignore momentum residual in stagnant regions.  Need an alternative
1025     ! criteria for residual scaling for such cases.
1026                 VEL = SQRT(U_M(IJK)**2+V_M(IJK)**2+W_M(IJK)**2)
1027                 IF (VEL > SMALL_NUMBER) THEN
1028                    NUM1 = ABS(NUM1)
1029                    DEN1 = ABS(A_M(IJK,0,M)*VEL)
1030     ! storing value of residual at each ijk location
1031                    RESID_IJK(IJK) = NUM1
1032     ! adding to terms that are accumulated
1033                    NCELLS = NCELLS + 1
1034                    NUM = NUM + NUM1
1035                    DEN = DEN + DEN1
1036                 ENDIF
1037              ENDIF
1038           ENDDO
1039     !$omp end parallel
1040     
1041           if(.not.debug_resid) return
1042     
1043     ! Collecting all the information among all the procesors -
1044     ! determining the global sum
1045           call global_all_sum(NUM)
1046           call global_all_sum(DEN)
1047           call global_all_sum(NCELLS)
1048     
1049           IJK_RESID = 1
1050           MAX_RESID = RESID_IJK( IJK_RESID )
1051           DO IJK = ijkstart3, ijkend3
1052           IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
1053               IF (RESID_IJK( IJK ) > MAX_RESID) THEN
1054                   IJK_RESID = IJK
1055                   MAX_RESID = RESID_IJK( IJK_RESID )
1056               ENDIF
1057           ENDDO
1058     
1059     ! Determining the max residual
1060           do nproc=0,NumPEs-1
1061              if(nproc.eq.myPE) then
1062                 MAX_RESID_L(nproc) = MAX_RESID
1063                 IJK_RESID_L(nproc) = FUNIJK_GL(I_OF(IJK_RESID), J_OF(IJK_RESID), K_OF(IJK_RESID))
1064              else
1065                 MAX_RESID_L(nproc) = 0.0
1066                 IJK_RESID_L(nproc) = 0
1067              endif
1068           enddo
1069     
1070     ! Determine the maximum among all the procesors
1071           call global_all_max(MAX_RESID)
1072     
1073     ! Collecting all the information among all the procesors
1074           call global_all_sum(MAX_RESID_L, MAX_RESID_GL)
1075           call global_all_sum(IJK_RESID_L, IJK_RESID_GL)
1076     
1077     ! Determining the global IJK location w.r.t. serial version
1078           IJK_RESID = IJKMAX2
1079           do nproc=0,NumPEs-1
1080             if(MAX_RESID_GL(nproc).eq.MAX_RESID.and.IJK_RESID_GL(nproc).lt.IJK_RESID) then
1081               IJK_RESID = IJK_RESID_GL(nproc)
1082             endif
1083           enddo
1084     
1085     ! Normalizing the residual
1086           IF (DEN > ZERO) THEN
1087              RESID = NUM/DEN
1088              MAX_RESID = NCELLS*MAX_RESID/DEN
1089           ELSEIF (NUM == ZERO) THEN
1090              RESID = ZERO
1091              MAX_RESID = ZERO
1092              IJK_RESID = 0
1093           ELSE
1094              RESID = LARGE_NUMBER
1095              MAX_RESID = LARGE_NUMBER
1096     !         WRITE (LINE, *) 'Warning: All center coefficients are zero.'
1097     !         CALL WRITE_ERROR ('CALC_RESID_V', LINE, 1)
1098           ENDIF
1099     
1100           RETURN
1101           END SUBROUTINE CALC_RESID_V
1102     
1103     
1104     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1105     !                                                                      C
1106     !  Subroutine: CALC_RESID_W                                            C
1107     !  Purpose: Calculate residuals for w-momentum equations               C
1108     !                                                                      C
1109     !  Author: M. Syamlal                                 Date: 21-MAY-96  C
1110     !  Reviewer:                                          Date:            C
1111     !                                                                      C
1112     !  Literature/Document References:                                     C
1113     !  Variables referenced:                                               C
1114     !  Variables modified:                                                 C
1115     !  Local variables:                                                    C
1116     !                                                                      C
1117     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1118     
1119           SUBROUTINE CALC_RESID_W(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, &
1120              RESID, MAX_RESID, IJK_RESID)
1121     
1122     !-----------------------------------------------
1123     ! Modules
1124     !-----------------------------------------------
1125           USE param
1126           USE param1
1127           USE matrix
1128           USE parallel
1129           USE geometry
1130           USE indices
1131           USE compar
1132           USE mpi_utility
1133           USE run
1134           USE fldvar
1135           USE physprop
1136           USE toleranc
1137           USE fun_avg
1138           USE functions
1139     
1140           IMPLICIT NONE
1141     !-----------------------------------------------
1142     ! Dummy arguments
1143     !-----------------------------------------------
1144     ! U velocity (x-dir)
1145           DOUBLE PRECISION, INTENT(IN) :: U_m(DIMENSION_3)
1146     ! V velocity (y-dir), used here for scaling
1147           DOUBLE PRECISION, INTENT(IN) :: V_m(DIMENSION_3)
1148     ! W velocity (z-dir), used here for scaling
1149           DOUBLE PRECISION, INTENT(IN) :: W_m(DIMENSION_3)
1150     ! Septadiagonal matrix A_m
1151           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1152     ! Vector b_m
1153           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3, 0:DIMENSION_M)
1154     ! Phase index
1155           INTEGER, INTENT(IN) :: M
1156     ! Numerator and denominator
1157           DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
1158     ! Average value of Residual
1159           DOUBLE PRECISION, INTENT(OUT) :: RESID
1160     ! Maximum value of Residual
1161           DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
1162     ! IJK of Maximum value of Residual
1163           INTEGER, INTENT(OUT) :: IJK_RESID
1164     !-----------------------------------------------
1165     !     Local variables
1166     !-----------------------------------------------
1167     ! Velocity magnitude
1168           DOUBLE PRECISION :: VEL
1169     ! Indices
1170           INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
1171     ! Numerators and denominators
1172           DOUBLE PRECISION :: NUM1, DEN1
1173     ! Number of fluid cells
1174           INTEGER :: NCELLS
1175     ! New local variables for DMP version
1176           DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
1177           DOUBLE PRECISION :: MAX_RESID_GL(0:numPEs-1), MAX_RESID_L(0:numPEs-1)
1178           INTEGER :: IJK_RESID_GL(0:numPEs-1), IJK_RESID_L(0:numPEs-1)
1179           INTEGER :: nproc
1180     ! Solids volume fraction at face
1181           DOUBLE PRECISION :: EPSA
1182     
1183     !-----------------------------------------------
1184     
1185     ! initializing
1186           NUM = ZERO
1187           DEN = ZERO
1188           MAX_RESID = -ONE
1189           NCELLS = 0
1190     
1191     !$omp  parallel default(none) &
1192     !$omp  private( IMJK, IPJK, IJMK, IJPK, IJKM, IJKP, &
1193     !$omp           VEL,  NUM1, DEN1,EPSA) &
1194     !$omp  shared (ijkstart3, ijkend3, resid_ijk, i_of, j_of, k_of,m,a_m,b_m,w_m,do_k,u_m,v_m,num,den,ncells)
1195     !$omp do reduction(+:num, DEN, NCELLS )
1196           DO IJK = ijkstart3, ijkend3
1197             RESID_IJK(IJK) = ZERO
1198     
1199             IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
1200     
1201     ! Skip walls where some values are undefined.
1202             IF(WALL_AT(IJK)) cycle
1203     
1204              if(m/=0) then
1205                 EPSA = AVG_Z(EP_S(IJK,M),EP_S(TOP_OF(IJK),M),K_OF(IJK))
1206                 if(EPSA <= DIL_EP_s) CYCLE
1207              endif
1208     
1209              IF (.NOT.IP_AT_T(IJK)) THEN
1210                 IMJK = IM_OF(IJK)
1211                 IJMK = JM_OF(IJK)
1212                 IPJK = IP_OF(IJK)
1213                 IJPK = JP_OF(IJK)
1214     
1215     ! evaluating the residual at cell ijk:
1216     !   RESp = B-sum(Anb*VARnb)-Ap*VARp
1217     !   (where nb = neighbor cells and p = center/0 cell)
1218                 NUM1 = B_M(IJK,M) - (A_M(IJK,0,M)*W_M(IJK)+&
1219                    A_M(IJK,E,M)*W_M(IPJK)+A_M(IJK,W,M)*W_M(IMJK)+&
1220                    A_M(IJK,N,M)*W_M(IJPK)+A_M(IJK,S,M)*W_M(IJMK))
1221                 IF (DO_K) THEN
1222                    IJKM = KM_OF(IJK)
1223                    IJKP = KP_OF(IJK)
1224                    NUM1 = NUM1 - (A_M(IJK,T,M)*W_M(IJKP)+A_M(IJK,B,M)*W_M(IJKM))
1225                 ENDIF
1226     
1227     ! Ignore momentum residual in stagnant regions.  Need an alternative
1228     ! criteria for residual scaling for such cases.
1229                 VEL = SQRT(U_M(IJK)**2+V_M(IJK)**2+W_M(IJK)**2)
1230                 IF (VEL > SMALL_NUMBER) THEN
1231                    NUM1 = ABS(NUM1)
1232                    DEN1 = ABS(A_M(IJK,0,M)*VEL)
1233     ! storing value of residual at each ijk location
1234                    RESID_IJK(IJK) = NUM1
1235     ! adding to terms that are accumulated
1236                    NCELLS = NCELLS + 1
1237                    NUM = NUM + NUM1
1238                    DEN = DEN + DEN1
1239                 ENDIF
1240              ENDIF
1241           ENDDO
1242     !$omp end parallel
1243     
1244           if(.not.debug_resid) return
1245     
1246     ! Collecting all the information among all the procesors -
1247     ! determining the global sum
1248           call global_all_sum(NUM)
1249           call global_all_sum(DEN)
1250           call global_all_sum(NCELLS)
1251     
1252           IJK_RESID = 1
1253           MAX_RESID = RESID_IJK( IJK_RESID )
1254           DO IJK = ijkstart3, ijkend3
1255           IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK),J_OF(IJK), K_OF(IJK))) CYCLE
1256     
1257               IF (RESID_IJK( IJK ) > MAX_RESID) THEN
1258                   IJK_RESID = IJK
1259                   MAX_RESID = RESID_IJK( IJK_RESID )
1260               ENDIF
1261           ENDDO
1262     
1263     ! Determining the max residual
1264           do nproc=0,NumPEs-1
1265              if(nproc.eq.myPE) then
1266                 MAX_RESID_L(nproc) = MAX_RESID
1267                 IJK_RESID_L(nproc) = FUNIJK_GL(I_OF(IJK_RESID), J_OF(IJK_RESID), K_OF(IJK_RESID))
1268              else
1269                 MAX_RESID_L(nproc) = 0.0
1270                 IJK_RESID_L(nproc) = 0
1271              endif
1272           enddo
1273     
1274     ! Determining the maximum among all the procesors
1275           call global_all_max(MAX_RESID)
1276     
1277     ! Collecting all the information among all the procesors
1278           call global_all_sum(MAX_RESID_L, MAX_RESID_GL)
1279           call global_all_sum(IJK_RESID_L, IJK_RESID_GL)
1280     
1281     ! Determining the global IJK location w.r.t. serial version
1282           IJK_RESID = IJKMAX2
1283     
1284           do nproc=0,NumPEs-1
1285             if(MAX_RESID_GL(nproc).eq.MAX_RESID.and.IJK_RESID_GL(nproc).lt.IJK_RESID) then
1286               IJK_RESID = IJK_RESID_GL(nproc)
1287             endif
1288           enddo
1289     
1290     ! Normalizing the residual
1291           IF (DEN > ZERO) THEN
1292              RESID = NUM/DEN
1293              MAX_RESID = NCELLS*MAX_RESID/DEN
1294           ELSE IF (NUM == ZERO) THEN
1295              RESID = ZERO
1296              MAX_RESID = ZERO
1297              IJK_RESID = 0
1298           ELSE
1299              RESID = LARGE_NUMBER
1300              MAX_RESID = LARGE_NUMBER
1301     !         WRITE (LINE, *) 'Warning: All center coefficients are zero.'
1302     !         CALL WRITE_ERROR ('CALC_RESID_W', LINE, 1)
1303           ENDIF
1304     
1305           RETURN
1306           END SUBROUTINE CALC_RESID_W
1307     
1308     
1309