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