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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CONV_ROP                                                C
4     !  Purpose: Calculate the face value of density used for calculating   C
5     !  convection fluxes. Master routine.                                  C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 31-MAY-05  C
8     !                                                                      C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11           SUBROUTINE CONV_ROP()
12     
13     ! Modules
14     !---------------------------------------------------------------------//
15           USE fldvar, only: rop_g, u_g, v_g, w_g
16           USE fldvar, only: rop_s, u_s, v_s, w_s
17           USE mflux, only: rop_ge, rop_gn, rop_gt
18           USE mflux, only: rop_se, rop_sn, rop_st
19           USE physprop, only: mmax
20           USE run, only: discretize
21           IMPLICIT NONE
22     
23     !---------------------------------------------------------------------//
24     ! Local variables
25     !---------------------------------------------------------------------//
26     ! solids phase index
27           INTEGER :: M
28     !---------------------------------------------------------------------//
29     
30           IF (DISCRETIZE(1) == 0) THEN               ! 0 & 1 => first order upwinding
31              CALL CONV_ROP0 (ROP_g, U_g, V_g, W_g, &
32                              ROP_gE, ROP_gN, ROP_gT)
33           ELSE
34              CALL CONV_ROP1 (DISCRETIZE(1), ROP_g, U_g, V_g, W_g, &
35                              ROP_gE, ROP_gN, ROP_gT)
36           ENDIF
37     
38           IF (DISCRETIZE(2) == 0) THEN               ! 0 & 1 => first order upwinding
39              DO M = 1, MMAX
40                 CALL CONV_ROP0 (ROP_s(1,M), U_s(1,M), V_s(1,M), W_s(1,M), &
41                                 ROP_sE(1,M), ROP_sN(1,M), ROP_sT(1,M))
42             ENDDO
43           ELSE
44              DO M = 1, MMAX
45                 CALL CONV_ROP1 (DISCRETIZE(2), ROP_s(1,M), &
46                                 U_s(1,M), V_s(1,M), W_s(1,M), &
47                                 ROP_sE(1,M), ROP_sN(1,M), ROP_sT(1,M))
48              ENDDO
49           ENDIF
50     
51           RETURN
52           END SUBROUTINE CONV_ROP
53     
54     
55     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
56     !                                                                      C
57     !  Subroutine: CONV_ROP                                                C
58     !  Purpose: Calculate the face value of density used for calculating   C
59     !  convection fluxes. FOU routine.                                     C
60     !                                                                      C
61     !  Author: M. Syamlal                                 Date: 31-MAY-05  C
62     !                                                                      C
63     !                                                                      C
64     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
65           SUBROUTINE CONV_ROP0(ROP, U, V, W, ROP_E, ROP_N, ROP_T)
66     
67     ! Modules
68     !---------------------------------------------------------------------//
69           USE compar, only: ijkstart3, ijkend3
70           USE functions, only: fluid_at
71           USE functions, only: east_of, north_of, top_of
72           USE functions, only: west_of, south_of, bottom_of
73           USE functions, only: im_of, jm_of, km_of
74           USE geometry, only: do_k
75           USE param, only: dimension_3
76           USE param1, only: zero
77           IMPLICIT NONE
78     
79     ! Dummy arguments
80     !---------------------------------------------------------------------//
81     ! macroscopic density (rho_prime)
82           DOUBLE PRECISION, INTENT(IN) :: ROP(DIMENSION_3)
83     ! Velocity components
84           DOUBLE PRECISION, INTENT(IN) :: U(DIMENSION_3)
85           DOUBLE PRECISION, INTENT(IN) :: V(DIMENSION_3)
86           DOUBLE PRECISION, INTENT(IN) :: W(DIMENSION_3)
87     ! Face value of density (for calculating convective fluxes)
88           DOUBLE PRECISION, INTENT(OUT) :: ROP_E(DIMENSION_3)
89           DOUBLE PRECISION, INTENT(OUT) :: ROP_N(DIMENSION_3)
90           DOUBLE PRECISION, INTENT(OUT) :: ROP_T(DIMENSION_3)
91     
92     ! Local variables
93     !---------------------------------------------------------------------//
94     ! Indices
95           INTEGER :: IJK
96           INTEGER :: IJKE, IJKN, IJKT
97           INTEGER :: IJKW, IJKS, IJKB
98           INTEGER :: IMJK, IJMK, IJKM
99     !---------------------------------------------------------------------//
100     
101     
102     !$omp  parallel do default(none) &
103     !$omp              private(IJK, IJKE, IJKN, IJKT, IJKW, &
104     !$omp                      IJKS, IJKB, IMJK, IJMK, IJKM) &
105     !$omp              shared(ijkstart3, ijkend3, u, v, w, do_k, rop, &
106     !$omp                     rop_e, rop_n, rop_t)
107           DO IJK = ijkstart3, ijkend3
108     
109              IF (FLUID_AT(IJK)) THEN
110                 IJKE = EAST_OF(IJK)
111                 IJKN = NORTH_OF(IJK)
112                 IJKT = TOP_OF(IJK)
113     
114                 IMJK = IM_OF(IJK)
115                 IJMK = JM_OF(IJK)
116     
117     ! East face (i+1/2, j, k)
118                 IF (U(IJK) >= ZERO) THEN
119                    ROP_E(IJK) = ROP(IJK)
120                 ELSE
121                    ROP_E(IJK) = ROP(IJKE)
122                 ENDIF
123     ! West face (i-1/2, j, k)
124                 IF (.NOT.FLUID_AT(IMJK)) THEN
125                    IJKW = WEST_OF(IJK)
126                    IF (U(IMJK) >= ZERO) THEN
127                       ROP_E(IMJK) = ROP(IJKW)
128                    ELSE
129                       ROP_E(IMJK) = ROP(IJK)
130                    ENDIF
131                 ENDIF
132     
133     
134     ! North face (i, j+1/2, k)
135                 IF (V(IJK) >= ZERO) THEN
136                    ROP_N(IJK) = ROP(IJK)
137                 ELSE
138                    ROP_N(IJK) = ROP(IJKN)
139                 ENDIF
140     ! South face (i, j-1/2, k)
141                 IF (.NOT.FLUID_AT(IJMK)) THEN
142                    IJKS = SOUTH_OF(IJK)
143                    IF (V(IJMK) >= ZERO) THEN
144                      ROP_N(IJMK) = ROP(IJKS)
145                    ELSE
146                      ROP_N(IJMK) = ROP(IJK)
147                    ENDIF
148                 ENDIF
149     
150     
151                 IF (DO_K) THEN
152                    IJKM = KM_OF(IJK)
153     ! Top face (i, j, k+1/2)
154                    IF (W(IJK) >= ZERO) THEN
155                       ROP_T(IJK) = ROP(IJK)
156                    ELSE
157                       ROP_T(IJK) = ROP(IJKT)
158                    ENDIF
159     ! Bottom face (i, j, k-1/2)
160                    IF (.NOT.FLUID_AT(IJKM)) THEN
161                       IJKB = BOTTOM_OF(IJK)
162                       IF (W(IJKM) >= ZERO) THEN
163                          ROP_T(IJKM) = ROP(IJKB)
164                       ELSE
165                          ROP_T(IJKM) = ROP(IJK)
166                       ENDIF
167                    ENDIF
168                 ENDIF   ! end if do_k
169     
170              ENDIF   ! end if fluid_at
171           ENDDO    ! end do ijk
172     
173           RETURN
174           END SUBROUTINE CONV_ROP0
175     
176     
177     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
178     !                                                                      C
179     !  Subroutine: CONV_ROP1                                               C
180     !  Purpose: Calculate the face value of density used for calculating   C
181     !  convection fluxes. HR routine.  Here interpolate the face value of  C
182     !  density.                                                            C
183     !                                                                      C
184     !  Author: M. Syamlal                                 Date: 31-MAY-05  C
185     !                                                                      C
186     !                                                                      C
187     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
188           SUBROUTINE CONV_ROP1(DISC, ROP, U, V, W, ROP_E, ROP_N, ROP_T)
189     
190     ! Modules
191     !---------------------------------------------------------------------//
192           USE compar, only: ijkstart3, ijkend3
193           USE functions, only: fluid_at
194           USE functions, only: east_of, north_of, top_of
195           USE functions, only: west_of, south_of, bottom_of
196           USE functions, only: im_of, jm_of, km_of
197           USE geometry, only: do_k
198           USE param, only: dimension_3
199           USE param1, only: one
200           USE xsi, only: calc_xsi
201           IMPLICIT NONE
202     
203     ! Dummy arguments
204     !---------------------------------------------------------------------//
205     ! Discretization scheme
206           INTEGER, INTENT(IN) :: DISC
207     ! macroscopic density (rho_prime)
208           DOUBLE PRECISION, INTENT(IN) :: ROP(DIMENSION_3)
209     ! Velocity components
210           DOUBLE PRECISION, INTENT(IN) :: U(DIMENSION_3)
211           DOUBLE PRECISION, INTENT(IN) :: V(DIMENSION_3)
212           DOUBLE PRECISION, INTENT(IN) :: W(DIMENSION_3)
213     ! Face value of density (for calculating convective fluxes)
214           DOUBLE PRECISION, INTENT(OUT) :: ROP_E(DIMENSION_3)
215           DOUBLE PRECISION, INTENT(OUT) :: ROP_N(DIMENSION_3)
216           DOUBLE PRECISION, INTENT(OUT) :: ROP_T(DIMENSION_3)
217     !
218     ! Local variables
219     !---------------------------------------------------------------------//
220           INTEGER :: IJK, IJKE, IJKN, IJKT
221           INTEGER :: IJKW, IJKS, IJKB, IMJK, IJMK, IJKM
222           Integer :: incr
223     
224           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
225     
226     !---------------------------------------------------------------------//
227     
228     ! Calculate factors
229           incr=0
230           CALL CALC_XSI (DISC, ROP, U, V, W, XSI_E, XSI_N, XSI_T, incr)
231     
232     !!!$omp  parallel do private(IJK, IJKE, IJKN, IJKT, IJKW, IJKS, IJKB, &
233     !!!$omp                      IMJK, IJMK, IJKM) &
234     !!!$omp&  schedule(static)
235           DO IJK = ijkstart3, ijkend3
236     
237              IF (FLUID_AT(IJK)) THEN
238                 IJKE = EAST_OF(IJK)
239                 IJKN = NORTH_OF(IJK)
240                 IJKT = TOP_OF(IJK)
241     
242                 IMJK = IM_OF(IJK)
243                 IJMK = JM_OF(IJK)
244     
245     ! East face (i+1/2, j, k)
246                 ROP_E(IJK) = ((ONE-XSI_E(IJK))*ROP(IJK)+&
247                              XSI_E(IJK)*ROP(IJKE))
248     ! West face (i-1/2, j, k)
249                 IF (.NOT.FLUID_AT(IMJK)) THEN
250                    IJKW = WEST_OF(IJK)
251                    ROP_E(IMJK) = ((ONE - XSI_E(IMJK))*ROP(IJKW)+&
252                                  XSI_E(IMJK)*ROP(IJK))
253                 ENDIF
254     
255     
256     ! North face (i, j+1/2, k)
257                 ROP_N(IJK) = ((ONE-XSI_N(IJK))*ROP(IJK)+&
258                              XSI_N(IJK)*ROP(IJKN))
259     ! South face (i, j-1/2, k)
260                 IF (.NOT.FLUID_AT(IJMK)) THEN
261                    IJKS = SOUTH_OF(IJK)
262                    ROP_N(IJMK) = ((ONE - XSI_N(IJMK))*ROP(IJKS)+&
263                                  XSI_N(IJMK)*ROP(IJK))
264                 ENDIF
265     
266     
267                 IF (DO_K) THEN
268                    IJKM = KM_OF(IJK)
269     
270     ! Top face (i, j, k+1/2)
271                    ROP_T(IJK) = ((ONE - XSI_T(IJK))*ROP(IJK)+&
272                                 XSI_T(IJK)*ROP(IJKT))
273     ! Bottom face (i, j, k-1/2)
274                    IF (.NOT.FLUID_AT(IJKM)) THEN
275                       IJKB = BOTTOM_OF(IJK)
276                       ROP_T(IJKM) = ((ONE - XSI_T(IJKM))*ROP(IJKB)+&
277                                     XSI_T(IJKM)*ROP(IJK))
278                    ENDIF
279                 ENDIF   ! end if do_k
280     
281              ENDIF   ! end if fluid_at
282           ENDDO    ! end do ijk
283     
284           RETURN
285           END SUBROUTINE CONV_ROP1
286