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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Subroutine: MPPIC_COMP_EULERIAN_VELS_NON_CG                         !
3     !  Author: R. Garg                                                     !
4     !                                                                      !
5     !  Purpose:                                                            !
6     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
7           SUBROUTINE MPPIC_COMP_EULERIAN_VELS_NON_CG
8             USE compar
9             USE constant
10             use desmpi
11             USE discretelement
12             use fldvar, only: u_s, v_s, w_s
13             USE functions
14             USE geometry
15             USE indices
16             USE mfix_pic
17             USE parallel
18             USE param
19             USE param1
20             USE physprop, only: mmax
21             IMPLICIT NONE
22     !-----------------------------------------------
23     ! Local variables
24     !-----------------------------------------------
25     ! general i, j, k indices
26           INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
27     
28     ! index of solid phase that particle NP belongs to
29           INTEGER :: M
30           INTEGER :: MMAX_TOT
31     !-----------------------------------------------
32     
33           MMAX_TOT = MMAX+DES_MMAX
34           DO IJK = ijkstart3, ijkend3
35                 I = I_OF(IJK)
36                 J = J_OF(IJK)
37                 K = K_OF(IJK)
38                 !U_so(IJK, :) = U_s(IJK, :)
39                 !V_so(IJK, :) = V_s(IJK, :)
40                 !W_so(IJK, :) = W_s(IJK, :)
41     ! could these be replaced by u_s, w_s, v_s?
42                 PIC_U_S(IJK, :) = ZERO
43                 PIC_V_S(IJK, :) = ZERO
44                 PIC_W_S(IJK, :) = ZERO
45                 IF (.NOT.IS_ON_myPE_plus1layer(I,J,K)) CYCLE
46                 IF(.NOT.FLUID_AT(IJK)) CYCLE
47     
48                 IF(I.GE.IMIN1.AND.I.LT.IMAX1) then !because we don't want to
49                    !calculate solids velocity at the wall cells.
50                    IPJK = IP_OF(IJK)
51                    DO M = MMAX+1, MMAX_TOT
52                       PIC_U_S(IJK,M) = 0.5d0*(U_S(IJK, M) + U_S(IPJK,M))
53                    ENDDO
54                 ENDIF
55     
56                 if(J.GE.JMIN1.AND.J.LT.JMAX1) then !because we don't want to
57                    !calculate solids velocity at the wall cells.
58                    IJPK = JP_OF(IJK)
59                    DO M = MMAX+1, MMAX_TOT
60                       PIC_V_S(IJK,M) = 0.5d0*(V_S(IJK, M) + V_S(IJPK,M))
61                    ENDDO
62                 ENDIF
63     
64     
65                 if(K.GE.KMIN1.AND.K.LT.KMAX1.AND.DO_K) then !because we don't want to
66                    !calculate solids velocity at the wall cells.
67                    IJKP = KP_OF(IJK)
68                    DO M = MMAX+1, MMAX_TOT
69                       PIC_W_S(IJK,M) = 0.5d0*(W_S(IJK, M) + W_S(IJKP,M))
70                    ENDDO
71                 ENDIF
72              ENDDO
73              !CALL SET_WALL_BC(IER)
74              !the above routine will apply noslip or free slip BC as per the mfix  convention.
75              !currently, this implies NSW or FSW wall BC's will be re-applied to gas-phase
76              !field as well. This can be changed later on to be more specific to MPPIC case
77              !CALL WRITE_MPPIC_VEL_S
78              CALL MPPIC_BC_U_S
79              CALL MPPIC_BC_V_S
80              IF(DO_K) CALL MPPIC_BC_W_S
81     
82           END SUBROUTINE MPPIC_COMP_EULERIAN_VELS_NON_CG
83     
84     
85     
86     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
87     !  Subroutine: MPPIC_COMP+EULERIAN_VELS_CG                             !
88     !  Author: R. Garg                                                     !
89     !                                                                      !
90     !  Puryyse:                                                            !
91     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
92           SUBROUTINE MPPIC_COMP_EULERIAN_VELS_CG
93             USE compar
94             USE constant
95             USE cutcell
96             USE desmpi
97             USE discretelement
98             use fldvar, only: u_s, v_s, w_s
99             USE functions
100             USE geometry
101             USE indices
102             USE mfix_pic
103             USE parallel
104             USE param
105             USE param1
106             use physprop, only: mmax
107             IMPLICIT NONE
108     !-----------------------------------------------
109     ! Local variables
110     !-----------------------------------------------
111     ! general i, j, k indices
112           INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
113     
114     ! index of solid phase that particle NP belongs to
115           INTEGER :: M
116           INTEGER :: MMAX_TOT
117     !-----------------------------------------------
118     
119           MMAX_TOT = MMAX+DES_MMAX
120           DO IJK = ijkstart3, ijkend3
121              I = I_OF(IJK)
122              J = J_OF(IJK)
123              K = K_OF(IJK)
124              !U_so(IJK, :) = U_s(IJK, :)
125              !V_so(IJK, :) = V_s(IJK, :)
126              !W_so(IJK, :) = W_s(IJK, :)
127              PIC_U_S(IJK, :) = ZERO
128              PIC_V_S(IJK, :) = ZERO
129              PIC_W_S(IJK, :) = ZERO
130     
131              IF (.NOT.IS_ON_myPE_plus1layer(I,J,K)) CYCLE
132     
133              IF(WALL_U_AT(IJK)) THEN
134                 PIC_U_S(IJK, :) = ZERO
135                 !currently only No slip BC is being set on this mean
136                 !solid's velocity field. Later this wall part can be
137                 !treated separately and U_S set only for scalar cells
138                 !where FLUID_AT(IJK) is true.
139              ELSE
140                 if(.not.FLUID_AT(IJK)) cycle
141     
142                 IPJK = IP_OF(IJK)
143                 IF(FLUID_AT(IPJK)) THEN
144                    DO M = MMAX+1, MMAX_TOT
145                       PIC_U_S(IJK,M) = 0.5d0*(U_S(IJK, M) + U_S(IPJK,M))
146                    ENDDO
147                 ELSE
148                    PIC_U_S(IJK,:) = U_S(IJK, :)
149                 ENDIF
150              ENDIF
151     
152              IF(WALL_V_AT(IJK)) THEN
153                 PIC_V_S(IJK, :) = ZERO
154              ELSE
155                 if(.not.FLUID_AT(IJK)) cycle
156                 IJPK = JP_OF(IJK)
157                 IF(FLUID_AT(IJPK)) THEN
158                    DO M = MMAX+1, MMAX_TOT
159                       PIC_V_S(IJK,M) = 0.5d0*(V_S(IJK, M) + V_S(IJPK,M))
160                    ENDDO
161                 ELSE
162                    PIC_V_S(IJK,:) = V_S(IJK, :)
163                 ENDIF
164              ENDIF
165     
166              IF(DO_K) THEN
167                 IF(WALL_W_AT(IJK)) THEN
168                    PIC_W_S(IJK, :) = ZERO
169                 ELSE
170                    if(.not.FLUID_AT(IJK)) cycle
171                    IJKP = KP_OF(IJK)
172                    IF(FLUID_AT(IJKP)) THEN
173                       DO M = MMAX+1, MMAX_TOT
174                          PIC_W_S(IJK,M) = 0.5d0*(W_S(IJK, M) + W_S(IJKP,M))
175                       ENDDO
176                    ELSE
177                       PIC_W_S(IJK,:) = W_S(IJK, :)
178                    ENDIF
179                 ENDIF
180              ENDIF
181           ENDDO
182     
183           END SUBROUTINE MPPIC_COMP_EULERIAN_VELS_CG
184     
185     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
186     !  Subroutine: MPPIC_APPLY_PS_GRAD_PART                                !
187     !  Author: R. Garg                                                     !
188     !                                                                      !
189     !  Purpose:                                                            !
190     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
191           SUBROUTINE MPPIC_APPLY_PS_GRAD_PART(L)
192     
193           USE param
194           USE param1
195           USE parallel
196           USE constant
197           USE discretelement
198           USE mpi_utility
199           USE mfix_pic
200           USE cutcell
201           USE fldvar, only: ep_g, u_s, v_s, w_s
202           USE fun_avg
203           USE functions
204           IMPLICIT NONE
205           INTEGER, INTENT(IN) :: L
206     !-----------------------------------------------
207     ! Local Variables
208     !-----------------------------------------------
209           INTEGER M, IDIM
210           INTEGER IJK, IJK_C
211     
212           DOUBLE PRECISION COEFF_EN, COEFF_EN2
213     
214           DOUBLE PRECISION DELUP(DIMN), PS_FORCE(DIMN), VEL_ORIG(DIMN)
215     
216     ! dt's in each direction  based on cfl_pic for the mppic case
217     
218           DOUBLE PRECISION :: MEANUS(DIMN, DIMENSION_M)
219           DOUBLE PRECISION :: RELVEL(DIMN)
220           DOUBLE PRECISION :: MEANVEL(DIMN)
221           DOUBLE PRECISION :: VEL_NEW(DIMN)
222     !      INTEGER :: TOT_CASE, case1_count, case2_count, case3_count, case4_count
223     
224           LOGICAL :: INSIDE_DOMAIN
225     !-----------------------------------------------
226     
227           M = PIJK(L,5)
228           IJK = PIJK(L,4)
229           COEFF_EN  = MPPIC_COEFF_EN1
230           COEFF_EN2 = MPPIC_COEFF_EN2
231     
232           VEL_ORIG(1:DIMN) = DES_VEL_NEW(L,:)
233           VEL_NEW (1:DIMN) = DES_VEL_NEW(L,:)
234           !IF(L.eq.1) WRITE(*,*) 'MPPIC COEFFS = ', COEFF_EN, COEFF_EN2
235           IF(L.EQ.FOCUS_PARTICLE) THEN
236     
237              WRITE(*,'(A20,2x,3(2x,i4))') 'CELL ID = ', PIJK(L,1:3)
238              WRITE(*,'(A20,2x,3(2x,g17.8))') 'EPS = ', 1.d0 - EP_g(PIJK(L,4))
239              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL ORIG = ', DES_VEL_NEW(L,:)
240     
241              WRITE(*,'(A20,2x,3(2x,g17.8))') 'FC = ', FC(L,:)
242           ENDIF
243     
244           MEANVEL(1) = U_S(IJK,M)
245           MEANVEL(2) = V_S(IJK,M)
246           IF(DO_K) MEANVEL(3) = W_S(IJK,M)
247     
248           PS_FORCE(:) = PS_GRAD(:,L)
249           !IF(ABS(PS_FORCE(2)).GT.ZERO)  WRITE(*,*) 'PS_FORCE = ', PS_FORCE
250           DELUP(:) = -PS_FORCE(:)
251     
252           MEANUS(:,M) =  AVGSOLVEL_P (:,L)
253           !MEANUS(:,M) = MEANVEL(:)
254           RELVEL(:) = DES_VEL_NEW(L,:) - MEANUS(:,M)
255     
256           !IF(EPg_P(L).gt.1.2d0*ep_star) RETURN
257     
258           DO IDIM = 1, DIMN
259     
260              IF(ABS(PS_FORCE(IDIM)).eq.zero) cycle
261     
262              IF(VEL_ORIG(IDIM)*MEANUS(IDIM,M).GT.ZERO) THEN
263     
264                 IF(VEL_ORIG(IDIM)*DELUP(IDIM).GT.ZERO) THEN
265     
266                    IF(ABS(MEANUS(IDIM,M)) .GT. ABS(VEL_ORIG(IDIM))) THEN
267                       IJK_C = IJK
268                       IF(IDIM.eq.1) then
269                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
270                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
271                       ELSEIF(IDIM.eq.2) then
272                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
273                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
274                       ELSEIF(IDIM.eq.3) then
275                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
276                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
277                       ENDIF
278                       INSIDE_DOMAIN = .false.
279                       INSIDE_DOMAIN = fluid_at(IJK_C)!and.(.not.cut_cell_at(IJK_C))
280     
281                       if(INSIDE_DOMAIN) then
282                          VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
283                       endif
284     !                  case4_count = case4_count + 1
285                    ELSE
286                       !do nothing
287                    ENDIF
288                 ELSE
289                    IF(ABS(VEL_ORIG(IDIM)).GT.ABS(MEANUS(IDIM,M))) then
290                       !VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
291                       IJK_C = IJK
292                       IF(IDIM.eq.1) then
293                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
294                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
295                       ELSEIF(IDIM.eq.2) then
296                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
297                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
298                       ELSEIF(IDIM.eq.3) then
299                          if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
300                          if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
301                       ENDIF
302                       !VEL_NEW(IDIM) = (MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM))
303     
304                       INSIDE_DOMAIN = .false.
305                       INSIDE_DOMAIN = fluid_at(IJK_C)!.and.(.not.cut_cell_at(IJK_C))
306     
307                       if(INSIDE_DOMAIN) then
308                          VEL_NEW(IDIM) = (MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM))
309                       else
310                          VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
311                       endif
312     
313                       IJK_C = IJK
314                       IF(IDIM.eq.1) then
315                          if(VEL_NEW(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
316                          if(VEL_NEW(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
317                       ELSEIF(IDIM.eq.2) then
318                          if(VEL_NEW(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
319                          if(VEL_NEW(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
320                       ELSEIF(IDIM.eq.3) then
321                          if(VEL_NEW(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
322                          if(VEL_NEW(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
323                       ENDIF
324                       INSIDE_DOMAIN = .false.
325                       INSIDE_DOMAIN = fluid_at(IJK_C) !.and.(.not.cut_cell_at(IJK_C))
326     
327                       !if(.not.INSIDE_DOMAIN) then
328                       !   VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
329                       !ENDIF
330     
331     !                  case1_count = case1_count + 1
332                    ELSE
333                       !do nothing
334                       VEL_NEW(IDIM) = COEFF_EN2 * VEL_ORIG(IDIM)
335                       !turning on the above would make the model uncondtionally stable
336     !                  case1_count = case1_count + 1
337     
338                    ENDIF
339                 ENDIF
340              ELSE
341                 IF(MEANUS(IDIM,M)*DELUP(IDIM).GT.ZERO) THEN
342                    VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
343     !               case2_count = case2_count + 1
344                 ELSE
345     !               case3_count = case3_count + 1
346                    !DO NOTHING
347                 ENDIF
348              ENDIF
349     
350              IF(MPPIC_GRAV_TREATMENT) THEN
351                 IF(DELUP(IDIM)*GRAV(IDIM).LT.ZERO.AND.VEL_ORIG(IDIM)*GRAV(IDIM).GT.ZERO) THEN
352                    VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
353                 ENDIF
354              ENDIF
355              DES_VEL_NEW( L,IDIM) = VEL_NEW(IDIM)
356           ENDDO
357     
358              !
359           if(L.eq.FOCUS_PARTICLE) THEN
360              !iF((IJK.eq.epg_min_loc(1).or.IJK_OLD.eq.epg_min_loc(1)).and.epg_min2.lt.0.38) then
361              !if(j.ne.2) cycle
362     
363              WRITE(*,'(A20,2x,3(2x,i5))') 'PIJK I, J, K =', I_OF(IJK),J_OF(IJK),K_OF(IJK)
364              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL ORIG = ', VEL_ORIG(:)
365              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL NEW = ', DES_VEL_NEW(L,:)
366              WRITE(*,'(A20,2x,3(2x,g17.8))') 'MEANUS = ', MEANUS(:,1)
367              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_POS_NEW = ', DES_POS_NEW(L,:)
368              WRITE(*,'(A20,2x,3(2x,g17.8))') 'GRAD PS = ', PS_FORCE(:)
369              WRITE(*,'(A20,2x,3(2x,g17.8))') 'DELUP =  ', DELUP(:)
370     
371              !WRITE(*,'(A20,2x,3(2x,g17.8))') 'UPRIMETAU_INT = ', UPRIMETAU_INT(:)
372              !WRITE(*,'(A20,2x,3(2x,g17.8))') 'UPRIMETAU = ', UPRIMETAU(:)
373              !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) WRITE(*,'(A20,2x,3(2x,g17.8))') 'U*DT, MFP =', UPRIMEMOD*DTSOLID, MEAN_FREE_PATH
374              read(*,*)
375           ENDIF
376     
377           !TOT_CASE = case1_count + case2_count + case3_count + case4_count
378           !IF(TOT_CASE.GT.0) THEN
379           !WRITE(*,'(A, 4(2x,i10))') 'CASE COUNT NUMBERS  = ', case1_count ,case2_count ,case3_count ,case4_count
380           !WRITE(*,'(A, 4(2x,g12.7))') 'CASE COUNT %AGE = ', real(case1_count)*100./real(tot_case), &
381           !     real(case2_count)*100./real(tot_case), real(case3_count)*100./real(tot_case), real(case4_count)*100./real(tot_case)
382           !ENDIF
383           RETURN
384     
385           !MEANUS(:,M) = MEANVEL(:)
386           !RELVEL(:) = DES_VEL_NEW(L,:) - MEANUS(:,M)
387           !DO IDIM = 1, DIMN
388           !    IF(ABS(PS_FORCE(IDIM)).eq.zero) cycle
389     
390           !   IF(RELVEL(IDIM)*DELUP(IDIM).GT.ZERO) THEN
391           !do nothing
392           !    ELSE
393           !       DES_VEL_NEW(L,IDIM) = MEANUS(IDIM,M) - 0.4d0*RELVEL(IDIM)
394     
395           !IF(DES_VEL_NEW(L,IDIM)*DELUP(IDIM).LT.ZERO) DES_VEL_NEW(L,IDIM) = -0.5d0*DES_VEL_NEW(L,IDIM)
396           !    ENDIF
397           ! ENDDO
398           ! CYCLE
399     
400     
401           END SUBROUTINE MPPIC_APPLY_PS_GRAD_PART
402     
403     
404     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
405     !                                                                      !
406     !  Subroutine: MPPIC_BC_U_S                                            !
407     !  Purpose:                                                            !
408     !                                                                      !
409     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
410           SUBROUTINE MPPIC_BC_U_S
411     
412     !-----------------------------------------------
413     ! Modules
414     !-----------------------------------------------
415           USE parallel
416           USE constant
417           USE geometry
418           USE indices
419           USE compar
420           USE mfix_pic
421           USE functions
422           IMPLICIT NONE
423     !-----------------------------------------------
424     ! Local variables
425     !-----------------------------------------------
426     ! Indices
427           INTEGER :: I1, J1, K1, IJK,&
428                      IJK_WALL
429     !-----------------------------------------------
430     
431     ! Set the default boundary conditions
432           IF (DO_K) THEN
433              K1 = 1
434              DO J1 = jmin3,jmax3
435                 DO I1 = imin3, imax3
436                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
437                    IJK_WALL = FUNIJK(I1,J1,K1)
438                    IJK = FUNIJK(I1,J1,K1+1)
439                    PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
440                 END DO
441              END DO
442     
443              K1 = KMAX2
444              DO J1 = jmin3,jmax3
445                 DO I1 = imin3, imax3
446                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
447                    IJK_WALL = FUNIJK(I1,J1,K1)
448                    IJK = FUNIJK(I1,J1,K1-1)
449                    PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
450                 END DO
451              END DO
452           ENDIF
453     
454           J1 = 1
455           DO K1 = kmin3, kmax3
456              DO I1 = imin3, imax3
457                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
458                 IJK_WALL = FUNIJK(I1,J1,K1)
459                 IJK = FUNIJK(I1,J1+1,K1)
460                 PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
461              END DO
462           END DO
463           J1 = JMAX2
464           DO K1 = kmin3, kmax3
465              DO I1 = imin3, imax3
466                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
467                 IJK_WALL = FUNIJK(I1,J1,K1)
468                 IJK = FUNIJK(I1,J1-1,K1)
469                 PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
470              END DO
471           END DO
472     
473           RETURN
474           END SUBROUTINE MPPIC_BC_U_S
475     
476     
477     
478     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
479     !                                                                      !
480     !  Subroutine: MPPIC_BC_V_S                                            !
481     !  Purpose:                                                            !
482     !                                                                      !
483     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
484           SUBROUTINE MPPIC_BC_V_S
485     
486     
487     !-----------------------------------------------
488     ! Modules
489     !-----------------------------------------------
490           USE parallel
491           USE constant
492           USE geometry
493           USE indices
494           USE compar
495           USE mfix_pic
496           USE functions
497           IMPLICIT NONE
498     !-----------------------------------------------
499     ! Local variables
500     !-----------------------------------------------
501     ! Indices
502           INTEGER          I1, J1, K1, IJK,&
503                            IJK_WALL
504     !-----------------------------------------------
505     
506     ! Set the default boundary conditions
507           IF (DO_K) THEN
508              K1 = 1
509              DO J1 = jmin3,jmax3
510                 DO I1 = imin3, imax3
511                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
512                    IJK_WALL = FUNIJK(I1,J1,K1)
513                    IJK = FUNIJK(I1,J1,K1+1)
514                    PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
515                 END DO
516              END DO
517              K1 = KMAX2
518              DO J1 = jmin3,jmax3
519                 DO I1 = imin3, imax3
520                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
521                    IJK_WALL = FUNIJK(I1,J1,K1)
522                    IJK = FUNIJK(I1,J1,K1-1)
523                    PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
524                 END DO
525              END DO
526           ENDIF
527     
528           I1 = 1
529           DO K1 = kmin3, kmax3
530              DO J1 = jmin3, jmax3
531                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
532                 IJK_WALL = FUNIJK(I1,J1,K1)
533                 IJK = FUNIJK(I1+1,J1,K1)
534                 PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
535     
536              END DO
537           END DO
538           I1 = IMAX2
539           DO K1 = kmin3, kmax3
540              DO J1 = jmin3, jmax3
541                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
542                 IJK_WALL = FUNIJK(I1,J1,K1)
543                 IJK = FUNIJK(I1-1,J1,K1)
544                 PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
545     
546              END DO
547           END DO
548           END SUBROUTINE MPPIC_BC_V_S
549     
550     
551     
552     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
553     !                                                                      !
554     !  Subroutine: MPPIC_BC_W_S                                            !
555     !  Purpose:                                                            !
556     !                                                                      !
557     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
558           SUBROUTINE MPPIC_BC_W_S
559     
560     !-----------------------------------------------
561     ! Modules
562     !-----------------------------------------------
563           USE parallel
564           USE constant
565           USE geometry
566           USE indices
567           USE compar
568           USE mfix_pic
569           USE functions
570           IMPLICIT NONE
571     !-----------------------------------------------
572     ! Local variables
573     !-----------------------------------------------
574     ! Indices
575           INTEGER :: I1, J1, K1, IJK,&
576                      IJK_WALL
577     !-----------------------------------------------
578     
579     ! Set the default boundary conditions
580           J1 = 1
581           DO K1 = kmin3,kmax3
582              DO I1 = imin3,imax3
583                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
584                 IJK = FUNIJK(I1,J1+1,K1)
585                 IJK_WALL = FUNIJK(I1,J1,K1)
586                 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
587              END DO
588           END DO
589           J1 = JMAX2
590           DO K1 = kmin3, kmax3
591              DO I1 = imin3, imax3
592                IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
593                 IJK = FUNIJK(I1,J1-1,K1)
594                 IJK_WALL = FUNIJK(I1,J1,K1)
595                 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
596              END DO
597           END DO
598           I1 = 1
599           DO K1 = kmin3, kmax3
600              DO J1 = jmin3, jmax3
601                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
602                 IJK = FUNIJK(I1+1,J1,K1)
603                 IJK_WALL = FUNIJK(I1,J1,K1)
604                 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
605              END DO
606           END DO
607           I1 = IMAX2
608           DO K1 = kmin3,kmax3
609              DO J1 = jmin3,jmax3
610                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
611                 IJK = FUNIJK(I1-1,J1,K1)
612                 IJK_WALL = FUNIJK(I1,J1,K1)
613                 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
614              END DO
615           END DO
616     
617     
618           END SUBROUTINE MPPIC_BC_W_S
619     
620     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
621     !  Subroutine: WRITE_NODEDATA                                          !
622     !  Author: R. Garg                                                     !
623     !                                                                      !
624     !  Purpose:                                                            !
625     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
626           SUBROUTINE WRITE_NODEDATA(funit)
627           USE param
628           USE param1
629           USE parallel
630           USE constant
631           USE run
632           USE geometry
633           USE indices
634           USE compar
635           USE discretelement
636           use desmpi
637           USE functions
638           USE fun_avg
639     
640           IMPLICIT NONE
641           integer, intent(in) :: funit
642     
643           integer :: ijk, i, j,k
644     
645           write(funit,*)'VARIABLES= ',' "I" ',' "J" ',' "K" ',&
646              ' "DES_ROPS_NODE" '
647     
648           write(funit, *)'ZONE F=POINT, I=', (IEND1-ISTART2)+1,  ', J=', &
649              JEND1-JSTART2+1, ', K=', KEND1-KSTART2 + 1
650     
651           DO K = KSTART2, KEND1
652              DO J = JSTART2, JEND1
653                 DO I = ISTART2, IEND1
654                    IJK = funijk(I, J, K)
655                    !WRITE(*,*) 'IJK = ', IJK, I, J, K , SIZE(BUFIN,1)
656     
657                    WRITE(funit, '(3(2x, i10), 3x, g17.8)') I, J, K,&
658                       DES_ROPS_NODE(IJK,1)
659                 ENDDO
660              ENDDO
661           ENDDO
662           CLOSE(funit, status = 'keep')
663           end SUBROUTINE WRITE_NODEDATA
664     
665     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
666     !  Subroutine: WRITE_MPPIC_VEL_S                                       !
667     !  Author: R. Garg                                                     !
668     !                                                                      !
669     !  Purpose:                                                            !
670     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
671           SUBROUTINE WRITE_MPPIC_VEL_S
672           USE param
673           USE param1
674           USE parallel
675           USE constant
676           USE run
677           USE geometry
678           USE indices
679           USE compar
680           USE fldvar, only : ep_g, u_s, v_s, w_s
681           use physprop, only: mmax
682           USE discretelement
683           USE mfix_pic
684           USE functions
685           implicit none
686           integer :: i, j, k, ijk, fluid_ind, LL, PC, IDIM
687           double precision :: zcor
688           character(LEN=255) :: filename
689     
690           WRITE(filename,'(A,"_",I5.5,".dat")') TRIM(RUN_NAME)//'_U_S_',myPE
691           OPEN(1000, file = TRIM(filename), form ='formatted', &
692                status='unknown',CONVERT='BIG_ENDIAN')
693           IF(DIMN.eq.2) then
694              write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
695                   ' "EP_s " ', ' "pU_S" ', ' "pV_S" ',' "dU_s" ',&
696                   ' "dV_s" '!, ' "P_S_FOR1" ', ' "P_S_FOR2" '
697              write(1000,*)'ZONE F=POINT, I=', (IEND3-ISTART3)+1,  ', J=',&
698                 JEND3-JSTART3+1, ', K=', KEND3-KSTART3 + 1
699           else
700              write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
701                   ' "EP_s " ', ' "pU_S" ', ' "pV_S" ', ' "pW_S" ',&
702                   ' "dU_s" ', ' "dV_s" ', ' "dW_s" '!, &
703             ! & ' "P_S_FOR1" ', ' "P_S_FOR2" ', ' "P_S_FOR3" '
704              write(1000,*)'ZONE F=POINT, I=', (IEND3-ISTART3)+1,  ', J=',&
705                 JEND3-JSTART3+1, ', K=', KEND3-KSTART3 + 1
706           ENDIF
707           DO K=KSTART3, KEND3
708              DO J=JSTART3, JEND3
709                 DO I=ISTART3, IEND3
710                    IJK  = FUNIJK(I,J,K)
711                    IF(FLUID_AT(IJK)) THEN
712                       FLUID_IND = 1
713                    ELSE
714                       FLUID_IND = 0
715                    END IF
716                    IF(DIMN.EQ.2) THEN
717                       ZCOR = ZT(K)
718                       WRITE(1000,'(3(2X,G17.8),4( 2X, G17.8))') &
719                          XE(I-1)+DX(I), YN(J-1)+DY(J), ZCOR, 1.D0-EP_G(IJK),&
720                          PIC_U_S(IJK,MMAX+1), PIC_V_S(IJK,MMAX+1),&
721                          U_S(IJK,MMAX+1), V_S(IJK,MMAX+1) !, PS_FORCE_PIC(1,IJK), PS_FORCE_PIC(2,IJK)
722                    ELSE
723                       ZCOR = ZT(K-1) + DZ(K)
724                       WRITE(1000,'(3(2X,G17.8),4( 2X, G17.8))') XE(I-1)+DX(I),&
725                          YN(J-1)+DY(J), ZCOR, 1.D0-EP_G(IJK), &
726                          PIC_U_S(IJK,MMAX+1), PIC_V_S(IJK,MMAX+1), PIC_W_S(IJK,MMAX+1),&
727                          U_S(IJK,MMAX+1), V_S(IJK,MMAX+1), W_S(IJK,MMAX+1)!, PS_FORCE_PIC(1,IJK), PS_FORCE_PIC(2,IJK),  PS_FORCE_PIC(3,IJK)
728                    ENDIF
729                 ENDDO
730              ENDDO
731           ENDDO
732           CLOSE(1000, STATUS='KEEP')
733     
734           return
735     
736           WRITE(FILENAME,'(A,"_",I5.5,".DAT")') &
737              TRIM(RUN_NAME)//'_PS_FORCE_',myPE
738           OPEN(1000, file = TRIM(filename), form ='formatted',&
739                status='unknown',CONVERT='BIG_ENDIAN')
740     
741           IF(DIMN.eq.3) write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ',&
742                ' "DELPX" ', '"DELPY"', '"DELPZ" ',&
743                ' "US_part" ', '"VS_part"' , '"WS_part"', '"EP_s_part"'
744           IF(DIMN.eq.2) write(1000,*)'VARIABLES= ',' "X" ',' "Y" ', &
745                ' "DELPX" ', '"DELPY"', ' "US_part" ', '"VS_part"' , &
746                '"EP_S_part"'
747     
748           PC = 1
749           DO LL = 1, MAX_PIP
750     
751              IF(PC .GT. PIP) EXIT
752              IF(IS_NONEXISTENT(LL)) CYCLE
753              pc = pc+1
754              IF(IS_GHOST(LL) .OR. IS_ENTERING_GHOST(LL) .OR. &
755                 IS_EXITING_GHOST(LL)) CYCLE
756     
757              WRITE(1000,'(10( 2x, g17.8))') &
758                 (DES_POS_NEW(LL,IDIM), IDIM=1,DIMN), &
759                 (PS_GRAD(IDIM,LL), IDIM=1,DIMN), &
760                 (AVGSOLVEL_P(IDIM,LL),IDIM=1,DIMN), 1-EPg_P(LL)
761           ENDDO
762           close(1000, status='keep')
763     
764           !write(*,*) 'want to quit ?', LL, mAX_PIP, PIP
765           !read(*,*) finish
766           !if(finish) STOp
767           END SUBROUTINE WRITE_MPPIC_VEL_S
768