pic_routines.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: MPPIC_COMP_EULERIAN_VELS_NON_CG                         !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose:                                                            !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE mppic_comp_eulerian_vels_non_cg
        USE compar
        USE constant
        use desmpi
        USE discretelement
        use fldvar, only: u_s, v_s, w_s
        USE functions
        USE geometry
        USE indices
        USE mfix_pic
        USE parallel
        USE param
        USE param1
        USE physprop, only: mmax
        IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! general i, j, k indices
      INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP

! index of solid phase that particle NP belongs to
      INTEGER :: M
      INTEGER :: MMAX_TOT
!-----------------------------------------------

      mmax_tot = mmax+des_mmax
      DO ijk = ijkstart3, ijkend3
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
            !U_so(IJK, :) = U_s(IJK, :)
            !V_so(IJK, :) = V_s(IJK, :)
            !W_so(IJK, :) = W_s(IJK, :)
! could these be replaced by u_s, w_s, v_s?
            pic_u_s(ijk, :) = zero
            pic_v_s(ijk, :) = zero
            pic_w_s(ijk, :) = zero
            IF (.NOT.is_on_mype_plus1layer(i,j,k)) cycle
            IF(.NOT.fluid_at(ijk)) cycle

            IF(i.GE.imin1.AND.i.LT.imax1) then !because we don't want to
               !calculate solids velocity at the wall cells.
               ipjk = ip_of(ijk)
               DO m = mmax+1, mmax_tot
                  pic_u_s(ijk,m) = 0.5d0*(u_s(ijk, m) + u_s(ipjk,m))
               ENDDO
            ENDIF

            if(j.GE.jmin1.AND.j.LT.jmax1) then !because we don't want to
               !calculate solids velocity at the wall cells.
               ijpk = jp_of(ijk)
               DO m = mmax+1, mmax_tot
                  pic_v_s(ijk,m) = 0.5d0*(v_s(ijk, m) + v_s(ijpk,m))
               ENDDO
            ENDIF


            if(k.GE.kmin1.AND.k.LT.kmax1.AND.do_k) then !because we don't want to
               !calculate solids velocity at the wall cells.
               ijkp = kp_of(ijk)
               DO m = mmax+1, mmax_tot
                  pic_w_s(ijk,m) = 0.5d0*(w_s(ijk, m) + w_s(ijkp,m))
               ENDDO
            ENDIF
         ENDDO
         !CALL SET_WALL_BC(IER)
         !the above routine will apply noslip or free slip BC as per the mfix  convention.
         !currently, this implies NSW or FSW wall BC's will be re-applied to gas-phase
         !field as well. This can be changed later on to be more specific to MPPIC case
         !CALL WRITE_MPPIC_VEL_S
         CALL mppic_bc_u_s
         CALL mppic_bc_v_s
         IF(do_k) CALL mppic_bc_w_s

      END SUBROUTINE mppic_comp_eulerian_vels_non_cg



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: MPPIC_COMP+EULERIAN_VELS_CG                             !
!  Author: R. Garg                                                     !
!                                                                      !
!  Puryyse:                                                            !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE mppic_comp_eulerian_vels_cg
        USE compar
        USE constant
        USE cutcell
        USE desmpi
        USE discretelement
        use fldvar, only: u_s, v_s, w_s
        USE functions
        USE geometry
        USE indices
        USE mfix_pic
        USE parallel
        USE param
        USE param1
        use physprop, only: mmax
        IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! general i, j, k indices
      INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP

! index of solid phase that particle NP belongs to
      INTEGER :: M
      INTEGER :: MMAX_TOT
