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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: CALC_GS_DES1                                            !
4     !  Author: J.Musser                                   Date: 21-NOV-14  !
5     !                                                                      !
6     !  Purpose: This routine is called from the DISCRETE side to calculate !
7     !  the gas-based forces acting on each particle using interpolated     !
8     !  values for gas velocity, gas pressure, and gas volume fraction.     !
9     !                                                                      !
10     !  Notes:                                                              !
11     !                                                                      !
12     !   F_gp is obtained from des_drag_gp subroutine is given as:          !
13     !    F_GP = beta*VOL_P/EP_s where VOL_P is the particle volume.        !
14     !                                                                      !
15     !  The drag force on each particle is equal to:                        !
16     !    D_FORCE = beta*VOL_P/EP_s*(Ug - Us) = F_GP *(Ug - Us)             !
17     !                                                                      !
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
19           SUBROUTINE DRAG_GS_DES1
20     
21     ! Gas phase volume fraction
22           use fldvar, only: EP_G
23     ! Gas phase velocities
24           use fldvar, only: U_G, V_G, W_G
25     ! Size of particle arrays on this processor.
26           use discretelement, only: MAX_PIP
27     ! Flag to use interpolation
28           use particle_filter, only: DES_INTERP_ON
29     ! Interpolation cells and weights
30           use particle_filter, only: FILTER_CELL, FILTER_WEIGHT
31     ! IJK of fluid cell containing particles center
32           use discretelement, only: PIJK
33     ! Flags indicating the state of particle
34           use discretelement, only: PEA
35     ! Drag force on each particle
36           use discretelement, only: F_GP
37     ! Particle velocity
38           use discretelement, only: DES_VEL_NEW
39     ! Total forces acting on particle
40           use discretelement, only: FC
41     ! Gas pressure force by fluid cell
42           use discretelement, only: P_FORCE
43     ! Particle volume.
44           use discretelement, only: PVOL
45     ! Model B momentum equation
46           use run, only: MODEL_B
47     ! Cell-center gas velocities.
48           use tmp_array, only: UGC => ARRAY1
49           use tmp_array, only: VGC => ARRAY2
50           use tmp_array, only: WGC => ARRAY3
51     ! Flag for MPPIC runs.
52           use mfix_pic, only: MPPIC
53     ! Flag to use implicit drag for MPPIC
54           use mfix_pic, only: MPPIC_PDRAG_IMPLICIT
55     ! Fluid grid loop bounds.
56           use compar, only: IJKStart3, IJKEnd3
57     ! Flag for 3D simulatoins.
58           use geometry, only: DO_K
59     ! Function to deterine if a cell contains fluid.
60           use functions, only: FLUID_AT
61     
62     ! Global Parameters:
63     !---------------------------------------------------------------------//
64     ! Double precision values.
65           use param1, only: ZERO
66     
67     ! Lock/Unlock the temp arrays to prevent double usage.
68           use tmp_array, only: LOCK_TMP_ARRAY
69           use tmp_array, only: UNLOCK_TMP_ARRAY
70     
71     
72           IMPLICIT NONE
73     
74     ! Loop counters: Particle, fluid cell, neighbor cells
75           INTEGER :: NP, IJK, LC
76     ! Interpolation weight
77           DOUBLE PRECISION :: WEIGHT
78     ! Interpolated gas phase quanties.
79           DOUBLE PRECISION :: lEPg, VELFP(3), lPF(3)
80     ! Drag force acting on each particle.
81           DOUBLE PRECISION :: D_FORCE(3)
82     ! Flag for Model A momentum equation
83           LOGICAL :: MODEL_A
84     ! Loop bound for filter
85           INTEGER :: LP_BND
86     
87     
88     ! Set flag for Model A momentum equation.
89           MODEL_A = .NOT.MODEL_B
90     ! Loop bounds for interpolation.
91           LP_BND = merge(27,9,DO_K)
92     
93     ! Lock the temp arrays.
94           CALL LOCK_TMP_ARRAY
95     
96     ! Calculate the cell center gas velocities.
97           CALL CALC_CELL_CENTER_GAS_VEL
98     
99     ! Calculate the gas phae forces acting on each particle.
100           DO NP=1,MAX_PIP
101              IF(.NOT.PEA(NP,1)) CYCLE
102              IF(any(PEA(NP,2:4))) CYCLE
103     
104              lEPG = ZERO
105              VELFP = ZERO
106              lPF = ZERO
107     
108     ! Calculate the gas volume fraction, velocity, and pressure force at
109     ! the particle's position.
110              IF(DES_INTERP_ON) THEN
111                 DO LC=1,LP_BND
112                    IJK = FILTER_CELL(LC,NP)
113                    WEIGHT = FILTER_WEIGHT(LC,NP)
114     ! Gas phase volume fraction.
115                    lEPG = lEPG + EP_G(IJK)*WEIGHT
116     ! Gas phase velocity.
117                    VELFP(1) = VELFP(1) + UGC(IJK)*WEIGHT
118                    VELFP(2) = VELFP(2) + VGC(IJK)*WEIGHT
119                    VELFP(3) = VELFP(3) + WGC(IJK)*WEIGHT
120     ! Gas pressure force.
121                    lPF = lPF + P_FORCE(:,IJK)*WEIGHT
122                 ENDDO
123              ELSE
124                 IJK = PIJK(NP,4)
125                 lEPG = EP_G(IJK)
126                 VELFP(1) = UGC(IJK)
127                 VELFP(2) = VGC(IJK)
128                 VELFP(3) = WGC(IJK)
129                 lPF = P_FORCE(:,IJK)
130              ENDIF
131     
132              CALL DES_DRAG_GP(NP, DES_VEL_NEW(:,NP), VELFP, lEPg)
133     
134     ! Calculate the gas-solids drag force on the particle
135              IF(MPPIC .AND. MPPIC_PDRAG_IMPLICIT) THEN
136     ! implicit treatment of the drag term for mppic
137                 D_FORCE = F_GP(NP)*VELFP
138              ELSE
139     ! default case
140                 D_FORCE = F_GP(NP)*(VELFP - DES_VEL_NEW(:,NP))
141              ENDIF
142     
143     ! Update the contact forces (FC) on the particle to include gas
144     ! pressure and gas-solids drag
145              FC(:,NP) = FC(:,NP) + D_FORCE(:)
146     
147              IF(MODEL_A) FC(:,NP) = FC(:,NP) + lPF*PVOL(NP)
148     
149           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
150     
151     ! Unlock the temp arrays.
152           CALL UNLOCK_TMP_ARRAY
153     
154           RETURN
155           END SUBROUTINE DRAG_GS_DES1
156     
157     
158     
159     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
160     !                                                                      !
161     !  Subroutine: DRAG_GS_GAS1                                            !
162     !  Author: J.Musser                                   Date: 21-NOV-14  !
163     !                                                                      !
164     !                                                                      !
165     !  Purpose: This routine is called from the CONTINUUM. It calculates   !
166     !  the scalar cell center drag force acting on the fluid using         !
167     !  interpolated values for the gas velocity and volume fraction. The   !
168     !  The resulting sources are interpolated back to the fluid grid.      !
169     !                                                                      !
170     !  NOTE: The loop over particles includes ghost particles so that MPI  !
171     !  communications are needed to distribute overlapping force between   !
172     !  neighboring grid cells. This is possible because only cells "owned" !
173     !  by the current process will have non-zero weights.                  !
174     !                                                                      !
175     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
176           SUBROUTINE DRAG_GS_GAS1
177     
178     ! Gas phase volume fraction
179           use fldvar, only: EP_G
180     ! Gas phase velocities
181           use fldvar, only: U_G, V_G, W_G
182     ! Size of particle array on this process.
183           use discretelement, only: MAX_PIP
184     ! Flag to use interpolation
185           use particle_filter, only: DES_INTERP_ON
186     ! Interpolation cells and weights
187           use particle_filter, only: FILTER_CELL, FILTER_WEIGHT
188     ! Flags indicating the state of particle
189           use discretelement, only: PEA
190     ! IJK of fluid cell containing particles center
191           use discretelement, only: PIJK
192     ! Drag force on each particle
193           use discretelement, only: F_GP
194     ! Particle velocity
195           use discretelement, only: DES_VEL_NEW
196     ! Contribution to gas momentum equation due to drag
197           use discretelement, only: DRAG_BM
198     ! Scalar cell center total drag force
199           use discretelement, only: F_GDS
200     ! Flag for MPPIC runs
201           use mfix_pic, only: MPPIC
202     ! Statical weight of each MPPIC parcel
203           use mfix_pic, only: DES_STAT_WT
204     ! Volume of scalar cell.
205           use geometry, only: VOL
206     ! Flag for 3D simulatoins.
207           use geometry, only: DO_K
208     ! Cell-center gas velocities.
209           use tmp_array, only: UGC => ARRAY1
210           use tmp_array, only: VGC => ARRAY2
211           use tmp_array, only: WGC => ARRAY3
212     ! Lock/Unlock the temp arrays to prevent double usage.
213           use tmp_array, only: LOCK_TMP_ARRAY
214           use tmp_array, only: UNLOCK_TMP_ARRAY
215     ! MPI wrapper for halo exchange.
216           use sendrecv, only: SEND_RECV
217     
218     
219     ! Global Parameters:
220     !---------------------------------------------------------------------//
221     ! Double precision values.
222           use param1, only: ZERO, ONE
223     
224           IMPLICIT NONE
225     
226     ! Loop counters: Particle, fluid cell, neighbor cells
227           INTEGER :: NP, IJK, LC
228     ! Interpolation weight
229           DOUBLE PRECISION :: WEIGHT
230     ! Interpolated gas phase quanties.
231           DOUBLE PRECISION :: lEPg, VELFP(3)
232     ! Loop bound for filter
233           INTEGER :: LP_BND
234     ! Drag force (intermediate calculation)
235           DOUBLE PRECISION :: lFORCE
236     ! Drag sources for fluid (intermediate calculation)
237           DOUBLE PRECISION :: lDRAG_BM(3)
238     
239     
240     ! Initialize fluid cell values.
241           F_GDS = ZERO
242           DRAG_BM = ZERO
243     
244     ! Loop bounds for interpolation.
245           LP_BND = merge(27,9,DO_K)
246     
247     ! Lock the temp arrays.
248           CALL LOCK_TMP_ARRAY
249     
250     ! Calculate the cell center gas velocities.
251           CALL CALC_CELL_CENTER_GAS_VEL
252     
253     ! Calculate the gas phae forces acting on each particle.
254           DO NP=1,MAX_PIP
255              IF(.NOT.PEA(NP,1)) CYCLE
256     
257     ! The drag force is not calculated on entering or exiting particles
258     ! as their velocities are fixed and may exist in 'non fluid' cells.
259             IF(any(PEA(NP,2:3))) CYCLE
260     
261              lEPG = ZERO
262              VELFP = ZERO
263     
264     ! Calculate the gas volume fraction, velocity, and at the
265     ! particle's position.
266              IF(DES_INTERP_ON) THEN
267                 DO LC=1,LP_BND
268                    IJK = FILTER_CELL(LC,NP)
269                    WEIGHT = FILTER_WEIGHT(LC,NP)
270     ! Gas phase volume fraction.
271                    lEPG = lEPG + EP_G(IJK)*WEIGHT
272     ! Gas phase velocity.
273                    VELFP(1) = VELFP(1) + UGC(IJK)*WEIGHT
274                    VELFP(2) = VELFP(2) + VGC(IJK)*WEIGHT
275                    VELFP(3) = VELFP(3) + WGC(IJK)*WEIGHT
276                 ENDDO
277              ELSE
278                 IJK = PIJK(NP,4)
279                 lEPG = EP_G(IJK)
280                 VELFP(1) = UGC(IJK)
281                 VELFP(2) = VGC(IJK)
282                 VELFP(3) = WGC(IJK)
283              ENDIF
284     
285     ! This avoids FP exceptions for some ghost particles.
286              IF(lEPg == ZERO) lEPG = EP_g(PIJK(NP,4))
287     
288              CALL DES_DRAG_GP(NP, DES_VEL_NEW(:,NP), VELFP, lEPg)
289     
290              lFORCE = F_GP(NP)
291              IF(MPPIC) lFORCE = lFORCE*DES_STAT_WT(NP)
292     
293              lDRAG_BM = lFORCE*DES_VEL_NEW(:,NP)
294     
295              IF(DES_INTERP_ON) THEN
296                 DO LC=1,LP_BND
297                    IJK = FILTER_CELL(LC,NP)
298                    WEIGHT = FILTER_WEIGHT(LC,NP)/VOL(IJK)
299     
300                    DRAG_BM(IJK,:) = DRAG_BM(IJK,:) + lDRAG_BM*WEIGHT
301                    F_GDS(IJK) = F_GDS(IJK) + lFORCE*WEIGHT
302                 ENDDO
303              ELSE
304                 IJK = PIJK(NP,4)
305                 WEIGHT = ONE/VOL(IJK)
306     
307                 DRAG_BM(IJK,:) = DRAG_BM(IJK,:) + lDRAG_BM*WEIGHT
308                 F_GDS(IJK) = F_GDS(IJK) + lFORCE*WEIGHT
309              ENDIF
310     
311           ENDDO
312     
313     ! Unlock the temp arrays.
314           CALL UNLOCK_TMP_ARRAY
315     
316     ! Update the drag force and sources in ghost layers.
317           CALL SEND_RECV(F_GDS, 2)
318           CALL SEND_RECV(DRAG_BM, 2)
319     
320     
321           RETURN
322           END SUBROUTINE DRAG_GS_GAS1
323     
324     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
325     !                                                                      !
326     !  Subroutine: DRAG_GS_EXPLICIT1                                       !
327     !  Author: J.Musser                                   Date: 21-NOV-14  !
328     !                                                                      !
329     !                                                                      !
330     !  Purpose: This routine is called from the CONTINUUM. It calculates   !
331     !  the scalar cell center drag force acting on the fluid using         !
332     !  interpolated values for the gas velocity and volume fraction. The   !
333     !  The resulting sources are interpolated back to the fluid grid.      !
334     !                                                                      !
335     !  NOTE: The loop over particles includes ghost particles so that MPI  !
336     !  communications are needed to distribute overlapping force between   !
337     !  neighboring grid cells. This is possible because only cells "owned" !
338     !  by the current process will have non-zero weights.                  !
339     !                                                                      !
340     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
341           SUBROUTINE DRAG_GS_EXPLICIT1
342     
343     ! Gas phase volume fraction
344           use fldvar, only: EP_G
345     ! Gas phase velocities
346           use fldvar, only: U_G, V_G, W_G
347     ! Size of particle array on this process.
348           use discretelement, only: MAX_PIP
349     ! Flag to use interpolation
350           use particle_filter, only: DES_INTERP_ON
351     ! Interpolation cells and weights
352           use particle_filter, only: FILTER_CELL, FILTER_WEIGHT
353     ! Flags indicating the state of particle
354           use discretelement, only: PEA
355     ! IJK of fluid cell containing particles center
356           use discretelement, only: PIJK
357     ! Drag force on each particle
358           use discretelement, only: F_GP
359     ! Particle velocity
360           use discretelement, only: DES_VEL_NEW
361     ! Particle drag force
362           use discretelement, only: DRAG_FC
363     ! Contribution to gas momentum equation due to drag
364           use discretelement, only: DRAG_BM
365     ! Scalar cell center total drag force
366           use discretelement, only: F_GDS
367     ! Flag for MPPIC runs
368           use mfix_pic, only: MPPIC
369     ! Statical weight of each MPPIC parcel
370           use mfix_pic, only: DES_STAT_WT
371     ! Volume of scalar cell.
372           use geometry, only: VOL
373     ! Flag for 3D simulatoins.
374           use geometry, only: DO_K
375     ! Cell-center gas velocities.
376           use tmp_array, only: UGC => ARRAY1
377           use tmp_array, only: VGC => ARRAY2
378           use tmp_array, only: WGC => ARRAY3
379     ! Lock/Unlock the temp arrays to prevent double usage.
380           use tmp_array, only: LOCK_TMP_ARRAY
381           use tmp_array, only: UNLOCK_TMP_ARRAY
382     ! MPI wrapper for halo exchange.
383           use sendrecv, only: SEND_RECV
384     
385     
386     
387     ! Global Parameters:
388     !---------------------------------------------------------------------//
389     ! Double precision values.
390           use param1, only: ZERO, ONE
391     
392           IMPLICIT NONE
393     
394     ! Loop counters: Particle, fluid cell, neighbor cells
395           INTEGER :: NP, IJK, LC
396     ! Interpolation weight
397           DOUBLE PRECISION :: WEIGHT
398     ! Interpolated gas phase quanties.
399           DOUBLE PRECISION :: lEPg, VELFP(3)
400     ! Loop bound for
401           INTEGER :: LP_BND
402     ! Drag force (intermediate calculation)
403           DOUBLE PRECISION :: lFORCE
404     ! Drag source for fluid (intermediate calculation)
405           DOUBLE PRECISION :: lDRAG_BM(3)
406     
407     
408     ! Initialize fluid cell values.
409           F_GDS = ZERO
410           DRAG_BM = ZERO
411     
412     ! Loop bounds for interpolation.
413           LP_BND = merge(27,9,DO_K)
414     
415     ! Lock the temp arrays.
416           CALL LOCK_TMP_ARRAY
417     
418     ! Calculate the cell center gas velocities.
419           CALL CALC_CELL_CENTER_GAS_VEL
420     
421     ! Calculate the gas phase forces acting on each particle.
422     !---------------------------------------------------------------------//
423     !$omp parallel default(none)                                           &
424     !$omp private(np, lepg, velfp, lc, ijk, weight, lforce, ldrag_bm)      &
425     !$omp shared(max_pip, pea, lp_bnd, drag_fc, f_gp, des_vel_new, ugc,    &
426     !$omp   vgc, wgc, mppic, drag_bm, vol, f_gds, ep_g, filter_weight,     &
427     !$omp   des_stat_wt, filter_cell, pijk, des_interp_on)
428     !$omp do
429           DO NP=1,MAX_PIP
430              IF(.NOT.PEA(NP,1)) CYCLE
431     
432     ! The drag force is not calculated on entering or exiting particles
433     ! as their velocities are fixed and may exist in 'non fluid' cells.
434             IF(any(PEA(NP,2:3))) THEN
435                 DRAG_FC(:,NP) = ZERO
436                 CYCLE
437              ENDIF
438     
439              lEPG = ZERO
440              VELFP = ZERO
441     
442     ! Calculate the gas volume fraction, velocity, and at the
443     ! particle's position.
444              IF(DES_INTERP_ON) THEN
445                 DO LC=1,LP_BND
446                    IJK = FILTER_CELL(LC,NP)
447                    WEIGHT = FILTER_WEIGHT(LC,NP)
448     ! Gas phase volume fraction.
449                    lEPG = lEPG + EP_G(IJK)*WEIGHT
450     ! Gas phase velocity.
451                    VELFP(1) = VELFP(1) + UGC(IJK)*WEIGHT
452                    VELFP(2) = VELFP(2) + VGC(IJK)*WEIGHT
453                    VELFP(3) = VELFP(3) + WGC(IJK)*WEIGHT
454                 ENDDO
455              ELSE
456                 IJK = PIJK(NP,4)
457                 lEPG = EP_G(IJK)
458                 VELFP(1) = UGC(IJK)
459                 VELFP(2) = VGC(IJK)
460                 VELFP(3) = WGC(IJK)
461              ENDIF
462     
463     ! This avoids FP exceptions for some ghost particles.
464              IF(lEPG == ZERO) lEPG = EP_G(PIJK(NP,4))
465     
466              CALL DES_DRAG_GP(NP, DES_VEL_NEW(:,NP), VELFP, lEPg)
467     
468     ! Evaluate the drag force acting on the particle.
469              DRAG_FC(:,NP) = F_GP(NP)*(VELFP - DES_VEL_NEW(:,NP))
470     
471     ! Calculate the force on the fluid.
472              lFORCE = F_GP(NP)
473              IF(MPPIC) lFORCE = lFORCE*DES_STAT_WT(NP)
474     
475              lDRAG_BM(:) = lFORCE*(DES_VEL_NEW(:,NP) - VELFP)
476     
477              IF(DES_INTERP_ON) THEN
478                 DO LC=1,LP_BND
479                    IJK = FILTER_CELL(LC,NP)
480                    WEIGHT = FILTER_WEIGHT(LC,NP)/VOL(IJK)
481                    !$omp atomic
482                    DRAG_BM(IJK,1) = DRAG_BM(IJK,1) + lDRAG_BM(1)*WEIGHT
483                    !$omp atomic
484                    DRAG_BM(IJK,2) = DRAG_BM(IJK,2) + lDRAG_BM(2)*WEIGHT
485                    !$omp atomic
486                    DRAG_BM(IJK,3) = DRAG_BM(IJK,3) + lDRAG_BM(3)*WEIGHT
487                    !$omp atomic
488                    F_GDS(IJK) = F_GDS(IJK) + lFORCE*WEIGHT
489                 ENDDO
490              ELSE
491                 IJK = PIJK(NP,4)
492                 WEIGHT = ONE/VOL(IJK)
493                 !$omp atomic
494                 DRAG_BM(IJK,1) = DRAG_BM(IJK,1) + lDRAG_BM(1)*WEIGHT
495                 !$omp atomic
496                 DRAG_BM(IJK,2) = DRAG_BM(IJK,2) + lDRAG_BM(2)*WEIGHT
497                 !$omp atomic
498                 DRAG_BM(IJK,3) = DRAG_BM(IJK,3) + lDRAG_BM(3)*WEIGHT
499                 !$omp atomic
500                 F_GDS(IJK) = F_GDS(IJK) + lFORCE*WEIGHT
501              ENDIF
502           ENDDO
503     !$omp end do
504     !$omp end parallel
505     
506     ! Unlock the temp arrays.
507           CALL UNLOCK_TMP_ARRAY
508     
509     ! Update the drag force and sources in ghost layers.
510           CALL SEND_RECV(F_GDS, 2)
511           CALL SEND_RECV(DRAG_BM, 2)
512     
513     
514           RETURN
515           END SUBROUTINE DRAG_GS_EXPLICIT1
516     
517     
518     
519     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
520     !                                                                      !
521     !  Subroutine: CALC_CELL_CENTER_GAS_VEL                                !
522     !  Author: J.Musser                                   Date: 07-NOV-14  !
523     !                                                                      !
524     !  Purpose: Calculate the scalar cell center gas velocity. This code   !
525     !  is common to the DEM and GAS calls for non-interpolated drag        !
526     !  routines.                                                           !
527     !                                                                      !
528     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
529           SUBROUTINE CALC_CELL_CENTER_GAS_VEL
530     
531     ! Global Variables:
532     !---------------------------------------------------------------------//
533     ! Fluid velocities.
534           use fldvar, only: U_G, V_G, W_G
535     ! Functions to average momentum to scalar cell center.
536           use fun_avg, only: AVG_X_E, AVG_Y_N, AVG_Z_T
537     ! Flags and correction factors for cut momentum cells.
538           use cutcell, only: CUT_U_TREATMENT_AT, THETA_UE, THETA_UE_BAR
539           use cutcell, only: CUT_V_TREATMENT_AT, THETA_VN, THETA_VN_BAR
540           use cutcell, only: CUT_W_TREATMENT_AT, THETA_WT, THETA_WT_BAR
541     ! Functions to lookup adjacent cells by index.
542           use functions, only: IM_OF, JM_OF, KM_OF
543           use indices, only: I_OF
544     ! Fluid grid loop bounds.
545           use compar, only: IJKStart3, IJKEnd3
546     ! Flag for 3D simulatoins.
547           use geometry, only: DO_K
548     ! Function to deterine if a cell contains fluid.
549           use functions, only: FLUID_AT
550     
551           use tmp_array, only: UGC => ARRAY1
552           use tmp_array, only: VGC => ARRAY2
553           use tmp_array, only: WGC => ARRAY3
554     
555     ! Global Parameters:
556     !---------------------------------------------------------------------//
557     ! Double precision parameters.
558           use param1, only: ZERO
559     
560           IMPLICIT NONE
561     
562     ! Local variables:
563     !---------------------------------------------------------------------//
564     ! Indices of adjacent cells
565           INTEGER :: IJK, IMJK, IJMK, IJKM
566     
567     ! Calculate the cell center gas velocity components.
568           DO IJK=IJKSTART3, IJKEND3
569              IF(FLUID_AT(IJK)) THEN
570                 IMJK = IM_OF(IJK)
571                 IF(CUT_U_TREATMENT_AT(IMJK)) THEN
572                    UGC(IJK) = (THETA_UE_BAR(IMJK)*U_G(IMJK) +              &
573                       THETA_UE(IMJK)*U_G(IJK))
574                 ELSE
575                    UGC(IJK) = AVG_X_E(U_G(IMJK),U_G(IJK),I_OF(IJK))
576                 ENDIF
577     
578                 IJMK = JM_OF(IJK)
579                 IF(CUT_V_TREATMENT_AT(IJMK)) THEN
580                    VGC(IJK) = (THETA_VN_BAR(IJMK)*V_G(IJMK) +              &
581                       THETA_VN(IJMK)*V_G(IJK))
582                 ELSE
583                    VGC(IJK) = AVG_Y_N(V_G(IJMK),V_G(IJK))
584                 ENDIF
585     
586                 IF(DO_K) THEN
587                    IJKM = KM_OF(IJK)
588                    IF(CUT_W_TREATMENT_AT(IJKM)) THEN
589                       WGC(IJK) = (THETA_WT_BAR(IJKM)*W_G(IJKM) +           &
590                          THETA_WT(IJKM)* W_G(IJK))
591                    ELSE
592                       WGC(IJK) = AVG_Z_T(W_G(IJKM),W_G(IJK))
593                    ENDIF
594                 ELSE
595                    WGC(IJK) = ZERO
596                 ENDIF
597              ELSE
598                 UGC(IJK) = ZERO
599                 VGC(IJK) = ZERO
600                 WGC(IJK) = ZERO
601              ENDIF
602           ENDDO
603     
604     
605           RETURN
606           END SUBROUTINE CALC_CELL_CENTER_GAS_VEL
607     
608