File: N:\mfix\model\check_mod.f

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