nsquare.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: NSQUARE                                                C
!                                                                      C
!                                                                      C
!  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
!  Reviewer: Sreekanth Pannala                        Date: 09-Nov-06  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE nsquare

      use des_allocate
      use des_bc
      use des_thermo
      use discretelement
      use functions
      use geometry, only: do_k, xlength, ylength, zlength
      use param1, only: zero

      IMPLICIT NONE

!-----------------------------------------------
! Local variables
!-----------------------------------------------
      INTEGER I, J, K
      INTEGER L, LL, CC
      DOUBLE PRECISION DISTVEC(3), DIST, R_LM
! Temporary variables to adjust particle position in event of periodic
! boundaries
      DOUBLE PRECISION XPOS(2), YPOS(2), ZPOS(2), TMPPOS(3)
! Max loop limit in each coordinate direction
      INTEGER II, JJ, KK
! Index to track accounted for particles
      INTEGER PC
! Index to track unaccounted for particles
      INTEGER PNPC
!-----------------------------------------------

      pc=1
      DO l=1, max_pip

         neighbor_index(l) = 1
         if (l .gt. 1) neighbor_index(l) = neighbor_index(l-1)

         IF(pc .GE. pip ) EXIT
         IF(is_nonexistent(l)) cycle

         pnpc = pip - pc
         DO ll = l+1, max_pip

            IF(pnpc .LE. 0) EXIT
            IF(is_nonexistent(ll)) cycle

            r_lm = des_radius(l) + des_radius(ll)
            r_lm = factor_rlm*r_lm

! the following section adjusts the neighbor check routine when
! any boundary is periodic
! ------------------------------
            IF (des_periodic_walls) THEN
               xpos(:) = des_pos_new(ll,1)
               ypos(:) = des_pos_new(ll,2)
               ii = 1
               jj = 1
               kk = 1

               IF(des_periodic_walls_x) THEN
                  IF (des_pos_new(l,1) + r_lm > xlength) THEN
                     ii = 2
                     xpos(ii) = des_pos_new(ll,1) + xlength
                  ELSEIF (des_pos_new(l,1) - r_lm < zero) THEN
                     ii = 2
                     xpos(ii) = des_pos_new(ll,1) - xlength
                  ENDIF
               ENDIF
               IF(des_periodic_walls_y) THEN
                  IF (des_pos_new(l,2) + r_lm > ylength) THEN
                     jj = 2
                     ypos(jj) = des_pos_new(ll,2) + ylength
                  ELSEIF (des_pos_new(l,2) - r_lm < ylength) THEN
                     jj = 2
                     ypos(jj) = des_pos_new(ll,2) - ylength
                  ENDIF
               ENDIF
               IF(do_k) THEN
                  zpos(:) = des_pos_new(ll,3)
                  IF(des_periodic_walls_z) THEN
                     IF (des_pos_new(l,3) + r_lm > zlength) THEN
                        kk = 2
                        zpos(kk) = des_pos_new(ll,3) + zlength
                     ELSEIF (des_pos_new(l,3) - r_lm < zero) THEN
                        kk = 2
                        zpos(kk) = des_pos_new(ll,3) - zlength
                     ENDIF
                  ENDIF
               ENDIF

! if particle L is within R_LM of a periodic boundary then check
! particles LL current position and its position shifted to the
! opposite boundary for neighbor contact
               outer: DO i = 1,ii
                  DO j = 1,jj
                     tmppos(1) = xpos(i)
                     tmppos(2) = ypos(j)
                     IF (do_k) THEN
                        DO k = 1,kk
                           tmppos(3) = zpos(k)
                           distvec(:) = tmppos(:) - des_pos_new(l,:)
                           dist = dot_product(distvec,distvec)
                           IF (dist.LE.r_lm) EXIT outer
                        ENDDO
                     ELSE
                        distvec(:) = tmppos(:) - des_pos_new(l,:)
                        dist = dot_product(distvec,distvec)
                        IF (dist.LE.r_lm) EXIT outer
                     ENDIF
                  ENDDO
               ENDDO outer

            ELSE   ! if .not.des_periodic_walls
               distvec(:) = des_pos_new(ll,:) - des_pos_new(l,:)
               dist = dot_product(distvec,distvec)
            ENDIF    ! endif des_periodic_walls
! ------------------------------

            IF (dist < r_lm**2) THEN
               cc = add_pair(l, ll)
            ENDIF
            pnpc = pnpc - 1
         ENDDO   ! end loop over LL

         pc = pc + 1
      ENDDO   ! end loop over L

      RETURN
      END SUBROUTINE nsquare