!----------------------------------------------------------------------- ! ! SUBROUTINE ilist(p,nterms,iterms,pos_comm,pmass, drdotdr,dx,dy,dz,pquad) ! ! !----------------------------------------------------------------------- ! ! ! Subroutine to form the interaction list of a particle. ! !----------------------------------------------------------------------- USE fly_h implicit none INCLUDE 'mpif.h' ! Declaration of local variables. ! ------------------------------- INTEGER(KIND=4) :: j,ind_in,ind_in3,indx,lmax1,ii INTEGER(KIND=4) ::p,i,nnodes,nsubdiv,nterms,nl INTEGER(KIND=4), DIMENSION (nsubcell) :: isub INTEGER(KIND=4), DIMENSION (ndim) ::dire INTEGER(KIND=4), DIMENSION (maxilf)::nodelist,twasubp,twisubset INTEGER :: il_sh LOGICAL ::tolcrit REAL(KIND=8) :: dp_x,dp_y,dp_z,cc_test REAL(KIND=8) :: ddx,ddy,ddz,dd2,cl_s REAL(KIND=8) :: ddr,dd_min REAL(KIND=8), DIMENSION (ndim) ::pos_l, pos_cl,d_c,dr2,dp REAL(KIND=8), DIMENSION (ndim,ndim) ::dd_r REAL(KIND=8), DIMENSION (ndim) ::sh INTEGER(KIND=4), DIMENSION (maxnterm) :: iterms REAL(KIND=8), DIMENSION (ndim) :: pos_comm REAL(KIND=8), DIMENSION (maxnterm) :: pmass REAL(KIND=8), DIMENSION (maxnterm) ::drdotdr ,dx ,dy ,dz REAL(KIND=8), DIMENSION (2*ndim-1,maxnterm) :: pquad !======================================================================= !----------------------------------------------------------------------- ! Initialize list of cells to examine. !----------------------------------------------------------------------- nterms=0 nnodes=1 nodelist(1)=root lmax1 = 0 IF(rmt_acc) THEN pos_l(1:ndim)=pos_rmt(1:ndim,p) ELSE pos_l(1:ndim)=pos(1:ndim,p) ENDIF pos_comm=pos_l !----------------------------------------------------------------------- ! Loop until no cells are left to examine. !----------------------------------------------------------------------- DO WHILE (nnodes.GT.0) lmax1 = lmax1 + 1 cl_s=size_level(lmax1) nsubdiv=0 !----------------------------------------------------------------------- ! Apply tolerance criterion to list of cells. !----------------------------------------------------------------------- DO i=1,nnodes nl=nodelist(i) IF(nl .le. nb_loc) THEN pos_cl(1:ndim)=pos(1:ndim,nl) nterms=nterms+1 iterms(nterms)=nl pmass(nterms)=mass_read dx(nterms)=pos_l(1)-pos_cl(1) dy(nterms)=pos_l(2)-pos_cl(2) dz(nterms)=pos_l(3)-pos_cl(3) drdotdr(nterms)=dx(nterms)**2+ & dy(nterms)**2+dz(nterms)**2 ELSE !IF(nl .le. nb_loc) il_sh=nl-nb_loc pos_cl(1:ndim)=pos_cell(1:ndim,il_sh) cc_test = cl_s**2*tol2inv ddx=pos_l(1)-pos_cl(1) ddy=pos_l(2)-pos_cl(2) ddz=pos_l(3)-pos_cl(3) dd2=ddx**2+ddy**2+ddz**2 tolcrit= dd2 .GE. cc_test !----------------------------------------------------------------------- ! Ewald boundary periodic conditions section !----------------------------------------------------------------------- IF(tolcrit) THEN !----------------------------------------------------------------------- ! Ewald: Search for "mirror" cells. ! ! For each simmetry plane checks whether the ! distance body[pos(p)]--cell[pos(pc)] is such ! that the reflexion w.r.t. direction i ! will make a nearest cell to body pos(p). !----------------------------------------------------------------------- dp(1)=0. dp(2)=0. dp(3)=0. ind_in=0 DO ii=1, ndim sh(ii)=pos_l(ii)-pos_cl(ii) ddr=abs(sh(ii)) if(ddr.gt.L2)then ind_in=ind_in+1 dire(ind_in)=ii ! d_c: + or -Lbox. ! ind_in3=INT(sh(ii)/ddr) ind_in3=1 IF(sh(ii).lt.0) ind_in3=-1 d_c(ind_in)=ind_in3*Lbox endif ENDDO ! If no reflection w.r.t. any plane can produce nearest cell, go to end IF(ind_in.eq.0) THEN nterms=nterms+1 pmass(nterms)=mass_cell(il_sh) pquad(1:2*ndim-1,nterms)=quad(1:2*ndim-1,il_sh) iterms(nterms)=nl dx(nterms)=ddx dy(nterms)=ddy dz(nterms)=ddz drdotdr(nterms)=dd2 CYCLE ! goto next nodelist(i) element ENDIF ! i=index of the mirror cell j=direction (x,y or z) do ii=1, ind_in do j=1, ndim dd_r(ii,j)=pos_cl(j) if(j.eq.dire(ii))dd_r(ii,j)=dd_r(ii,j)+d_c(ii) enddo enddo ! Search the nearest cell among the 3 reflexions if(ind_in.eq.1)then indx=ind_in go to 200 endif DO ii=1, ind_in dr2(ii)=(pos_l(1)-dd_r(ii,1))* & (pos_l(1)-dd_r(ii,1)) + & (pos_l(2)-dd_r(ii,2))*(pos_l(2)-dd_r(ii,2))+ & (pos_l(3)-dd_r(ii,3))*(pos_l(3)-dd_r(ii,3)) ENDDO indx=1 dd_min=dr2(1) if(ind_in.eq.2)then if(dr2(2).le.dd_min)indx=2 endif if(ind_in.eq.3)then if(dr2(2).le.dd_min)indx=2 if(dr2(3).le.dd_min.and.dr2(3).le.dr2(2))indx=3 endif ! Assigns to dp(i=1,ndim) the coordinate shift to reach ! from cell "pc" the nearest (to particle "p") cell outside the domain. 200 continue DO ii=1, ndim dp(ii)=dd_r(indx,ii)-pos_cl(ii) ENDDO dp_x=pos_l(1)-(pos_cl(1)+dp(1)) dp_y=pos_l(2)-(pos_cl(2)+dp(2)) dp_z=pos_l(3)-(pos_cl(3)+dp(3)) tolcrit=(dp_x*dp_x+dp_y*dp_y+dp_z*dp_z).GE.cc_test ENDIF !----------------------------------------------------------------------- ! End of Ewald section !----------------------------------------------------------------------- 120 IF(tolcrit) THEN nterms=nterms+1 pmass(nterms)=mass_cell(il_sh) pquad(1:2*ndim-1,nterms)=quad(1:2*ndim-1,il_sh) iterms(nterms)=nl dx(nterms)=ddx dy(nterms)=ddy dz(nterms)=ddz drdotdr(nterms)=dd2 ELSE !IF(tolcrit) nsubdiv=nsubdiv+1 twisubset(nsubdiv)=i !nodelist index: cell to be opened ENDIF !IF(tolcrit) ENDIF !IF(nl .le. nb_loc) END DO ! DO i=1,nnodes !----------------------------------------------------------------------- ! Add cells which satisfy criterion to interaction list. Note that, ! depending on theta, self-interaction term will be included. !----------------------------------------------------------------------- IF(nterms.GT.maxnterm) THEN write(uterm,*)'PE=',me,' nterms=',nterms, & ' p=',p,' nnodes=',nnodes call flush(uterm) CALL error('ilist error 1: overflow') ENDIF !----------------------------------------------------------------------- ! Add subcells of cells which fail tolerance criterion to list of ! cells to examine. !----------------------------------------------------------------------- IF(8*nsubdiv.GT.maxilf.OR.8*nsubdiv.GT.ncells) THEN CALL error('ilist error 2: overflow') ENDIF DO 40 i=1,nsubdiv il_sh=nodelist(twisubset(i))-nb_loc isub(1:nsubcell)=subp(1:nsubcell,il_sh) twasubp(i)=isub(1) twasubp(i+nsubdiv)=isub(2) twasubp(i+2*nsubdiv)=isub(3) twasubp(i+3*nsubdiv)=isub(4) twasubp(i+4*nsubdiv)=isub(5) twasubp(i+5*nsubdiv)=isub(6) twasubp(i+6*nsubdiv)=isub(7) twasubp(i+7*nsubdiv)=isub(8) 40 CONTINUE nnodes=0 DO i=1,8*nsubdiv IF(twasubp(i).NE.0) THEN nnodes=nnodes+1 twisubset(nnodes)=i ENDIF ENDDO DO 60 i=1,nnodes nodelist(i)=twasubp(twisubset(i)) 60 CONTINUE ENDDO ! DO WHILE (nnodes.GT.0) ! start new cycle with new nodelist(nnodes) RETURN END