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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: B_m_P_star_e(IJK, IER)                                 C
4     !  Purpose: Determine p_star when there is an interface at east        C
5     !                                                                      C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 6-AUG-96   C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !                                                                      C
11     !  Literature/Document References:                                     C
12     !                                                                      C
13     !  Variables referenced:                                               C
14     !  Variables modified:                                                 C
15     !                                                                      C
16     !  Local variables:                                                    C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     !
20           SUBROUTINE B_M_P_STAR_E(B_M, IJK)
21     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
22     !...Switches: -xf
23     !
24     !  Include param.inc file to specify parameter values
25     !
26     !-----------------------------------------------
27     !   M o d u l e s
28     !-----------------------------------------------
29           USE param
30           USE param1
31           USE scales
32           USE constant
33           USE physprop
34           USE fldvar
35           USE run
36           USE rxns
37           USE toleranc
38           USE geometry
39           USE indices
40           USE compar
41           USE bodyforce
42           USE fun_avg
43           USE functions
44           IMPLICIT NONE
45     !-----------------------------------------------
46     !   G l o b a l   P a r a m e t e r s
47     !-----------------------------------------------
48     !-----------------------------------------------
49     !   D u m m y   A r g u m e n t s
50     !-----------------------------------------------
51     !
52     !                      Indices
53           INTEGER          I, J, K, IJK, IJKE
54     !
55     !                      Phase index
56           INTEGER          M
57     !
58     !                      Average volume fraction
59           DOUBLE PRECISION EPGA
60     !
61     !                      Average density
62           DOUBLE PRECISION ROPGA
63     !
64     !                      LHS (similar to A_m)
65           DOUBLE PRECISION A
66     !
67     !                      RHS (similar to b_m)
68           DOUBLE PRECISION B
69     !
70     !                      b_m
71           DOUBLE PRECISION b_m
72     !
73     !                      sum of ep_s
74           DOUBLE PRECISION Eps
75     !
76     !                      Source terms (Surface)
77           DOUBLE PRECISION Sdp, Sdps
78     !
79     !                      Source terms (Volumetric)
80           DOUBLE PRECISION V0, Vmt, Vbf
81     !-----------------------------------------------
82     
83           I = I_OF(IJK)
84           J = J_OF(IJK)
85           K = K_OF(IJK)
86           IJKE = EAST_OF(IJK)
87     !
88           A = ZERO
89           B = ZERO
90           EPS = ZERO
91     !
92           DO M = 1, MMAX
93              IF (CLOSE_PACKED(M)) THEN
94                 EPGA = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
95     !
96     !         Surface forces
97     !
98     !         Pressure term
99                 SDP = -P_SCALE*EPGA*(P_G(IJKE)-P_G(IJK))*AYZ(IJK)
100                 SDPS = -EPGA*(P_S(IJKE,M)-P_S(IJK,M))*AYZ(IJK)
101     !
102     !           Shear stress terms
103     !
104     !         Volumetric forces
105                 ROPGA = AVG_X(ROP_S(IJK,M),ROP_S(IJKE,M),I)
106     !
107     !         Previous time step
108                 V0 = AVG_X(ROP_SO(IJK,M),ROP_SO(IJKE,M),I)*ODT
109     !
110     !         Interphase mass transfer
111                 VMT = AVG_X(SUM_R_S(IJK,M),SUM_R_S(IJKE,M),I)
112     !
113     !         Body force
114                 VBF = ROPGA*BFX_S(IJK,M)
115     !
116     !         Collect the terms
117                 A = A - ((V0 + ZMAX(VMT))*VOL_U(IJK))
118                 B=B-(SDP+SDPS+((V0+ZMAX((-VMT)))*U_SO(IJK,M)+VBF)*VOL_U(IJK))
119                 EPS = EPS + EPGA
120              ENDIF
121           END DO
122           B_M = -((B - A)/(EPS*AYZ(IJK))-P_STAR(IJK))
123     !
124           RETURN
125           END SUBROUTINE B_M_P_STAR_E
126     !
127     !
128     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
129     !                                                                      C
130     !  Module name: B_m_P_star_n(IJK, IER)                                 C
131     !  Purpose: Determine p_star when there is an interface at north       C
132     !                                                                      C
133     !                                                                      C
134     !  Author: M. Syamlal                                 Date: 6-AUG-96   C
135     !  Reviewer:                                          Date:            C
136     !                                                                      C
137     !                                                                      C
138     !  Literature/Document References:                                     C
139     !                                                                      C
140     !  Variables referenced:                                               C
141     !  Variables modified:                                                 C
142     !                                                                      C
143     !  Local variables:                                                    C
144     !                                                                      C
145     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
146     !
147           SUBROUTINE B_M_P_STAR_N(B_M, IJK)
148     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
149     !...Switches: -xf
150     !
151     !  Include param.inc file to specify parameter values
152     !
153     !-----------------------------------------------
154     !   M o d u l e s
155     !-----------------------------------------------
156           USE param
157           USE param1
158           USE scales
159           USE constant
160           USE physprop
161           USE fldvar
162           USE run
163           USE rxns
164           USE toleranc
165           USE geometry
166           USE indices
167           USE compar
168           USE bodyforce
169           USE fun_avg
170           USE functions
171           IMPLICIT NONE
172     !-----------------------------------------------
173     !   G l o b a l   P a r a m e t e r s
174     !-----------------------------------------------
175     !-----------------------------------------------
176     !   D u m m y   A r g u m e n t s
177     !-----------------------------------------------
178     !
179     !                      Indices
180           INTEGER          I, J, K, IJK, IJKN
181     !
182     !                      Phase index
183           INTEGER          M
184     !
185     !                      Average volume fraction
186           DOUBLE PRECISION EPGA
187     !
188     !                      Average density
189           DOUBLE PRECISION ROPGA
190     !
191     !                      LHS (similar to A_m)
192           DOUBLE PRECISION A
193     !
194     !                      RHS (similar to b_m)
195           DOUBLE PRECISION B
196     !
197     !                      b_m
198           DOUBLE PRECISION b_m
199     !
200     !                      sum of ep_s
201           DOUBLE PRECISION Eps
202     !
203     !                      Source terms (Surface)
204           DOUBLE PRECISION Sdp, Sdps
205     !
206     !                      Source terms (Volumetric)
207           DOUBLE PRECISION V0, Vmt, Vbf
208     !
209     !-----------------------------------------------
210     
211           I = I_OF(IJK)
212           J = J_OF(IJK)
213           K = K_OF(IJK)
214           IJKN = NORTH_OF(IJK)
215     !
216           A = ZERO
217           B = ZERO
218           EPS = ZERO
219     !
220           DO M = 1, MMAX
221              IF (CLOSE_PACKED(M)) THEN
222                 EPGA = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
223     !
224     !         Surface forces
225     !
226     !         Pressure term
227                 SDP = -P_SCALE*EPGA*(P_G(IJKN)-P_G(IJK))*AXZ(IJK)
228                 SDPS = -EPGA*(P_S(IJKN,M)-P_S(IJK,M))*AXZ(IJK)
229     !
230     !           Shear stress terms
231     !
232     !         Volumetric forces
233                 ROPGA = AVG_Y(ROP_S(IJK,M),ROP_S(IJKN,M),J)
234     !
235     !         Previous time step
236                 V0 = AVG_Y(ROP_SO(IJK,M),ROP_SO(IJKN,M),J)*ODT
237     !
238     !         Interphase mass transfer
239                 VMT = AVG_Y(SUM_R_S(IJK,M),SUM_R_S(IJKN,M),J)
240     !
241     !         Body force
242                 VBF = ROPGA*BFY_S(IJK,M)
243     !
244     !         Collect the terms
245                 A = A - ((V0 + ZMAX(VMT))*VOL_V(IJK))
246                 B=B-(SDP+SDPS+((V0+ZMAX((-VMT)))*V_SO(IJK,M)+VBF)*VOL_V(IJK))
247                 EPS = EPS + EPGA
248              ENDIF
249           END DO
250           B_M = -((B - A)/(EPS*AXZ(IJK))-P_STAR(IJK))
251     !
252           RETURN
253           END SUBROUTINE B_M_P_STAR_N
254     !
255     !
256     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
257     !                                                                      C
258     !  Module name: B_m_P_star_t(IJK, IER)                                 C
259     !  Purpose: Determine p_star when there is an interface at top         C
260     !                                                                      C
261     !                                                                      C
262     !  Author: M. Syamlal                                 Date: 6-AUG-96   C
263     !  Reviewer:                                          Date:            C
264     !                                                                      C
265     !                                                                      C
266     !  Literature/Document References:                                     C
267     !                                                                      C
268     !  Variables referenced:                                               C
269     !  Variables modified:                                                 C
270     !                                                                      C
271     !  Local variables:                                                    C
272     !                                                                      C
273     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
274     !
275           SUBROUTINE B_M_P_STAR_T(B_M, IJK)
276     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
277     !...Switches: -xf
278     !
279     !  Include param.inc file to specify parameter values
280     !
281     !-----------------------------------------------
282     !   M o d u l e s
283     !-----------------------------------------------
284           USE param
285           USE param1
286           USE scales
287           USE constant
288           USE physprop
289           USE fldvar
290           USE run
291           USE rxns
292           USE toleranc
293           USE geometry
294           USE indices
295           USE compar
296           USE bodyforce
297           USE fun_avg
298           USE functions
299           IMPLICIT NONE
300     !-----------------------------------------------
301     !   G l o b a l   P a r a m e t e r s
302     !-----------------------------------------------
303     !-----------------------------------------------
304     !   D u m m y   A r g u m e n t s
305     !-----------------------------------------------
306     !
307     !                      Indices
308           INTEGER          I, J, K, IJK, IJKT
309     !
310     !                      Phase index
311           INTEGER          M
312     !
313     !                      Average volume fraction
314           DOUBLE PRECISION EPGA
315     !
316     !                      Average density
317           DOUBLE PRECISION ROPGA
318     !
319     !                      LHS (similar to A_m)
320           DOUBLE PRECISION A
321     !
322     !                      RHS (similar to b_m)
323           DOUBLE PRECISION B
324     !
325     !                      b_m
326           DOUBLE PRECISION b_m
327     !
328     !                      sum of ep_s
329           DOUBLE PRECISION Eps
330     !
331     !                      Source terms (Surface)
332           DOUBLE PRECISION Sdp, Sdps
333     !
334     !                      Source terms (Volumetric)
335           DOUBLE PRECISION V0, Vmt, Vbf
336     !
337     !-----------------------------------------------
338     
339           I = I_OF(IJK)
340           J = J_OF(IJK)
341           K = K_OF(IJK)
342           IJKT = TOP_OF(IJK)
343     !
344           A = ZERO
345           B = ZERO
346           EPS = ZERO
347     !
348           DO M = 1, MMAX
349              IF (CLOSE_PACKED(M)) THEN
350                 EPGA = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
351     !
352     !         Surface forces
353     !
354     !         Pressure term
355                 SDP = -P_SCALE*EPGA*(P_G(IJKT)-P_G(IJK))*AXY(IJK)
356                 SDPS = -EPGA*(P_S(IJKT,M)-P_S(IJK,M))*AXY(IJK)
357     !
358     !           Shear stress terms
359     !
360     !         Volumetric forces
361                 ROPGA = AVG_Z(ROP_S(IJK,M),ROP_S(IJKT,M),K)
362     !
363     !         Previous time step
364                 V0 = AVG_Z(ROP_SO(IJK,M),ROP_SO(IJKT,M),K)*ODT
365     !
366     !         Interphase mass transfer
367                 VMT = AVG_Z(SUM_R_S(IJK,M),SUM_R_S(IJKT,M),K)
368     !
369     !         Body force
370                 VBF = ROPGA*BFZ_S(IJK,M)
371     !
372     !         Collect the terms
373                 A = A - ((V0 + ZMAX(VMT))*VOL_W(IJK))
374                 B=B-(SDP+SDPS+((V0+ZMAX((-VMT)))*W_SO(IJK,M)+VBF)*VOL_W(IJK))
375                 EPS = EPS + EPGA
376              ENDIF
377           END DO
378           B_M = -((B - A)/(EPS*AXY(IJK))-P_STAR(IJK))
379     !
380           RETURN
381           END SUBROUTINE B_M_P_STAR_T
382     !
383     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
384     !                                                                      C
385     !  Module name: B_m_P_star_w(b_m, IJK, IER)                            C
386     !  Purpose: Determine p_star when there is an interface at west        C
387     !                                                                      C
388     !                                                                      C
389     !  Author: M. Syamlal                                 Date: 6-AUG-96   C
390     !  Reviewer:                                          Date:            C
391     !                                                                      C
392     !                                                                      C
393     !  Literature/Document References:                                     C
394     !                                                                      C
395     !  Variables referenced:                                               C
396     !  Variables modified:                                                 C
397     !                                                                      C
398     !  Local variables:                                                    C
399     !                                                                      C
400     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
401     !
402           SUBROUTINE B_M_P_STAR_W(B_M, IJK)
403     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
404     !...Switches: -xf
405     !
406     !  Include param.inc file to specify parameter values
407     !
408     !-----------------------------------------------
409     !   M o d u l e s
410     !-----------------------------------------------
411           USE param
412           USE param1
413           USE scales
414           USE constant
415           USE physprop
416           USE fldvar
417           USE run
418           USE rxns
419           USE toleranc
420           USE geometry
421           USE indices
422           USE compar
423           USE bodyforce
424           USE fun_avg
425           USE functions
426           IMPLICIT NONE
427     !-----------------------------------------------
428     !   G l o b a l   P a r a m e t e r s
429     !-----------------------------------------------
430     !-----------------------------------------------
431     !   D u m m y   A r g u m e n t s
432     !-----------------------------------------------
433     !
434     !                      Indices
435           INTEGER          IM, J, K, IJK, IJKW, IMJK
436     !
437     !                      Phase index
438           INTEGER          M
439     !
440     !                      Average volume fraction
441           DOUBLE PRECISION EPGA
442     !
443     !                      Average density
444           DOUBLE PRECISION ROPGA
445     !
446     !                      LHS (similar to A_m)
447           DOUBLE PRECISION A
448     !
449     !                      RHS (similar to b_m)
450           DOUBLE PRECISION B
451     !
452     !                      b_m
453           DOUBLE PRECISION b_m
454     !
455     !                      sum of ep_s
456           DOUBLE PRECISION Eps
457     !
458     !                      Source terms (Surface)
459           DOUBLE PRECISION Sdp, Sdps
460     !
461     !                      Source terms (Volumetric)
462           DOUBLE PRECISION V0, Vmt, Vbf
463     !
464     !-----------------------------------------------
465     
466           IM = IM1(I_OF(IJK))
467           J = J_OF(IJK)
468           K = K_OF(IJK)
469           IJKW = WEST_OF(IJK)
470           IMJK = IM_OF(IJK)
471     !
472           A = ZERO
473           B = ZERO
474           EPS = ZERO
475     !
476           DO M = 1, MMAX
477              IF (CLOSE_PACKED(M)) THEN
478                 EPGA = AVG_X(EP_S(IJKW,M),EP_S(IJK,M),IM)
479     !
480     !         Surface forces
481     !
482     !         Pressure term
483                 SDP = -P_SCALE*EPGA*(P_G(IJK)-P_G(IJKW))*AYZ(IMJK)
484                 SDPS = -EPGA*(P_S(IJK,M)-P_S(IJKW,M))*AYZ(IMJK)
485     !
486     !           Shear stress terms
487     !
488     !         Volumetric forces
489                 ROPGA = AVG_X(ROP_S(IJKW,M),ROP_S(IJK,M),IM)
490     !
491     !         Previous time step
492                 V0 = AVG_X(ROP_SO(IJKW,M),ROP_SO(IJK,M),IM)*ODT
493     !
494     !         Interphase mass transfer
495                 VMT = AVG_X(SUM_R_S(IJKW,M),SUM_R_S(IJK,M),IM)
496     !
497     !         Body force
498                 VBF = ROPGA*BFX_S(IJKW,M)
499     !
500     !         Collect the terms
501                 A = A - ((V0 + ZMAX(VMT))*VOL_U(IJKW))
502                 B=B-(SDP+SDPS+((V0+ZMAX((-VMT)))*U_SO(IMJK,M)+VBF)*VOL_U(IJKW))
503                 EPS = EPS + EPGA
504              ENDIF
505           END DO
506           B_M = -(((-(B - A)/(EPS*AYZ(IMJK))))-P_STAR(IJK))
507     !
508           RETURN
509           END SUBROUTINE B_M_P_STAR_W
510     !
511     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
512     !                                                                      C
513     !  Module name: B_m_P_star_s(b_m, IJK, IER)                            C
514     !  Purpose: Determine p_star when there is an interface at south       C
515     !                                                                      C
516     !                                                                      C
517     !  Author: M. Syamlal                                 Date: 6-AUG-96   C
518     !  Reviewer:                                          Date:            C
519     !                                                                      C
520     !                                                                      C
521     !  Literature/Document References:                                     C
522     !                                                                      C
523     !  Variables referenced:                                               C
524     !  Variables modified:                                                 C
525     !                                                                      C
526     !  Local variables:                                                    C
527     !                                                                      C
528     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
529     !
530           SUBROUTINE B_M_P_STAR_S(B_M, IJK)
531     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
532     !...Switches: -xf
533     !
534     !  Include param.inc file to specify parameter values
535     !
536     !-----------------------------------------------
537     !   M o d u l e s
538     !-----------------------------------------------
539           USE param
540           USE param1
541           USE scales
542           USE constant
543           USE physprop
544           USE fldvar
545           USE run
546           USE rxns
547           USE toleranc
548           USE geometry
549           USE indices
550           USE compar
551           USE bodyforce
552           USE fun_avg
553           USE functions
554           IMPLICIT NONE
555     !-----------------------------------------------
556     !   G l o b a l   P a r a m e t e r s
557     !-----------------------------------------------
558     !-----------------------------------------------
559     !   D u m m y   A r g u m e n t s
560     !-----------------------------------------------
561     !
562     !
563     !                      Indices
564           INTEGER          I, JM, K, IJK, IJKS, IJMK
565     !
566     !                      Phase index
567           INTEGER          M
568     !
569     !                      Average volume fraction
570           DOUBLE PRECISION EPGA
571     !
572     !                      Average density
573           DOUBLE PRECISION ROPGA
574     !
575     !                      LHS (similar to A_m)
576           DOUBLE PRECISION A
577     !
578     !                      RHS (similar to b_m)
579           DOUBLE PRECISION B
580     !
581     !                      b_m
582           DOUBLE PRECISION b_m
583     !
584     !                      sum of ep_s
585           DOUBLE PRECISION Eps
586     !
587     !                      Source terms (Surface)
588           DOUBLE PRECISION Sdp, Sdps
589     !
590     !                      Source terms (Volumetric)
591           DOUBLE PRECISION V0, Vmt, Vbf
592     !-----------------------------------------------
593     
594           I = I_OF(IJK)
595           JM = JM1(J_OF(IJK))
596           K = K_OF(IJK)
597           IJKS = SOUTH_OF(IJK)
598           IJMK = JM_OF(IJK)
599     !
600           A = ZERO
601           B = ZERO
602           EPS = ZERO
603     !
604           DO M = 1, MMAX
605              IF (CLOSE_PACKED(M)) THEN
606                 EPGA = AVG_Y(EP_S(IJKS,M),EP_S(IJK,M),JM)
607     !
608     !         Surface forces
609     !
610     !         Pressure term
611                 SDP = -P_SCALE*EPGA*(P_G(IJK)-P_G(IJKS))*AXZ(IJK)
612                 SDPS = -EPGA*(P_S(IJK,M)-P_S(IJKS,M))*AXZ(IJK)
613     !
614     !           Shear stress terms
615     !
616     !         Volumetric forces
617                 ROPGA = AVG_Y(ROP_S(IJKS,M),ROP_S(IJK,M),JM)
618     !
619     !         Previous time step
620                 V0 = AVG_Y(ROP_SO(IJKS,M),ROP_SO(IJK,M),JM)*ODT
621     !
622     !         Interphase mass transfer
623                 VMT = AVG_Y(SUM_R_S(IJKS,M),SUM_R_S(IJK,M),JM)
624     !
625     !         Body force
626                 VBF = ROPGA*BFY_S(IJKS,M)
627     !
628     !         Collect the terms
629                 A = A - ((V0 + ZMAX(VMT))*VOL_V(IJKS))
630                 B=B-(SDP+SDPS+((V0+ZMAX((-VMT)))*V_SO(IJMK,M)+VBF)*VOL_V(IJKS))
631                 EPS = EPS + EPGA
632              ENDIF
633           END DO
634           B_M = -(((-(B - A)/(EPS*AXZ(IJK))))-P_STAR(IJK))
635     !
636           RETURN
637           END SUBROUTINE B_M_P_STAR_S
638     !
639     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
640     !                                                                      C
641     !  Module name: B_m_P_star_b(b_m, IJK, IER)                            C
642     !  Purpose: Determine p_star when there is an interface at bottom      C
643     !                                                                      C
644     !                                                                      C
645     !  Author: M. Syamlal                                 Date: 6-AUG-96   C
646     !  Reviewer:                                          Date:            C
647     !                                                                      C
648     !                                                                      C
649     !  Literature/Document References:                                     C
650     !                                                                      C
651     !  Variables referenced:                                               C
652     !  Variables modified:                                                 C
653     !                                                                      C
654     !  Local variables:                                                    C
655     !                                                                      C
656     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
657     !
658           SUBROUTINE B_M_P_STAR_B(B_M, IJK)
659     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
660     !...Switches: -xf
661     !
662     !  Include param.inc file to specify parameter values
663     !
664     !-----------------------------------------------
665     !   M o d u l e s
666     !-----------------------------------------------
667           USE param
668           USE param1
669           USE scales
670           USE constant
671           USE physprop
672           USE fldvar
673           USE run
674           USE rxns
675           USE toleranc
676           USE geometry
677           USE indices
678           USE compar
679           USE bodyforce
680           USE fun_avg
681           USE functions
682           IMPLICIT NONE
683     !-----------------------------------------------
684     !   G l o b a l   P a r a m e t e r s
685     !-----------------------------------------------
686     !-----------------------------------------------
687     !   D u m m y   A r g u m e n t s
688     !-----------------------------------------------
689     !
690     !                      Indices
691           INTEGER          I, J, KM, IJK, IJKB, IJKM
692     !
693     !                      Phase index
694           INTEGER          M
695     !
696     !                      Average volume fraction
697           DOUBLE PRECISION EPGA
698     !
699     !                      Average density
700           DOUBLE PRECISION ROPGA
701     !
702     !                      LHS (similar to A_m)
703           DOUBLE PRECISION A
704     !
705     !                      RHS (similar to b_m)
706           DOUBLE PRECISION B
707     !
708     !                      b_m
709           DOUBLE PRECISION b_m
710     !
711     !                      sum of ep_s
712           DOUBLE PRECISION Eps
713     !
714     !                      Source terms (Surface)
715           DOUBLE PRECISION Sdp, Sdps
716     !
717     !                      Source terms (Volumetric)
718           DOUBLE PRECISION V0, Vmt, Vbf
719     !-----------------------------------------------
720     
721           I = I_OF(IJK)
722           J = J_OF(IJK)
723           KM = KM1(K_OF(IJK))
724           IJKB = BOTTOM_OF(IJK)
725           IJKM = KM_OF(IJK)
726     !
727           A = ZERO
728           B = ZERO
729           EPS = ZERO
730     !
731           DO M = 1, MMAX
732              IF (CLOSE_PACKED(M)) THEN
733                 EPGA = AVG_Z(EP_S(IJKB,M),EP_S(IJK,M),KM)
734     !
735     !         Surface forces
736     !
737     !         Pressure term
738                 SDP = -P_SCALE*EPGA*(P_G(IJK)-P_G(IJKB))*AXY(IJK)
739                 SDPS = -EPGA*(P_S(IJK,M)-P_S(IJKB,M))*AXY(IJK)
740     !
741     !           Shear stress terms
742     !
743     !         Volumetric forces
744                 ROPGA = AVG_Z(ROP_S(IJKB,M),ROP_S(IJK,M),KM)
745     !
746     !         Previous time step
747                 V0 = AVG_Z(ROP_SO(IJKB,M),ROP_SO(IJK,M),KM)*ODT
748     !
749     !         Interphase mass transfer
750                 VMT = AVG_Z(SUM_R_S(IJKB,M),SUM_R_S(IJK,M),KM)
751     !
752     !         Body force
753                 VBF = ROPGA*BFZ_S(IJKB,M)
754     !
755     !         Collect the terms
756                 A = A - ((V0 + ZMAX(VMT))*VOL_W(IJKB))
757                 B=B-(SDP+SDPS+((V0+ZMAX((-VMT)))*W_SO(IJKM,M)+VBF)*VOL_W(IJKB))
758                 EPS = EPS + EPGA
759              ENDIF
760           END DO
761           B_M = -(((-(B - A)/(EPS*AXY(IJK))))-P_STAR(IJK))
762     !
763           RETURN
764           END SUBROUTINE B_M_P_STAR_B
765     
766