!-----------------------------------------------

      mmax_tot = mmax+des_mmax
      DO ijk = ijkstart3, ijkend3
         i = i_of(ijk)
         j = j_of(ijk)
         k = k_of(ijk)
         !U_so(IJK, :) = U_s(IJK, :)
         !V_so(IJK, :) = V_s(IJK, :)
         !W_so(IJK, :) = W_s(IJK, :)
         pic_u_s(ijk, :) = zero
         pic_v_s(ijk, :) = zero
         pic_w_s(ijk, :) = zero

         IF (.NOT.is_on_mype_plus1layer(i,j,k)) cycle

         IF(wall_u_at(ijk)) THEN
            pic_u_s(ijk, :) = zero
            !currently only No slip BC is being set on this mean
            !solid's velocity field. Later this wall part can be
            !treated separately and U_S set only for scalar cells
            !where FLUID_AT(IJK) is true.
         ELSE
            if(.not.fluid_at(ijk)) cycle

            ipjk = ip_of(ijk)
            IF(fluid_at(ipjk)) THEN
               DO m = mmax+1, mmax_tot
                  pic_u_s(ijk,m) = 0.5d0*(u_s(ijk, m) + u_s(ipjk,m))
               ENDDO
            ELSE
               pic_u_s(ijk,:) = u_s(ijk, :)
            ENDIF
         ENDIF

         IF(wall_v_at(ijk)) THEN
            pic_v_s(ijk, :) = zero
         ELSE
            if(.not.fluid_at(ijk)) cycle
            ijpk = jp_of(ijk)
            IF(fluid_at(ijpk)) THEN
               DO m = mmax+1, mmax_tot
                  pic_v_s(ijk,m) = 0.5d0*(v_s(ijk, m) + v_s(ijpk,m))
               ENDDO
            ELSE
               pic_v_s(ijk,:) = v_s(ijk, :)
            ENDIF
         ENDIF

         IF(do_k) THEN
            IF(wall_w_at(ijk)) THEN
               pic_w_s(ijk, :) = zero
            ELSE
               if(.not.fluid_at(ijk)) cycle
               ijkp = kp_of(ijk)
               IF(fluid_at(ijkp)) THEN
                  DO m = mmax+1, mmax_tot
                     pic_w_s(ijk,m) = 0.5d0*(w_s(ijk, m) + w_s(ijkp,m))
                  ENDDO
               ELSE
                  pic_w_s(ijk,:) = w_s(ijk, :)
               ENDIF
            ENDIF
         ENDIF
      ENDDO

      END SUBROUTINE mppic_comp_eulerian_vels_cg

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: MPPIC_APPLY_PS_GRAD_PART                                !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose:                                                            !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE mppic_apply_ps_grad_part(L)

      USE param
      USE param1
      USE parallel
      USE constant
      USE discretelement
      USE mpi_utility
      USE mfix_pic
      USE cutcell
      USE fldvar, only: ep_g, u_s, v_s, w_s
      USE fun_avg
      USE functions
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: L
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
      INTEGER M, IDIM
      INTEGER IJK, IJK_C

      DOUBLE PRECISION COEFF_EN, COEFF_EN2

      DOUBLE PRECISION DELUP(dimn), PS_FORCE(dimn), VEL_ORIG(dimn)

! dt's in each direction  based on cfl_pic for the mppic case

      DOUBLE PRECISION :: MEANUS(dimn, dimension_m)
      DOUBLE PRECISION :: RELVEL(dimn)
      DOUBLE PRECISION :: MEANVEL(dimn)
      DOUBLE PRECISION :: VEL_NEW(dimn)
!      INTEGER :: TOT_CASE, case1_count, case2_count, case3_count, case4_count

      LOGICAL :: INSIDE_DOMAIN
