40
41
42
45
46
47
48#include "implicit_f.inc"
49#include "comlock.inc"
50
51
52
53#include "com04_c.inc"
54#include "sphcom.inc"
55#include "param_c.inc"
56#include "task_c.inc"
57
58
59
60 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),ITAB(NUMNOD),
61 . ISPCOND(NISPCOND,*),ISPSYM(NSPCOND,*),
62 . IREDUCE, WSP2SORT(*), ITASK, KREDUCE(*),
63 . LGAUGE(3,NBGAUGE)
65 . x(3,numnod) ,v(3,numnod) ,ms(*) ,spbuf(nspbuf,*) ,
66 . xframe(nxframe,*) ,
67 . xspsym(3,*) ,vspsym(3,*),
68 . myspatrue, dmax, myspatrue2, gauge(llgauge,nbgauge)
69
70
71
72 INTEGER K,N,INOD,JNOD,J,NVOIS,M,SM,JS,I,NN,IG,
73 . NVOIS1,NVOIS2,NVOISS,NVOISS1,NVOISS2,
74 . IS,IC,NC,L,ISLIDE, IAUX, NSPHSYM_L,
75 . JK,JL,IERROR,NS,NB,JJ1,JJ2, JJ, JJJ,
76 . JVOIS(NSPHSYM+KVOISPH+1), JSTOR(NSPHSYM+KVOISPH+1),
77 . JPERM(NSPHSYM+KVOISPH+1), JCOND(KVOISPH)
79 . xi,yi,zi,di,rhoi,xj,yj,zj,dj,rhoj,dij,
80 . vxi,vyi,vzi,vxj,vyj,vzj,
81 . ox,oy,oz,nx,ny,nz,
82 . xs,ys,zs,vxs,vys,vzs,vn,dd,dm,dk,dl,
83 . xisort,yisort,zisort,disort,
84 . xjsort,yjsort,zjsort,djsort,
85 . spalinr, dvois(nsphsym+kvoisph+1)
86
87
88 spalinr=sqrt(one + myspatrue)
89
90
91 myspatrue2 = myspatrue
92
93
94 nvois1 = 0
95 nvoiss1 = 0
96
97 DO nc=1,nspcond
98 is=ispcond(3,nc)
99 ic=ispcond(2,nc)
100 islide=ispcond(5,nc)
101 ox=xframe(10,is)
102 oy=xframe(11,is)
103 oz=xframe(12,is)
104 nx=xframe(3*(ic-1)+1,is)
105 ny=xframe(3*(ic-1)+2,is)
106 nz=xframe(3*(ic-1)+3,is)
107 DO ns=1+itask,nsp2sort,nthread
108 n=wsp2sort(ns)
109 inod =kxsp(3,n)
110 xi =x(1,inod)
111 yi =x(2,inod)
112 zi =x(3,inod)
113 vxi=v(1,inod)
114 vyi=v(2,inod)
115 vzi=v(3,inod)
116 dd=(xi-ox)*nx+(yi-oy)*ny+(zi-oz)*nz
117 IF (ispsym(nc,n)/=-1) THEN
118 nsphsym_l = ispsym(nc,n)
119 xs=xi - two*dd*nx
120 ys=yi - two*dd*ny
121 zs=zi - two*dd*nz
122 IF(islide==0)THEN
123 vxs=-vxi
124 vys=-vyi
125 vzs=-vzi
126 ELSE
127 vn=vxi*nx+vyi*ny+vzi*nz
128 vxs=vxi - two*vn*nx
129 vys=vyi - two*vn*ny
130 vzs=vzi - two*vn*nz
131 ENDIF
132 xspsym(1,nsphsym_l)= xs
133 xspsym(2,nsphsym_l)= ys
134 xspsym(3,nsphsym_l)= zs
135 vspsym(1,nsphsym_l)=vxs
136 vspsym(2,nsphsym_l)=vys
137 vspsym(3,nsphsym_l)=vzs
138 ENDIF
139 ENDDO
140
141
142
143 DO ns = itask+1,
nsphr,nthread
144 xi =xsphr(3,ns)
145 yi =xsphr(4,ns)
146 zi =xsphr(5,ns)
147 vxi=xsphr(9,ns)
148 vyi=xsphr(10,ns)
149 vzi=xsphr(11,ns)
150 dd=(xi-ox)*nx+(yi-oy)*ny+(zi-oz)*nz
153 xs=xi - two*dd*nx
154 ys=yi - two*dd*ny
155 zs=zi - two*dd*nz
156 IF(islide==0)THEN
157 vxs=-vxi
158 vys=-vyi
159 vzs=-vzi
160 ELSE
161 vn=vxi*nx+vyi*ny+vzi*nz
162 vxs=vxi - two*vn*nx
163 vys=vyi - two*vn*ny
164 vzs=vzi - two*vn*nz
165 ENDIF
166 xspsym(1,nsphsym_l)= xs
167 xspsym(2,nsphsym_l)= ys
168 xspsym(3,nsphsym_l)= zs
169 vspsym(1,nsphsym_l)=vxs
170 vspsym(2,nsphsym_l)=vys
171 vspsym(3,nsphsym_l)=vzs
172 END IF
173 END DO
174 END DO
175
176
177
179
180
181
182
183
184 IF (nspcond/=0)THEN
185
186 DO ns=itask+1,nsp2sort,nthread
187 n=wsp2sort(ns)
188 inod=kxsp(3,n)
189 xi =x(1,inod)
190 yi =x(2,inod)
191 zi =x(3,inod)
192 di =spbuf(1,n)
193
194 nvois2 =kxsp(5,n)
195 kxsp(7,n)=0
196
197 DO nc=1,nspcond
198 js=ispsym(nc,n)
199 IF(js>0)THEN
200 xj =xspsym(1,js)
201 yj =xspsym(2,js)
202 zj =xspsym(3,js)
203 dij=di+di
204 dij=dij*dij
205 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
206 IF(dd<=(one+myspatrue2)*dij) THEN
207 nvoiss=kxsp(7,n)
208 nvoiss=nvoiss+1
209 jvois(nvois2+nvoiss)=nc+n*(nspcond+1)
210 dvois(nvois2+nvoiss)=dd/dij
211 kxsp(7,n)=nvoiss
212 ENDIF
213 ENDIF
214 ENDDO
215
216 DO i=1,nvois2
217 jnod=ixsp(i,n)
218 IF(jnod>0)THEN
219 m=nod2sp(jnod)
220 dj =spbuf(1,m)
221 dij=(di+dj)
222 dij=dij*dij
223 DO nc=1,nspcond
224 js=ispsym(nc,m)
225 IF(js>0)THEN
226 xj =xspsym(1,js)
227 yj =xspsym(2,js)
228 zj =xspsym(3,js)
229 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
230 IF(dd<=(one+myspatrue2)*dij) THEN
231 nvoiss=kxsp(7,n)
232 nvoiss=nvoiss+1
233 jvois(nvois2+nvoiss)=nc+m*(nspcond+1)
234 dvois(nvois2+nvoiss)=dd/dij
235 kxsp(7,n)=nvoiss
236 END IF
237 END IF
238 END DO
239 ELSE
240 nn = -jnod
241 dj =xsphr(2,nn)
242 dij=(di+dj)
243 dij=dij*dij
244 DO nc=1,nspcond
246 IF(js>0)THEN
247 xj =xspsym(1,js)
248 yj =xspsym(2,js)
249 zj =xspsym(3,js)
250 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
251 IF(dd<=(one+myspatrue2)*dij) THEN
252 nvoiss=kxsp(7,n)
253 nvoiss=nvoiss+1
254 jvois(nvois2+nvoiss)=-nc-nn*(nspcond+1)
255
256
257
258 dvois(nvois2+nvoiss)=dd/dij
259 kxsp(7,n)=nvoiss
260 END IF
261 END IF
262 END DO
263 END IF
264 END DO
265 nvoiss=kxsp(7,n)
266 IF(nvois2+nvoiss<=kvoisph)THEN
267
268
269
270 nvoiss1=0
271 nvoiss2=nvoiss
272 DO k=nvois2+1,nvois2+nvoiss
273 dk=dvois(k)
274 jk=jvois(k)
275 IF(dk<one)THEN
276 nvoiss1=nvoiss1+1
277 ixsp(nvois2+nvoiss1,n)=jk
278 ELSE
279 ixsp(nvois2+nvoiss2,n)=jk
280 nvoiss2=nvoiss2-1
281 ENDIF
282 ENDDO
283 kxsp(6,n)=nvoiss1
284 ELSE
285
286
287 ireduce=1
288 kreduce(n)=1
289
290
291 DO j=1,nvois2
292 jnod=ixsp(j,n)
293 IF(jnod>0)THEN
294 xj =x(1,jnod)
295 yj =x(2,jnod)
296 zj =x(3,jnod)
297 m =nod2sp(jnod)
298 dj=spbuf(1,m)
299 dij=di+dj
300 dij=dij*dij
301 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
302 dvois(j)=dd/dij
303 jvois(j)=m*(nspcond+1)
304 ELSE
305 nn = -jnod
306 xj = xsphr(3,nn)
307 yj = xsphr(4,nn)
308 zj = xsphr(5,nn)
309 dj = xsphr(2,nn)
310 dij=di+dj
311 dij=dij*dij
312 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
313 dvois(j)=dd/dij
314 jvois(j)=-nn*(nspcond+1)
315 END IF
316 ENDDO
317
318 CALL myqsort(nvois2+nvoiss,dvois,jperm,ierror)
319 DO k=1,nvois2+nvoiss
320 jstor(k)=jvois(k)
321 ENDDO
322 DO k=1,kvoisph
323 jvois(k)=jstor(jperm(k))
324 ENDDO
325
326
327
328 dk=dvois(kvoisph)
329
330 nvois=0
331 DO k=1,kvoisph
332 IF(dvois(k)<dk)THEN
333 nvois=nvois+1
334 END IF
335 END DO
336
337
338 nvois1=0
339 DO k=1,nvois
340 jk=jvois(k)
341 dk=dvois(k)
342 IF(jk>0)THEN
343 nc=mod(jk,nspcond+1)
344 IF(nc==0.AND.dk<one)THEN
345 nvois1=nvois1+1
346 ixsp(nvois1,n)=kxsp(3,jk/(nspcond+1))
347 ENDIF
348 ELSE
349 jk = -jk
350 nc=mod(jk,nspcond+1)
351 IF(nc==0.AND.dk<one)THEN
352 nvois1=nvois1+1
353 ixsp(nvois1,n)=-jk/(nspcond+1)
354 END IF
355 END IF
356 END DO
357 kxsp(4,n)=nvois1
358 nvois2=nvois1
359
360 DO k=1,nvois
361 jk=jvois(k)
362 dk=dvois(k)
363 IF(jk>0)THEN
364 nc=mod(jk,nspcond+1)
365 IF(nc==0.AND.dk>=one)THEN
366 nvois2=nvois2+1
367 ixsp(nvois2,n)=kxsp(3,jk/(nspcond+1))
368 END IF
369 ELSE
370 jk = -jk
371 nc=mod(jk,nspcond+1)
372 IF(nc==0.AND.dk>=one)THEN
373 nvois2=nvois2+1
374 ixsp(nvois2,n)=-jk/(nspcond+1)
375 END IF
376 END IF
377 END DO
378 kxsp(5,n)=nvois2
379
380
381 nvoiss1=0
382 DO k=1,nvois
383 jk=jvois(k)
384 dk=dvois(k)
385 IF(jk>0)THEN
386 nc=mod(jk,nspcond+1)
387 IF(nc/=0.AND.dk<one)THEN
388 nvoiss1=nvoiss1+1
389 ixsp(nvois2+nvoiss1,n)=jk
390 END IF
391 ELSE
392 jk = -jk
393 nc=mod(jk,nspcond+1)
394 IF(nc/=0.AND.dk<one)THEN
395 nvoiss1=nvoiss1+1
396 ixsp(nvois2+nvoiss1,n)=-jk
397 END IF
398 END IF
399 END DO
400 kxsp(6,n)=nvoiss1
401 nvoiss2=nvoiss1
402
403 DO k=1,nvois
404 jk=jvois(k)
405 dk=dvois(k)
406 IF(jk>0)THEN
407 nc=mod(jk,nspcond+1)
408 IF(nc/=0.AND.dk>=one)THEN
409 nvoiss2=nvoiss2+1
410 ixsp(nvois2+nvoiss2,n)=jk
411 END IF
412 ELSE
413 jk = -jk
414 nc=mod(jk,nspcond+1)
415 IF(nc/=0.AND.dk>=one)THEN
416 nvoiss2=nvoiss2+1
417 ixsp(nvois2+nvoiss2,n)=-jk
418 END IF
419 END IF
420 END DO
421 kxsp(7,n)=nvoiss2
422
423 END IF
424
425 IF(nvois1+nvoiss1 > lvoisph)ireduce=1
426
427 END DO
428
429
430
431
432
433 DO ig=1,nbgauge
434 IF(lgauge(1,ig) <= -(numels+1))THEN
435 n=numsph+ig
436
437 xi =gauge(2,ig)
438 yi =gauge(3,ig)
439 zi =gauge(4,ig)
440
441 nvois2 =kxsp(5,n)
442 kxsp(7,n)=0
443
444 DO i=1,nvois2
445 jnod=ixsp(i,n)
446 IF(jnod>0)THEN
447 m=nod2sp(jnod)
448 dj =spbuf(1,m)
449 dij=two*dj
450 dij=dij*dij
451 DO nc=1,nspcond
452 js=ispsym(nc,m)
453 IF(js>0)THEN
454 xj =xspsym(1,js)
455 yj =xspsym(2,js)
456 zj =xspsym(3,js)
457 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
458 IF(dd<=(one+myspatrue2)*dij) THEN
459 nvoiss=kxsp(7,n)
460 nvoiss=nvoiss+1
461 jvois(nvois2+nvoiss)=nc+m*(nspcond+1)
462 dvois(nvois2+nvoiss)=dd/dij
463 kxsp(7,n)=nvoiss
464 END IF
465 END IF
466 END DO
467 ELSE
468 nn = -jnod
469 dj =xsphr(2,nn)
470 dij=two*dj
471 dij=dij*dij
472 DO nc=1,nspcond
474 IF(js>0)THEN
475 xj =xspsym(1,js)
476 yj =xspsym(2,js)
477 zj =xspsym(3,js)
478 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
479 IF(dd<=(one+myspatrue2)*dij) THEN
480 nvoiss=kxsp(7,n)
481 nvoiss=nvoiss+1
482 jvois(nvois2+nvoiss)=-nc-nn*(nspcond+1)
483
484
485
486 dvois(nvois2+nvoiss)=dd/dij
487 kxsp(7,n)=nvoiss
488 END IF
489 END IF
490 END DO
491 END IF
492 END DO
493 nvoiss=kxsp(7,n)
494 IF(nvois2+nvoiss<=kvoisph)THEN
495
496
497
498 nvoiss1=0
499 nvoiss2=nvoiss
500 DO k=nvois2+1,nvois2+nvoiss
501 dk=dvois(k)
502 jk=jvois(k)
503 IF(dk<one)THEN
504 nvoiss1=nvoiss1+1
505 ixsp(nvois2+nvoiss1,n)=jk
506 ELSE
507 ixsp(nvois2+nvoiss2,n)=jk
508 nvoiss2=nvoiss2-1
509 ENDIF
510 ENDDO
511 kxsp(6,n)=nvoiss1
512 ELSE
513
514
515 DO j=1,nvois2
516 jnod=ixsp(j,n)
517 IF(jnod>0)THEN
518 xj =x(1,jnod)
519 yj =x(2,jnod)
520 zj =x(3,jnod)
521 m =nod2sp(jnod)
522 dj=spbuf(1,m)
523 dij=two*dj
524 dij=dij*dij
525 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
526 dvois(j)=dd/dij
527 jvois(j)=m*(nspcond+1)
528 ELSE
529 nn = -jnod
530 xj = xsphr(3,nn)
531 yj = xsphr(4,nn)
532 zj = xsphr(5,nn)
533 dj = xsphr(2,nn)
534 dij=two*dj
535 dij=dij*dij
536 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
537 dvois(j)=dd/dij
538 jvois(j)=-nn*(nspcond+1)
539 END IF
540 ENDDO
541
542 CALL myqsort(nvois2+nvoiss,dvois,jperm,ierror)
543 DO k=1,nvois2+nvoiss
544 jstor(k)=jvois(k)
545 ENDDO
546 DO k=1,kvoisph
547 jvois(k)=jstor(jperm(k))
548 ENDDO
549
550
551
552 dk=dvois(kvoisph)
553
554 nvois=0
555 DO k=1,kvoisph
556 IF(dvois(k)<dk)THEN
557 nvois=nvois+1
558 END IF
559 END DO
560
561
562 nvois1=0
563 DO k=1,nvois
564 jk=jvois(k)
565 dk=dvois(k)
566 IF(jk>0)THEN
567 nc=mod(jk,nspcond+1)
568 IF(nc==0.AND.dk<one)THEN
569 nvois1=nvois1+1
570 ixsp(nvois1,n)=kxsp(3,jk/(nspcond+1))
571 ENDIF
572 ELSE
573 jk = -jk
574 nc=mod(jk,nspcond+1)
575 IF(nc==0.AND.dk<one)THEN
576 nvois1=nvois1+1
577 ixsp(nvois1,n)=-jk/(nspcond+1)
578 END IF
579 END IF
580 END DO
581 kxsp(4,n)=nvois1
582 nvois2=nvois1
583
584 DO k=1,nvois
585 jk=jvois(k)
586 dk=dvois(k)
587 IF(jk>0)THEN
588 nc=mod(jk,nspcond+1)
589 IF(nc==0.AND.dk>=one)THEN
590 nvois2=nvois2+1
591 ixsp(nvois2,n)=kxsp(3,jk/(nspcond+1))
592 END IF
593 ELSE
594 jk = -jk
595 nc=mod(jk,nspcond+1)
596 IF(nc==0.AND.dk>=one)THEN
597 nvois2=nvois2+1
598 ixsp(nvois2,n)=-jk/(nspcond+1)
599 END IF
600 END IF
601 END DO
602 kxsp(5,n)=nvois2
603
604
605 nvoiss1=0
606 DO k=1,nvois
607 jk=jvois(k)
608 dk=dvois(k)
609 IF(jk>0)THEN
610 nc=mod(jk,nspcond+1)
611 IF(nc/=0.AND.dk<one)THEN
612 nvoiss1=nvoiss1+1
613 ixsp(nvois2+nvoiss1,n)=jk
614 END IF
615 ELSE
616 jk = -jk
617 nc=mod(jk,nspcond+1)
618 IF(nc/=0.AND.dk<one)THEN
619 nvoiss1=nvoiss1+1
620 ixsp(nvois2+nvoiss1,n)=-jk
621 END IF
622 END IF
623 END DO
624 kxsp(6,n)=nvoiss1
625 nvoiss2=nvoiss1
626
627 DO k=1,nvois
628 jk=jvois(k)
629 dk=dvois(k)
630 IF(jk>0)THEN
631 nc=mod(jk,nspcond+1)
632 IF(nc/=0.AND.dk>=one)THEN
633 nvoiss2=nvoiss2+1
634 ixsp(nvois2+nvoiss2,n)=jk
635 END IF
636 ELSE
637 jk = -jk
638 nc=mod(jk,nspcond+1)
639 IF(nc/=0.AND.dk>=one)THEN
640 nvoiss2=nvoiss2+1
641 ixsp(nvois2+nvoiss2,n)=-jk
642 END IF
643 END IF
644 END DO
645 kxsp(7,n)=nvoiss2
646
647 END IF
648
649 END IF
650 END DO
651
652 END IF
653
654 RETURN
subroutine myqsort(n, a, perm, error)
integer, dimension(:,:), allocatable ispsymr