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