!-----------------------------------------------

      m = pijk(l,5)
      ijk = pijk(l,4)
      coeff_en  = mppic_coeff_en1
      coeff_en2 = mppic_coeff_en2

      vel_orig(1:dimn) = des_vel_new(l,:)
      vel_new(1:dimn) = des_vel_new(l,:)
      !IF(L.eq.1) WRITE(*,*) 'MPPIC COEFFS = ', COEFF_EN, COEFF_EN2
      IF(l.EQ.focus_particle) THEN

         WRITE(*,'(A20,2x,3(2x,i4))') 'CELL ID = ', pijk(l,1:3)
         WRITE(*,'(A20,2x,3(2x,g17.8))') 'EPS = ', 1.d0 - ep_g(pijk(l,4))
         WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL ORIG = ', des_vel_new(l,:)

         WRITE(*,'(A20,2x,3(2x,g17.8))') 'FC = ', fc(l,:)
      ENDIF

      meanvel(1) = u_s(ijk,m)
      meanvel(2) = v_s(ijk,m)
      IF(do_k) meanvel(3) = w_s(ijk,m)

      ps_force(:) = ps_grad(:,l)
      !IF(ABS(PS_FORCE(2)).GT.ZERO)  WRITE(*,*) 'PS_FORCE = ', PS_FORCE
      delup(:) = -ps_force(:)

      meanus(:,m) =  avgsolvel_p(:,l)
      !MEANUS(:,M) = MEANVEL(:)
      relvel(:) = des_vel_new(l,:) - meanus(:,m)

      !IF(EPg_P(L).gt.1.2d0*ep_star) RETURN

      DO idim = 1, dimn

         IF(abs(ps_force(idim)).eq.zero) cycle

         IF(vel_orig(idim)*meanus(idim,m).GT.zero) THEN

            IF(vel_orig(idim)*delup(idim).GT.zero) THEN

               IF(abs(meanus(idim,m)) .GT. abs(vel_orig(idim))) THEN
                  ijk_c = ijk
                  IF(idim.eq.1) then
                     if(vel_orig(idim).LT.zero) ijk_c = im_of(ijk)
                     if(vel_orig(idim).GT.zero) ijk_c = ip_of(ijk)
                  ELSEIF(idim.eq.2) then
                     if(vel_orig(idim).LT.zero) ijk_c = jm_of(ijk)
                     if(vel_orig(idim).GT.zero) ijk_c = jp_of(ijk)
                  ELSEIF(idim.eq.3) then
                     if(vel_orig(idim).LT.zero) ijk_c = km_of(ijk)
                     if(vel_orig(idim).GT.zero) ijk_c = kp_of(ijk)
                  ENDIF
                  inside_domain = .false.
                  inside_domain = fluid_at(ijk_c)!and.(.not.cut_cell_at(IJK_C))

                  if(inside_domain) then
                     vel_new(idim) = meanus(idim,m) - coeff_en*relvel(idim)
                  endif
!                  case4_count = case4_count + 1
               ELSE
                  !do nothing
               ENDIF
            ELSE
               IF(abs(vel_orig(idim)).GT.abs(meanus(idim,m))) then
                  !VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
                  ijk_c = ijk
                  IF(idim.eq.1) then
                     if(vel_orig(idim).LT.zero) ijk_c = im_of(ijk)
                     if(vel_orig(idim).GT.zero) ijk_c = ip_of(ijk)
                  ELSEIF(idim.eq.2) then
                     if(vel_orig(idim).LT.zero) ijk_c = jm_of(ijk)
                     if(vel_orig(idim).GT.zero) ijk_c = jp_of(ijk)
                  ELSEIF(idim.eq.3) then
                     if(vel_orig(idim).LT.zero) ijk_c = km_of(ijk)
                     if(vel_orig(idim).GT.zero) ijk_c = kp_of(ijk)
                  ENDIF
                  !VEL_NEW(IDIM) = (MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM))

                  inside_domain = .false.
                  inside_domain = fluid_at(ijk_c)!.and.(.not.cut_cell_at(IJK_C))

                  if(inside_domain) then
                     vel_new(idim) = (meanus(idim,m) - coeff_en*relvel(idim))
                  else
                     vel_new(idim) = -coeff_en*vel_orig(idim)
                  endif

                  ijk_c = ijk
                  IF(idim.eq.1) then
                     if(vel_new(idim).LT.zero) ijk_c = im_of(ijk)
                     if(vel_new(idim).GT.zero) ijk_c = ip_of(ijk)
                  ELSEIF(idim.eq.2) then
                     if(vel_new(idim).LT.zero) ijk_c = jm_of(ijk)
                     if(vel_new(idim).GT.zero) ijk_c = jp_of(ijk)
                  ELSEIF(idim.eq.3) then
                     if(vel_new(idim).LT.zero) ijk_c = km_of(ijk)
                     if(vel_new(idim).GT.zero) ijk_c = kp_of(ijk)
                  ENDIF
                  inside_domain = .false.
                  inside_domain = fluid_at(ijk_c) !.and.(.not.cut_cell_at(IJK_C))

                  !if(.not.INSIDE_DOMAIN) then
                  !   VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
                  !ENDIF

