integrate_time

File: N:\mfix\model\des\pic\integrate_time_pic.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module name: CFNEWVALUES                                            !
4     !                                                                      !
5     !  Purpose: DES - Calculate the new values of particle velocity,       !
6     !           position, angular velocity etc                             !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE INTEGRATE_TIME_PIC
10     
11           USE discretelement
12           USE mfix_pic
13     
14           IMPLICIT NONE
15     
16     ! call the coloring function like approach
17           IF(MPPIC_SOLID_STRESS_SNIDER) THEN
18              CALL INTEGRATE_TIME_PIC_SNIDER
19           ELSE
20              CALL INTEGRATE_TIME_PIC_GARG
21           ENDIF
22     
23           RETURN
24           END SUBROUTINE INTEGRATE_TIME_PIC
25     
26     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
27     !                                                                      !
28     !  Module name: INTEGRATE_TIME_PIC_SNIDER                              !
29     !                                                                      !
30     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
31           SUBROUTINE INTEGRATE_TIME_PIC_SNIDER
32     
33           USE param
34           USE param1
35           USE parallel
36           USE physprop
37           USE constant
38           USE fldvar
39           USE discretelement
40           USE des_bc
41           USE mpi_utility
42           USE fldvar
43           USE cutcell
44           USE mfix_pic
45           USE randomno
46           USE geometry, only: DO_K, NO_K
47           USE fun_avg
48           USE functions
49     
50           IMPLICIT NONE
51     
52     ! Local Variables
53     !---------------------------------------------------------------------//
54     ! Loop counters
55           INTEGER :: NP, LC, PC
56     ! Drag coefficient divided by parcel mass
57           DOUBLE PRECISION :: DP_BAR
58     ! Restitution coefficient plus one
59           DOUBLE PRECISION :: ENp1
60     ! Integrated particle velocity without particle normal stress
61           DOUBLE PRECISION :: VEL(3)
62     ! Estimated parcel velocity from normal particle stress
63           DOUBLE PRECISION :: DELUP(3)
64     ! Actual parcel velocity from particle normal stress.
65           DOUBLE PRECISION :: UPRIMETAU(3)
66     ! Slip velocity between parcel and bulk solids
67           DOUBLE PRECISION :: SLIPVEL(3)
68     ! Distance traveled by the parcel (and magnitude)
69           DOUBLE PRECISION :: DIST(3), DIST_MAG
70     ! Max velocity magnitude of all parcels
71           DOUBLE PRECISION :: MAX_VEL
72     ! Mininum fluid cell dimenstion
73           DOUBLE PRECISION :: DXYZ_MIN
74     ! Three times the square root of two.
75           DOUBLE PRECISION :: l3SQT2
76     ! Volume fraction used to limit parcel movement in close pack regions.
77           DOUBLE PRECISION :: EPg_MFP
78     ! Mean free path of parcels based on volume fraction.
79           DOUBLE PRECISION :: MFP
80     ! Solids volume fraction at the parcel
81           DOUBLE PRECISION :: EPs
82     ! Loop bound
83           INTEGER :: LC_BND
84     ! Gravity normal and magnitude
85           DOUBLE PRECISION :: GRAV_NORM(3)
86     !......................................................................!
87     
88     ! Volume fraction used to limit parcel movement in close pack regions.
89           EPg_MFP = 0.95d0*EP_STAR
90     ! Initialize max parcel velocity
91           MAX_VEL = -UNDEFINED
92     ! Restitution coefficient plus one.
93           ENp1 = MPPIC_COEFF_EN1 + 1.0d0
94     ! Three times the square root of two.
95           l3SQT2 = merge(3.0d0, 2.0d0, DO_K)*sqrt(2.0d0)
96     ! Dimenion loop bound
97           LC_BND = merge(2,3,NO_K)
98     
99     ! Calculate the normal of gravity
100           IF(GRAV_MAG > SMALL_NUMBER) GRAV_NORM = GRAV/GRAV_MAG
101     
102           PC = 1
103           DO NP = 1, MAX_PIP
104     
105              IF(PC.GT.PIP) EXIT
106              IF(IS_NONEXISTENT(NP)) CYCLE
107     
108              PC = PC+1
109     
110     ! If a particle is classified as new, then forces are ignored.
111     ! Classification from new to existing is performed in routine
112     ! des_check_new_particle.f
113              IF(IS_ENTERING(NP))THEN
114                 FC(NP,:) = ZERO
115              ELSE
116                 FC(NP,:) = FC(NP,:)/PMASS(NP) + GRAV(:)
117              ENDIF
118     
119     ! DP_BAR is the D_p in the Snider paper (JCP, vol 170, 523-549, 2001)
120     ! By comparing MFIX equations and equations in those papers,
121     ! D_p = Beta/(EP_S*RHOP)
122     ! F_gp in drag_fgs.f  = Beta*PVOL/EP_S
123     ! Therefore, D_p = F_gp/(PVOL*RHOP) = F_gp/PMASS
124              IF(DES_CONTINUUM_COUPLED .AND. MPPIC_PDRAG_IMPLICIT)&
125                 DP_BAR = F_gp(NP)/PMASS(NP)
126     
127     ! Solids volume fraction
128     !         EPs = max(ONE-EP_G(PIJK(NP,4)), 1.0d-4)
129              EPs = max(ONE-EPG_P(NP), 1.0d-4)
130     
131     ! Numerically integrated particle velocity without particle normal
132     ! stress force :: Snider (Eq 38)
133              VEL(:) = (DES_VEL_NEW(NP,:) + &
134                 FC(NP,:)*DTSOLID)/(1.d0+DP_BAR*DTSOLID)
135     
136     ! Estimated discrete particle velocity from continuum particle
137     ! normal stress graident :: Snider (Eq 39)
138              DELUP(:) = -(DTSOLID*PS_GRAD(:,NP)) / &
139                 (EPs*RO_S0(1)*(1.d0+DP_BAR*DTSOLID))
140     
141     ! Slip velocity between parcel and average solids velocity adjusted
142     ! for gravity driven flows.
143                 SLIPVEL = AVG_VEL_wGRAV(NP,VEL) - VEL
144     
145     ! Calculate particle velocity from particle normal stress.
146              DO LC = 1, LC_BND
147     ! Snider (40)
148                 IF(PS_GRAD(LC,NP) <= ZERO) THEN
149                    UPRIMETAU(LC) = min(DELUP(LC), ENp1*SLIPVEL(LC))
150                    UPRIMETAU(LC) = max(UPRIMETAU(LC), ZERO)
151     ! Snider (41)
152                 ELSE
153                    UPRIMETAU(LC) = max(DELUP(LC), ENp1*SLIPVEL(LC))
154                    UPRIMETAU(LC) = min(UPRIMETAU(LC), ZERO)
155                 ENDIF
156              ENDDO
157     
158     ! Total particle velocity is the sum of the particle velocity from
159     ! normal stress and velocity from all other foces :: Sinder (Eq 37)
160              DES_VEL_NEW(NP,:) = VEL + UPRIMETAU
161     ! Total distance traveled over the current time step
162              DIST(:) = DES_VEL_NEW(NP,:)*DTSOLID
163     
164     ! Limit parcel movement in close pack regions to the mean free path. 
165              IF(EPG_P(NP) < EPg_MFP) THEN
166                 IF(dot_product(DES_VEL_NEW(NP,:),PS_GRAD(:,NP)) > ZERO) THEN
167     ! Total distance traveled given current velocity.
168                    DIST_MAG = dot_product(DIST, DIST)
169     ! Mean free path based on volume fraction :: Snider (43)
170                    MFP = 3.7d0*DES_RADIUS(NP)/(l3SQT2*EPs)
171     ! Impose the distance limit and calculate the corresponding velocity.
172                    IF(MFP**2 < DIST_MAG) THEN
173                       DIST = MFP*DIST/sqrt(DIST_MAG)
174                       DES_VEL_NEW(NP,:) = DIST/DTSOLID
175                    ENDIF
176                 ENDIF
177              ENDIF
178     
179     
180     ! Update the parcel position.
181              DES_POS_NEW(NP,:) = DES_POS_NEW(NP,:) + DIST
182     ! Clear out the force array.
183              FC(NP,:) = ZERO
184     
185     ! Determine the maximum distance traveled by any single parcel.
186              DIST_MAG = dot_product(DES_VEL_NEW(NP,:),DES_VEL_NEW(NP,:)) 
187              MAX_VEL = max(MAX_VEL, DIST_MAG)
188     
189           ENDDO
190     
191     ! Get the minimum fluid cell dimension
192           DXYZ_MIN = min(minval(DX(IMIN1:IMAX1)),minval(DY(JMIN1:JMAX1)))
193           IF(DO_K) DXYZ_MIN = min(DXYZ_MIN,minval(DZ(KMIN1:KMAX1)))
194     
195     ! Calculate the maximum time step to prevent a parcel from moving more
196     ! than a fluid cell in one time step.  The global_all_max could get
197     ! expensive in large MPI applications.
198           DTPIC_CFL = CFL_PIC*DXYZ_MIN/(MAX_VEL+SMALL_NUMBER)
199           CALL global_all_max(DTPIC_CFL)
200     
201           RETURN
202           contains
203     
204     
205     !......................................................................!
206     !                                                                      !
207     ! Function: AVG_SOLIDS_vGRAV                                           !
208     ! Author: J.Musser                                                     !
209     !                                                                      !
210     ! Purpose: Return the average solids velocity around the parcel with   !
211     ! a correction for gravity dominated flows.                            !
212     !                                                                      !
213     ! Notes: The average solids velocity is reduced by 1/2 when it is      !
214     ! moving within 5 degrees of parall of the direction of gravity. This  !
215     ! is helps prevent overpacking while solids are settling. The value is !
216     ! blended from AVG to 0.5AVG from 0.1 to 0.0 radians.                  !
217     !                                                                      !
218     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
219           FUNCTION AVG_VEL_wGRAV(lNP,lVEL) RESULT(aVEL)
220     
221           IMPLICIT NONE
222     
223     ! Particle index
224           INTEGER, INTENT(IN) :: lNP
225     ! Intermediate parcel velocity
226           DOUBLE PRECISION, INTENT(IN) :: lVEL(3)
227     ! Average solids velocity modified for gravity
228           DOUBLE PRECISION :: aVEL(3)
229     ! Temp double precision variables.
230           DOUBLE PRECISION :: tDP1, tDP2
231     
232           aVEL = AVGSOLVEL_P(:,lNP)
233     
234     ! Reduce the avreage solids velocity for gravity driven flows.
235           IF(GRAV_MAG > SMALL_NUMBER) THEN
236              tDP1 = dot_product(GRAV_NORM,aVEL)
237              IF(tDP1 > ZERO) THEN
238                 IF(dot_product(aVEL,lVEL) > ZERO) THEN
239                    tDP2 = dot_product(aVEL, aVEL)
240                    IF(tDP1**2 > 0.99*tDP2) aVEL = &
241                       (50.75d0- 50.0d0*tDP1/sqrt(tDP2))*aVEL   ! 3/4
242     
243     !               IF(tDP1**2 > 0.99*tDP2) aVEL = aVEL * &
244     !                  (100.50d0-100.0d0*tDP1/sqrt(tDP2))*aVEL  ! 1/2
245     
246     !               IF(tDP1**2 > 0.99*tDP2) aVEL = 0.75d0* aVEL
247                 ENDIF
248              ENDIF
249           ENDIF
250     
251           RETURN
252           END FUNCTION AVG_VEL_wGRAV
253     
254           END SUBROUTINE INTEGRATE_TIME_PIC_SNIDER
255     
256     
257     
258     
259     
260     
261     
262     
263     
264     
265     
266     
267     
268     
269     
270     
271     
272     
273     
274     
275     
276     
277     
278     
279     
280     
281     
282     
283     
284     
285           SUBROUTINE INTEGRATE_TIME_PIC_GARG
286     
287           USE param
288           USE param1
289           USE parallel
290           USE physprop
291           USE constant
292           USE fldvar
293           USE discretelement
294           USE des_bc
295           USE mpi_utility
296           USE mfix_pic
297           USE randomno
298           USE cutcell
299           USE fun_avg
300           USE functions
301           IMPLICIT NONE
302     !-----------------------------------------------
303     ! Local Variables
304     !-----------------------------------------------
305           INTEGER NP, M, LC
306           INTEGER I, J, K, IJK, IJK_OLD
307     
308           DOUBLE PRECISION DD(3), DIST, &
309                            DP_BAR, MEANVEL(3), D_GRIDUNITS(3)
310     
311           DOUBLE PRECISION MEAN_FREE_PATH
312     ! index to track accounted for particles
313           INTEGER PC
314     
315     ! Logical for local debug warnings
316           LOGICAL DES_LOC_DEBUG
317     
318     ! maximum distance particles can move in MPPIC
319           DOUBLE PRECISION UPRIMEMOD
320     
321     ! dt's in each direction  based on cfl_pic for the mppic case
322     
323           DOUBLE PRECISION :: DTPIC_TMPX, DTPIC_TMPY , DTPIC_TMPZ
324           DOUBLE PRECISION :: THREEINTOSQRT2, RAD_EFF, POS_Z
325           DOUBLE PRECISION :: VELP_INT(3)
326     
327           LOGICAL :: DELETE_PART, INSIDE_DOMAIN
328           INTEGER :: PIP_DEL_COUNT
329     
330           INTEGER :: LPIP_DEL_COUNT_ALL(0:numPEs-1), PIJK_OLD(5)
331     
332           double precision  sig_u, mean_u
333           DOUBLE PRECISION :: DTPIC_MIN_X, DTPIC_MIN_Y, DTPIC_MIN_Z
334           double precision, allocatable, dimension(:,:) ::  rand_vel
335     
336           double precision :: norm1, norm2, norm3
337           Logical :: OUTER_STABILITY_COND, DES_FIXED_BED
338     !-----------------------------------------------
339     
340           OUTER_STABILITY_COND = .false.
341           DES_FIXED_BED = .false.
342           PC = 1
343           DTPIC_CFL = LARGE_NUMBER
344           DTPIC_MIN_X =  LARGE_NUMBER
345           DTPIC_MIN_Y =  LARGE_NUMBER
346           DTPIC_MIN_Z =  LARGE_NUMBER
347     
348           if(NO_K) THREEINTOSQRT2 = 2.D0*SQRT(2.D0)
349           if(DO_K) THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
350           THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
351           PIP_DEL_COUNT = 0
352     
353           !EPG_MIN2 = MINVAL(EP_G(:))
354           !epg_min_loc = MINLOC(EP_G)
355           !IJK = epg_min_loc(1)
356           allocate(rand_vel(3, MAX_PIP))
357           do LC = 1, merge(2,3,NO_K)
358              mean_u = zero
359              sig_u = 1.d0
360              CALL NOR_RNO(RAND_VEL(LC, 1:MAX_PIP), MEAN_U, SIG_U)
361           enddo
362           !WRITE(*, '(A,2x,3(g17.8))') 'FC MIN AND MAX IN Z = ', MINVAL(FC(3,:)), MAXVAL(FC(3,:))
363           !WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'FC MIN AND MAX IN Z = ', MINVAL(FC(3,:)), MAXVAL(FC(3,:))
364     
365           DO NP = 1, MAX_PIP
366              DELETE_PART = .false.
367     ! pradeep skip ghost particles
368              if(pc.gt.pip) exit
369              if(is_nonexistent(NP)) cycle
370              pc = pc+1
371              if(is_ghost(NP) .or. is_entering_ghost(NP) .or. is_exiting_ghost(NP)) cycle
372     
373              DES_LOC_DEBUG = .FALSE.
374     
375              IF(.NOT.IS_ENTERING(NP))THEN
376                 FC(NP,:) = FC(NP,:)/PMASS(NP) + GRAV(:)
377              ELSE
378                 FC(NP,:) = ZERO
379              ENDIF
380     
381              IF(DES_FIXED_BED) FC(NP,:) = ZERO
382     
383              !DP_BAR is the D_p in the Snider paper (JCP, vol 170, 523-549, 2001)
384              !By comparing the MFIX and equations in those papers,
385              !D_p = Beta/(EP_S*RHOP)
386              !F_gp in drag_fgs.f  = Beta*PVOL/EP_S
387              !Therefore, D_p = F_gp/(PVOL*RHOP) = F_gp/PMASS
388              DP_BAR = F_gp(NP)/(PMASS(NP))
389              IF(.NOT.DES_CONTINUUM_COUPLED) DP_BAR = ZERO
390     
391              if(.not.MPPIC_PDRAG_IMPLICIT) DP_BAR = ZERO
392              M = PIJK(NP,5)
393              IJK = PIJK(NP,4)
394              IJK_OLD = IJK
395              PIJK_OLD(:) = PIJK(NP,:)
396     
397              I = I_OF(IJK)
398              J = J_OF(IJK)
399              K = K_OF(IJK)
400     
401              DES_VEL_NEW(NP,:) = (DES_VEL_NEW(NP,:) + &
402              & FC(NP,:)*DTSOLID)/(1.d0+DP_BAR*DTSOLID)
403     
404              IF(DES_FIXED_BED) DES_VEL_NEW(NP,:) = ZERO
405     
406              !MPPIC_VPTAU(NP,:) = DES_VEL_NEW(NP,:)
407     
408              VELP_INT(:) = DES_VEL_NEW(NP,:)
409     
410              MEANVEL(1) = U_S(IJK_OLD,M)
411              MEANVEL(2) = V_S(IJK_OLD,M)
412              IF(DO_K) MEANVEL(3) = W_S(IJK_OLD,M)
413     
414              RAD_EFF = DES_RADIUS(NP)
415              !RAD_EFF = (DES_STAT_WT(NP)**(1.d0/3.d0))*DES_RADIUS(NP)
416              IF(.not.DES_ONEWAY_COUPLED) then
417                 MEAN_FREE_PATH  = MAX(1.d0/(1.d0-EP_STAR), 1.d0/(1.D0-EP_G(IJK_OLD)))
418                 MEAN_FREE_PATH  = MEAN_FREE_PATH*RAD_EFF/THREEINTOSQRT2
419              endif
420     
421              DO LC = 1, merge(2,3,NO_K)
422                 !SIG_U = 0.05D0*MEANVEL(LC)
423                 !DES_VEL_NEW(NP,LC) = DES_VEL_NEW(NP,LC) + SIG_U*RAND_VEL(LC, L )
424                 !PART_TAUP = RO_Sol(NP)*((2.d0*DES_RADIUS(NP))**2.d0)/(18.d0* MU_G(IJK))
425                 SIG_U = 0.005D0      !*MEANVEL(LC)
426                 RAND_VEL(LC, NP)  = SIG_U*RAND_VEL(LC, NP)*DES_VEL_NEW(NP,LC)
427                 IF(DES_FIXED_BED) RAND_VEL(LC,NP) = ZERO
428                 !rand_vel(LC, NP)  = sig_u*rand_vel(LC, NP)/part_taup
429                 !rand_vel(LC, NP)  = sig_u* mean_free_path*rand_vel(LC, NP)/part_taup
430                 DES_VEL_NEW(NP,LC) = DES_VEL_NEW(NP,LC) + rand_vel(LC, NP)
431              enddo
432     
433              IF(.not.DES_ONEWAY_COUPLED.and.(.not.des_fixed_bed)) CALL MPPIC_APPLY_PS_GRAD_PART(NP)
434     
435              UPRIMEMOD = SQRT(DOT_PRODUCT(DES_VEL_NEW(NP,1:), DES_VEL_NEW(NP,1:)))
436     
437              !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) then
438              !   DES_VEL_NEW(NP,:) = (DES_VEL_NEW(NP,:)/UPRIMEMOD)*MEAN_FREE_PATH/DTSOLID
439              !ENDIF
440     
441              IF(DES_FIXED_BED) THEN
442                 DD(:) = ZERO
443              ELSE
444                 DD(:) = DES_VEL_NEW(NP,:)*DTSOLID !+ rand_vel(:, NP)*dtsolid
445              ENDIF
446     
447              CALL PIC_FIND_NEW_CELL(NP)
448     
449              IJK = PIJK(NP,4)
450     
451              !IF((EP_G(IJK).LT.0.35.and.fluid_at(ijk)).or.(ep_g(ijk_old).lt.0.35)) then !.and.(ijk.ne.ijk_old)) then
452              !IF((EP_G(IJK).LT.EP_STAR.and.fluid_at(ijk)).and.(ijk.ne.ijk_old)) then
453              INSIDE_DOMAIN = .true.
454              INSIDE_DOMAIN = FLUID_AT(IJK)
455     
456              IF(CUT_CELL_AT(IJK)) THEN
457                 POS_Z = zero
458                 IF(DO_K) POS_Z = DES_POS_NEW(NP,3)
459                 CALL GET_DEL_H_DES(IJK,'SCALAR', &
460                 & DES_POS_NEW(NP,1),  DES_POS_NEW(NP,2), &
461                 & POS_Z, &
462                 & DIST, NORM1, NORM2, NORM3, .true.)
463     
464                 IF(DIST.LE.ZERO) INSIDE_DOMAIN = .false.
465              ENDIF
466     
467              !IF((EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN)) then
468              !IF(EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN) then
469              IF(1.d0 - EP_G(IJK).GT. 1.3d0*(1.d0 - EP_STAR).and.INSIDE_DOMAIN.and.IJK.NE.IJK_OLD.and.OUTER_STABILITY_COND) then
470                 DD(:) = ZERO
471                 DES_VEL_NEW(NP,:) = 0.8d0*DES_VEL_NEW(NP,:)
472              ENDIF
473     
474              PIJK(NP,:) = PIJK_OLD(:)
475     
476              DES_POS_NEW(NP,:) = DES_POS_NEW(NP,:) + DD(:)
477     
478              DIST = SQRT(DOT_PRODUCT(DD,DD))
479     
480              D_GRIDUNITS(1) = ABS(DD(1))/DX(PIJK(NP,1))
481              D_GRIDUNITS(2) = ABS(DD(2))/DY(PIJK(NP,2))
482              D_GRIDUNITS(3) = 1.d0
483              IF(DO_K) D_GRIDUNITS(3) = ABS(DD(3))/DZ(PIJK(NP,3))
484     
485              DIST = SQRT(DOT_PRODUCT(DD,DD))
486              DTPIC_TMPX = (CFL_PIC*DX(PIJK(NP,1)))/(ABS(DES_VEL_NEW(NP,1))+SMALL_NUMBER)
487              DTPIC_TMPY = (CFL_PIC*DY(PIJK(NP,2)))/(ABS(DES_VEL_NEW(NP,2))+SMALL_NUMBER)
488              DTPIC_TMPZ = LARGE_NUMBER
489              IF(DO_K) DTPIC_TMPZ = (CFL_PIC*DZ(PIJK(NP,3)))/(ABS(DES_VEL_NEW(NP,3))+SMALL_NUMBER)
490     
491     
492     ! Check if the particle has moved a distance greater than or equal to grid spacing
493     ! if so, then exit
494     
495              DO LC = 1, merge(2,3,NO_K)
496                 IF(D_GRIDUNITS(LC).GT.ONE) THEN
497                    IF(DMP_LOG.OR.myPE.eq.pe_IO) THEN
498     
499                       WRITE(UNIT_LOG, 2001) NP, D_GRIDUNITS(:), DES_VEL_NEW(NP,:)
500                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'rand_vel = ', rand_vel(:,NP)
501                       IF (DO_OLD) WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'des_pos_old = ', des_pos_old(NP,:)
502                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'des_pos_new = ', des_pos_new(NP,:)
503                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'FC          = ', FC(NP,:)
504     
505                       WRITE(*, 2001) NP, D_GRIDUNITS(:), DES_VEL_NEW(NP,:)
506     
507                       WRITE(*, '(A,2x,3(g17.8))') 'rand_vel = ', rand_vel(:,NP)
508                       IF (DO_OLD) WRITE(*, '(A,2x,3(g17.8))') 'des_pos_old = ', des_pos_old(NP,:)
509                       WRITE(*, '(A,2x,3(g17.8))') 'des_pos_new = ', des_pos_new(NP,:)
510                       WRITE(*, '(A,2x,3(g17.8))') 'FC          = ', FC(NP,:)
511     !                  read(*,*)
512                       DELETE_PART = .true.
513     
514                    ENDIF
515                    !CALL mfix_exit(myPE)
516                 END IF
517     
518              END DO
519     
520              IF(.not.DELETE_PART) then
521                 DTPIC_CFL = MIN(DTPIC_CFL, DTPIC_TMPX, DTPIC_TMPY, DTPIC_TMPZ)
522                 DTPIC_MIN_X = MIN(DTPIC_MIN_X, DTPIC_TMPX)
523                 DTPIC_MIN_Y = MIN(DTPIC_MIN_Y, DTPIC_TMPY)
524                 DTPIC_MIN_Z = MIN(DTPIC_MIN_Z, DTPIC_TMPZ)
525              ENDIF
526              FC(NP,:) = ZERO
527     
528              IF(DELETE_PART) THEN
529                 CALL SET_NONEXISTENT(NP)
530                 PIP_DEL_COUNT = PIP_DEL_COUNT + 1
531              ENDIF
532              IF (DES_LOC_DEBUG) WRITE(*,1001)
533           ENDDO
534     
535     
536           DEALLOCATE(RAND_VEL)
537     
538           CALL global_all_max(DTPIC_CFL)
539           PIP = PIP - PIP_DEL_COUNT
540     
541           LPIP_DEL_COUNT_ALL(:) = 0
542           LPIP_DEL_COUNT_ALL(mype) = PIP_DEL_COUNT
543           CALL GLOBAL_ALL_SUM(LPIP_DEL_COUNT_ALL)
544           IF((DMP_LOG).AND.SUM(LPIP_DEL_COUNT_ALL(:)).GT.0) THEN
545              IF(PRINT_DES_SCREEN) THEN
546                 WRITE(*,'(/,2x,A,2x,i10,/,A)') &
547                 'TOTAL NUMBER OF PARTS STEPPING MORE THAN ONE GRID SPACE = ',&
548                  SUM(LPIP_DEL_COUNT_ALL(:)), &
549                 'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
550              ENDIF
551     
552              WRITE(UNIT_LOG,'(/,2x,A,2x,i10,/,A)') &
553                 'TOTAL NUMBER OF PARTS  STEPPING MORE THAN ONE GRID SPACEC=',&
554                 SUM(LPIP_DEL_COUNT_ALL(:)), &
555                 'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
556              !DO IPROC = 0, NUMPES-1
557              !   WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
558              !     'PARTICLES OUTSIDE DOMAIN (PIC)  ON PROC:', IPROC,' EQUAL TO', LPIP_DEL_COUNT_ALL(IPROC)
559              !ENDDO
560     
561           ENDIF
562     
563           DTPIC_MAX = MIN(DTPIC_CFL, DTPIC_TAUP)
564     
565           RETURN
566           IF(DTSOLID.GT.DTPIC_MAX) THEN
567              !IF(DMP_LOG) WRITE(UNIT_LOG, 2001) DTSOLID
568              IF(myPE.eq.pe_IO) WRITE(*, 2004) DTSOLID
569           ELSEIF(DTSOLID.LT.DTPIC_MAX) THEN
570     
571              !IF(DMP_LOG) WRITE(UNIT_LOG, 2002) DTSOLID
572              IF(myPE.eq.pe_IO) WRITE(*, 2002) DTPIC_MAX
573           ELSE
574             !WRITE(*,'(A40,2x,g17.8)') 'DT
575             !IF(mype.eq.pe_IO) WRITE(*,2003) DTSOLID
576           END IF
577     
578           WRITE(UNIT_LOG, '(10x,A,2x,3(g17.8))')&
579              'DTPIC MINS IN EACH DIRECTION = ', DTPIC_MIN_X, &
580              DTPIC_MIN_Y, DTPIC_MIN_Z
581           WRITE(*, '(10x,A,2x,3(g17.8))')'DTPIC MINS IN EACH DIRECTION = ',&
582             DTPIC_MIN_X, DTPIC_MIN_Y, DTPIC_MIN_Z
583     
584      1001 FORMAT(3X,'<---------- END INTEGRATE_TIME_PIC ----------')
585     
586     2001  FORMAT(/1X,70('*'),//,10X,  &
587                & 'MOVEMENT UNDESIRED IN INTEGRATE_TIME_PIC: PARTICLE', i5, /,10X, &
588                & ' MOVED MORE THAN A GRID SPACING IN ONE TIME STEP', /,10X, &
589                & 'MOVEMENT IN GRID UNITS = ', 3(2x, g17.8),/,10X,  &
590                & 'TERMINAL ERROR: NOT STOPPING, BUT DELETING THE PARTICLE', &
591                & /1X,70('*'), /10X, &
592                & 'DES_VEL_NEW = ',  3(2x, g17.8))
593     
594      2004 FORMAT(/10x, &
595           & 'DTSOLID SHUD BE REDUCED TO', g17.8)
596     
597      2002 FORMAT(/10x, &
598           & 'DTSOLID CAN BE INCREASED TO', g17.8)
599     
600           RETURN
601           END SUBROUTINE INTEGRATE_TIME_PIC_GARG
602     
603     
604     !------------------------------------------------------------------------
605     ! subroutine       : des_dbgpic
606     ! Author           : Pradeep G.
607     ! Purpose          : For debugging the pic values
608     ! Parameters       : pstart - start indices of the particle
609     !                    pend - end indices of the particle
610     !                    pfreq - optional frequency (when the local count matches the
611     !                    frequency the filw will be written)
612     !                    if not send then it prints the file
613     !------------------------------------------------------------------------
614     
615           subroutine des_dbgpic (pstart,pend,pfreq)
616     
617     !-----------------------------------------------
618     ! Modules
619     !-----------------------------------------------
620           use discretelement
621           USE fldvar
622           use physprop, only: mmax
623           use functions
624           implicit none
625     !-----------------------------------------------
626     ! Dummy arguments
627     !-----------------------------------------------
628           INTEGER, INTENT(IN) :: pstart,pend
629           INTEGER, INTENT(IN), OPTIONAL :: pfreq
630     !-----------------------------------------------
631     ! Local variables
632     !-----------------------------------------------
633           integer lp,lijk
634           integer, save :: lfcount = 0 ,lfreq =0
635           character(255) :: filename
636     !-----------------------------------------------
637           if (present(pfreq)) then
638              lfreq = lfreq+1
639              if (lfreq .ne. pfreq) return
640              lfreq =0
641           end if
642           lfcount = lfcount + 1
643           write(filename,'("debug",I3.3)') lfcount
644           open (unit=100,file=filename)
645           do lp = pstart,pend
646              if (is_normal(lp) .or. is_entering(lp) .or. is_exiting(lp)) then
647                 lijk = pijk(lp,4)
648                 write(100,*)"positon =",lijk,pijk(lp,1),pijk(lp,2), &
649                    pijk(lp,3),ep_g(lijk),U_s(lijk,MMAX+1)
650                 write(100,*)"forces =", FC(lp,2),tow(lp,1)
651              end if
652           end do
653           close (100)
654     
655           RETURN
656           END SUBROUTINE des_dbgpic
657     
658     
659     !------------------------------------------------------------------------
660     ! subroutine       : des_dbgtecplot
661     ! Author           : Pradeep G.
662     ! Purpose          : prints the tecplot file for particle location
663     ! Parameters       : pstart - start indices of the particle
664     !                    pend - end indices of the particle
665     !                    pfreq - optional frequency (when the local count matches the
666     !                    frequency the filw will be written)
667     !                    if not send then it prints the file
668     !------------------------------------------------------------------------
669     
670           subroutine des_dbgtecplot (pstart,pend,pfreq)
671     
672     !-----------------------------------------------
673     ! Modules
674     !-----------------------------------------------
675           use discretelement
676           USE fldvar
677           use functions
678           implicit none
679     !-----------------------------------------------
680     ! Dummy arguments
681     !-----------------------------------------------
682           INTEGER, INTENT(IN) :: pstart,pend
683           INTEGER, INTENT(IN), OPTIONAL :: pfreq
684     !-----------------------------------------------
685     ! Local variables
686     !-----------------------------------------------
687           integer lp,lijk
688           integer, save :: lfcount = 0 ,lfreq =0
689           character(255) :: filename
690     !-----------------------------------------------
691     
692           if (present(pfreq)) then
693              lfreq = lfreq+1
694              if (lfreq .ne. pfreq) return
695              lfreq =0
696           end if
697           lfcount = lfcount + 1
698           write(filename,'("new_tec",I3.3,".dat")') lfcount
699           open (unit=100,file=filename)
700           write(100,'(9(A,3X),A)') &
701              'VARIABLES = ', '"ijk"', '"x"', '"y"', '"vx"', '"vy"', &
702              '"ep_g"', '"FCX"' ,'"FCY"', '"TOW"'
703           write(100,'(A,F14.7,A)') 'zone T = "' , s_time , '"'
704           do lp = pstart,pend
705              if (.not.is_nonexistent(lp)) then
706                 lijk = pijk(lp,4)
707                 write(100,*)lijk,des_pos_new(lp,1),des_pos_new(lp,2), &
708                    des_vel_new(lp,1),des_vel_new(lp,2),ep_g(lijk),&
709                    fc(lp,1),fc(lp,2),tow(1,lp)
710              endif
711           enddo
712           close (100)
713           RETURN
714           END SUBROUTINE DES_DBGTECPLOT
715