File: N:\mfix\model\des\nsquare.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: NSQUARE                                                C
4     !>  Purpose: DES - N-Square neighbor search
5     !                                                                      C
6     !                                                                      C
7     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
8     !  Reviewer: Sreekanth Pannala                        Date: 09-Nov-06  C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11     
12           SUBROUTINE NSQUARE
13     
14           use des_allocate
15           use des_bc
16           use des_thermo
17           use discretelement
18           use functions
19           use geometry, only: DO_K, xlength, ylength, zlength
20           use param1, only: zero
21     
22           IMPLICIT NONE
23     
24     !-----------------------------------------------
25     ! Local variables
26     !-----------------------------------------------
27           INTEGER I, J, K
28           INTEGER L, LL, CC
29           DOUBLE PRECISION DISTVEC(3), DIST, R_LM
30     ! Temporary variables to adjust particle position in event of periodic
31     ! boundaries
32           DOUBLE PRECISION XPOS(2), YPOS(2), ZPOS(2), TMPPOS(3)
33     ! Max loop limit in each coordinate direction
34           INTEGER II, JJ, KK
35     ! Index to track accounted for particles
36           INTEGER PC
37     ! Index to track unaccounted for particles
38           INTEGER PNPC
39     !-----------------------------------------------
40     
41           PC=1
42           DO L=1, MAX_PIP
43     
44              NEIGHBOR_INDEX(L) = 1
45              if (L .gt. 1) NEIGHBOR_INDEX(L) = NEIGHBOR_INDEX(L-1)
46     
47              IF(PC .GE. PIP ) EXIT
48              IF(IS_NONEXISTENT(L)) CYCLE
49     
50              PNPC = PIP - PC
51              DO LL = L+1, MAX_PIP
52     
53                 IF(PNPC .LE. 0) EXIT
54                 IF(IS_NONEXISTENT(LL)) CYCLE
55     
56                 R_LM = DES_RADIUS(L) + DES_RADIUS(LL)
57                 R_LM = FACTOR_RLM*R_LM
58     
59     ! the following section adjusts the neighbor check routine when
60     ! any boundary is periodic
61     ! ------------------------------
62                 IF (DES_PERIODIC_WALLS) THEN
63                    XPOS(:) = DES_POS_NEW(LL,1)
64                    YPOS(:) = DES_POS_NEW(LL,2)
65                    II = 1
66                    JJ = 1
67                    KK = 1
68     
69                    IF(DES_PERIODIC_WALLS_X) THEN
70                       IF (DES_POS_NEW(L,1) + R_LM > XLENGTH) THEN
71                          II = 2
72                          XPOS(II) = DES_POS_NEW(LL,1) + XLENGTH
73                       ELSEIF (DES_POS_NEW(L,1) - R_LM < ZERO) THEN
74                          II = 2
75                          XPOS(II) = DES_POS_NEW(LL,1) - XLENGTH
76                       ENDIF
77                    ENDIF
78                    IF(DES_PERIODIC_WALLS_Y) THEN
79                       IF (DES_POS_NEW(L,2) + R_LM > YLENGTH) THEN
80                          JJ = 2
81                          YPOS(JJ) = DES_POS_NEW(LL,2) + YLENGTH
82                       ELSEIF (DES_POS_NEW(L,2) - R_LM < YLENGTH) THEN
83                          JJ = 2
84                          YPOS(JJ) = DES_POS_NEW(LL,2) - YLENGTH
85                       ENDIF
86                    ENDIF
87                    IF(DO_K) THEN
88                       ZPOS(:) = DES_POS_NEW(LL,3)
89                       IF(DES_PERIODIC_WALLS_Z) THEN
90                          IF (DES_POS_NEW(L,3) + R_LM > ZLENGTH) THEN
91                             KK = 2
92                             ZPOS(KK) = DES_POS_NEW(LL,3) + ZLENGTH
93                          ELSEIF (DES_POS_NEW(L,3) - R_LM < ZERO) THEN
94                             KK = 2
95                             ZPOS(KK) = DES_POS_NEW(LL,3) - ZLENGTH
96                          ENDIF
97                       ENDIF
98                    ENDIF
99     
100     ! if particle L is within R_LM of a periodic boundary then check
101     ! particles LL current position and its position shifted to the
102     ! opposite boundary for neighbor contact
103                    OUTER: DO I = 1,II
104                       DO J = 1,JJ
105                          TMPPOS(1) = XPOS(I)
106                          TMPPOS(2) = YPOS(J)
107                          IF (DO_K) THEN
108                             DO K = 1,KK
109                                TMPPOS(3) = ZPOS(K)
110                                DISTVEC(:) = TMPPOS(:) - DES_POS_NEW(L,:)
111                                DIST = dot_product(DISTVEC,DISTVEC)
112                                IF (DIST.LE.R_LM) EXIT OUTER
113                             ENDDO
114                          ELSE
115                             DISTVEC(:) = TMPPOS(:) - DES_POS_NEW(L,:)
116                             DIST = dot_product(DISTVEC,DISTVEC)
117                             IF (DIST.LE.R_LM) EXIT OUTER
118                          ENDIF
119                       ENDDO
120                    ENDDO OUTER
121     
122                 ELSE   ! if .not.des_periodic_walls
123                    DISTVEC(:) = DES_POS_NEW(LL,:) - DES_POS_NEW(L,:)
124                    DIST = dot_product(DISTVEC,DISTVEC)
125                 ENDIF    ! endif des_periodic_walls
126     ! ------------------------------
127     
128                 IF (DIST < R_LM**2) THEN
129                    cc = add_pair(L, LL)
130                 ENDIF
131                 PNPC = PNPC - 1
132              ENDDO   ! end loop over LL
133     
134              PC = PC + 1
135           ENDDO   ! end loop over L
136     
137           RETURN
138           END SUBROUTINE NSQUARE
139