!                  case1_count = case1_count + 1
               ELSE
                  !do nothing
                  vel_new(idim) = coeff_en2 * vel_orig(idim)
                  !turning on the above would make the model uncondtionally stable
!                  case1_count = case1_count + 1

               ENDIF
            ENDIF
         ELSE
            IF(meanus(idim,m)*delup(idim).GT.zero) THEN
               vel_new(idim) = meanus(idim,m) - coeff_en*relvel(idim)
!               case2_count = case2_count + 1
            ELSE
!               case3_count = case3_count + 1
               !DO NOTHING
            ENDIF
         ENDIF

         IF(mppic_grav_treatment) THEN
            IF(delup(idim)*grav(idim).LT.zero.AND.vel_orig(idim)*grav(idim).GT.zero) THEN
               vel_new(idim) = -coeff_en*vel_orig(idim)
            ENDIF
         ENDIF
         des_vel_new( l,idim) = vel_new(idim)
      ENDDO

         !
      if(l.eq.focus_particle) THEN
         !iF((IJK.eq.epg_min_loc(1).or.IJK_OLD.eq.epg_min_loc(1)).and.epg_min2.lt.0.38) then
         !if(j.ne.2) cycle

         WRITE(*,'(A20,2x,3(2x,i5))') 'PIJK I, J, K =', i_of(ijk),j_of(ijk),k_of(ijk)
         WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL ORIG = ', vel_orig(:)
         WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL NEW = ', des_vel_new(l,:)
         WRITE(*,'(A20,2x,3(2x,g17.8))') 'MEANUS = ', meanus(:,1)
         WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_POS_NEW = ', des_pos_new(l,:)
         WRITE(*,'(A20,2x,3(2x,g17.8))') 'GRAD PS = ', ps_force(:)
         WRITE(*,'(A20,2x,3(2x,g17.8))') 'DELUP =  ', delup(:)

         !WRITE(*,'(A20,2x,3(2x,g17.8))') 'UPRIMETAU_INT = ', UPRIMETAU_INT(:)
         !WRITE(*,'(A20,2x,3(2x,g17.8))') 'UPRIMETAU = ', UPRIMETAU(:)
         !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) WRITE(*,'(A20,2x,3(2x,g17.8))') 'U*DT, MFP =', UPRIMEMOD*DTSOLID, MEAN_FREE_PATH
         read(*,*)
      ENDIF

      !TOT_CASE = case1_count + case2_count + case3_count + case4_count
      !IF(TOT_CASE.GT.0) THEN
      !WRITE(*,'(A, 4(2x,i10))') 'CASE COUNT NUMBERS  = ', case1_count ,case2_count ,case3_count ,case4_count
      !WRITE(*,'(A, 4(2x,g12.7))') 'CASE COUNT %AGE = ', real(case1_count)*100./real(tot_case), &
      !     real(case2_count)*100./real(tot_case), real(case3_count)*100./real(tot_case), real(case4_count)*100./real(tot_case)
      !ENDIF
      RETURN

      !MEANUS(:,M) = MEANVEL(:)
      !RELVEL(:) = DES_VEL_NEW(L,:) - MEANUS(:,M)
      !DO IDIM = 1, DIMN
      !    IF(ABS(PS_FORCE(IDIM)).eq.zero) cycle

      !   IF(RELVEL(IDIM)*DELUP(IDIM).GT.ZERO) THEN
      !do nothing
      !    ELSE
      !       DES_VEL_NEW(L,IDIM) = MEANUS(IDIM,M) - 0.4d0*RELVEL(IDIM)

      !IF(DES_VEL_NEW(L,IDIM)*DELUP(IDIM).LT.ZERO) DES_VEL_NEW(L,IDIM) = -0.5d0*DES_VEL_NEW(L,IDIM)
      !    ENDIF
      ! ENDDO
      ! CYCLE


      END SUBROUTINE mppic_apply_ps_grad_part


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: MPPIC_BC_U_S                                            !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE mppic_bc_u_s

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE parallel
      USE constant
      USE geometry
      USE indices
      USE compar
      USE mfix_pic
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I1, J1, K1, IJK,&
                 IJK_WALL
