File: N:\mfix\model\qmomk\qmomk_time_march.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                         C
3     !     Module name: QMOMK_TIME_MARCH                                       C
4     !     Purpose: Called in time_march.f to do QMOMK calculations            C
5     !                                                                         C
6     !     Author: Alberto Passalacqua                        Date:            C
7     !     Reviewer:                                          Date:            C
8     !                                                                         C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10     !
11     SUBROUTINE QMOMK_TIME_MARCH
12     
13       USE param
14       USE param1
15       USE constant
16       USE run
17       USE output
18       USE physprop
19       USE fldvar
20       USE geometry
21       USE cont
22       USE tau_g
23       USE tau_s
24       USE visc_g
25       USE visc_s
26       USE funits
27       USE vshear
28       USE scalars
29       USE drag
30       USE rxns
31       USE compar
32       USE time_cpu
33       USE is
34       USE indices
35       USE sendrecv
36       USE qmom_kinetic_equation
37       USE qmomk_fluxes
38       USE qmomk_quadrature
39       USE qmomk_collision
40       USE qmomk_parameters
41       USE drag
42       USE ur_facs
43       USE fun_avg
44       USE functions
45     
46       IMPLICIT NONE
47     
48       INTEGER :: I,J,K,M,M2,IN
49       INTEGER :: IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
50       DOUBLE PRECISION :: Vmax, min_space_delta
51       DOUBLE PRECISION :: vrel, Rep, CD, beta_drag, min_tau_drag
52       DOUBLE PRECISION :: UGC, VGC, WGC, mom0, QMOMK_TCOL, drag_exp, QMOMK_OMEGA
53       DOUBLE PRECISION QMOMK_TIME, QMOMK_DT_TMP, TMP_DTS
54       INTEGER :: TIME_FACTOR, TIME_I
55       LOGICAL, SAVE ::  FIRST_PASS = .TRUE.
56     
57       DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Nlminus, Nlplus, Nrminus, Nrplus
58       DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Ulminus, Ulplus, Urminus, Urplus
59       DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Vlminus, Vlplus, Vrminus, Vrplus
60       DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Wlminus, Wlplus, Wrminus, Wrplus
61       DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: F_x_left, F_x_right, F_y_left
62       DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: F_y_right, F_z_left, F_z_right
63     
64       IF (FIRST_PASS) THEN
65          FIRST_PASS = .FALSE.
66          ! Setting initial guess for drag time
67          QMOMK_TIME = 0.D0
68          QMOMK_TAU_DRAG = 1.0D-5
69     
70          ! A.P. For restart runs, IC's are already set and would be overwritten
71          IF (RUN_TYPE == 'NEW') THEN
72            CALL QMOMK_INITIAL_CONDITIONS
73            CALL QMOMK_INIT_BC
74          ELSE
75            CALL QMOMK_INIT_BC
76          ENDIF
77     
78          QMOMK_M1 = QMOMK_M0
79          QMOMK_N1 = QMOMK_N0
80          QMOMK_U1 = QMOMK_U0
81          QMOMK_V1 = QMOMK_V0
82          QMOMK_W1 = QMOMK_W0
83       END IF
84     
85       ! Finding initial QMOMK time step
86       ! CFL condition
87       Vmax = 0.d0
88       DO M = 1, MMAX
89         DO IJK = ijkstart3, ijkend3
90             DO I = 1, QMOMK_NN
91               IF (Vmax < MAX(ABS(QMOMK_U0(I,IJK,M)), ABS(QMOMK_V0(I,IJK,M)), ABS(QMOMK_W0(I,IJK,M)))) THEN
92                   Vmax = MAX(ABS(QMOMK_U0(I,IJK,M)), ABS(QMOMK_V0(I,IJK,M)), ABS(QMOMK_W0(I,IJK,M)))
93               END IF
94             END DO
95             CALL COMPUTE_COLLISION_TIME(QMOMK_M0(:,IJK,M), D_p0(M), THETA_M(IJK,M), QMOMK_COLLISION_TIME(IJK,M), DT)
96         END DO
97       END DO
98     
99       min_space_delta = MIN (MINVAL(DX), MINVAL(DY), MINVAL(DZ))
100     
101       QMOMK_DT = QMOMK_CFL*min_space_delta/Vmax
102     
103       ! Drag time and collision time
104       min_tau_drag = 1.D0
105       QMOMK_TCOL = 1.D10
106       DO IJK = ijkstart3, ijkend3
107         IF (FLUID_AT(IJK)) THEN
108             IF (min_tau_drag >  MINVAL(QMOMK_TAU_DRAG(:,IJK,:))) THEN
109               min_tau_drag = MINVAL(QMOMK_TAU_DRAG(:,IJK,:))
110             END IF
111             IF (QMOMK_TCOL > MINVAL(QMOMK_COLLISION_TIME(IJK,:))) THEN
112                QMOMK_TCOL = MINVAL(QMOMK_COLLISION_TIME(IJK,:))
113             END IF
114         END IF
115       END DO
116     
117       IF ( QMOMK_COLLISIONS /= 'NONE') THEN
118          QMOMK_DT = MIN(QMOMK_TCOL,QMOMK_DT)
119       END IF
120     
121       QMOMK_DT = MIN(QMOMK_DT,min_tau_drag/10.D0)
122     
123       ! Preparing for QMOMK internal time stepping, if required
124       QMOMK_TIME = TIME
125       IF (DT >= QMOMK_DT) THEN
126          TIME_FACTOR = CEILING(real(DT/QMOMK_DT)) + 1
127          QMOMK_DT_TMP = ZERO
128       ELSE
129          TIME_FACTOR = 1
130          QMOMK_DT_TMP = QMOMK_DT
131          QMOMK_DT = DT
132       END IF
133     
134       TMP_DTS = ZERO
135     
136       DO TIME_I = 1, TIME_FACTOR  ! Starting QMOMK internal time stepping
137     
138          ! Exit if time is greater than flow time
139          IF (QMOMK_TIME .GE. (TIME + DT)) EXIT
140     
141          IF ((QMOMK_TIME + QMOMK_DT) .GT. (TIME + DT)) THEN
142           TMP_DTS = QMOMK_DT
143           QMOMK_DT = TIME + DT - QMOMK_TIME
144          END IF
145     
146          PRINT *,'Flow time step: ',DT
147          PRINT *,'QMOM DT = ',QMOMK_DT
148          PRINT *,'QMOMK internal time step: ',TIME_I
149          PRINT *,'Time factor', TIME_FACTOR
150     
151          QMOMK_OMEGA = (1.D0 + C_e)/2.D0
152     
153          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
154          ! First step of Runge-Kutta time integration !
155          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
156          ! Calculating kinetic based fluxes (Half time step)
157          DO M = 1, MMAX
158             DO IJK = ijkstart3, ijkend3
159                IF (FLUID_AT(IJK)) THEN
160                   I = I_OF(IJK)
161                   J = J_OF(IJK)
162                   K = K_OF(IJK)
163     
164                   IMJK = IM_OF(IJK)
165                   IJMK = JM_OF(IJK)
166                   IJKM = KM_OF(IJK)
167     
168                   IPJK = IP_OF(IJK)
169                   IJPK = JP_OF(IJK)
170                   IJKP = KP_OF(IJK)
171     
172                   ! 2D Case
173                   IF (NO_K) THEN
174                      ! x - direction
175                      Nlminus = QMOMK_N0(:,IMJK,M)
176                      Ulminus = QMOMK_U0(:,IMJK,M)
177                      Vlminus = QMOMK_V0(:,IMJK,M)
178                      Wlminus = QMOMK_W0(:,IMJK,M)
179     
180                      Nlplus  = QMOMK_N0(:,IJK,M)
181                      Ulplus  = QMOMK_U0(:,IJK,M)
182                      Vlplus  = QMOMK_V0(:,IJK,M)
183                      Wlplus  = QMOMK_W0(:,IJK,M)
184     
185                      Nrminus = QMOMK_N0(:,IJK,M)
186                      Urminus = QMOMK_U0(:,IJK,M)
187                      Vrminus = QMOMK_V0(:,IJK,M)
188                      Wrminus = QMOMK_W0(:,IJK,M)
189     
190                      Nrplus  = QMOMK_N0(:,IPJK,M)
191                      Urplus  = QMOMK_U0(:,IPJK,M)
192                      Vrplus  = QMOMK_V0(:,IPJK,M)
193                      Wrplus  = QMOMK_W0(:,IPJK,M)
194     
195                      CALL KINETIC_FLUX_X_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
196                           Nlplus, Ulplus, Vlplus, Wlplus, F_x_left)
197                      CALL KINETIC_FLUX_X_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
198                           Nrplus, Urplus, Vrplus, Wrplus, F_x_right)
199     
200                      ! y - direction
201                      Nlminus = QMOMK_N0(:,IJMK,M)
202                      Ulminus = QMOMK_U0(:,IJMK,M)
203                      Vlminus = QMOMK_V0(:,IJMK,M)
204                      Wlminus = QMOMK_W0(:,IJMK,M)
205     
206                      Nlplus  = QMOMK_N0(:,IJK,M)
207                      Ulplus  = QMOMK_U0(:,IJK,M)
208                      Vlplus  = QMOMK_V0(:,IJK,M)
209                      Wlplus  = QMOMK_W0(:,IJK,M)
210     
211                      Nrminus = QMOMK_N0(:,IJK,M)
212                      Urminus = QMOMK_U0(:,IJK,M)
213                      Vrminus = QMOMK_V0(:,IJK,M)
214                      Wrminus = QMOMK_W0(:,IJK,M)
215     
216                      Nrplus  = QMOMK_N0(:,IJPK,M)
217                      Urplus  = QMOMK_U0(:,IJPK,M)
218                      Vrplus  = QMOMK_V0(:,IJPK,M)
219                      Wrplus  = QMOMK_W0(:,IJPK,M)
220     
221                      CALL KINETIC_FLUX_Y_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
222                           Nlplus, Ulplus, Vlplus, Wlplus, F_y_left)
223                      CALL KINETIC_FLUX_Y_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
224                           Nrplus, Urplus, Vrplus, Wrplus, F_y_right)
225     
226                      QMOMK_M1(:,IJK,M) = QMOMK_M0(:,IJK,M) &
227                           - 0.5*(QMOMK_DT/MINVAL(DX))*(F_x_right - F_x_left) - 0.5*(QMOMK_DT/MINVAL(DY))*(F_y_right - F_y_left)
228     
229     
230                   ELSE ! 3D case
231                      ! x - direction
232                      Nlminus = QMOMK_N0(:,IMJK,M)
233                      Ulminus = QMOMK_U0(:,IMJK,M)
234                      Vlminus = QMOMK_V0(:,IMJK,M)
235                      Wlminus = QMOMK_W0(:,IMJK,M)
236     
237                      Nlplus  = QMOMK_N0(:,IJK,M)
238                      Ulplus  = QMOMK_U0(:,IJK,M)
239                      Vlplus  = QMOMK_V0(:,IJK,M)
240                      Wlplus  = QMOMK_W0(:,IJK,M)
241     
242                      Nrminus = QMOMK_N0(:,IJK,M)
243                      Urminus = QMOMK_U0(:,IJK,M)
244                      Vrminus = QMOMK_V0(:,IJK,M)
245                      Wrminus = QMOMK_W0(:,IJK,M)
246     
247                      Nrplus  = QMOMK_N0(:,IPJK,M)
248                      Urplus  = QMOMK_U0(:,IPJK,M)
249                      Vrplus  = QMOMK_V0(:,IPJK,M)
250                      Wrplus  = QMOMK_W0(:,IPJK,M)
251     
252                      CALL KINETIC_FLUX_X_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
253                           Nlplus, Ulplus, Vlplus, Wlplus, F_x_left)
254                      CALL KINETIC_FLUX_X_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
255                           Nrplus, Urplus, Vrplus, Wrplus, F_x_right)
256     
257                      ! y - direction
258                      Nlminus = QMOMK_N0(:,IJMK,M)
259                      Ulminus = QMOMK_U0(:,IJMK,M)
260                      Vlminus = QMOMK_V0(:,IJMK,M)
261                      Wlminus = QMOMK_W0(:,IJMK,M)
262     
263                      Nlplus  = QMOMK_N0(:,IJK,M)
264                      Ulplus  = QMOMK_U0(:,IJK,M)
265                      Vlplus  = QMOMK_V0(:,IJK,M)
266                      Wlplus  = QMOMK_W0(:,IJK,M)
267     
268                      Nrminus = QMOMK_N0(:,IJK,M)
269                      Urminus = QMOMK_U0(:,IJK,M)
270                      Vrminus = QMOMK_V0(:,IJK,M)
271                      Wrminus = QMOMK_W0(:,IJK,M)
272     
273                      Nrplus  = QMOMK_N0(:,IJPK,M)
274                      Urplus  = QMOMK_U0(:,IJPK,M)
275                      Vrplus  = QMOMK_V0(:,IJPK,M)
276                      Wrplus  = QMOMK_W0(:,IJPK,M)
277     
278                      CALL KINETIC_FLUX_Y_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
279                           Nlplus, Ulplus, Vlplus, Wlplus, F_y_left)
280                      CALL KINETIC_FLUX_Y_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
281                           Nrplus, Urplus, Vrplus, Wrplus, F_y_right)
282     
283                      ! z - direction
284                      Nlminus = QMOMK_N0(:,IJKM,M)
285                      Ulminus = QMOMK_U0(:,IJKM,M)
286                      Vlminus = QMOMK_V0(:,IJKM,M)
287                      Wlminus = QMOMK_W0(:,IJKM,M)
288     
289                      Nlplus  = QMOMK_N0(:,IJK,M)
290                      Ulplus  = QMOMK_U0(:,IJK,M)
291                      Vlplus  = QMOMK_V0(:,IJK,M)
292                      Wlplus  = QMOMK_W0(:,IJK,M)
293     
294                      Nrminus = QMOMK_N0(:,IJK,M)
295                      Urminus = QMOMK_U0(:,IJK,M)
296                      Vrminus = QMOMK_V0(:,IJK,M)
297                      Wrminus = QMOMK_W0(:,IJK,M)
298     
299                      Nrplus  = QMOMK_N0(:,IJKP,M)
300                      Urplus  = QMOMK_U0(:,IJKP,M)
301                      Vrplus  = QMOMK_V0(:,IJKP,M)
302                      Wrplus  = QMOMK_W0(:,IJKP,M)
303     
304                      CALL KINETIC_FLUX_Z_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
305                           Nlplus, Ulplus, Vlplus, Wlplus, F_z_left)
306                      CALL KINETIC_FLUX_Z_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
307                           Nrplus, Urplus, Vrplus, Wrplus, F_z_right)
308     
309                      QMOMK_M1(:,IJK,M) = QMOMK_M0(:,IJK,M) - 0.5*(QMOMK_DT/MINVAL(DX))*(F_x_right - F_x_left) - &
310                           0.5*(QMOMK_DT/MINVAL(DY))*(F_y_right - F_y_left) - 0.5*(QMOMK_DT/MINVAL(DZ))*(F_z_right - F_z_left)
311                   END IF
312                END IF
313             END DO
314          END DO
315          ! End of fluxes calcuation
316     
317          ! Update weights and abscissas
318          DO M = 1, MMAX
319             DO IJK = ijkstart3, ijkend3
320                IF (FLUID_AT(IJK))  THEN
321                   CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), &
322                        QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
323                   CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
324                        QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
325                END IF
326             END DO
327          END DO
328     
329          ! RHS of moments equations
330          ! Gravity and drag contribution (Half time step)
331          DO M = 1, MMAX
332             DO IJK = ijkstart3, ijkend3
333              IF (FLUID_AT(IJK)) THEN
334                I = I_OF(IJK)
335                J = J_OF(IJK)
336                K = K_OF(IJK)
337                IMJK = IM_OF(IJK)
338                IJMK = JM_OF(IJK)
339                IJKM = KM_OF(IJK)
340     
341                UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
342                VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
343                WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
344     
345                DO IN = 1, QMOMK_NN
346                   vrel = SQRT(((QMOMK_U1(IN,IJK,M) - UGC)**2) + ((QMOMK_V1(IN,IJK,M) - VGC)**2) + ((QMOMK_W1(IN,IJK,M) - WGC)**2))
347     
348                   ! Particle Reynolds number
349                   Rep = RO_g0*D_p0(M)*vrel/MU_g0
350                   Cd = 24.D0*(EP_G(IJK)**(-2.65))*(1.D0 + 0.15D0*((EP_G(IJK)*Rep)**0.687))/(Rep + SMALL_NUMBER)
351                   beta_drag = 3.D0*RO_g0*Cd*vrel/(4.D0*D_p0(M)*RO_s(IJK,M))
352     
353                   drag_exp = EXP(-0.5*QMOMK_DT*beta_drag)
354     
355                   IF (beta_drag > SMALL_NUMBER) THEN   !CHECK HERE, it was /= 0
356                     QMOMK_TAU_DRAG(IN,IJK,M) = 1.D0/beta_drag
357                   ELSE
358                     QMOMK_TAU_DRAG(IN,IJK,M) = LARGE_NUMBER
359                   END IF
360     
361                   ! Drag calculation for QMOMK
362                   QMOMK_U1 (IN,IJK,M) = drag_exp*QMOMK_U1 (IN,IJK,M) + (1.D0 - drag_exp)*UGC
363                   QMOMK_V1 (IN,IJK,M) = drag_exp*QMOMK_V1 (IN,IJK,M) + (1.D0 - drag_exp)*(VGC - GRAVITY*QMOMK_TAU_DRAG(IN,IJK,M))
364                   QMOMK_W1 (IN,IJK,M) = drag_exp*QMOMK_W1 (IN,IJK,M) + (1.D0 - drag_exp)*WGC
365                END DO
366                CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
367                     QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
368                CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), &
369                     QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
370                CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
371                     QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
372               END IF
373             END DO
374          END DO
375     
376          DO M = 1, MMAX
377             DO IJK = ijkstart3, ijkend3
378               IF (FLUID_AT(IJK)) THEN
379                 CALL BIND_THETA(QMOMK_M1(:,IJK,M), MINVAL(QMOMK_TAU_DRAG(:,IJK,M)))
380                 CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), &
381                      QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
382                CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
383                     QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
384               END IF
385             END DO
386          END DO
387     
388          ! End of gravity and drag contribution
389     
390          PRINT *,'Before collisions'
391     
392          ! Collisions contribution (Half time step)
393          IF (QMOMK_COLLISIONS == 'BGK' .AND. MMAX == 1) THEN
394             DO M = 1, MMAX
395                DO IJK = ijkstart3, ijkend3
396                   IF (FLUID_AT(IJK))  THEN
397                      mom0 = QMOMK_M1(1,IJK,M)
398                      IF (mom0 > epsn) THEN
399                         CALL COLLISIONS_BGK(QMOMK_M1(:,IJK,M), 0.5*QMOMK_DT, QMOMK_TCOL, D_p0 (M), C_e)
400                         CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), &
401                              QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
402                         CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
403                              QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
404                      END IF
405                   END IF
406                END DO
407             END DO
408          ELSE IF (QMOMK_COLLISIONS == 'BOLTZMANN' .AND. MMAX == 1) THEN
409             DO IJK = ijkstart3, ijkend3
410                IF (FLUID_AT(IJK)) THEN
411                   CALL SOLVE_BOLTZMANN_COLLISIONS_ONE_SPECIE(QMOMK_M1(:,IJK,1), QMOMK_N1(:,IJK,1), &
412                        QMOMK_U1(:,IJK,1), QMOMK_V1(:,IJK,1), QMOMK_W1(:,IJK,1), &
413                        0.5*QMOMK_DT, C_e, D_p0 (1), QMOMK_COLLISIONS_ORDER)
414                   CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,1), QMOMK_N1(:,IJK,1), QMOMK_U1(:,IJK,1), QMOMK_V1(:,IJK,1), QMOMK_W1(:,IJK,1))
415                   CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,1), &
416                        QMOMK_U1(:,IJK,1), QMOMK_V1(:,IJK,1), QMOMK_W1(:,IJK,1), QMOMK_M1(:,IJK,1))
417                END IF
418             END DO
419          ELSE
420             DO M = 1, MMAX
421                DO IJK = ijkstart3, ijkend3
422                   IF (FLUID_AT(IJK)) THEN
423                      DO M2 = 1, MMAX
424                         CALL SOLVE_BOLTZMANN_COLLISIONS_TWO_SPECIES(QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), &
425                              QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), &
426                              QMOMK_M1(:,IJK,M2), QMOMK_N1(:,IJK,M2), QMOMK_U1(:,IJK,M2), QMOMK_V1(:,IJK,M2), QMOMK_W1(:,IJK,M2), &
427                              0.5*QMOMK_DT, 1.D0/6.D0*Pi*(D_p0(M)**3)*RO_s(IJK,M), 1.D0/6.D0*Pi*(D_p0(M2)**3)*RO_s(IJK,M2), &
428                              D_p0(M), D_p0(M2), C_e, C_e, QMOMK_COLLISIONS_ORDER)
429                      END DO
430                      CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
431                      CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
432                           QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
433                   END IF
434                END DO
435             END DO
436          END IF
437     
438          PRINT *,'After collisions'
439     
440          ! End of collisions contribution
441     
442          ! Boundary conditions update
443           CALL QMOMK_SET_BC
444     
445          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446          ! Second step of Runge-Kutta time integration !
447          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448     
449          ! M0 -> old values
450          ! M1 -> first step RK values
451     
452          ! Calculating kinetic based fluxes (Full time step)
453          DO M = 1, MMAX
454             DO IJK = ijkstart3, ijkend3
455                IF (FLUID_AT(IJK)) THEN
456                   I = I_OF(IJK)
457                   J = J_OF(IJK)
458                   K = K_OF(IJK)
459     
460                   IMJK = IM_OF(IJK)
461                   IJMK = JM_OF(IJK)
462                   IJKM = KM_OF(IJK)
463     
464                   IPJK = IP_OF(IJK)
465                   IJPK = JP_OF(IJK)
466                   IJKP = KP_OF(IJK)
467     
468                   ! 2D Case
469     
470                   IF (NO_K) THEN
471                      ! x - direction
472                      Nlminus = QMOMK_N1(:,IMJK,M)
473                      Ulminus = QMOMK_U1(:,IMJK,M)
474                      Vlminus = QMOMK_V1(:,IMJK,M)
475                      Wlminus = QMOMK_W1(:,IMJK,M)
476     
477                      Nlplus  = QMOMK_N1(:,IJK,M)
478                      Ulplus  = QMOMK_U1(:,IJK,M)
479                      Vlplus  = QMOMK_V1(:,IJK,M)
480                      Wlplus  = QMOMK_W1(:,IJK,M)
481     
482                      Nrminus = QMOMK_N1(:,IJK,M)
483                      Urminus = QMOMK_U1(:,IJK,M)
484                      Vrminus = QMOMK_V1(:,IJK,M)
485                      Wrminus = QMOMK_W1(:,IJK,M)
486     
487                      Nrplus  = QMOMK_N1(:,IPJK,M)
488                      Urplus  = QMOMK_U1(:,IPJK,M)
489                      Vrplus  = QMOMK_V1(:,IPJK,M)
490                      Wrplus  = QMOMK_W1(:,IPJK,M)
491     
492                      CALL KINETIC_FLUX_X_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
493                           Nlplus, Ulplus, Vlplus, Wlplus, F_x_left)
494                      CALL KINETIC_FLUX_X_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
495                           Nrplus, Urplus, Vrplus, Wrplus, F_x_right)
496     
497                      ! y - direction
498                      Nlminus = QMOMK_N1(:,IJMK,M)
499                      Ulminus = QMOMK_U1(:,IJMK,M)
500                      Vlminus = QMOMK_V1(:,IJMK,M)
501                      Wlminus = QMOMK_W1(:,IJMK,M)
502     
503                      Nlplus  = QMOMK_N1(:,IJK,M)
504                      Ulplus  = QMOMK_U1(:,IJK,M)
505                      Vlplus  = QMOMK_V1(:,IJK,M)
506                      Wlplus  = QMOMK_W1(:,IJK,M)
507     
508                      Nrminus = QMOMK_N1(:,IJK,M)
509                      Urminus = QMOMK_U1(:,IJK,M)
510                      Vrminus = QMOMK_V1(:,IJK,M)
511                      Wrminus = QMOMK_W1(:,IJK,M)
512     
513                      Nrplus  = QMOMK_N1(:,IJPK,M)
514                      Urplus  = QMOMK_U1(:,IJPK,M)
515                      Vrplus  = QMOMK_V1(:,IJPK,M)
516                      Wrplus  = QMOMK_W1(:,IJPK,M)
517     
518                      CALL KINETIC_FLUX_Y_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
519                           Nlplus, Ulplus, Vlplus, Wlplus, F_y_left)
520                      CALL KINETIC_FLUX_Y_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
521                           Nrplus, Urplus, Vrplus, Wrplus, F_y_right)
522     
523                      QMOMK_M1(:,IJK,M) = QMOMK_M0(:,IJK,M) - (QMOMK_DT/MINVAL(DX))*(F_x_right - F_x_left) &
524                           - (QMOMK_DT/MINVAL(DY))*(F_y_right - F_y_left)
525     
526                      ! 3D case
527                   ELSE
528                      ! x - direction
529                      Nlminus = QMOMK_N1(:,IMJK,M)
530                      Ulminus = QMOMK_U1(:,IMJK,M)
531                      Vlminus = QMOMK_V1(:,IMJK,M)
532                      Wlminus = QMOMK_W1(:,IMJK,M)
533     
534                      Nlplus  = QMOMK_N1(:,IJK,M)
535                      Ulplus  = QMOMK_U1(:,IJK,M)
536                      Vlplus  = QMOMK_V1(:,IJK,M)
537                      Wlplus  = QMOMK_W1(:,IJK,M)
538     
539                      Nrminus = QMOMK_N1(:,IJK,M)
540                      Urminus = QMOMK_U1(:,IJK,M)
541                      Vrminus = QMOMK_V1(:,IJK,M)
542                      Wrminus = QMOMK_W1(:,IJK,M)
543     
544                      Nrplus  = QMOMK_N1(:,IPJK,M)
545                      Urplus  = QMOMK_U1(:,IPJK,M)
546                      Vrplus  = QMOMK_V1(:,IPJK,M)
547                      Wrplus  = QMOMK_W1(:,IPJK,M)
548     
549                      CALL KINETIC_FLUX_X_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
550                           Nlplus, Ulplus, Vlplus, Wlplus, F_x_left)
551                      CALL KINETIC_FLUX_X_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
552                           Nrplus, Urplus, Vrplus, Wrplus, F_x_right)
553     
554                      ! y - direction
555                      Nlminus = QMOMK_N1(:,IJMK,M)
556                      Ulminus = QMOMK_U1(:,IJMK,M)
557                      Vlminus = QMOMK_V1(:,IJMK,M)
558                      Wlminus = QMOMK_W1(:,IJMK,M)
559     
560                      Nlplus  = QMOMK_N1(:,IJK,M)
561                      Ulplus  = QMOMK_U1(:,IJK,M)
562                      Vlplus  = QMOMK_V1(:,IJK,M)
563                      Wlplus  = QMOMK_W1(:,IJK,M)
564     
565                      Nrminus = QMOMK_N1(:,IJK,M)
566                      Urminus = QMOMK_U1(:,IJK,M)
567                      Vrminus = QMOMK_V1(:,IJK,M)
568                      Wrminus = QMOMK_W1(:,IJK,M)
569     
570                      Nrplus  = QMOMK_N1(:,IJPK,M)
571                      Urplus  = QMOMK_U1(:,IJPK,M)
572                      Vrplus  = QMOMK_V1(:,IJPK,M)
573                      Wrplus  = QMOMK_W1(:,IJPK,M)
574     
575                      CALL KINETIC_FLUX_Y_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
576                           Nlplus, Ulplus, Vlplus, Wlplus, F_y_left)
577                      CALL KINETIC_FLUX_Y_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
578                           Nrplus, Urplus, Vrplus, Wrplus, F_y_right)
579     
580                      ! z - direction
581                      Nlminus = QMOMK_N1(:,IJKM,M)
582                      Ulminus = QMOMK_U1(:,IJKM,M)
583                      Vlminus = QMOMK_V1(:,IJKM,M)
584                      Wlminus = QMOMK_W1(:,IJKM,M)
585     
586                      Nlplus  = QMOMK_N1(:,IJK,M)
587                      Ulplus  = QMOMK_U1(:,IJK,M)
588                      Vlplus  = QMOMK_V1(:,IJK,M)
589                      Wlplus  = QMOMK_W1(:,IJK,M)
590     
591                      Nrminus = QMOMK_N1(:,IJK,M)
592                      Urminus = QMOMK_U1(:,IJK,M)
593                      Vrminus = QMOMK_V1(:,IJK,M)
594                      Wrminus = QMOMK_W1(:,IJK,M)
595     
596                      Nrplus  = QMOMK_N1(:,IJKP,M)
597                      Urplus  = QMOMK_U1(:,IJKP,M)
598                      Vrplus  = QMOMK_V1(:,IJKP,M)
599                      Wrplus  = QMOMK_W1(:,IJKP,M)
600     
601                      CALL KINETIC_FLUX_Z_TWENTY_EIGHT_NODES (Nlminus, Ulminus, Vlminus, Wlminus, &
602                           Nlplus, Ulplus, Vlplus, Wlplus, F_z_left)
603                      CALL KINETIC_FLUX_Z_TWENTY_EIGHT_NODES (Nrminus, Urminus, Vrminus, Wrminus, &
604                           Nrplus, Urplus, Vrplus, Wrplus, F_z_right)
605     
606                      QMOMK_M1(:,IJK,M) = QMOMK_M0(:,IJK,M) - (QMOMK_DT/MINVAL(DX))*(F_x_right - F_x_left) - &
607                           (QMOMK_DT/MINVAL(DY))*(F_y_right - F_y_left) - (QMOMK_DT/MINVAL(DZ))*(F_z_right - F_z_left)
608                   END IF
609                END IF
610             END DO
611          END DO
612     
613          ! Update weights and abscissas
614          DO M = 1, MMAX
615             DO IJK = ijkstart3, ijkend3
616                IF (FLUID_AT(IJK))  THEN
617                   CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
618                   CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
619                        QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
620                END IF
621             END DO
622          END DO
623     
624          ! RHS of moments equations
625          ! Collisions contribution (Full time step)
626          IF (QMOMK_COLLISIONS == 'BGK' .AND. MMAX == 1) THEN
627             DO M = 1, MMAX
628                DO IJK = ijkstart3, ijkend3
629                   IF (FLUID_AT(IJK))  THEN
630                      mom0 = QMOMK_M1(1,IJK,M)
631                      IF (mom0 > epsn) THEN
632                         CALL COLLISIONS_BGK(QMOMK_M1(:,IJK,M), QMOMK_DT, QMOMK_TCOL, D_p0 (M), C_e)
633                         CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), &
634                              QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
635                         CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
636                              QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
637                      END IF
638                   END IF
639                END DO
640             END DO
641          ELSE IF (QMOMK_COLLISIONS == 'BOLTZMANN' .AND. MMAX == 1) THEN
642             DO IJK = ijkstart3, ijkend3
643                IF (FLUID_AT(IJK)) THEN
644                   CALL SOLVE_BOLTZMANN_COLLISIONS_ONE_SPECIE(QMOMK_M1(:,IJK,1), QMOMK_N1(:,IJK,1), &
645                        QMOMK_U1(:,IJK,1), QMOMK_V1(:,IJK,1), QMOMK_W1(:,IJK,1), &
646                        QMOMK_DT, C_e, D_p0 (1), QMOMK_COLLISIONS_ORDER)
647                   CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,1), QMOMK_N1(:,IJK,1), QMOMK_U1(:,IJK,1), QMOMK_V1(:,IJK,1), QMOMK_W1(:,IJK,1))
648                   CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,1), &
649                        QMOMK_U1(:,IJK,1), QMOMK_V1(:,IJK,1), QMOMK_W1(:,IJK,1), QMOMK_M1(:,IJK,1))
650                END IF
651             END DO
652          ELSE
653             DO M = 1, MMAX
654                DO IJK = ijkstart3, ijkend3
655                   IF (FLUID_AT(IJK)) THEN
656                      DO M2 = 1, MMAX
657                         CALL SOLVE_BOLTZMANN_COLLISIONS_TWO_SPECIES(QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), &
658                              QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), &
659                              QMOMK_M1(:,IJK,M2), QMOMK_N1(:,IJK,M2), &
660                              QMOMK_U1(:,IJK,M2), QMOMK_V1(:,IJK,M2), QMOMK_W1(:,IJK,M2), &
661                              QMOMK_DT, 1.D0/6.D0*Pi*(D_p0(M)**3)*RO_s(IJK,M), &
662                              1.D0/6.D0*Pi*(D_p0(M2)**3)*RO_s(IJK,M2),  D_p0(M), D_p0(M2), &
663                              C_e, C_e, QMOMK_COLLISIONS_ORDER)
664                      END DO
665                      CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
666                      CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
667                           QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
668                   END IF
669                END DO
670             END DO
671          END IF
672          ! End of collisions contribution
673     
674          ! Gravity and drag contribution (Full time step)
675          DO M = 1, MMAX
676             DO IJK = ijkstart3, ijkend3
677                IF (FLUID_AT(IJK)) THEN
678                   I = I_OF(IJK)
679     
680                   IMJK = IM_OF(IJK)
681                   IJMK = JM_OF(IJK)
682                   IJKM = KM_OF(IJK)
683     
684                   UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
685                   VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
686                   WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
687     
688                   DO IN = 1, QMOMK_NN
689                      vrel = SQRT(((QMOMK_U1(IN,IJK,M) - UGC)**2) + ((QMOMK_V1(IN,IJK,M) - VGC)**2) + ((QMOMK_W1(IN,IJK,M) - WGC)**2))
690     
691                      ! Particle Reynolds number
692                      Rep = RO_g0*D_p0(M)*vrel/MU_g0
693                      Cd = 24.D0*(EP_G(IJK)**(-2.65))*(1.D0 + 0.15*((EP_G(IJK)*Rep)**0.687))/(Rep+SMALL_NUMBER)
694                      beta_drag = 3.D0*RO_g0*Cd*vrel/(4.D0*D_p0(M)*RO_s(IJK,M))
695     
696                      IF (beta_drag > SMALL_NUMBER) THEN ! CHECK HERE, it was /= 0.
697                         QMOMK_TAU_DRAG(IN,IJK,M) = 1.D0/beta_drag
698                      ELSE
699                         QMOMK_TAU_DRAG(IN,IJK,M) = LARGE_NUMBER
700                      END IF
701     
702                      drag_exp = EXP(-QMOMK_DT*beta_drag)
703     
704                      ! Drag calculation for QMOMK
705                      QMOMK_U1 (IN,IJK,M) = drag_exp*QMOMK_U1 (IN,IJK,M) + (1.D0 - drag_exp)*UGC
706                      QMOMK_V1 (IN,IJK,M) = drag_exp*QMOMK_V1 (IN,IJK,M) + (1.D0 - drag_exp)*(VGC - GRAVITY*QMOMK_TAU_DRAG(IN,IJK,M))
707                      QMOMK_W1 (IN,IJK,M) = drag_exp*QMOMK_W1 (IN,IJK,M) + (1.D0 - drag_exp)*WGC
708                   END DO
709                   CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
710                        QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
711                   CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
712                   CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
713                        QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
714                END IF
715             END DO
716          END DO
717          ! End of gravity and drag contribution
718     
719          DO M = 1, MMAX
720            DO IJK = ijkstart3, ijkend3
721              IF (FLUID_AT(IJK)) THEN
722                CALL BIND_THETA(QMOMK_M1(:,IJK,M), MINVAL(QMOMK_TAU_DRAG(:,IJK,M)))
723                CALL EIGHT_NODE_3D (QMOMK_M1(:,IJK,M), QMOMK_N1(:,IJK,M), QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M))
724                CALL MOMENTS_TWENTY_EIGHT_NODES (QMOMK_N1(:,IJK,M), &
725                     QMOMK_U1(:,IJK,M), QMOMK_V1(:,IJK,M), QMOMK_W1(:,IJK,M), QMOMK_M1(:,IJK,M))
726              END IF
727            END DO
728          END DO
729     
730          ! Boundary conditions update
731          CALL QMOMK_SET_BC
732     
733          ! Preparing data for next time step
734          QMOMK_M0 = QMOMK_M1
735          QMOMK_N0 = QMOMK_N1
736          QMOMK_U0 = QMOMK_U1
737          QMOMK_V0 = QMOMK_V1
738          QMOMK_W0 = QMOMK_W1
739     
740          PRINT *,'QMOM DT 5 = ',QMOMK_DT
741     
742          QMOMK_TIME = QMOMK_TIME + QMOMK_DT
743     
744          PRINT *,'QMOM DT 6 = ',QMOMK_DT
745       END DO  ! QMOMK internal time stepping
746     
747       IF (DT .LT. QMOMK_DT_TMP) THEN
748         QMOMK_DT = QMOMK_DT_TMP
749       ENDIF
750     
751       IF (TMP_DTS .NE. ZERO) THEN
752         QMOMK_DT = TMP_DTS
753         TMP_DTS = ZERO
754       ENDIF
755     
756       PRINT *,'Time = ',TIME
757       PRINT *,'Time QMOMK = ',QMOMK_TIME
758     
759       ! Resetting gas volume fraction to one, in order to update it
760       DO IJK = ijkstart3, ijkend3
761          IF (FLUID_AT(IJK)) THEN
762             EP_G(IJK) = ONE
763          END IF
764       END DO
765     
766       ! Passing values to the fluid solver and for storage
767       DO M = 1, MMAX
768          DO IJK = ijkstart3, ijkend3
769             IF (FLUID_AT(IJK)) THEN
770                ! Mean particle velocities
771                U_S(IJK, M) = QMOMK_M1(2,IJK,M)/QMOMK_M1(1,IJK,M)
772                V_S(IJK, M) = QMOMK_M1(3,IJK,M)/QMOMK_M1(1,IJK,M)
773                W_S(IJK, M) = QMOMK_M1(4,IJK,M)/QMOMK_M1(1,IJK,M)
774     
775                ! Gas-phase volume fraction
776                EP_G(IJK) = EP_G(IJK) - QMOMK_M1(1,IJK,M)
777     
778                ! Particle phases "densities"
779                ROP_S(IJK, M) = QMOMK_M1(1,IJK,M)*RO_s(IJK,M)
780     
781                ! Particle phases granular temperatures
782                THETA_M(IJK,M) = ((QMOMK_M1(5,IJK,M)/QMOMK_M1(1,IJK,M) - &
783                                 (QMOMK_M1(2,IJK,M)/QMOMK_M1(1,IJK,M))*(QMOMK_M1(2,IJK,M)/QMOMK_M1(1,IJK,M))) + &
784                                 (QMOMK_M1(8,IJK,M)/QMOMK_M1(1,IJK,M) - &
785                                    (QMOMK_M1(3,IJK,M)/QMOMK_M1(1,IJK,M))*(QMOMK_M1(3,IJK,M)/QMOMK_M1(1,IJK,M))) + &
786                                 (QMOMK_M1(10,IJK,M)/QMOMK_M1(1,IJK,M) - &
787                                    (QMOMK_M1(4,IJK,M)/QMOMK_M1(1,IJK,M))*(QMOMK_M1(4,IJK,M)/QMOMK_M1(1,IJK,M))))/3.D0
788     
789             END IF
790          END DO
791       END DO
792     
793       ! Calculating momentum exchange terms
794       IF (QMOMK_COUPLED) THEN
795          DO M = 1, MMAX
796             DO IJK = ijkstart3, ijkend3
797                IF (FLUID_AT(IJK)) THEN
798                   I = I_OF(IJK)
799     
800                   IMJK = IM_OF(IJK)
801                   IJMK = JM_OF(IJK)
802                   IJKM = KM_OF(IJK)
803     
804                   UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
805                   VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
806                   WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
807     
808                   DO IN = 1, QMOMK_NN
809                      vrel = SQRT(((QMOMK_U1(IN,IJK,M) - UGC)**2) + ((QMOMK_V1(IN,IJK,M) - VGC)**2) + ((QMOMK_W1(IN,IJK,M) - WGC)**2))
810     
811                      Rep = RO_g0*D_p0(M)*vrel/MU_g0
812                      Cd = 24.D0*(EP_G(IJK)**(-2.65))*(1.D0 + 0.15*((EP_G(IJK)*Rep)**0.687))/(Rep+SMALL_NUMBER)
813                      beta_drag = 3.D0*QMOMK_N1(IN,IJK,M)*RO_g0*Cd*vrel/(4.D0*D_p0(M))
814     
815                      QMOMK_F_GS(IN,IJK,M) = beta_drag
816                   END DO
817                END IF
818             END DO
819          END DO
820       END IF
821     END SUBROUTINE QMOMK_TIME_MARCH
822