File: RELATIVE:/../../../mfix.git/model/calc_mflux.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_MFLUX                                              C
4     !  Purpose: Calculate the convection fluxes. Master routine.           C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 31-MAY-05  C
7     !                                                                      C
8     !                                                                      C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10           SUBROUTINE CALC_MFLUX()
11     
12     ! Modules
13     !---------------------------------------------------------------------//
14           USE fldvar, only: u_g, v_g, w_g
15           USE fldvar, only: u_s, v_s, w_s
16           USE mflux, only: rop_ge, rop_gn, rop_gt
17           USE mflux, only: rop_se, rop_sn, rop_st
18           USE mflux, only: flux_ge, flux_gn, flux_gt
19           USE mflux, only: flux_se, flux_sn, flux_st
20           USE mflux, only: flux_gse, flux_gsn, flux_gst
21           USE mflux, only: flux_sse, flux_ssn, flux_sst
22           USE physprop, only: mmax
23           USE run, only: added_mass, m_am
24           IMPLICIT NONE
25     
26     ! Local variables
27     !---------------------------------------------------------------------//
28     ! solids phase index
29           INTEGER ::  M
30     !---------------------------------------------------------------------//
31     
32           IF(.NOT.Added_Mass) THEN
33              CALL CALC_MFLUX0 (U_g, V_g, W_g, ROP_gE, ROP_gN, ROP_gT, &
34                                Flux_gE, Flux_gN, Flux_gT)
35              DO M = 1, MMAX
36                CALL CALC_MFLUX0 (U_s(1,M), V_s(1,M), W_s(1,M), &
37                                  ROP_sE(1,M), ROP_sN(1,M), ROP_sT(1,M), &
38                                  Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M))
39              ENDDO
40     
41           ELSE
42     ! New fluxes are defined for gas based on added mass
43              CALL CALC_MFLUX_AM (U_g, V_g, W_g, ROP_gE, ROP_gN, ROP_gT, &
44                                ROP_sE(1,M_AM), ROP_sN(1,M_AM), ROP_sT(1,M_AM), &
45                                Flux_gE, Flux_gN, Flux_gT, M_AM)
46     
47              CALL CALC_MFLUX0 (U_g, V_g, W_g, ROP_gE, ROP_gN, ROP_gT, &
48                                Flux_gSE, Flux_gSN, Flux_gST)
49              DO M = 1, MMAX
50     ! New fluxes are defined for M = M_am where virtual mass force is added.
51                IF(M==M_AM) THEN
52                   CALL CALC_MFLUX_AM (U_s(1,M), V_s(1,M), W_s(1,M), &
53                                  ROP_sE(1,M), ROP_sN(1,M), ROP_sT(1,M), &
54                                  ROP_gE, ROP_gN, ROP_gT, &
55                                  Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M),&
56                                  M_AM)
57                   CALL CALC_MFLUX0 (U_s(1,M), V_s(1,M), W_s(1,M), &
58                                  ROP_sE(1,M), ROP_sN(1,M), ROP_sT(1,M), &
59                                  Flux_sSE, Flux_sSN, Flux_sST)
60                ELSE
61                   CALL CALC_MFLUX0 (U_s(1,M), V_s(1,M), W_s(1,M), &
62                                  ROP_sE(1,M), ROP_sN(1,M), ROP_sT(1,M), &
63                                  Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M))
64                ENDIF
65              ENDDO
66           ENDIF
67     
68           RETURN
69           END SUBROUTINE CALC_MFLUX
70     
71     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
72     !                                                                      C
73     !  Subroutine:: CALC_MFLUX0                                            C
74     !  Purpose: Calculate the convection fluxes.                           C
75     !                                                                      C
76     !  Author: M. Syamlal                                 Date: 31-MAY-05  C
77     !                                                                      C
78     !                                                                      C
79     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
80           SUBROUTINE CALC_MFLUX0(U, V, W, ROP_E, ROP_N, ROP_T,&
81                                  Flux_E, Flux_N, Flux_T)
82     
83     ! Modules
84     !---------------------------------------------------------------------//
85           USE compar, only: ijkstart3, ijkend3
86           USE functions, only: fluid_at
87           USE functions, only: im_of, jm_of, km_of
88           USE geometry, only: do_k
89           USE geometry, only: ayz, axz, axy
90           USE param, only: dimension_3
91           IMPLICIT NONE
92     
93     ! Dummy arguments
94     !---------------------------------------------------------------------//
95     ! Velocity components
96           DOUBLE PRECISION, INTENT(IN) :: U(DIMENSION_3)
97           DOUBLE PRECISION, INTENT(IN) :: V(DIMENSION_3)
98           DOUBLE PRECISION, INTENT(IN) :: W(DIMENSION_3)
99     ! Face value of density (for calculating convective fluxes)
100           DOUBLE PRECISION, INTENT(IN) :: ROP_E(DIMENSION_3)
101           DOUBLE PRECISION, INTENT(IN) :: ROP_N(DIMENSION_3)
102           DOUBLE PRECISION, INTENT(IN) :: ROP_T(DIMENSION_3)
103     ! Convective mass fluxes
104           DOUBLE PRECISION, INTENT(OUT) :: Flux_E(DIMENSION_3)
105           DOUBLE PRECISION, INTENT(OUT) :: Flux_N(DIMENSION_3)
106           DOUBLE PRECISION, INTENT(OUT) :: Flux_T(DIMENSION_3)
107     
108     ! Local variables
109     !---------------------------------------------------------------------//
110     ! Indices
111           INTEGER :: IJK, IMJK, IJMK, IJKM
112     !---------------------------------------------------------------------//
113     
114     !$omp  parallel do default(none) private( IJK, IMJK, IJMK, IJKM) &
115     !$omp              shared(ijkstart3, ijkend3, flux_e, flux_n, flux_t,&
116     !$omp                     rop_e, rop_n, rop_t, axz, axy, ayz, u, v, w, do_k)
117     
118           DO IJK = ijkstart3, ijkend3
119              IF (FLUID_AT(IJK)) THEN
120     
121                 IMJK = IM_OF(IJK)
122                 IJMK = JM_OF(IJK)
123     
124     ! East face (i+1/2, j, k)
125                 Flux_E(IJK) = ROP_E(IJK)*AYZ(IJK)*U(IJK)
126     ! West face (i-1/2, j, k)
127                 IF (.NOT.FLUID_AT(IMJK)) THEN
128                    Flux_E(IMJK) = ROP_E(IMJK)*AYZ(IMJK)*U(IMJK)
129                 ENDIF
130     
131     ! North face (i, j+1/2, k)
132                 Flux_N(IJK) = ROP_N(IJK)*AXZ(IJK)*V(IJK)
133     ! South face (i, j-1/2, k)
134                 IF (.NOT.FLUID_AT(IJMK)) THEN
135                   Flux_N(IJMK) = ROP_N(IJMK)*AXZ(IJMK)*V(IJMK)
136                 ENDIF
137     
138                 IF (DO_K) THEN
139                    IJKM = KM_OF(IJK)
140     ! Top face (i, j, k+1/2)
141                    Flux_T(IJK) = ROP_T(IJK)*AXY(IJK)*W(IJK)
142     ! Bottom face (i, j, k-1/2)
143                    IF (.NOT.FLUID_AT(IJKM)) THEN
144                      Flux_T(IJKM) = ROP_T(IJKM)*AXY(IJKM)*W(IJKM)
145                    ENDIF
146                 ENDIF   ! end if do_k
147              ENDIF   ! end if fluid_at
148           ENDDO   ! end do ijk
149     
150           RETURN
151           END SUBROUTINE CALC_MFLUX0
152     
153     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
154     !                                                                      C
155     !                                                                      C
156     !  Subroutine: CALC_MFLUX_AM                                           C
157     !  Purpose: Calculate the convection fluxes for case with added_mass.  C
158     !                                                                      C
159     !  Author: S. Benyahia                                 Date:           C
160     !                                                                      C
161     !                                                                      C
162     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
163           SUBROUTINE CALC_MFLUX_AM(U, V, W, ROP_E, ROP_N, ROP_T, &
164                                  ROPa_E, ROPa_N, ROPa_T, Flux_E, &
165                                  Flux_N, Flux_T, M_AM)
166     
167     ! Modules
168     !---------------------------------------------------------------------//
169           USE compar, only: ijkstart3, ijkend3
170           USE functions, only: fluid_at
171           USE functions, only: im_of, jm_of, km_of
172           USE fldvar, only: ro_s
173           USE geometry, only: do_k
174           USE geometry, only: ayz, axz, axy
175           USE param, only: dimension_3
176           USE param1, only: one
177           USE physprop, only: cv
178           IMPLICIT NONE
179     
180     ! Dummy arguments
181     !---------------------------------------------------------------------//
182     ! Velocity components
183           DOUBLE PRECISION, INTENT(IN) :: U(DIMENSION_3)
184           DOUBLE PRECISION, INTENT(IN) :: V(DIMENSION_3)
185           DOUBLE PRECISION, INTENT(IN) :: W(DIMENSION_3)
186     ! Face value of density (for calculating convective fluxes)
187           DOUBLE PRECISION, INTENT(IN) :: ROP_E(DIMENSION_3)
188           DOUBLE PRECISION, INTENT(IN) :: ROP_N(DIMENSION_3)
189           DOUBLE PRECISION, INTENT(IN) :: ROP_T(DIMENSION_3)
190     ! Face value of density of added mass phase
191           DOUBLE PRECISION, INTENT(IN) :: ROPa_E(DIMENSION_3)
192           DOUBLE PRECISION, INTENT(IN) :: ROPa_N(DIMENSION_3)
193           DOUBLE PRECISION, INTENT(IN) :: ROPa_T(DIMENSION_3)
194     ! Convective mass fluxes
195           DOUBLE PRECISION, INTENT(OUT) :: Flux_E(DIMENSION_3)
196           DOUBLE PRECISION, INTENT(OUT) :: Flux_N(DIMENSION_3)
197           DOUBLE PRECISION, INTENT(OUT) :: Flux_T(DIMENSION_3)
198     ! phase index for phase of added_mass
199           INTEGER, INTENT(IN) :: M_AM
200     
201     ! Local variables
202     !---------------------------------------------------------------------//
203     ! Indices
204           INTEGER :: IJK, IMJK, IJMK, IJKM
205     !---------------------------------------------------------------------//
206     
207     !!!$omp  parallel do private( IJK, IMJK, IJMK, IJKM) &
208     !!!$omp&  schedule(static)
209           DO IJK = ijkstart3, ijkend3
210     
211              IF (FLUID_AT(IJK)) THEN
212     
213                 IMJK = IM_OF(IJK)
214                 IJMK = JM_OF(IJK)
215     
216     ! East face (i+1/2, j, k)
217                 Flux_E(IJK) = ROP_E(IJK)* &
218                               (ONE + Cv*ROPa_E(IJK)/RO_S(IJK,M_AM))*&
219                               AYZ(IJK)*U(IJK)
220     ! West face (i-1/2, j, k)
221                 IF (.NOT.FLUID_AT(IMJK)) THEN
222                    Flux_E(IMJK) = ROP_E(IMJK)*&
223                                   (ONE + Cv*ROPa_E(IMJK)/RO_S(IJK,M_AM))*&
224                                   AYZ(IMJK)*U(IMJK)
225                 ENDIF
226     
227     ! North face (i, j+1/2, k)
228                 Flux_N(IJK) = ROP_N(IJK)*&
229                               (ONE + Cv*ROPa_N(IJK)/RO_S(IJK,M_AM))*&
230                               AXZ(IJK)*V(IJK)
231     ! South face (i, j-1/2, k)
232                 IF (.NOT.FLUID_AT(IJMK)) THEN
233                   Flux_N(IJMK) = ROP_N(IJMK)*&
234                                  (ONE + Cv*ROPa_N(IJMK)/RO_S(IJK,M_AM))*&
235                                  AXZ(IJMK)*V(IJMK)
236                 ENDIF
237     
238                 IF (DO_K) THEN
239                    IJKM = KM_OF(IJK)
240     
241     ! Top face (i, j, k+1/2)
242                    Flux_T(IJK) = ROP_T(IJK)*&
243                                  (ONE + Cv*ROPa_T(IJK)/RO_S(IJK,M_AM))*&
244                                  AXY(IJK)*W(IJK)
245     ! Bottom face (i, j, k-1/2)
246                    IF (.NOT.FLUID_AT(IJKM)) THEN
247                      Flux_T(IJKM) = ROP_T(IJKM)*&
248                                     (ONE + Cv*ROPa_T(IJKM)/RO_S(IJK,M_AM))*&
249                                     AXY(IJKM)*W(IJKM)
250                    ENDIF
251                 ENDIF   ! end if do_k
252              ENDIF   ! end if fluid_at
253           ENDDO    ! end do ijk
254     
255           RETURN
256           END SUBROUTINE CALC_MFLUX_AM
257