!-----------------------------------------------

! Set the default boundary conditions
      IF (do_k) THEN
         k1 = 1
         DO j1 = jmin3,jmax3
            DO i1 = imin3, imax3
               IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
               ijk_wall = funijk(i1,j1,k1)
               ijk = funijk(i1,j1,k1+1)
               pic_u_s(ijk_wall, :) = pic_u_s(ijk,:)
            END DO
         END DO

         k1 = kmax2
         DO j1 = jmin3,jmax3
            DO i1 = imin3, imax3
               IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
               ijk_wall = funijk(i1,j1,k1)
               ijk = funijk(i1,j1,k1-1)
               pic_u_s(ijk_wall, :) = pic_u_s(ijk,:)
            END DO
         END DO
      ENDIF

      j1 = 1
      DO k1 = kmin3, kmax3
         DO i1 = imin3, imax3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk_wall = funijk(i1,j1,k1)
            ijk = funijk(i1,j1+1,k1)
            pic_u_s(ijk_wall, :) = pic_u_s(ijk,:)
         END DO
      END DO
      j1 = jmax2
      DO k1 = kmin3, kmax3
         DO i1 = imin3, imax3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk_wall = funijk(i1,j1,k1)
            ijk = funijk(i1,j1-1,k1)
            pic_u_s(ijk_wall, :) = pic_u_s(ijk,:)
         END DO
      END DO

      RETURN
      END SUBROUTINE mppic_bc_u_s



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: MPPIC_BC_V_S                                            !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE mppic_bc_v_s


!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE parallel
      USE constant
      USE geometry
      USE indices
      USE compar
      USE mfix_pic
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER          I1, J1, K1, IJK,&
                       IJK_WALL
!-----------------------------------------------

! Set the default boundary conditions
      IF (do_k) THEN
         k1 = 1
         DO j1 = jmin3,jmax3
            DO i1 = imin3, imax3
               IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
               ijk_wall = funijk(i1,j1,k1)
               ijk = funijk(i1,j1,k1+1)
               pic_v_s(ijk_wall, :) = pic_v_s(ijk,:)
            END DO
         END DO
         k1 = kmax2
         DO j1 = jmin3,jmax3
            DO i1 = imin3, imax3
               IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
               ijk_wall = funijk(i1,j1,k1)
               ijk = funijk(i1,j1,k1-1)
               pic_v_s(ijk_wall, :) = pic_v_s(ijk,:)
            END DO
         END DO
      ENDIF

      i1 = 1
      DO k1 = kmin3, kmax3
         DO j1 = jmin3, jmax3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk_wall = funijk(i1,j1,k1)
            ijk = funijk(i1+1,j1,k1)
            pic_v_s(ijk_wall, :) = pic_v_s(ijk,:)

         END DO
      END DO
      i1 = imax2
      DO k1 = kmin3, kmax3
         DO j1 = jmin3, jmax3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk_wall = funijk(i1,j1,k1)
            ijk = funijk(i1-1,j1,k1)
            pic_v_s(ijk_wall, :) = pic_v_s(ijk,:)

         END DO
      END DO
      END SUBROUTINE mppic_bc_v_s



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: MPPIC_BC_W_S                                            !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE mppic_bc_w_s

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE parallel
      USE constant
      USE geometry
      USE indices
      USE compar
      USE mfix_pic
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I1, J1, K1, IJK,&
                 IJK_WALL
