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