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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Subroutine: PIC_TIME_MARCH                                          !
3     !  Author: R. Garg                                                     !
4     !                                                                      !
5     !  Purpose: Main PIC driver routine.                                   !
6     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
7           SUBROUTINE PIC_TIME_MARCH
8     
9           USE param
10           USE param1
11           USE run
12           USE output
13           USE funits
14           USE discretelement
15           USE constant
16           USE sendrecv
17           USE des_bc
18           USE cutcell
19           USE mfix_pic
20           USE pic_bc
21           USE error_manager
22           USE fldvar, only: P_g
23     
24           use mpi_utility, only: GLOBAL_ALL_SUM
25           use mpi_funs_des, only: DES_PAR_EXCHANGE
26     
27           USE fun_avg
28           USE functions
29     
30           IMPLICIT NONE
31     
32     ! Local variables
33     !---------------------------------------------------------------------//
34           INTEGER IJK, I, J, K
35           INTEGER NN, FACTOR, NP, BCV_I, LL, PC, TIME_LOOP_COUNT
36     
37     ! Local variables to keep track of time when dem restart and des
38     ! write data need to be written when des_continuum_coupled is F
39           DOUBLE PRECISION DES_RES_TIME, DES_SPX_TIME
40     
41     ! Temporary variables when des_continuum_coupled is T to track
42     ! changes in solid time step
43           DOUBLE PRECISION TMP_DTS, DTSOLID_TMP
44           CHARACTER(LEN=5) :: FILENAME
45     
46     
47     ! Logical to see whether this is the first entry to this routine
48           LOGICAL,SAVE:: FIRST_PASS = .TRUE.
49     
50     ! temp
51           DOUBLE PRECISION xpos,ypos, zpos, NORM_CF(3), DIST
52     
53     
54           DOUBLE PRECISION :: DES_KE_VEC(DIMN)
55     
56     ! time till which the PIC loop will be run
57           double precision :: TEND_PIC_LOOP
58     ! number of PIC time steps
59           Integer :: PIC_ITERS
60     ! Identifies that the indicated particle is of interest for debugging
61           LOGICAL FOCUS
62     ! Global number of parcels.
63           INTEGER :: gPIP
64     
65     
66     
67     
68           CALL INIT_ERR_MSG("PIC_TIME_MARCH")
69     
70           S_TIME = TIME
71           TIME_LOOP_COUNT = 0
72     
73     ! compute the gas-phase pressure gradient at the beginning of the
74     ! des loop as the gas-phase pressure field will not change during
75     ! des calls
76           IF(DES_CONTINUUM_COUPLED) CALL CALC_PG_GRAD
77     
78           TEND_PIC_LOOP = MERGE(TIME+DT, TSTOP, DES_CONTINUUM_COUPLED)
79           PIC_ITERS = 0
80     
81     ! If the current time in the discrete loop exceeds the current time in
82     ! the continuum simulation, exit the lagrangian loop
83           DO WHILE(S_TIME.LT.TEND_PIC_LOOP)
84     
85              !DTPIC_MAX = MIN( 1e-04, DTPIC_MAX)
86              DTSOLID = MERGE(MIN(DTPIC_MAX, DT), DTPIC_MAX, DES_CONTINUUM_COUPLED)
87              DTSOLID_ORIG  = DTSOLID
88              IF(MOD(PIC_ITERS, 10).eq.0) then
89                 IF(DES_CONTINUUM_COUPLED) then
90                    WRITE(ERR_MSG, 2000) DTSOLID, DTPIC_CFL, DTPIC_TAUP, DT
91                 ELSE
92                    WRITE(ERR_MSG, 2001) S_TIME, DTSOLID, DTPIC_CFL, DTPIC_TAUP, DT
93                 ENDIF
94                 CALL FLUSH_ERR_MSG(HEADER = .FALSE., FOOTER = .FALSE.)
95              ENDIF
96     
97      2000 FORMAT(/5x,'DTSOLID CURRENT  = ',g17.8,/5x,'DTPIC_CFL',8x,'= ',  &
98              g17.8, /5x,'DTPIC TAUP',7x,'= ',g17.8,/5x,'DT FLOW',10x,'= ', &
99              g17.8)
100     
101      2001 FORMAT(/5x,'TIME',13X,'= ',g17.8,/5x,'DTSOLID CURRENT  = ',g17.8,&
102              /5x,'DTPIC_CFL',8X,'= ', g17.8,/5x,'DTPIC TAUP',7x,'= ',g17.8,&
103              /5x,'DT FLOW',10X,'= ', g17.8)
104     
105     
106              PIC_ITERS  = PIC_ITERS + 1
107     
108     
109     ! If next time step in the discrete loop will exceed the current time
110     ! in the continuum simulation, modify the discrete time step so final
111     ! time will match
112              IF(S_TIME + DTSOLID > TEND_PIC_LOOP) THEN
113                 WRITE(ERR_MSG, 2010)  DTSOLID, TIME + DT - S_TIME
114                 CALL FLUSH_ERR_MSG(HEADER = .FALSE., FOOTER = .FALSE.)
115                 DTSOLID = TEND_PIC_LOOP - S_TIME
116              ENDIF
117     
118      2010 FORMAT(/5X,'REDUCING DTSOLID TO ENSURE STIME + DTSOLID LE ',     &
119           'TIME + DT', /5x,'DTSOLID ORIG         = ', g17.8, /5x,          &
120           'DTSOLID ACTUAL       = ', g17.8)
121     
122     ! exchange particle crossing boundaries
123              CALL DES_PAR_EXCHANGE
124     
125              CALL PARTICLES_IN_CELL
126     ! Calculate mean fields using either interpolation or cell averaging.
127              CALL COMP_MEAN_FIELDS
128     
129     ! This was moved from particles in cell and the passed variables should be
130     ! added to particles in cell or made global.
131              CALL REPORT_PIC_STATS
132     
133              CALL MPPIC_COMPUTE_PS_GRAD
134              IF(DES_CONTINUUM_COUPLED)   CALL CALC_DRAG_DES
135              IF (DO_OLD) CALL CFUPDATEOLD
136     
137              CALL CFNEWVALUES
138     
139     ! Impose the wall-particle boundary condition for pic case
140     ! The same routine also applies the mass inflow/outflow BCs as well
141              CALL PIC_APPLY_WALLBC_STL
142     
143     
144     ! Update time to reflect changes
145              S_TIME = S_TIME + DTSOLID
146     
147     
148              DTPIC_MAX = MIN(DTPIC_CFL, DTPIC_TAUP)
149     
150     ! When coupled, all write calls are made in time_march (the continuum
151     ! portion) according to user settings for spx_time and res_time.
152     ! The following section targets data writes for DEM only cases:
153              IF(.NOT.DES_CONTINUUM_COUPLED) THEN
154     ! Keep track of TIME for DEM simulations
155                 TIME = S_TIME
156     
157     ! Write data using des_spx_time and des_res_time; note the time will
158     ! reflect current position of particles
159                 IF(PRINT_DES_DATA) THEN
160                    IF ( (S_TIME+0.1d0*DTSOLID >= DES_SPX_TIME) .OR. &
161                         (S_TIME+0.1d0*DTSOLID >= TSTOP)) then
162                       DES_SPX_TIME = &
163                          ( INT((S_TIME+0.1d0*DTSOLID)/DES_SPX_DT) &
164                          + 1 )*DES_SPX_DT
165                       CALL WRITE_DES_DATA
166                       IF(DMP_LOG) WRITE(UNIT_LOG,'(3X,A,1X,ES15.5)') &
167                          'DES data file written at time =', S_TIME
168                    ENDIF
169                 ENDIF
170     
171                 IF ( (S_TIME+0.1d0*DTSOLID >= DES_RES_TIME) .OR. &
172                      (S_TIME+0.1d0*DTSOLID >= TSTOP) .OR. &
173                      (NN == FACTOR) ) THEN
174                    DES_RES_TIME = &
175                       ( INT((S_TIME+0.1d0*DTSOLID)/DES_RES_DT) &
176                       + 1 )*DES_RES_DT
177                       CALL WRITE_RES0_DES
178     ! Write RES1 here since it won't be called in time_march.  This will
179     ! also keep track of TIME
180                    CALL WRITE_RES1
181                    IF(DMP_LOG) WRITE(UNIT_LOG,'(3X,A,1X,ES15.5)') &
182                    'DES.RES and .RES files written at time =', S_TIME
183                 ENDIF
184              ENDIF  ! end if (.not.des_continuum_coupled)
185     
186     
187     
188     
189           ENDDO
190     
191           CALL GLOBAL_ALL_SUM(PIP, gPIP)
192           WRITE(ERR_MSG, 3000) trim(iVal(PIC_ITERS)), trim(iVal(gPIP))
193           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
194     
195      3000 FORMAT(/'PIC NITs: ',A,3x,'Total PIP: ', A)
196     
197     !      IJK_BOT = funijk(imin1, 2,kmin1)
198     !      IJK_TOP = funijk(imin1, jmax1, kmin1)
199     !      WRITE(*,'(/5X, A, 3(2x,g17.8))') &
200     !         'MPPIC: PRES BOTTOM, TOP, AND DIFF KPA', P_G(IJK_BOT)/10000.d0, &
201     !         P_G(IJK_TOP)/10000.d0, (P_G(IJK_BOT) -  P_G(IJK_TOP))/10000.d0
202           !WRITE(*,'(/5X, A, 3(2x,g17.8))') &
203           !'PRES BOTTOM, TOP, AND DIFF ', P_G(IJK_BOT), P_G(IJK_TOP), P_G(IJK_BOT) -  P_G(IJK_TOP)
204           IF(.NOT.DES_CONTINUUM_COUPLED)then
205              if(dmp_log)write(unit_log,'(1X,A)')&
206              '<---------- END MPPIC_TIME_MARCH ----------'
207     !Pradeep call send recv for variables
208           else
209              !call send_recv(ep_g,2)
210              !call send_recv(rop_g,2)
211              !call send_recv(des_u_s,2)
212              !call send_recv(des_v_s,2)
213              !if(dimn.eq.3) call send_recv(des_w_s,2)
214              !call send_recv(rop_s,2)
215     ! now the above are communicated in comp_mean_fields_interp itself.
216     ! so no need to communicate them here.
217           END IF
218     
219     
220           CALL FINL_ERR_MSG
221         end SUBROUTINE PIC_TIME_MARCH
222     
223     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
224     !  Subroutine: MPPIC_COMP_EULERIAN_VELS_NON_CG                         !
225     !  Author: R. Garg                                                     !
226     !                                                                      !
227     !  Purpose:                                                            !
228     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
229           SUBROUTINE MPPIC_COMP_EULERIAN_VELS_NON_CG
230             USE compar
231             USE constant
232             USE discretelement
233             USE functions
234             USE geometry
235             USE indices
236             USE mfix_pic
237             USE parallel
238             USE param
239             USE param1
240             use desmpi
241             IMPLICIT NONE
242     !-----------------------------------------------
243     ! Local variables
244     !-----------------------------------------------
245     ! local variable used for debugging
246           LOGICAL :: FOCUS
247     ! general i, j, k indices
248           INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,&
249                      IPJPK, IPJKP, IJPKP, IPJPKP, II, JJ, KK
250     
251     ! index of solid phase that particle NP belongs to
252           INTEGER :: M
253     
254           DO IJK = ijkstart3, ijkend3
255                 I = I_OF(IJK)
256                 J = J_OF(IJK)
257                 K = K_OF(IJK)
258                 !U_so(IJK, :) = U_s(IJK, :)
259                 !V_so(IJK, :) = V_s(IJK, :)
260                 !W_so(IJK, :) = W_s(IJK, :)
261                 PIC_U_S(IJK, :) = ZERO
262                 PIC_V_S(IJK, :) = ZERO
263                 PIC_W_S(IJK, :) = ZERO
264                 IF (.NOT.IS_ON_myPE_plus1layer(I,J,K)) CYCLE
265                 IF(.NOT.FLUID_AT(IJK)) CYCLE
266     
267                 IF(I.GE.IMIN1.AND.I.LT.IMAX1) then !because we don't want to
268                    !calculate solids velocity at the wall cells.
269                    IPJK = IP_OF(IJK)
270                    DO M = 1, DES_MMAX
271                       PIC_U_S(IJK,M) = 0.5d0*(DES_U_S(IJK, M) + DES_U_S(IPJK,M))
272                    ENDDO
273                 ENDIF
274     
275                 if(J.GE.JMIN1.AND.J.LT.JMAX1) then !because we don't want to
276                    !calculate solids velocity at the wall cells.
277                    IJPK = JP_OF(IJK)
278                    DO M = 1, DES_MMAX
279                       PIC_V_S(IJK,M) = 0.5d0*(DES_V_S(IJK, M) + DES_V_S(IJPK,M))
280                    ENDDO
281                 ENDIF
282     
283     
284                 if(K.GE.KMIN1.AND.K.LT.KMAX1.AND.DO_K) then !because we don't want to
285                    !calculate solids velocity at the wall cells.
286                    IJKP = KP_OF(IJK)
287                    DO M = 1, DES_MMAX
288                       PIC_W_S(IJK,M) = 0.5d0*(DES_W_S(IJK, M) + DES_W_S(IJKP,M))
289                    ENDDO
290                 ENDIF
291              ENDDO
292              !CALL SET_WALL_BC(IER)
293              !the above routine will apply noslip or free slip BC as per the mfix  convention.
294              !currently, this implies NSW or FSW wall BC's will be re-applied to gas-phase
295              !field as well. This can be changed later on to be more specific to MPPIC case
296              !CALL WRITE_MPPIC_VEL_S
297              CALL MPPIC_BC_U_S
298              CALL MPPIC_BC_V_S
299              IF(DO_K) CALL MPPIC_BC_W_S
300     
301           END SUBROUTINE MPPIC_COMP_EULERIAN_VELS_NON_CG
302     
303     
304     
305     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
306     !  Subroutine: MPPIC_COMP+EULERIAN_VELS_CG                             !
307     !  Author: R. Garg                                                     !
308     !                                                                      !
309     !  Puryyse:                                                            !
310     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
311           SUBROUTINE MPPIC_COMP_EULERIAN_VELS_CG
312             USE compar
313             USE constant
314             USE cutcell
315             USE desmpi
316             USE discretelement
317             USE functions
318             USE geometry
319             USE indices
320             USE mfix_pic
321             USE parallel
322             USE param
323             USE param1
324             IMPLICIT NONE
325     !-----------------------------------------------
326     ! Local variables
327     !-----------------------------------------------
328     ! local variable used for debugging
329           LOGICAL :: FOCUS
330     ! general i, j, k indices
331           INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,&
332                      IPJPK, IPJKP, IJPKP, IPJPKP, II, JJ, KK
333     
334     ! index of solid phase that particle NP belongs to
335           INTEGER :: M
336     
337           DO IJK = ijkstart3, ijkend3
338              I = I_OF(IJK)
339              J = J_OF(IJK)
340              K = K_OF(IJK)
341              !U_so(IJK, :) = U_s(IJK, :)
342              !V_so(IJK, :) = V_s(IJK, :)
343              !W_so(IJK, :) = W_s(IJK, :)
344              PIC_U_S(IJK, :) = ZERO
345              PIC_V_S(IJK, :) = ZERO
346              PIC_W_S(IJK, :) = ZERO
347     
348              IF (.NOT.IS_ON_myPE_plus1layer(I,J,K)) CYCLE
349     
350              IF(WALL_U_AT(IJK)) THEN
351                 PIC_U_S(IJK, :) = ZERO
352                 !currently only No slip BC is being set on this mean
353                 !solid's velocity field. Later this wall part can be
354                 !treated separately and U_S set only for scalar cells
355                 !where FLUID_AT(IJK) is true.
356              ELSE
357                 if(.not.FLUID_AT(IJK)) cycle
358     
359                 IPJK = IP_OF(IJK)
360                 IF(FLUID_AT(IPJK)) THEN
361                    DO M = 1, DES_MMAX
362                       PIC_U_S(IJK,M) = 0.5d0*(DES_U_S(IJK, M) + DES_U_S(IPJK,M))
363                    ENDDO
364                 ELSE
365                    PIC_U_S(IJK,:) = DES_U_S(IJK, :)
366                 ENDIF
367              ENDIF
368     
369              IF(WALL_V_AT(IJK)) THEN
370                 PIC_V_S(IJK, :) = ZERO
371              ELSE
372                 if(.not.FLUID_AT(IJK)) cycle
373                 IJPK = JP_OF(IJK)
374                 IF(FLUID_AT(IJPK)) THEN
375                    DO M = 1, DES_MMAX
376                       PIC_V_S(IJK,M) = 0.5d0*(DES_V_S(IJK, M) + DES_V_S(IJPK,M))
377                    ENDDO
378                 ELSE
379                    PIC_V_S(IJK,:) = DES_V_S(IJK, :)
380                 ENDIF
381              ENDIF
382     
383              IF(DO_K) THEN
384                 IF(WALL_W_AT(IJK)) THEN
385                    PIC_W_S(IJK, :) = ZERO
386                 ELSE
387                    if(.not.FLUID_AT(IJK)) cycle
388                    IJKP = KP_OF(IJK)
389                    IF(FLUID_AT(IJKP)) THEN
390                       DO M = 1, DES_MMAX
391                          PIC_W_S(IJK,M) = 0.5d0*(DES_W_S(IJK, M) + DES_W_S(IJKP,M))
392                       ENDDO
393                    ELSE
394                       PIC_W_S(IJK,:) = DES_W_S(IJK, :)
395                    ENDIF
396                 ENDIF
397              ENDIF
398           ENDDO
399     
400           END SUBROUTINE MPPIC_COMP_EULERIAN_VELS_CG
401     
402     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
403     !  Subroutine: MPPIC_APPLY_PS_GRAD_PART                                !
404     !  Author: R. Garg                                                     !
405     !                                                                      !
406     !  Purpose:                                                            !
407     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
408           SUBROUTINE MPPIC_APPLY_PS_GRAD_PART(L)
409     
410           USE param
411           USE param1
412           USE parallel
413           USE constant
414           USE discretelement
415           USE mpi_utility
416           USE mfix_pic
417           USE cutcell
418           USE fldvar, only: ep_g
419           USE fun_avg
420           USE functions
421           IMPLICIT NONE
422           INTEGER, INTENT(IN) :: L
423     !-----------------------------------------------
424     ! Local Variables
425     !-----------------------------------------------
426           INTEGER M, IDIM
427           INTEGER I, J, K, IJK, IJK_C
428     
429           DOUBLE PRECISION D(DIMN), DIST,  DP_BAR, COEFF_EN, COEFF_EN2
430     
431           DOUBLE PRECISION DELUP(DIMN), UPRIMETAU(DIMN), UPRIMETAU_INT(DIMN), PS_FORCE(DIMN), VEL_ORIG(DIMN)
432     ! index to track accounted for particles
433           INTEGER PC , epg_min_loc(1)
434     
435     ! Logical for local debug warnings
436           LOGICAL DES_LOC_DEBUG
437     
438     ! maximum distance particles can move in MPPIC
439           DOUBLE PRECISION MAXDIST_PIC, UPRIMEMOD, UPRIMEMODNEW, signvel
440     
441     ! dt's in each direction  based on cfl_pic for the mppic case
442     
443           DOUBLE PRECISION  MEANUS(DIMN, DES_MMAX), RELVEL(DIMN), MEANVEL(DIMN), VEL_NEW(DIMN)
444     !      INTEGER :: TOT_CASE, case1_count, case2_count, case3_count, case4_count
445     
446           LOGICAL :: INSIDE_DOMAIN
447     !-----------------------------------------------
448     
449           M = PIJK(L,5)
450           IJK = PIJK(L,4)
451           COEFF_EN  = MPPIC_COEFF_EN1
452           COEFF_EN2 = MPPIC_COEFF_EN2
453     
454           VEL_ORIG(1:DIMN) = DES_VEL_NEW(:,L)
455           VEL_NEW (1:DIMN) = DES_VEL_NEW(:,L)
456           !IF(L.eq.1) WRITE(*,*) 'MPPIC COEFFS = ', COEFF_EN, COEFF_EN2
457           IF(L.EQ.FOCUS_PARTICLE) THEN
458     
459              WRITE(*,'(A20,2x,3(2x,i4))') 'CELL ID = ', PIJK(L,1:3)
460              WRITE(*,'(A20,2x,3(2x,g17.8))') 'EPS = ', 1.d0 - EP_g(PIJK(L,4))
461              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL ORIG = ', DES_VEL_NEW(:,L)
462     
463              WRITE(*,'(A20,2x,3(2x,g17.8))') 'FC = ', FC(:,L)
464           ENDIF
465     
466           MEANVEL(1) = DES_U_S(IJK,M)
467           MEANVEL(2) = DES_V_S(IJK,M)
468           IF(DO_K) MEANVEL(3) = DES_W_S(IJK,M)
469     
470           PS_FORCE(:) = PS_GRAD(L, :)
471           !IF(ABS(PS_FORCE(2)).GT.ZERO)  WRITE(*,*) 'PS_FORCE = ', PS_FORCE
472           DELUP(:) = -PS_FORCE(:)
473     
474           MEANUS(:,M) =  AVGSOLVEL_P (:,L)
475           !MEANUS(:,M) = MEANVEL(:)
476           RELVEL(:) = DES_VEL_NEW(:,L) - MEANUS(:,M)
477     
478           !IF(EPg_P(L).gt.1.2d0*ep_star) RETURN
479     
480           DO IDIM = 1, DIMN
481     
482              IF(ABS(PS_FORCE(IDIM)).eq.zero) cycle
483     
484              IF(VEL_ORIG(IDIM)*MEANUS(IDIM,M).GT.ZERO) THEN
485     
486                 IF(VEL_ORIG(IDIM)*DELUP(IDIM).GT.ZERO) THEN
487     
488                    IF(ABS(MEANUS(IDIM,M)) .GT. ABS(VEL_ORIG(IDIM))) THEN
489                       IJK_C = IJK
490                       IF(IDIM.eq.1) then
491                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
492                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
493                       ELSEIF(IDIM.eq.2) then
494                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
495                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
496                       ELSEIF(IDIM.eq.3) then
497                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
498                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
499                       ENDIF
500                       INSIDE_DOMAIN = .false.
501                       INSIDE_DOMAIN = fluid_at(IJK_C)!and.(.not.cut_cell_at(IJK_C))
502     
503                       if(INSIDE_DOMAIN) then
504                          VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
505                       endif
506     !                  case4_count = case4_count + 1
507                    ELSE
508                       !do nothing
509                    ENDIF
510                 ELSE
511                    IF(ABS(VEL_ORIG(IDIM)).GT.ABS(MEANUS(IDIM,M))) then
512                       !VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
513                       IJK_C = IJK
514                       IF(IDIM.eq.1) then
515                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
516                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
517                       ELSEIF(IDIM.eq.2) then
518                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
519                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
520                       ELSEIF(IDIM.eq.3) then
521                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
522                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
523                       ENDIF
524                       !VEL_NEW(IDIM) = (MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM))
525     
526                       INSIDE_DOMAIN = .false.
527                       INSIDE_DOMAIN = fluid_at(IJK_C)!.and.(.not.cut_cell_at(IJK_C))
528     
529                       if(INSIDE_DOMAIN) then
530                          VEL_NEW(IDIM) = (MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM))
531                       else
532                          VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
533                       endif
534     
535                       IJK_C = IJK
536                       IF(IDIM.eq.1) then
537                          if(VEL_NEW(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
538                          if(VEL_NEW(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
539                       ELSEIF(IDIM.eq.2) then
540                          if(VEL_NEW(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
541                          if(VEL_NEW(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
542                       ELSEIF(IDIM.eq.3) then
543                          if(VEL_NEW(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
544                          if(VEL_NEW(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
545                       ENDIF
546                       INSIDE_DOMAIN = .false.
547                       INSIDE_DOMAIN = fluid_at(IJK_C) !.and.(.not.cut_cell_at(IJK_C))
548     
549                       !if(.not.INSIDE_DOMAIN) then
550                       !   VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
551                       !ENDIF
552     
553     !                  case1_count = case1_count + 1
554                    ELSE
555                       !do nothing
556                       VEL_NEW(IDIM) = COEFF_EN2 * VEL_ORIG(IDIM)
557                       !turning on the above would make the model uncondtionally stable
558     !                  case1_count = case1_count + 1
559     
560                    ENDIF
561                 ENDIF
562              ELSE
563                 IF(MEANUS(IDIM,M)*DELUP(IDIM).GT.ZERO) THEN
564                    VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
565     !               case2_count = case2_count + 1
566                 ELSE
567     !               case3_count = case3_count + 1
568                    !DO NOTHING
569                 ENDIF
570              ENDIF
571     
572              IF(MPPIC_GRAV_TREATMENT) THEN
573                 IF(DELUP(IDIM)*GRAV(IDIM).LT.ZERO.AND.VEL_ORIG(IDIM)*GRAV(IDIM).GT.ZERO) THEN
574                    VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
575                    !MPPIC_VPTAU(L, IDIM) = VEL_NEW(IDIM) -  DES_VEL_NEW(IDIM, L)
576                 ENDIF
577              ENDIF
578              !MPPIC_VPTAU(L, IDIM) = VEL_NEW(IDIM) -  DES_VEL_NEW(IDIM, L)
579              DES_VEL_NEW(IDIM, L) = VEL_NEW(IDIM)
580           ENDDO
581     
582              !
583           if(L.eq.FOCUS_PARTICLE) THEN
584              !iF((IJK.eq.epg_min_loc(1).or.IJK_OLD.eq.epg_min_loc(1)).and.epg_min2.lt.0.38) then
585              !if(j.ne.2) cycle
586     
587              WRITE(*,'(A20,2x,3(2x,i5))') 'PIJK I, J, K =', I_OF(IJK),J_OF(IJK),K_OF(IJK)
588              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL ORIG = ', VEL_ORIG(:)
589              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL NEW = ', DES_VEL_NEW(:,L)
590              WRITE(*,'(A20,2x,3(2x,g17.8))') 'MEANUS = ', MEANUS(:,1)
591              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_POS_NEW = ', DES_POS_NEW(:,L)
592              WRITE(*,'(A20,2x,3(2x,g17.8))') 'GRAD PS = ', PS_FORCE(:)
593              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DELUP =  ', DELUP(:)
594     
595              !WRITE(*,'(A20,2x,3(2x,g17.8))') 'UPRIMETAU_INT = ', UPRIMETAU_INT(:)
596              !WRITE(*,'(A20,2x,3(2x,g17.8))') 'UPRIMETAU = ', UPRIMETAU(:)
597              !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) WRITE(*,'(A20,2x,3(2x,g17.8))') 'U*DT, MFP =', UPRIMEMOD*DTSOLID, MEAN_FREE_PATH
598              read(*,*)
599           ENDIF
600     
601           !TOT_CASE = case1_count + case2_count + case3_count + case4_count
602           !IF(TOT_CASE.GT.0) THEN
603           !WRITE(*,'(A, 4(2x,i10))') 'CASE COUNT NUMBERS  = ', case1_count ,case2_count ,case3_count ,case4_count
604           !WRITE(*,'(A, 4(2x,g12.7))') 'CASE COUNT %AGE = ', real(case1_count)*100./real(tot_case), &
605           !     real(case2_count)*100./real(tot_case), real(case3_count)*100./real(tot_case), real(case4_count)*100./real(tot_case)
606           !ENDIF
607           RETURN
608     
609           !MEANUS(:,M) = MEANVEL(:)
610           !RELVEL(:) = DES_VEL_NEW(:,L) - MEANUS(:,M)
611           !DO IDIM = 1, DIMN
612           !    IF(ABS(PS_FORCE(IDIM)).eq.zero) cycle
613     
614           !   IF(RELVEL(IDIM)*DELUP(IDIM).GT.ZERO) THEN
615           !do nothing
616           !    ELSE
617           !       DES_VEL_NEW(IDIM,L) = MEANUS(IDIM,M) - 0.4d0*RELVEL(IDIM)
618     
619           !IF(DES_VEL_NEW(IDIM,L)*DELUP(IDIM).LT.ZERO) DES_VEL_NEW(IDIM,L) = -0.5d0*DES_VEL_NEW(IDIM,L)
620           !    ENDIF
621           ! ENDDO
622           ! CYCLE
623     
624     
625           END SUBROUTINE MPPIC_APPLY_PS_GRAD_PART
626     
627     
628     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
629     !  Subroutine: MPPIC_COMPUTE_PS_GRAD                                   !
630     !  Author: R. Garg                                                     !
631     !                                                                      !
632     !  Purpose:                                                            !
633     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
634           SUBROUTINE MPPIC_COMPUTE_PS_GRAD
635           USE param
636           USE param1
637           USE parallel
638           USE physprop
639           USE fldvar, only: ep_g, u_g, v_g, w_g
640           USE run
641           USE geometry
642           USE indices
643           USE bc
644           USE compar
645           USE sendrecv
646           USE discretelement
647           USE constant
648           USE cutcell
649           USE interpolation
650           USE mfix_pic
651           USE fun_avg
652           USE functions
653           implicit none
654     
655           ! general i, j, k indices
656           INTEGER I, J, K, IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM, IDIM, M
657     
658           ! temporary variables used to calculate pressure at scalar cell edge
659           DOUBLE PRECISION TEMP1, TEMP2, avg_factor, VOL_TOT_VEC(3), VOL_TOT_SCAL
660     
661           integer :: korder, iw,ie,js,jn,kb,ktp, onew, pcell(3), cur_ijk, NP, nindx
662     
663           integer :: ii,jj,kk, ipjpk, ijpkp, ipjkp, ipjpkp, I1, J1, K1
664     
665           double precision :: vol_ijk, vol_ipjk, vol_ijpk, vol_ipjpk
666           double precision :: vol_ijkp, vol_ipjkp, vol_ijpkp, vol_ipjpkp
667     
668           if(MPPIC_SOLID_STRESS_SNIDER) then
669     
670              DO IJK = IJKSTART3, IJKEND3
671                 IF(FLUID_AT(IJK)) THEN
672                    PIC_P_S(IJK,1) = PSFAC_FRIC_PIC * ((1.d0 - EP_G(IJK))**FRIC_EXP_PIC)/ MAX(EP_G(IJK) &
673                         - EP_STAR, FRIC_NON_SING_FAC*EP_G(IJK))
674                    !write(102,'(2(2x,i4),2(2x,g17.8))') I_OF(IJK), J_OF(IJK), EP_S(IJK,1), P_STAR(IJK)
675                 ELSE
676                    !So that ghost cells have higher pressure
677                    PIC_P_S(IJK,1) = PSFAC_FRIC_PIC * ((1.d0 - 0.1d0)**FRIC_EXP_PIC)/ MAX(0.1d0 - EP_STAR, FRIC_NON_SING_FAC*0.1d0)
678     
679                 ENDIF
680              ENDDO
681     
682     
683           ELSE
684              DO IJK = IJKSTART3, IJKEND3
685                 PS_FORCE_PIC(IJK,:) = ZERO
686     
687                 IF(FLUID_AT(IJK)) THEN
688                    if(EP_G(IJK).lt.ep_star) then
689                       PIC_P_S(IJK,1) = one*(one-ep_g(ijk))
690                    else
691     
692                       PIC_P_S(IJK,1) = ZERO
693                    endif
694     
695                 ELSE
696                    PIC_P_S(IJK,1) = 1.!\*0.d0
697                 ENDIF
698              ENDDO
699     
700           ENDIF
701     
702           CALL SEND_RECV(PIC_P_S,1)
703     
704           !Since EP_G is already shared across the processors, the pressure gradient calculation
705           !can be made a function call so that the extra communication of P_S can be avoided.
706     
707           !DO k = kstart2, kend1
708           !do j = jstart2, jend1
709           !do i = istart2, iend1
710           DO IJK = IJKSTART3, IJKEND3
711              PS_FORCE_PIC(IJK,:) = ZERO
712              IF(.NOT.FLUID_AT(IJK)) CYCLE
713     
714              I = I_OF(IJK)
715              J = J_OF(IJK)
716              K = K_OF(IJK)
717              IPJK = IP_OF(IJK)
718              IJPK = JP_OF(IJK)
719              IJKP = KP_OF(IJK)
720     
721              if(fluid_at(ipjk)) then
722                 PS_FORCE_PIC(IJK,1) = 2.d0*(PIC_P_S(IPJK,1) - PIC_P_S(IJK,1))/(DX(I)+DX(I_of(ipjk)))
723              else
724                 if(PIC_P_S(IJK,1).GT.ZERO) then
725                    !this will always be true for Snider's case
726                    PS_FORCE_PIC(IJK,1) = 2.d0*(PIC_P_S(IPJK,1) - PIC_P_S(IJK,1))/(DX(I)+DX(I_of(ipjk)))
727                 ELSE
728                    PS_FORCE_PIC(IJK,1) = zero
729                 endif
730              ENDIF
731     
732              IF(FLUID_AT(IJPK)) THEN
733                 PS_FORCE_PIC(IJK,2) = 2.d0*(PIC_P_S(IJPK,1) - PIC_P_S(IJK,1))/(DY(j)+Dy(j_of(ijpk)))
734              ELSE
735                 IF(PIC_P_S(IJK,1).GT.ZERO) then
736                    PS_FORCE_PIC(IJK,2) = 2.d0*(PIC_P_S(IJPK,1) - PIC_P_S(IJK,1))/(DY(j)+Dy(j_of(ijpk)))
737                 ELSE
738                    PS_FORCE_PIC(IJK,2) = zero
739                 ENDIF
740              ENDIF
741     
742              IF(DO_K) THEN
743                 IF(FLUID_AT(IJKP)) then
744                    PS_FORCE_PIC(IJK,3) = 2.d0*(PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/(Dz(k)+Dz(k_of(ijkp)))
745                 ELSE
746                    IF(PIC_P_S(IJK,1).GT.ZERO) then
747                       PS_FORCE_PIC(IJK,3) = 2.d0*(PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/(Dz(k)+Dz(k_of(ijkp)))
748                    ELSE
749                       PS_FORCE_PIC(IJK,3) = ZERO
750                    ENDIF
751                 ENDIF
752              ENDIF
753           ENDDO
754           !Rahul:
755           !the above will not compute pressure gradients normal to the east, south and bottom faces
756           !which are very important
757           I1 = 1
758           DO K1 = KSTART3, KEND3
759              DO J1 = JSTART3, JEND3
760     !//
761                 IF(I1.NE.ISTART2)   EXIT
762                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
763                 IJK = FUNIJK(I1,J1,K1)
764                 I = I_OF(IJK)
765                 J = J_OF(IJK)
766                 K = K_OF(IJK)
767     
768                 IPJK = IP_OF(IJK)
769                 IF(PIC_P_S(IPJK,1).GT.ZERO) then
770                    PS_FORCE_PIC(IJK,1) = 2.d0*(PIC_P_S(IPJK,1) - PIC_P_S(IJK,1))/(DX(I)+DX(I_of(ipjk)))
771                 else
772                    PS_FORCE_PIC(IJK,1) = zero
773                 endif
774     
775              ENDDO
776           ENDDO
777           J1 = 1
778           DO K1 = KSTART3, KEND3
779              DO I1 = ISTART3, IEND3
780     !//
781                 IF(J1.NE.JSTART2)   EXIT
782                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
783                 IJK = FUNIJK(I1,J1,K1)
784                 I = I_OF(IJK)
785                 J = J_OF(IJK)
786                 K = K_OF(IJK)
787     
788     
789                 IJPK = JP_OF(IJK)
790     
791                 IF(PIC_P_S(IJPK,1).GT.ZERO) then
792     
793                    PS_FORCE_PIC(IJK,2) = 2.d0*(PIC_P_S(IJPK,1) - PIC_P_S(IJK,1))/(DY(j)+Dy(j_of(ijpk)))
794                 ELSE
795                    PS_FORCE_PIC(IJK,2) = ZERO
796                 ENDIF
797     
798              END DO
799           END DO
800     
801           IF(DO_K) then
802              K1 = 1
803              DO J1 = JSTART3, JEND3
804                 DO I1 = ISTART3, IEND3
805                    IF(K1.NE.KSTART2)   EXIT
806                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
807                    IJK = FUNIJK(I1,J1,K1)
808                    I = I_OF(IJK)
809                    J = J_OF(IJK)
810                    K = K_OF(IJK)
811     
812     
813                    IJKP = KP_OF(IJK)
814     
815                    IF(PIC_P_S(IJKP,1).GT.ZERO) then
816     
817                       PS_FORCE_PIC(IJK,3) = 2.d0*(PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/(Dz(k)+Dz(k_of(ijkp)))
818                    ELSE
819                       PS_FORCE_PIC(IJK, 3) = ZERO
820                    ENDIF
821     
822                 END DO
823              END DO
824           ENDIF
825     
826           DO IDIM = 1, merge(2,3,NO_K)
827              CALL SEND_RECV(PS_FORCE_PIC(:,IDIM),1)
828           ENDDO
829     
830           CALL SET_INTERPOLATION_SCHEME(2)
831     
832           KORDER = merge ( 1, 2, no_k) !1+(DIMN-2)
833     
834     ! avg_factor=0.25 (in 3D) or =0.5 (in 2D)
835           !avg_factor = 0.25d0*(dimn-2) + 0.5d0*(3-dimn)
836           AVG_FACTOR = merge(0.5d0, 0.25D0, NO_K)
837     
838           do ijk = ijkstart3,ijkend3
839     
840              if(.not.fluid_at(ijk) .or. pinc(ijk).eq.0) cycle
841              i = i_of(ijk)
842              j = j_of(ijk)
843              k = k_of(ijk)
844              pcell(1) = i-1
845              pcell(2) = j-1
846              pcell(3) = merge(1, k-1, no_k) ! =k-1 (in 3d) or =1 (in 2d)
847              call set_interpolation_stencil(pcell,iw,ie,js,jn,kb,&
848              ktp,interp_scheme,merge(2,3,no_k),ordernew = onew)
849     
850     !Compute velocity at grid nodes and set the geometric stencil
851              do k = 1, merge(1, onew, no_K)!  (3-dimn)*1+(dimn-2)*onew
852                 do j = 1,onew
853                    do i = 1,onew
854                       ii = iw + i-1
855                       jj = js + j-1
856                       kk = kb + k-1
857                       cur_ijk = funijk_map_c(ii,jj,kk)
858                       ipjk    = funijk_map_c(ii+1,jj,kk)
859                       ijpk    = funijk_map_c(ii,jj+1,kk)
860                       ipjpk   = funijk_map_c(ii+1,jj+1,kk)
861     
862                       vol_ijk = zero
863                       vol_ipjk = zero
864                       vol_ijpk = zero
865                       vol_ipjpk = zero
866     
867                       vol_ijkp = zero
868                       vol_ipjkp = zero
869                       vol_ijpkp = zero
870                       vol_ipjpkp = zero
871     
872                       if(fluid_at(cur_ijk)) vol_ijk   = vol(cur_ijk)
873                       if(fluid_at(ipjk))    vol_ipjk  = vol(ipjk)
874                       if(fluid_at(ijpk))    vol_ijpk  = vol(ijpk)
875                       if(fluid_at(ipjpk))   vol_ipjpk = vol(ipjpk)
876     
877     
878                       if(DO_K) then
879                          ijkp    = funijk_map_c(ii,jj,kk+1)
880                          ijpkp   = funijk_map_c(ii,jj+1,kk+1)
881                          ipjkp   = funijk_map_c(ii+1,jj,kk+1)
882                          ipjpkp  = funijk_map_c(ii+1,jj+1,kk+1)
883     
884     
885                          if(fluid_at(ijkp))     vol_ijkp   = vol(ijkp)
886                          if(fluid_at(ipjkp))    vol_ipjkp  = vol(ipjkp)
887                          if(fluid_at(ijpkp))    vol_ijpkp  = vol(ijpkp)
888                          if(fluid_at(ipjpkp))   vol_ipjpkp = vol(ipjpkp)
889     
890                       endif
891                       gstencil(i,j,k,1) = xe(ii)
892                       gstencil(i,j,k,2) = yn(jj)
893                       gstencil(i,j,k,3) = merge(DZ(1), ZT(KK), NO_K)
894     
895                       VOL_TOT_SCAL = ZERO
896     
897                       VOL_TOT_SCAL = vol_ijk + vol_ipjk + vol_ijpk + vol_ipjpk + &
898                       & vol_ijkp + vol_ipjkp + vol_ijpkp + vol_ipjpkp
899     
900                       VOL_TOT_VEC = ZERO
901     
902                       VOL_TOT_VEC(1) = VOL(CUR_IJK) + VOL(IJPK)
903                       VOL_TOT_VEC(2) = VOL(CUR_IJK) + VOL(IPJK)
904     
905                       DO M = 1, DES_MMAX
906                          VEL_SOL_STENCIL(i,j,k,1,M) = pic_u_s(cur_ijk,M)*vol(cur_ijk) + pic_u_s(ijpk,M)*vol(ijpk)
907     
908                          VEL_SOL_STENCIL(i,j,k,2,M) = pic_v_s(cur_ijk,M)*vol(cur_ijk) + pic_v_s(ipjk,M)*vol(ipjk)
909                       ENDDO
910                       !ep_g(cur_ijk)*vol_ijk+ ep_g(ipjk)*vol_ipjk+ &
911                       !   &  ep_g(ijpk)*vol_ijpk + ep_g(ipjpk)*vol_ipjpk,
912                       sstencil(i,j,k) = ep_g(cur_ijk)*vol_ijk+ ep_g(ipjk)*vol_ipjk+ &
913                       &  ep_g(ijpk)*vol_ijpk + ep_g(ipjpk)*vol_ipjpk
914     
915                       psgradstencil(i,j,k,1) = PS_FORCE_PIC(cur_ijk,1)*VOL(CUR_IJK) + PS_FORCE_PIC(ijpk,1)*VOL(IJPK)
916     
917                       psgradstencil(i,j,k,2) = PS_FORCE_PIC(cur_ijk,2)*VOL(CUR_IJK) + PS_FORCE_PIC(ipjk,2)*VOL(IPJK)
918     
919                       vstencil(i,j,k,1) = u_g(cur_ijk)*vol(cur_ijk) + u_g(ijpk)*vol(ijpk)
920                       vstencil(i,j,k,2) = v_g(cur_ijk)*vol(cur_ijk) + v_g(ipjk)*vol(ipjk)
921                       if(DO_K) then
922                          VOL_TOT_VEC(1) = VOL_TOT_VEC(1) + VOL(IJKP) + VOL(IJPKP)
923                          VOL_TOT_VEC(2) = VOL_TOT_VEC(2) + VOL(IJKP) + VOL(IPJKP)
924                          VOL_TOT_VEC(3) = VOL(CUR_IJK) + VOL(IPJK) + VOL(IJPK) + VOL(IPJPK)
925                          sstencil(i,j,k) = sstencil(i,j,k) + ep_g(ijkp)*vol_ijkp + ep_g(ipjkp)*vol_ipjkp &
926                          & + ep_g(ijpkp)*vol_ijpkp + ep_g(ipjpkp)*vol_ipjpkp
927     
928                          psgradstencil(i,j,k,1) = psgradstencil(i,j,k,1) &
929                               + PS_FORCE_PIC(ijkp,1)*VOL(ijkp) + PS_FORCE_PIC(ijpkp,1)*vol(ijpkp)
930                          psgradstencil(i,j,k,2) = psgradstencil(i,j,k,2) &
931                               + PS_FORCE_PIC(ijkp,2)*vol(ijkp) + PS_FORCE_PIC(ipjkp,2)*vol(ipjkp)
932                          psgradstencil(i,j,k,3) = PS_FORCE_PIC(cur_ijk,3)*vol(cur_ijk)+&
933                          & PS_FORCE_PIC(ijpk,3)*vol(ijpk)+PS_FORCE_PIC(ipjk,3)*vol(ipjk)+&
934                          & PS_FORCE_PIC(ipjpk,3)*vol(ipjpk)
935     
936                          vstencil(i,j,k,1) = vstencil(i,j,k,1) + u_g(ijkp)*vol(ijkp) + u_g(ijpkp)*vol(ijpkp)
937     
938                          vstencil(i,j,k,2) = vstencil(i,j,k,2) + v_g(ijkp)*vol(ijkp) + v_g(ipjkp)*vol(ipjkp)
939                          vstencil(i,j,k,3) = w_g(cur_ijk)*vol(cur_ijk)+&
940                          & w_g(ijpk)*vol(ijpk)+w_g(ipjk)*vol(ipjk)+w_g(ipjpk)*vol(ipjpk)
941     
942                          DO M = 1, DES_MMAX
943                             VEL_SOL_STENCIL(i,j,k,1, M) = VEL_SOL_STENCIL(i,j,k,1,M) &
944                             & + PIC_u_s(ijkp,M)*vol(ijkp) + PIC_u_s(ijpkp,M)*vol(ijpkp)
945                             VEL_SOL_STENCIL(i,j,k,2, M) = VEL_SOL_STENCIL(i,j,k,2,M) &
946                             & + PIC_v_s(ijkp,M)*vol(ijkp) + PIC_v_s(ipjkp,M)*vol(ipjkp)
947                             VEL_SOL_STENCIL(i,j,k,3, M) = PIC_w_s(cur_ijk,M)*vol(cur_ijk) +&
948                             PIC_w_s(ijpk,M)*vol(ijpk)+PIC_w_s(ipjk,M)*vol(ipjk)+PIC_w_s(ipjpk,M)*vol(ipjpk)
949                          ENDDO
950                       else
951                          psgradstencil(i,j,k,3) = 0.d0
952                          VEL_SOL_STENCIL(i,j,k,3, 1:DES_MMAX) = 0.d0
953                          vstencil(i,j,k,3) = 0.d0
954     
955                       endif
956                       DO IDIM = 1, merge(2,3,NO_K)
957                          IF(VOL_TOT_VEC(IDIM).GT.ZERO)  THEN
958                             psgradstencil(i,j,k,idim) = psgradstencil(i,j,k,idim)/VOL_TOT_VEC(idim)
959     
960                             VEL_SOL_STENCIL(i,j,k,idim, 1:DES_MMAX) = VEL_SOL_STENCIL(i,j,k,idim, 1:DES_MMAX)/VOL_TOT_VEC(idim)
961     
962                             vstencil(i,j,k,idim) = vstencil(i,j,k,idim)/VOL_TOT_VEC(idim)
963     
964                             !no need for if as sum of positive numbers can only be zero
965                             !if and only if each one of them are zero
966                          ENDIF
967     
968                       ENDDO
969     
970     
971                       !write(*,*) 'epg*vol     = ', ep_g(cur_ijk)*vol_ijk, ep_g(ipjk)*vol_ipjk, &
972                       !&  ep_g(ijpk)*vol_ijpk , ep_g(ipjpk)*vol_ipjpk,  ep_g(cur_ijk)*vol_ijk+ ep_g(ipjk)*vol_ipjk+ &
973                       !&  ep_g(ijpk)*vol_ijpk + ep_g(ipjpk)*vol_ipjpk,sstencil(i,j,k)
974     
975     
976                       if(VOL_TOT_SCAL.gt.zero) sstencil(i,j,k) = sstencil(i,j,k)/VOL_TOT_SCAL
977     
978     
979                    enddo
980                 enddo
981              enddo
982     
983                       !loop through particles in the cell
984              do nindx = 1,pinc(ijk)
985                 np = pic(ijk)%p(nindx)
986                 m = pijk(np,5)
987     
988                 if(NO_K) then !2-D
989                    call interpolator(gstencil(1:onew,1:onew,1,1:dimn), &
990                    psgradstencil(1:onew,1:onew,1,1:dimn), &
991                    des_pos_new(1:dimn,np),PS_GRAD(np,1:dimn),  &
992                    onew,interp_scheme,weightp)
993                 else !3-D, diff in psgradstencil size
994                    call interpolator(gstencil(1:onew,1:onew,1:onew,1:dimn), &
995                    psgradstencil(1:onew,1:onew,1:onew,1:dimn), &
996                    des_pos_new(1:dimn,np),PS_GRAD(np,1:dimn),  &
997                    onew,interp_scheme,weightp)
998                 endif
999     
1000                 do idim = 1,  merge(2,3,NO_K)
1001                    AVGSOLVEL_P(IDIM,NP) = ARRAY_DOT_PRODUCT(               &
1002                       VEL_SOL_STENCIL(:,:,:,IDIM,M),WEIGHTP(:,:,:))
1003                 ENDDO
1004     
1005                 EPG_P(NP) = ARRAY_DOT_PRODUCT(SSTENCIL(:,:,:),WEIGHTP(:,:,:))
1006     
1007              END DO
1008           END DO
1009     
1010     
1011           END SUBROUTINE MPPIC_COMPUTE_PS_GRAD
1012     
1013     
1014     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1015     !                                                                      !
1016     !  Subroutine: MPPIC_BC_U_S                                            !
1017     !  Purpose:                                                            !
1018     !                                                                      !
1019     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1020           SUBROUTINE MPPIC_BC_U_S
1021     
1022     !-----------------------------------------------
1023     ! Modules
1024     !-----------------------------------------------
1025           USE parallel
1026           USE constant
1027           USE geometry
1028           USE indices
1029           USE compar
1030           USE mfix_pic
1031           USE functions
1032           IMPLICIT NONE
1033     !-----------------------------------------------
1034     ! Local variables
1035     !-----------------------------------------------
1036     ! Error index
1037           INTEGER :: IER
1038     ! Boundary condition
1039           INTEGER ::  L
1040     ! Indices
1041           INTEGER :: I, J, K, IM, I1, I2, J1, J2, K1, K2, IJK,&
1042                      JM, KM, IJKW, IMJK, IPJK, IP, IJK_WALL
1043     !-----------------------------------------------
1044     
1045     ! Set the default boundary conditions
1046           IF (DO_K) THEN
1047              K1 = 1
1048              DO J1 = jmin3,jmax3
1049                 DO I1 = imin3, imax3
1050                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1051                    IJK_WALL = FUNIJK(I1,J1,K1)
1052                    IJK = FUNIJK(I1,J1,K1+1)
1053                    PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
1054                 END DO
1055              END DO
1056     
1057              K1 = KMAX2
1058              DO J1 = jmin3,jmax3
1059                 DO I1 = imin3, imax3
1060                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1061                    IJK_WALL = FUNIJK(I1,J1,K1)
1062                    IJK = FUNIJK(I1,J1,K1-1)
1063                    PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
1064                 END DO
1065              END DO
1066           ENDIF
1067     
1068           J1 = 1
1069           DO K1 = kmin3, kmax3
1070              DO I1 = imin3, imax3
1071                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1072                 IJK_WALL = FUNIJK(I1,J1,K1)
1073                 IJK = FUNIJK(I1,J1+1,K1)
1074                 PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
1075              END DO
1076           END DO
1077           J1 = JMAX2
1078           DO K1 = kmin3, kmax3
1079              DO I1 = imin3, imax3
1080                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1081                 IJK_WALL = FUNIJK(I1,J1,K1)
1082                 IJK = FUNIJK(I1,J1-1,K1)
1083                 PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
1084              END DO
1085           END DO
1086     
1087           RETURN
1088           END SUBROUTINE MPPIC_BC_U_S
1089     
1090     
1091     
1092     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1093     !                                                                      !
1094     !  Subroutine: MPPIC_BC_V_S                                            !
1095     !  Purpose:                                                            !
1096     !                                                                      !
1097     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1098           SUBROUTINE MPPIC_BC_V_S
1099     
1100     
1101     !-----------------------------------------------
1102     ! Modules
1103     !-----------------------------------------------
1104           USE parallel
1105           USE constant
1106           USE geometry
1107           USE indices
1108           USE compar
1109           USE mfix_pic
1110           USE functions
1111           IMPLICIT NONE
1112     !-----------------------------------------------
1113     ! Local variables
1114     !-----------------------------------------------
1115     ! Error index
1116           INTEGER          IER
1117     ! Boundary condition
1118           INTEGER          L
1119     ! Indices
1120           INTEGER          I,  J, K, JM, I1, I2, J1, J2, K1, K2, IJK,&
1121                            IM, KM, IJKS, IJMK, IJPK, IJK_WALL
1122     !-----------------------------------------------
1123     
1124     ! Set the default boundary conditions
1125           IF (DO_K) THEN
1126              K1 = 1
1127              DO J1 = jmin3,jmax3
1128                 DO I1 = imin3, imax3
1129                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1130                    IJK_WALL = FUNIJK(I1,J1,K1)
1131                    IJK = FUNIJK(I1,J1,K1+1)
1132                    PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
1133                 END DO
1134              END DO
1135              K1 = KMAX2
1136              DO J1 = jmin3,jmax3
1137                 DO I1 = imin3, imax3
1138                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1139                    IJK_WALL = FUNIJK(I1,J1,K1)
1140                    IJK = FUNIJK(I1,J1,K1-1)
1141                    PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
1142                 END DO
1143              END DO
1144           ENDIF
1145     
1146           I1 = 1
1147           DO K1 = kmin3, kmax3
1148              DO J1 = jmin3, jmax3
1149                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1150                 IJK_WALL = FUNIJK(I1,J1,K1)
1151                 IJK = FUNIJK(I1+1,J1,K1)
1152                 PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
1153     
1154              END DO
1155           END DO
1156           I1 = IMAX2
1157           DO K1 = kmin3, kmax3
1158              DO J1 = jmin3, jmax3
1159                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1160                 IJK_WALL = FUNIJK(I1,J1,K1)
1161                 IJK = FUNIJK(I1-1,J1,K1)
1162                 PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
1163     
1164              END DO
1165           END DO
1166           END SUBROUTINE MPPIC_BC_V_S
1167     
1168     
1169     
1170     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1171     !                                                                      !
1172     !  Subroutine: MPPIC_BC_W_S                                            !
1173     !  Purpose:                                                            !
1174     !                                                                      !
1175     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1176           SUBROUTINE MPPIC_BC_W_S
1177     
1178     !-----------------------------------------------
1179     ! Modules
1180     !-----------------------------------------------
1181           USE parallel
1182           USE constant
1183           USE geometry
1184           USE indices
1185           USE compar
1186           USE mfix_pic
1187           USE functions
1188           IMPLICIT NONE
1189     !-----------------------------------------------
1190     ! Local variables
1191     !-----------------------------------------------
1192     ! Error index
1193           INTEGER :: IER
1194     ! Boundary condition
1195           INTEGER :: L
1196     ! Indices
1197           INTEGER :: I, J, K, KM, I1, I2, J1, J2, K1, K2, IJK,&
1198                      IM, JM, IJKB, IJKM, IJKP, IJK_WALL
1199     !-----------------------------------------------
1200     
1201     ! Set the default boundary conditions
1202           J1 = 1
1203           DO K1 = kmin3,kmax3
1204              DO I1 = imin3,imax3
1205                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1206                 IJK = FUNIJK(I1,J1+1,K1)
1207                 IJK_WALL = FUNIJK(I1,J1,K1)
1208                 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
1209              END DO
1210           END DO
1211           J1 = JMAX2
1212           DO K1 = kmin3, kmax3
1213              DO I1 = imin3, imax3
1214                IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1215                 IJK = FUNIJK(I1,J1-1,K1)
1216                 IJK_WALL = FUNIJK(I1,J1,K1)
1217                 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
1218              END DO
1219           END DO
1220           I1 = 1
1221           DO K1 = kmin3, kmax3
1222              DO J1 = jmin3, jmax3
1223                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1224                 IJK = FUNIJK(I1+1,J1,K1)
1225                 IJK_WALL = FUNIJK(I1,J1,K1)
1226                 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
1227              END DO
1228           END DO
1229           I1 = IMAX2
1230           DO K1 = kmin3,kmax3
1231              DO J1 = jmin3,jmax3
1232                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1233                 IJK = FUNIJK(I1-1,J1,K1)
1234                 IJK_WALL = FUNIJK(I1,J1,K1)
1235                 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
1236              END DO
1237           END DO
1238     
1239     
1240           END SUBROUTINE MPPIC_BC_W_S
1241     
1242     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1243     !  Subroutine: WRITE_NODEDATA                                          !
1244     !  Author: R. Garg                                                     !
1245     !                                                                      !
1246     !  Purpose:                                                            !
1247     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1248           SUBROUTINE WRITE_NODEDATA(bufin, funit)
1249           USE param
1250           USE param1
1251           USE parallel
1252           USE constant
1253           USE run
1254           USE geometry
1255           USE indices
1256           USE compar
1257           USE discretelement
1258           use desmpi
1259           USE functions
1260           USE fun_avg
1261     
1262           IMPLICIT NONE
1263           integer, intent(in) :: funit
1264           double precision, dimension(:), intent(in)  :: bufin
1265     
1266           integer :: ijk, i, j,k
1267     
1268           write(funit,*)'VARIABLES= ',' "I" ',' "J" ',' "K" ',' "DES_ROPS_NODE" '
1269     
1270           write(funit, *)'ZONE F=POINT, I=', (IEND1-ISTART2)+1,  ', J=', JEND1-JSTART2+1, ', K=', KEND1-KSTART2 + 1
1271     
1272           DO K = KSTART2, KEND1
1273              DO J = JSTART2, JEND1
1274                 DO I = ISTART2, IEND1
1275                    IJK = funijk(I, J, K)
1276                    !WRITE(*,*) 'IJK = ', IJK, I, J, K , SIZE(BUFIN,1)
1277     
1278                    WRITE(funit, '(3(2x, i10), 3x, g17.8)') I, J, K , DES_ROPS_NODE(IJK,1)
1279                 ENDDO
1280              ENDDO
1281           ENDDO
1282           CLOSE(funit, status = 'keep')
1283           end SUBROUTINE WRITE_NODEDATA
1284     
1285     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1286     !  Subroutine: WRITE_MPPIC_VEL_S                                       !
1287     !  Author: R. Garg                                                     !
1288     !                                                                      !
1289     !  Purpose:                                                            !
1290     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1291           SUBROUTINE WRITE_MPPIC_VEL_S
1292           USE param
1293           USE param1
1294           USE parallel
1295           USE constant
1296           USE run
1297           USE geometry
1298           USE indices
1299           USE compar
1300           USE fldvar, only : ep_g
1301           USE discretelement
1302           USE mfix_pic
1303           USE functions
1304           implicit none
1305           integer :: i, j, k, ijk, fluid_ind, LL, PC, IDIM
1306           double precision :: zcor
1307           character(LEN=100) :: filename
1308           logical finish
1309     
1310           WRITE(filename,'(A,"_",I5.5,".dat")') TRIM(RUN_NAME)//'_U_S_',myPE
1311           OPEN(1000, file = TRIM(filename), form ='formatted', status='unknown')
1312           IF(DIMN.eq.2) then
1313              write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
1314                   ' "EP_s " ', ' "U_S" ', ' "V_S" ',' "DES_U_s" ', ' "DES_V_s" '!, ' "P_S_FOR1" ', ' "P_S_FOR2" '
1315              write(1000,*)'ZONE F=POINT, I=', (IEND3-ISTART3)+1,  ', J=', JEND3-JSTART3+1, ', K=', KEND3-KSTART3 + 1
1316           else
1317              write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
1318                   ' "EP_s " ', ' "U_S" ', ' "V_S" ', ' "W_S" ',' "DES_U_s" ', ' "DES_V_s" ', ' "DES_W_s" '!, &
1319             ! & ' "P_S_FOR1" ', ' "P_S_FOR2" ', ' "P_S_FOR3" '
1320              write(1000,*)'ZONE F=POINT, I=', (IEND3-ISTART3)+1,  ', J=', JEND3-JSTART3+1, ', K=', KEND3-KSTART3 + 1
1321           ENDIF
1322           DO K=KSTART3, KEND3
1323              DO J=JSTART3, JEND3
1324                 DO I=ISTART3, IEND3
1325                    IJK  = FUNIJK(I,J,K)
1326                    IF(FLUID_AT(IJK)) THEN
1327                       FLUID_IND = 1
1328                    ELSE
1329                       FLUID_IND = 0
1330                    END IF
1331                    IF(DIMN.EQ.2) THEN
1332                       ZCOR = ZT(K)
1333                       WRITE(1000,'(3(2X,G17.8),4( 2X, G17.8))')  XE(I-1)+DX(I), YN(J-1)+DY(J),ZCOR, &
1334                          1.D0 - EP_G(IJK), PIC_U_S(IJK,1), PIC_V_S(IJK,1), DES_U_S(IJK,1), &
1335                          DES_V_S(IJK,1) !, PS_FORCE_PIC(IJK,1), PS_FORCE_PIC(IJK,2)
1336                    ELSE
1337                       ZCOR = ZT(K-1) + DZ(K)
1338                       WRITE(1000,'(3(2X,G17.8),4( 2X, G17.8))')  XE(I-1)+DX(I), YN(J-1)+DY(J),ZCOR, &
1339                          1.D0 - EP_G(IJK), PIC_U_S(IJK,1), PIC_V_S(IJK,1), PIC_W_S(IJK,1), DES_U_S(IJK,1), &
1340                         DES_V_S(IJK,1), DES_W_S(IJK,1)!, PS_FORCE_PIC(IJK,1), PS_FORCE_PIC(IJK,2),  PS_FORCE_PIC(IJK,3)
1341                    ENDIF
1342                 ENDDO
1343              ENDDO
1344           ENDDO
1345           CLOSE(1000, STATUS='KEEP')
1346     
1347           return
1348     
1349           WRITE(FILENAME,'(A,"_",I5.5,".DAT")') TRIM(RUN_NAME)//'_PS_FORCE_',myPE
1350           OPEN(1000, file = TRIM(filename), form ='formatted', status='unknown')
1351     
1352           IF(DIMN.eq.3) write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
1353                ' "DELPX" ', '"DELPY"', '"DELPZ" ',' "US_part" ', '"VS_part"' , '"WS_part"', '"EP_s_part"'
1354           IF(DIMN.eq.2) write(1000,*)'VARIABLES= ',' "X" ',' "Y" ', &
1355                ' "DELPX" ', '"DELPY"', ' "US_part" ', '"VS_part"' , '"EP_S_part"'
1356     
1357           PC = 1
1358           DO LL = 1, MAX_PIP
1359     
1360              IF(PC .GT. PIP) EXIT
1361              IF(.NOT.PEA(LL,1)) CYCLE
1362              pc = pc+1
1363              IF(PEA(LL,4)) CYCLE
1364     
1365              WRITE(1000,'(10( 2x, g17.8))') (DES_POS_NEW(IDIM, LL), IDIM = 1, DIMN), &
1366                   (PS_GRAD(LL, IDIM) , IDIM = 1, DIMN), (AVGSOLVEL_P (IDIM, LL) , IDIM = 1, DIMN), 1-EPg_P(LL)
1367           ENDDO
1368           close(1000, status='keep')
1369     
1370           !write(*,*) 'want to quit ?', LL, mAX_PIP, PIP
1371           !read(*,*) finish
1372           !if(finish) STOp
1373           END SUBROUTINE WRITE_MPPIC_VEL_S
1374     
1375     
1376     
1377     
1378     
1379