!-----------------------------------------------

! Set the default boundary conditions
      j1 = 1
      DO k1 = kmin3,kmax3
         DO i1 = imin3,imax3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk = funijk(i1,j1+1,k1)
            ijk_wall = funijk(i1,j1,k1)
            pic_w_s(ijk_wall,:) = pic_w_s(ijk,:)
         END DO
      END DO
      j1 = jmax2
      DO k1 = kmin3, kmax3
         DO i1 = imin3, imax3
           IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk = funijk(i1,j1-1,k1)
            ijk_wall = funijk(i1,j1,k1)
            pic_w_s(ijk_wall,:) = pic_w_s(ijk,:)
         END DO
      END DO
      i1 = 1
      DO k1 = kmin3, kmax3
         DO j1 = jmin3, jmax3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk = funijk(i1+1,j1,k1)
            ijk_wall = funijk(i1,j1,k1)
            pic_w_s(ijk_wall,:) = pic_w_s(ijk,:)
         END DO
      END DO
      i1 = imax2
      DO k1 = kmin3,kmax3
         DO j1 = jmin3,jmax3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk = funijk(i1-1,j1,k1)
            ijk_wall = funijk(i1,j1,k1)
            pic_w_s(ijk_wall,:) = pic_w_s(ijk,:)
         END DO
      END DO


      END SUBROUTINE mppic_bc_w_s

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: WRITE_NODEDATA                                          !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose:                                                            !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE write_nodedata(funit)
      USE param
      USE param1
      USE parallel
      USE constant
      USE run
      USE geometry
      USE indices
      USE compar
      USE discretelement
      use desmpi
      USE functions
      USE fun_avg

      IMPLICIT NONE
      integer, intent(in) :: funit

      integer :: ijk, i, j,k

      write(funit,*)'VARIABLES= ',' "I" ',' "J" ',' "K" ',&
         ' "DES_ROPS_NODE" '

      write(funit, *)'ZONE F=POINT, I=', (iend1-istart2)+1,  ', J=', &
         jend1-jstart2+1, ', K=', kend1-kstart2 + 1

      DO k = kstart2, kend1
         DO j = jstart2, jend1
            DO i = istart2, iend1
               ijk = funijk(i, j, k)
               !WRITE(*,*) 'IJK = ', IJK, I, J, K , SIZE(BUFIN,1)

               WRITE(funit, '(3(2x, i10), 3x, g17.8)') i, j, k,&
                  des_rops_node(ijk,1)
            ENDDO
         ENDDO
      ENDDO
      CLOSE(funit, status = 'keep')
      end SUBROUTINE write_nodedata

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: WRITE_MPPIC_VEL_S                                       !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose:                                                            !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE write_mppic_vel_s
      USE param
      USE param1
      USE parallel
      USE constant
      USE run
      USE geometry
      USE indices
      USE compar
      USE fldvar, only : ep_g, u_s, v_s, w_s
      use physprop, only: mmax
      USE discretelement
      USE mfix_pic
      USE functions
      implicit none
      integer :: i, j, k, ijk, fluid_ind, LL, PC, IDIM
      double precision :: zcor
      character(LEN=255) :: filename

      WRITE(filename,'(A,"_",I5.5,".dat")') trim(run_name)//'_U_S_',mype
      OPEN(1000, file = trim(filename), form ='formatted', &
           status='unknown',convert='BIG_ENDIAN')
      IF(dimn.eq.2) then
         write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
              ' "EP_s " ', ' "pU_S" ', ' "pV_S" ',' "dU_s" ',&
              ' "dV_s" '!, ' "P_S_FOR1" ', ' "P_S_FOR2" '
         write(1000,*)'ZONE F=POINT, I=', (iend3-istart3)+1,  ', J=',&
            jend3-jstart3+1, ', K=', kend3-kstart3 + 1
      else
         write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
              ' "EP_s " ', ' "pU_S" ', ' "pV_S" ', ' "pW_S" ',&
              ' "dU_s" ', ' "dV_s" ', ' "dW_s" '!, &
        ! & ' "P_S_FOR1" ', ' "P_S_FOR2" ', ' "P_S_FOR3" '
         write(1000,*)'ZONE F=POINT, I=', (iend3-istart3)+1,  ', J=',&
            jend3-jstart3+1, ', K=', kend3-kstart3 + 1
      ENDIF
      DO k=kstart3, kend3
         DO j=jstart3, jend3
            DO i=istart3, iend3
               ijk  = funijk(i,j,k)
               IF(fluid_at(ijk)) THEN
                  fluid_ind = 1
               ELSE
                  fluid_ind = 0
               END IF
               IF(dimn.EQ.2) THEN
                  zcor = zt(k)
                  WRITE(1000,'(3(2X,G17.8),4( 2X, G17.8))') &
                     xe(i-1)+dx(i), yn(j-1)+dy(j), zcor, 1.d0-ep_g(ijk),&
                     pic_u_s(ijk,mmax+1), pic_v_s(ijk,mmax+1),&
                     u_s(ijk,mmax+1), v_s(ijk,mmax+1) !, PS_FORCE_PIC(1,IJK), PS_FORCE_PIC(2,IJK)
               ELSE
                  zcor = zt(k-1) + dz(k)
                  WRITE(1000,'(3(2X,G17.8),4( 2X, G17.8))') xe(i-1)+dx(i),&
                     yn(j-1)+dy(j), zcor, 1.d0-ep_g(ijk), &
                     pic_u_s(ijk,mmax+1), pic_v_s(ijk,mmax+1), pic_w_s(ijk,mmax+1),&
                     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)
               ENDIF
            ENDDO
         ENDDO
      ENDDO
      CLOSE(1000, status='KEEP')

      return

      WRITE(filename,'(A,"_",I5.5,".DAT")') &
         trim(run_name)//'_PS_FORCE_',mype
      OPEN(1000, file = trim(filename), form ='formatted',&
           status='unknown',convert='BIG_ENDIAN')

      IF(dimn.eq.3) write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ',&
           ' "DELPX" ', '"DELPY"', '"DELPZ" ',&
           ' "US_part" ', '"VS_part"' , '"WS_part"', '"EP_s_part"'
      IF(dimn.eq.2) write(1000,*)'VARIABLES= ',' "X" ',' "Y" ', &
           ' "DELPX" ', '"DELPY"', ' "US_part" ', '"VS_part"' , &
           '"EP_S_part"'

      pc = 1
      DO ll = 1, max_pip

         IF(pc .GT. pip) EXIT
         IF(is_nonexistent(ll)) cycle
         pc = pc+1
         IF(is_ghost(ll) .OR. is_entering_ghost(ll) .OR. &
            is_exiting_ghost(ll)) cycle

         WRITE(1000,'(10( 2x, g17.8))') &
            (des_pos_new(ll,idim), idim=1,dimn), &
            (ps_grad(idim,ll), idim=1,dimn), &
            (avgsolvel_p(idim,ll),idim=1,dimn), 1-epg_p(ll)
      ENDDO
      close(1000, status='keep')

      !write(*,*) 'want to quit ?', LL, mAX_PIP, PIP
      !read(*,*) finish
      !if(finish) STOp
      END SUBROUTINE write_mppic_vel_s