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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CHECK_Mass_balance(init)                               C
4     !  Purpose: Check global species and elemental balances                C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 27-Nov-02  C
7     !  Reviewer:                                          Date:            C
8     !                                                                      C
9     !  Literature/Document References:                                     C
10     !                                                                      C
11     !  Variables referenced:                                               C
12     !  Variables modified:                                                 C
13     !                                                                      C
14     !  Local variables: ABORT, SUM                                         C
15     !                                                                      C
16     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
17     !
18           SUBROUTINE CHECK_Mass_balance (init)
19     !
20     !  Include param.inc file to specify parameter values
21     !
22     !-----------------------------------------------
23     !   M o d u l e s
24     !-----------------------------------------------
25           USE param
26           USE param1
27           USE toleranc
28           USE fldvar
29           USE rxns
30           USE geometry
31           USE run
32           USE bc
33           USE constant
34           USE physprop
35           USE indices
36           USE funits
37           USE compar
38           USE mpi_utility
39           USE output
40           USE check
41           USE mflux
42           USE functions
43           IMPLICIT NONE
44     !-----------------------------------------------
45     !   G l o b a l   P a r a m e t e r s
46     !-----------------------------------------------
47     !-----------------------------------------------
48     !   L o c a l   P a r a m e t e r s
49     !-----------------------------------------------
50     !-----------------------------------------------
51     !   L o c a l   V a r i a b l e s
52     !-----------------------------------------------
53     !
54     !                      Flag to check whether this is an initialization call
55     !                      0 -> initialization call; 1 -> integration call
56           INTEGER          init
57     !
58           INTEGER          IER
59     !
60     !                      Solids phase
61           INTEGER          M
62     !
63     !                      Species index
64           INTEGER          N
65     !
66     !                      Do-loop counter
67           INTEGER          L
68     !
69     !                      Starting and ending indices (Make sure these represent a plane and NOT a volume)
70           INTEGER          I1, I2, J1, J2, K1, K2
71     !
72     !
73           DOUBLE PRECISION Accumulation_old, Accumulation_delta, flux, flux_in_tot, flux_out_tot, error_percent
74           DOUBLE PRECISION fin, fout
75     
76     !     functions
77           DOUBLE PRECISION Accumulation, Accumulation_sp
78     !-----------------------------------------------
79     
80           if(report_mass_balance_dt == UNDEFINED) return
81     
82           if(init == 0) then
83     !       allocate arrays
84     !       Initilaize this routine
85             start_time = time
86             report_time = time + report_mass_balance_dt
87     
88             !initialize flux and reaction rate arrays
89             DO L = 1, DIMENSION_BC
90               flux_out_g(L) = ZERO
91               flux_in_g(L) = ZERO
92               DO N = 1, NMAX(0)
93                 flux_out_X_g(L, N) = ZERO
94                 flux_in_X_g(L, N) = ZERO
95               END DO
96     
97               DO M = 1, MMAX
98                 flux_out_s(L, M) = ZERO
99                 flux_in_s(L, M) = ZERO
100                 DO N = 1, NMAX(M)
101                   flux_out_X_s(L, M, N) = ZERO
102                   flux_in_X_s(L, M, N) = ZERO
103                 END DO
104               END DO
105     
106             END DO
107     
108             Integral_SUM_R_g = ZERO
109             DO N = 1, NMAX(0)
110               Integral_R_g(N) = ZERO
111             END DO
112     
113             DO M = 1, MMAX
114               Integral_SUM_R_s(M) = ZERO
115               DO N = 1, NMAX(M)
116                 Integral_R_s(M, N) = ZERO
117               END DO
118             END DO
119     
120             !Accumulation
121             Accumulation_g = Accumulation(ROP_g)
122             DO N = 1, NMAX(0)
123               Accumulation_X_g(N) = Accumulation_sp(ROP_g, X_g(1, N))
124             END DO
125     
126             DO M = 1, MMAX
127               Accumulation_s(M) = Accumulation(ROP_s(1,M))
128               DO N = 1, NMAX(M)
129                 Accumulation_X_s(M, N) = Accumulation_sp(ROP_s(1, M), X_s(1, M, N))
130               END DO
131             END DO
132     
133             return
134     
135     
136           else
137     !       Flow in and out of the boundaries
138     !       Use dt_prev for these integrations because adjust_dt may have changed the dt
139             DO L = 1, DIMENSION_BC
140               IF (BC_DEFINED(L)) THEN
141                 I1 = BC_I_W(L)
142                 I2 = BC_I_E(L)
143                 J1 = BC_J_S(L)
144                 J2 = BC_J_N(L)
145                 K1 = BC_K_B(L)
146                 K2 = BC_K_T(L)
147     
148     !            call Calc_mass_flux(I1, I2, J1, J2, K1, K2, BC_PLANE(L), U_g, V_g, W_g, ROP_g, fin, fout, IER)
149                 IF(.NOT.Added_Mass) THEN
150                   call Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_gE, Flux_gN, Flux_gT, fin, fout, IER)
151                 ELSE
152                   call Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_gSE, Flux_gSN, Flux_gST, fin, fout, IER)
153                 ENDIF
154                 flux_out_g(L) = flux_out_g(L) + fout  * dt_prev
155                 flux_in_g(L) = flux_in_g(L) + fin * dt_prev
156     
157                 DO M = 1, MMAX
158     !              call Calc_mass_flux(I1, I2, J1, J2, K1, K2, BC_PLANE(L), U_s(1,m), V_s(1,m), W_s(1,m), ROP_s(1,m), fin, fout, IER)
159                   IF(.NOT.Added_Mass .OR. M /= M_AM) THEN
160                     call Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_sE(1,m), Flux_sN(1,m), Flux_sT(1,m), fin, fout, IER)
161                   ELSE
162                     call Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_sSE, Flux_sSN, Flux_sST, fin, fout, IER)
163                   ENDIF
164                   flux_out_s(L, M) = flux_out_s(L, M) + fout  * dt_prev
165                   flux_in_s(L, M) = flux_in_s(L, M) + fin * dt_prev
166                 END DO
167     
168               ENDIF
169             END DO
170     
171             Integral_SUM_R_g = Integral_SUM_R_g + Accumulation(SUM_R_g) * dt_prev
172             IF(SPECIES_EQ(0))THEN
173               DO N = 1, NMAX(0)
174                 DO L = 1, DIMENSION_BC
175                   IF (BC_DEFINED(L)) THEN
176                     I1 = BC_I_W(L)
177                     I2 = BC_I_E(L)
178                     J1 = BC_J_S(L)
179                     J2 = BC_J_N(L)
180                     K1 = BC_K_B(L)
181                     K2 = BC_K_T(L)
182     !                call Calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, BC_PLANE(L), U_g, V_g, W_g, ROP_g, X_g(1, N), fin, fout, IER)
183                     IF(.NOT.Added_Mass) THEN
184                       call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_gE, &
185                       Flux_gN, Flux_gT, X_g(1, N), DISCRETIZE(7), fin, fout, IER)
186                     ELSE
187                       call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_gSE, &
188                       Flux_gSN, Flux_gST, X_g(1, N), DISCRETIZE(7), fin, fout, IER)
189                     ENDIF
190                     flux_out_X_g(L, N) = flux_out_X_g(L, N) + fout  * dt_prev
191                     flux_in_X_g(L, N) = flux_in_X_g(L, N) + fin * dt_prev
192                   ENDIF
193                 END DO
194     
195                  Integral_R_g(N) = Integral_R_g(N) + (Accumulation(R_gp(1,N)) - &
196                     Accumulation_sp(ROX_gc(1,N), X_g(1,N)) )* dt_prev
197               END DO
198             ENDIF
199     
200             DO M = 1, MMAX
201               IF(SPECIES_EQ(M))THEN
202                 Integral_SUM_R_s(M) = Integral_SUM_R_s(M) + Accumulation(SUM_R_s(1,M)) * dt_prev
203                 DO N = 1, NMAX(M)
204                   DO L = 1, DIMENSION_BC
205                     IF (BC_DEFINED(L)) THEN
206                       I1 = BC_I_W(L)
207                       I2 = BC_I_E(L)
208                       J1 = BC_J_S(L)
209                       J2 = BC_J_N(L)
210                       K1 = BC_K_B(L)
211                       K2 = BC_K_T(L)
212     !                  call Calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, BC_PLANE(L), &
213     !                  U_s(1,M), V_s(1,M), W_s(1,M), ROP_s(1,M), X_s(1, M, N), fin, fout, &
214                       IF(.NOT.Added_Mass .OR. M /= M_AM) THEN
215                         call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), &
216                         Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M), X_s(1, M, N), DISCRETIZE(7), fin, fout, &
217                         IER)
218                       ELSE
219                         call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), &
220                         Flux_sSE, Flux_sSN, Flux_sST, X_s(1, M, N), DISCRETIZE(7), fin, fout, IER)
221                       ENDIF
222                       flux_out_X_s(L, M, N) = flux_out_X_s(L, M, N) + fout  * dt_prev
223                       flux_in_X_s(L, M, N) = flux_in_X_s(L, M, N) + fin * dt_prev
224                     ENDIF
225                   END DO
226     
227                    Integral_R_s(M, N) = Integral_R_s(M, N) + (Accumulation(R_sp(1,M,N)) - &
228                       Accumulation_sp(ROX_sc(1,M,N), X_s(1,M,N)) )* dt_prev
229                 END DO
230               ENDIF
231             END DO
232     
233           endif
234     
235     
236           if (time >= report_time)then
237     
238             CALL START_LOG
239             IF(DMP_LOG)WRITE(UNIT_LOG, '(/A,G12.5,A,G12.5)') 'Mass balance for interval ', start_time, ' to ', time
240     
241             Accumulation_old = Accumulation_g
242             Accumulation_g = Accumulation(ROP_g)
243             Accumulation_delta = Accumulation_g - Accumulation_old - Integral_SUM_R_g
244             IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Total Fluid Accumulation (g)'
245             IF(DMP_LOG)WRITE(UNIT_LOG, '(4(A,G12.5))') '  Old = ', Accumulation_old, ', New = ', &
246               Accumulation_g, ', Production = ', Integral_SUM_R_g, ', net accu(New - Old - Production) = ', Accumulation_delta
247     
248             IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Integral of boundary flux (g)'
249             IF(DMP_LOG)WRITE(Unit_log, '(A, T8, A, T21, A, T34, A)')'  BC#', 'in', 'out', '(in - out)'
250             flux = zero
251             flux_in_tot = zero
252             flux_out_tot = zero
253             DO L = 1, DIMENSION_BC
254               if(flux_out_g(L) /= ZERO .OR. flux_in_g(L) /= ZERO) then
255                 IF(DMP_LOG)WRITE(Unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')L, flux_in_g(L), flux_out_g(L), (flux_in_g(L)-flux_out_g(L))
256               endif
257               flux = flux + flux_in_g(L) - flux_out_g(L)
258               flux_in_tot = flux_in_tot + flux_in_g(L)
259               flux_out_tot = flux_out_tot + flux_out_g(L)
260             END DO
261             if((flux - Accumulation_delta) /= zero) then
262               error_percent = undefined
263               if(flux_in_tot /= zero) error_percent = (flux - Accumulation_delta)*100./flux_in_tot
264             else
265               error_percent = zero
266             endif
267             IF(DMP_LOG)WRITE(Unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
268             IF(DMP_LOG)WRITE(Unit_log, '(A, G12.5, A, G12.5)')'Error (net influx - net accu) = ', &
269               (flux - Accumulation_delta), ' %Error_g = ', error_percent
270     
271             DO M = 1, MMAX
272               Accumulation_old = Accumulation_s(M)
273               Accumulation_s(M) = Accumulation(ROP_s(1, M))
274               Accumulation_delta = Accumulation_s(M) - Accumulation_old - Integral_SUM_R_s(M)
275               IF(DMP_LOG)WRITE(UNIT_LOG, '(/A, I1, A)') 'Total Solids-', M, ' Accumulation (g)'
276               IF(DMP_LOG)WRITE(UNIT_LOG, '(4(A,G12.5))') '  Old = ', Accumulation_old, ', New = ', &
277                 Accumulation_s(M), ', Production = ', Integral_SUM_R_s(M), ', net accu(New - Old - Production) = ', Accumulation_delta
278     
279               IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Integral of boundary flux (g)'
280               IF(DMP_LOG)WRITE(Unit_log, '(A, T8, A, T21, A, T34, A)')'  BC#', 'in', 'out', '(in - out)'
281               flux = zero
282               flux_in_tot = zero
283               flux_out_tot = zero
284               DO L = 1, DIMENSION_BC
285                 if(flux_out_s(L,M) /= ZERO .OR. flux_in_s(L,M) /= ZERO) then
286                   IF(DMP_LOG) THEN
287                      WRITE(Unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')L, &
288                           flux_in_s(L,M), flux_out_s(L,M), (flux_in_s(L,M)-flux_out_s(L,M))
289                   ENDIF
290                 endif
291                 flux = flux + flux_in_s(L,M) - flux_out_s(L,M)
292                 flux_in_tot = flux_in_tot + flux_in_s(L,M)
293                 flux_out_tot = flux_out_tot + flux_out_s(L,M)
294               END DO
295               if((flux - Accumulation_delta) /= zero) then
296                 if(flux_in_tot /= zero) then
297                   error_percent = (flux - Accumulation_delta)*100./flux_in_tot
298                 else
299                   error_percent = (flux - Accumulation_delta)*100./Accumulation_old
300                 endif
301               else
302                 error_percent = zero
303               endif
304               IF(DMP_LOG)WRITE(Unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
305               IF(DMP_LOG)WRITE(Unit_log, '(A, G12.5, A, I1, A, G12.5)')'Error (net influx - net accu) = ', &
306                 (flux - Accumulation_delta), ' %Error_',M,' = ', error_percent
307             END DO
308     
309     
310             IF(SPECIES_EQ(0) )THEN
311               DO N = 1, NMAX(0)
312     
313                 IF(DMP_LOG)WRITE(UNIT_LOG, '(/A,I2)') 'Gas species - ', N
314     
315                 Accumulation_old = Accumulation_X_g(N)
316                 Accumulation_X_g(N) = Accumulation_sp(ROP_g, X_g(1, N))
317                 Accumulation_delta = Accumulation_X_g(N) - Accumulation_old - Integral_R_g(N)
318                 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Species Accumulation (g)'
319                 IF(DMP_LOG)WRITE(UNIT_LOG, '(4(A,G12.5))') '  Old = ', Accumulation_old, ', New = ', &
320                   Accumulation_X_g(N), ', Production = ', Integral_R_g(N), ', net accu(New - Old - Production) = ', Accumulation_delta
321     
322                 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Integral of boundary flux (g)'
323                 IF(DMP_LOG)WRITE(Unit_log, '(A, T8, A, T21, A, T34, A)')'  BC#', 'in', 'out', '(in - out)'
324                 flux = zero
325                 flux_in_tot = zero
326                 flux_out_tot = zero
327                 DO L = 1, DIMENSION_BC
328                   if(flux_out_X_g(L, N) /= ZERO .OR. flux_in_X_g(L, N) /= ZERO) then
329                     IF(DMP_LOG)WRITE(Unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')L, flux_in_X_g(L, N), flux_out_X_g(L, N), &
330                       (flux_in_X_g(L, N)-flux_out_X_g(L, N))
331                   endif
332                   flux = flux + flux_in_X_g(L, N) - flux_out_X_g(L, N)
333                   flux_in_tot = flux_in_tot + flux_in_X_g(L, N)
334                   flux_out_tot = flux_out_tot + flux_out_X_g(L, N)
335                 END DO
336                 if((flux - Accumulation_delta) /= zero) then
337                   error_percent = undefined
338                   if(flux_in_tot /= zero) then
339                     error_percent = (flux - Accumulation_delta)*100./flux_in_tot
340                   elseif (Integral_R_g(N) /= zero) then
341                     error_percent = (flux - Accumulation_delta)*100./Integral_R_g(N)
342                   endif
343                 else
344                   error_percent = zero
345                 endif
346                 IF(DMP_LOG)WRITE(Unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
347                 IF(DMP_LOG)WRITE(Unit_log, '(A, G12.5, A, I2, A, G12.5)')'Error (net influx - net accu) = ', &
348                  (flux - Accumulation_delta), ' %Error_g(',N,') = ', error_percent
349     
350               END DO
351             ENDIF
352     
353             DO M = 1, MMAX
354               IF(SPECIES_EQ(M) )THEN
355                 DO N = 1, NMAX(M)
356     
357                   IF(DMP_LOG)WRITE(UNIT_LOG, '(/A,I1, A, I2)') 'Solids-', M, ' species - ', N
358     
359                   Accumulation_old = Accumulation_X_s(M,N)
360                   Accumulation_X_s(M,N) = Accumulation_sp(ROP_s(1,M), X_s(1, M, N))
361                   Accumulation_delta = Accumulation_X_s(M,N) - Accumulation_old - Integral_R_s(M,N)
362                   IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Species Accumulation (g)'
363                   IF(DMP_LOG)WRITE(UNIT_LOG, '(4(A,G12.5))') '  Old = ', Accumulation_old, ', New = ', &
364                     Accumulation_X_s(M,N), ', Production = ', Integral_R_s(M,N), &
365                     ', net accu(New - Old - Production) = ', Accumulation_delta
366     
367                   IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Integral of boundary flux (g)'
368                   IF(DMP_LOG)WRITE(Unit_log, '(A, T8, A, T21, A, T34, A)')'  BC#', 'in', 'out', '(in - out)'
369                   flux = zero
370                   flux_in_tot = zero
371                   flux_out_tot = zero
372                   DO L = 1, DIMENSION_BC
373                     if(flux_out_X_s(L, M, N) /= ZERO .OR. flux_in_X_s(L, M, N) /= ZERO) then
374                       IF(DMP_LOG)WRITE(Unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')L, flux_in_X_s(L, M, N), flux_out_X_s(L, M, N), &
375                        (flux_in_X_s(L, M, N)-flux_out_X_s(L, M, N))
376                     endif
377                     flux = flux + flux_in_X_s(L, M, N) - flux_out_X_s(L, M, N)
378                     flux_in_tot = flux_in_tot + flux_in_X_s(L, M, N)
379                     flux_out_tot = flux_out_tot + flux_out_X_s(L, M, N)
380                   END DO
381                   if((flux - Accumulation_delta) /= zero) then
382                     error_percent = undefined
383                     if(flux_in_tot /= zero) then
384                       error_percent = (flux - Accumulation_delta)*100./flux_in_tot
385                     elseif (Integral_R_s(M,N) /= zero)then
386                       error_percent = (flux - Accumulation_delta)*100./Integral_R_s(M,N)
387                     endif
388                   else
389                     error_percent = zero
390                   endif
391                   IF(DMP_LOG)WRITE(Unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
392                   IF(DMP_LOG)WRITE(Unit_log, '(A, G12.5, A, I1, A, I2, A, G12.5)')'Error (net influx - net accu) = ', &
393                     (flux - Accumulation_delta), ' %Error_',M,'(',N,') = ', error_percent
394     
395                 END DO
396               ENDIF
397             END DO
398             IF(DMP_LOG)WRITE(UNIT_LOG, '(/)')
399             CALL END_LOG
400     
401             start_time = time
402             report_time = time + report_mass_balance_dt
403     
404             !initialize flux and reaction rate arrays
405             DO L = 1, DIMENSION_BC
406               flux_out_g(L) = ZERO
407               flux_in_g(L) = ZERO
408               DO N = 1, NMAX(0)
409                 flux_out_X_g(L, N) = ZERO
410                 flux_in_X_g(L, N) = ZERO
411               END DO
412     
413               DO M = 1, MMAX
414                 flux_out_s(L, M) = ZERO
415                 flux_in_s(L, M) = ZERO
416                 DO N = 1, NMAX(M)
417                   flux_out_X_s(L, M, N) = ZERO
418                   flux_in_X_s(L, M, N) = ZERO
419                 END DO
420               END DO
421     
422             END DO
423     
424             Integral_SUM_R_g = ZERO
425             DO N = 1, NMAX(0)
426               Integral_R_g(N) = ZERO
427             END DO
428     
429             DO M = 1, MMAX
430               Integral_SUM_R_s(M) = ZERO
431               DO N = 1, NMAX(M)
432                 Integral_R_s(M, N) = ZERO
433               END DO
434             END DO
435     
436     
437           endif
438     
439     
440           RETURN
441           END SUBROUTINE CHECK_Mass_balance
442     
443     
444     
445     
446     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
447     !                                                                      C
448     !  Module name: Calc_mass_flux                                         C
449     !  Purpose: Calculate the mass flux (g/s) across a plane               C
450     !                                                                      C
451     !  Author: M. Syamlal                                 Date: 9-Dec-02   C
452     !  Reviewer:                                          Date:            C
453     !                                                                      C
454     !  Literature/Document References:                                     C
455     !                                                                      C
456     !  Variables referenced:                                               C
457     !  Variables modified:                                                 C
458     !                                                                      C
459     !                                                                      C
460     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
461     !
462           SUBROUTINE  Calc_mass_flux(I1, I2, J1, J2, K1, K2, Plane, U, V, W, ROP, flux_in, flux_out, IER)
463           USE param
464           USE param1
465           IMPLICIT NONE
466     
467     !
468     !                      Starting and ending indices (Make sure these represent a plane and NOT a volume)
469           INTEGER ::           I1, I2, J1, J2, K1, K2
470     
471     !                      For each (scalar) cell the plane (W, E, S, N, B, T) to be used for flux calc
472           Character ::             Plane
473     
474     !                      Components of Velocity
475           DOUBLE PRECISION ::  U(DIMENSION_3), V(DIMENSION_3), W(DIMENSION_3)
476     
477     !                      Macroscopic density
478           DOUBLE PRECISION ::  ROP(DIMENSION_3)
479     
480     !                      Species Mass fraction
481           DOUBLE PRECISION :: Xn(DIMENSION_3)
482     
483     !                      flux across the plane (composed of many cell faces)
484     !                      flux_in: flux from the scalar cell into the rest of the domain
485     !                      flux_out: flux into the scalar cell from the rest of the domain
486           DOUBLE PRECISION ::  flux_in, flux_out
487     
488           INTEGER ::           IER
489     
490           Xn = one
491           call Calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, PLANE, U, V, W, ROP, Xn, flux_in, flux_out, IER)
492           return
493           end SUBROUTINE  Calc_mass_flux
494     
495     
496     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
497     !                                                                      C
498     !  Module name: Calc_mass_flux_sp                                      C
499     !  Purpose: Calculate the species mass flux (g/s) across a plane       C
500     !                                                                      C
501     !  Author: M. Syamlal                                 Date: 9-Dec-02   C
502     !  Reviewer:                                          Date:            C
503     !                                                                      C
504     !  Literature/Document References:                                     C
505     !                                                                      C
506     !  Variables referenced:                                               C
507     !  Variables modified:                                                 C
508     !                                                                      C
509     !                                                                      C
510     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
511     !
512           SUBROUTINE  Calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, Plane, U, V, W, ROP, Xn, flux_in, flux_out, IER)
513           USE param
514           USE param1
515           USE geometry
516           USE physprop
517           USE indices
518           USE compar
519           USE mpi_utility
520           USE functions
521           IMPLICIT NONE
522     
523     !
524     !                      Starting and ending indices (Make sure these represent a plane and NOT a volume)
525           INTEGER ::           I1, I2, J1, J2, K1, K2
526     
527     !                      For each (scalar) cell the plane (W, E, S, N, B, T) to be used for flux calc
528           Character ::             Plane
529     
530     !                      Components of Velocity
531           DOUBLE PRECISION ::  U(DIMENSION_3), V(DIMENSION_3), W(DIMENSION_3)
532     
533     !                      Macroscopic density
534           DOUBLE PRECISION ::  ROP(DIMENSION_3)
535     
536     !                      Species Mass fraction
537           DOUBLE PRECISION :: Xn(DIMENSION_3)
538     
539     !                      flux across the plane (composed of many cell faces)
540     !                      flux_in: flux from the scalar cell into the rest of the domain
541     !                      flux_out: flux into the scalar cell from the rest of the domain
542           DOUBLE PRECISION ::  flux_in, flux_out
543     
544           INTEGER ::           IER
545     !
546     !                      Indices
547           INTEGER          I, J, K, IJK
548     
549           flux_in = zero
550           flux_out = zero
551           IER = 0
552     
553     
554                 DO K = K1, K2
555                 DO J = J1, J2
556                 DO I = I1, I2
557                   IF(.NOT.IS_ON_myPE_owns(I, J, K)) cycle
558                   IF(DEAD_CELL_AT(I,J,K)) cycle
559     
560                   SELECT CASE (TRIM(PLANE))
561                   CASE ('W')
562                     IJK = IM_OF(FUNIJK(I,J,K))
563                     IF(U(IJK) > ZERO)THEN
564                       flux_out = flux_out + AYZ(IJK) * U(IJK) * ROP(IJK)  * Xn(IJK)
565                     ELSE
566                       flux_in = flux_in - AYZ(IJK) * U(IJK) * ROP(IP_OF(IJK))  * Xn(IP_OF(IJK))
567                     ENDIF
568     
569                   CASE ('E')
570                     IJK = FUNIJK(I,J,K)
571                     IF(U(IJK) > ZERO)THEN
572                       flux_in = flux_in + AYZ(IJK) * U(IJK) * ROP(IJK)  * Xn(IJK)
573                     ELSE
574                       flux_out = flux_out - AYZ(IJK) * U(IJK) * ROP(IP_OF(IJK))  * Xn(IP_OF(IJK))
575                     ENDIF
576     
577                   CASE ('S')
578                     IJK = JM_OF(FUNIJK(I,J,K))
579                     IF(V(IJK) > ZERO)THEN
580                       flux_out = flux_out + AXZ(IJK) * V(IJK) * ROP(IJK)  * Xn(IJK)
581                     ELSE
582                       flux_in = flux_in - AXZ(IJK) * V(IJK) * ROP(JP_OF(IJK))  * Xn(JP_OF(IJK))
583                     ENDIF
584     
585                   CASE ('N')
586                     IJK = FUNIJK(I,J,K)
587                     IF(V(IJK) > ZERO)THEN
588                       flux_in = flux_in + AXZ(IJK) * V(IJK) * ROP(IJK)  * Xn(IJK)
589                     ELSE
590                       flux_out = flux_out - AXZ(IJK) * V(IJK) * ROP(JP_OF(IJK))  * Xn(JP_OF(IJK))
591                     ENDIF
592     
593                   CASE ('B')
594                     IJK = KM_OF(FUNIJK(I,J,K))
595                     IF(W(IJK) > ZERO)THEN
596                       flux_out = flux_out + AXY(IJK) * W(IJK) * ROP(IJK)  * Xn(IJK)
597                     ELSE
598                       flux_in = flux_in - AXY(IJK) * W(IJK) * ROP(KP_OF(IJK))  * Xn(KP_OF(IJK))
599                     ENDIF
600     
601                   CASE ('T')
602                     IJK = FUNIJK(I,J,K)
603                     IF(W(IJK) > ZERO)THEN
604                       flux_in = flux_in + AXY(IJK) * W(IJK) * ROP(IJK)  * Xn(IJK)
605                     ELSE
606                       flux_out = flux_out - AXY(IJK) * W(IJK) * ROP(KP_OF(IJK))  *  Xn(KP_OF(IJK))
607                     ENDIF
608     
609                   CASE DEFAULT
610     !                IER = 1
611     !                CALL START_LOG
612     !                IF(DMP_LOG)WRITE (UNIT_LOG, '(A, A1)' ) 'From: Calc_mass_flux, Unknown Plane: ', Plane
613     !                CALL MFIX_EXIT(myPE)
614     
615                   END SELECT
616                 ENDDO
617                 ENDDO
618                 ENDDO
619     
620                 call global_all_sum(flux_in)
621                 call global_all_sum(flux_out)
622     
623           return
624           end SUBROUTINE  Calc_mass_flux_sp
625     
626     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
627     !                                                                      C
628     !  Module name: Calc_mass_fluxHR                                       C
629     !  Purpose: Calculate the mass flux (g/s) across a plane. HR version.  C
630     !                                                                      C
631     !  Author: M. Syamlal                                 Date: 26-Jul-05   C
632     !  Reviewer:                                          Date:            C
633     !                                                                      C
634     !  Literature/Document References:                                     C
635     !                                                                      C
636     !  Variables referenced:                                               C
637     !  Variables modified:                                                 C
638     !                                                                      C
639     !                                                                      C
640     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
641     !
642           SUBROUTINE  Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, Plane, Flux_E, Flux_N, Flux_T, flux_in, flux_out, IER)
643           USE param
644           USE param1
645           IMPLICIT NONE
646     
647     !
648     !                      Starting and ending indices (Make sure these represent a plane and NOT a volume)
649           INTEGER ::           I1, I2, J1, J2, K1, K2
650     
651     !                      For each (scalar) cell the plane (W, E, S, N, B, T) to be used for flux calc
652           Character ::             Plane
653     
654     !                      Components of flux
655           DOUBLE PRECISION ::  Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
656     
657     !                      Species Mass fraction
658           DOUBLE PRECISION :: Xn(DIMENSION_3)
659     
660     !                      flux across the plane (composed of many cell faces)
661     !                      flux_in: flux from the scalar cell into the rest of the domain
662     !                      flux_out: flux into the scalar cell from the rest of the domain
663           DOUBLE PRECISION ::  flux_in, flux_out
664     
665           INTEGER ::           IER
666     
667           Xn = one
668           call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, PLANE, Flux_E, Flux_N, Flux_T, Xn, 0, flux_in, flux_out, IER)
669           return
670           end SUBROUTINE  Calc_mass_fluxHR
671     
672     
673     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
674     !                                                                      C
675     !  Module name: Calc_mass_flux_spHR                                    C
676     !  Purpose: Calculate the species mass flux (g/s) across a plane. HR version       C
677     !                                                                      C
678     !  Author: M. Syamlal                                 Date: 26-Jul-05   C
679     !  Reviewer:                                          Date:            C
680     !                                                                      C
681     !  Literature/Document References:                                     C
682     !                                                                      C
683     !  Variables referenced:                                               C
684     !  Variables modified:                                                 C
685     !                                                                      C
686     !                                                                      C
687     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
688     !
689           SUBROUTINE  Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, Plane, Flux_E, Flux_N, Flux_T, Xn, disc, flux_in, flux_out, IER)
690           USE param
691           USE param1
692           USE geometry
693           USE physprop
694           USE indices
695           USE compar
696           Use xsi_array
697           USE mpi_utility
698           USE xsi
699           USE functions
700           IMPLICIT NONE
701     
702     !
703     !                      Starting and ending indices (Make sure these represent a plane and NOT a volume)
704           INTEGER ::           I1, I2, J1, J2, K1, K2
705     
706     !                      For each (scalar) cell the plane (W, E, S, N, B, T) to be used for flux calc
707           Character ::             Plane
708     
709     !                      Components of flux
710           DOUBLE PRECISION ::  Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
711     
712     !                      Species Mass fraction
713           DOUBLE PRECISION :: Xn(DIMENSION_3)
714     
715     !
716     !                      Discretizationindex
717           INTEGER          Disc
718     
719     !                      flux across the plane (composed of many cell faces)
720     !                      flux_in: flux from the scalar cell into the rest of the domain
721     !                      flux_out: flux into the scalar cell from the rest of the domain
722           DOUBLE PRECISION ::  flux_in, flux_out
723     
724     !                       Flux and Xn at the boundary plane
725           DOUBLE PRECISION ::  Flux, PHI_HO
726     
727           INTEGER ::           IER
728     !
729     !                      Indices
730           INTEGER          I, J, K, IJK
731     
732           call lock_xsi_array
733     
734     !     Unlike other calls to calc_xsi, the fluxes rather than the velocities are passed.
735     !     This should work fine because the velocities are used to determine the upwind bias,
736     !     for which fluxes should suffice. The flag incr is not needed and is set to 0.
737           CALL CALC_XSI (DISC, Xn, Flux_E, Flux_N, Flux_T, XSI_E, XSI_N, XSI_T, 0)
738     
739           flux_in = zero
740           flux_out = zero
741           IER = 0
742     
743                 DO K = K1, K2
744                 DO J = J1, J2
745                 DO I = I1, I2
746                   IF(.NOT.IS_ON_myPE_owns(I, J, K)) cycle
747                   IF(DEAD_CELL_AT(I,J,K)) cycle
748     
749                   SELECT CASE (TRIM(PLANE))
750                   CASE ('W')
751                     IJK = IM_OF(FUNIJK(I,J,K))
752                     PHI_HO = XSI_E(IJK)*Xn(IP_OF(IJK))+(1.0-XSI_E(IJK))*Xn(IJK)
753                     Flux = FLUX_E(IJK)
754     
755                   CASE ('E')
756                     IJK = FUNIJK(I,J,K)
757                     PHI_HO = XSI_E(IJK)*Xn(IP_OF(IJK))+(1.0-XSI_E(IJK))*Xn(IJK)
758                     Flux = -FLUX_E(IJK)
759     
760                   CASE ('S')
761                     IJK = JM_OF(FUNIJK(I,J,K))
762                     PHI_HO = XSI_N(IJK)*Xn(JP_OF(IJK))+(1.0-XSI_N(IJK))*Xn(IJK)
763                     Flux = FLUX_N(IJK)
764     
765                   CASE ('N')
766                     IJK = FUNIJK(I,J,K)
767                     PHI_HO = XSI_N(IJK)*Xn(JP_OF(IJK))+(1.0-XSI_N(IJK))*Xn(IJK)
768                     Flux = -FLUX_N(IJK)
769     
770                   CASE ('B')
771                     IJK = KM_OF(FUNIJK(I,J,K))
772                     PHI_HO = XSI_T(IJK)*Xn(KP_OF(IJK))+(1.0-XSI_T(IJK))*Xn(IJK)
773                     Flux = FLUX_T(IJK)
774     
775                   CASE ('T')
776                     IJK = FUNIJK(I,J,K)
777                     PHI_HO = XSI_T(IJK)*Xn(KP_OF(IJK))+(1.0-XSI_T(IJK))*Xn(IJK)
778                     Flux = -FLUX_T(IJK)
779     
780                   CASE DEFAULT
781                     PHI_HO = ZERO
782                     Flux = ZERO
783     
784                   END SELECT
785     
786                   IF(FLUX > ZERO)THEN
787                     flux_out = flux_out + FLUX * PHI_HO
788                   ELSE
789                     flux_in = flux_in - FLUX * PHI_HO
790                   ENDIF
791     
792                 ENDDO
793                 ENDDO
794                 ENDDO
795     
796                 call global_all_sum(flux_in)
797                 call global_all_sum(flux_out)
798           call unlock_xsi_array
799     
800           return
801           end SUBROUTINE  Calc_mass_flux_spHR
802     
803     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
804     !                                                                      C
805     !  Module name: Accumulation(ro)                                       C
806     !  Purpose: Intergrate density accumulation over the entire domain     C
807     !                                                                      C
808     !  Author: M. Syamlal                                 Date: 30-DEC-02  C
809     !  Local variables:  None                                              C
810     !                                                                      C
811     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
812     !
813           DOUBLE PRECISION FUNCTION Accumulation(ro)
814     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
815     !...Switches: -xf
816     !
817     !  Include param.inc file to specify parameter values
818     !
819     !-----------------------------------------------
820     !   M o d u l e s
821     !-----------------------------------------------
822           USE param
823           USE param1
824           IMPLICIT NONE
825     !-----------------------------------------------
826     !   G l o b a l   P a r a m e t e r s
827     !-----------------------------------------------
828     !-----------------------------------------------
829     !   D u m m y   A r g u m e n t s
830     !-----------------------------------------------
831     !                      density distributiom
832           DOUBLE PRECISION ::  RO(DIMENSION_3)
833     !                      Species Mass fraction
834           DOUBLE PRECISION :: Xn(DIMENSION_3)
835     
836     !     functions
837           DOUBLE PRECISION Accumulation_sp
838     !
839     !-----------------------------------------------
840     
841           Xn = ONE
842           Accumulation = Accumulation_sp(ro, xn)
843     
844           RETURN
845           END FUNCTION Accumulation
846     
847     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
848     !                                                                      C
849     !  Module name: Accumulation_sp(ro, Xn)                                C
850     !  Purpose: Intergrate density accumulation over the entire domain     C
851     !                                                                      C
852     !  Author: M. Syamlal                                 Date: 30-DEC-02  C
853     !  Local variables:  None                                              C
854     !                                                                      C
855     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
856     !
857           DOUBLE PRECISION FUNCTION Accumulation_sp(ro, Xn)
858     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
859     !...Switches: -xf
860     !
861     !  Include param.inc file to specify parameter values
862     !
863     !-----------------------------------------------
864     !   M o d u l e s
865     !-----------------------------------------------
866           USE param
867           USE param1
868           USE parallel
869           USE physprop
870           USE geometry
871           USE indices
872           USE compar
873           USE mpi_utility
874           USE functions
875           IMPLICIT NONE
876     !-----------------------------------------------
877     !   G l o b a l   P a r a m e t e r s
878     !-----------------------------------------------
879     !-----------------------------------------------
880     !   D u m m y   A r g u m e n t s
881     !-----------------------------------------------
882     !                      Indices
883           INTEGER          IJK
884     
885     !                      density distributiom
886           DOUBLE PRECISION ::  RO(DIMENSION_3)
887     !                      Species Mass fraction
888           DOUBLE PRECISION :: Xn(DIMENSION_3)
889     
890           DOUBLE PRECISION SUM
891     !
892     !-----------------------------------------------
893     
894           SUM = ZERO
895     !
896     !$omp   parallel do default(none) private(IJK) shared(ijkstart3, ijkend3, i_of, j_of, k_of, ro, xn, vol) &
897     !$omp   reduction(+:SUM)
898           DO IJK = ijkstart3, ijkend3
899             IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK), J_OF(IJK), K_OF(IJK))) cycle
900               IF (FLUID_AT(IJK)) SUM = SUM + RO(IJK) * Xn(IJK) * VOL(IJK)
901           END DO
902     
903           call global_all_sum(sum)
904           Accumulation_sp = sum
905     
906           RETURN
907           END FUNCTION Accumulation_sp
908     
909     
910     
911           DOUBLE PRECISION FUNCTION check_conservation (Phi, A_m, B_m, IJK)
912           USE param
913           USE param1
914           USE geometry
915           USE indices
916           USE matrix
917           USE compar
918           USE functions
919           IMPLICIT NONE
920     !
921     !                      Septadiagonal matrix A_m
922           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
923     !
924     !                      Vector b_m
925           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
926     
927           DOUBLE PRECISION Phi(DIMENSION_3)
928     
929     !                      Indices
930           INTEGER          IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP
931     
932               IMJK = IM_OF(IJK)
933               IJMK = JM_OF(IJK)
934               IJKM = KM_OF(IJK)
935               IPJK = IP_OF(IJK)
936               IJPK = JP_OF(IJK)
937               IJKP = KP_OF(IJK)
938     
939                check_conservation =  Phi(IJK) * A_m(IJK, 0, 0)  &
940                        + Phi(IPJK) * A_m(IJK, E, 0) &
941                        + Phi(IMJK) * A_m(IJK, W, 0) &
942                        + Phi(IJPK) * A_m(IJK, N, 0) &
943                        + Phi(IJMK) * A_m(IJK, S, 0) &
944                        + Phi(IJKP) * A_m(IJK, T, 0) &
945                        + Phi(IJKM) * A_m(IJK, B, 0) &
946                        - B_m(IJK,0)
947     
948           END FUNCTION check_conservation
949     
950     
951     
952