c***************************************************************************c c***************************************************************************c c subroutine caldist2d c c***************************************************************************c c***************************************************************************c subroutine caldist2d (drst,drmin,dset,istop,iwhere) implicit double precision (a-h,o-z) include 'header.h' parameter (c0=0.188030699d0,c1=1.48359853d0, & c2=-1.0979059d0,c3=0.430357353d0) parameter (one36=1.d0/36.d0) c*** c idifx=0 rcutnfc=rcellnfc rcutnf=rcellnf rcutrs=rcellrs ylo = -yc12 yhi = yc12 drtest=10.d0 c Loop over each particle pair do 10 kk=1,nsq1 ii = id(1,kk) jj = id(2,kk) CB****don't calculate distance btwn two wall particles! if (jj.le.npwall)then if (ii.le.npwall) then goto 10 endif endif CB****don't calculate distance btwn two stopped agglomerates c$$$ if (iwhat(jj).eq.2) then c$$$ if (iwhat(ii).eq.2) then c$$$c dx=2.01d0 c$$$ goto 10 c$$$ endif c$$$ endif CB********* iz = ii*6 jz = jj*6 dx = x(jz-5)-x(iz-5) dy = x(jz-4)-x(iz-4) dz = x(jz-3)-x(iz-3) c*** correct particle minimum image for free particles zlo = -zc12 zhi = zc12 if (dz.gt.zhi) then dz = dz - zcell c write(6,*)'here1' elseif (dz.lt.zlo) then dz = dz + zcell c write(6,*)'here2' endif if (dy.gt.yhi) then dy = dy - ycell c dx = dx - hgt c write(6,*)'here3' elseif (dy.lt.ylo) then dy = dy + ycell c dx = dx + hgt c write(6,*)'here4' endif c xlo =-xc12 !+ dy*gamt*gamma xhi = xlo + xcell if (dx.gt.xhi) then dx = dx - xcell c write(6,*)'here5' elseif (dx.lt.xlo) then dx = dx + xcell c write(6,*)'here6' endif drsq = dx*dx + dy*dy + dz*dz dr = dsqrt(drsq) c write (6,*) dr c particles in the pair overlapped (drset<2.0) !CB if (dr.lt.drset) then if (iwhat(jj).eq.0) then if (iwhat(ii).eq.0) then write (15,201)ii,jj,dr,rtime,iwhere,iwhat(ii),iwhat(jj) write (15,'(6(1x,e12.4))') x(iz-5),x(iz-4),x(iz-3), & x(jz-5),x(jz-4),x(jz-3) write(*,*)"OVERLAP! line77",rtime,ii,jj,dr c$$$ commented out, particles are overlapping, they might get apart again c$$$ c iwhat(jj)=2 c$$$ c if (iwhat(ii).eq.0) then c$$$ c iwhat(ii)=2 c$$$ c endif c WRITE(6,*) 'drltdrset,ii,jj',ii,jj,iwhat(ii),iwhat(jj) c istop = 0 c STOP c return endif else !(ii is 'wall' or agglomerate) write (15,201)ii,jj,dr,rtime,iwhere,iwhat(ii),iwhat(jj) write (15,*) ii,jj,'ii is not free, iwhat/=0' c write (6,*) ii,jj,'ii is not free' endif endif !CB c particles in the pair to be overlapp (->2.0) c for St=\0 stop running;for St=0, set the distance to be dset0=2+epsilon c note: the physical positions of particles are not moved, lubrication forces c would separate them eventually. 33 if (dr+tiny.lt.2.d0) then if(stk .gt. tiny) then !for St bigger than zero write(15,201) ii,jj,dr,rtime,iwhere,iwhat(ii),iwhat(jj) write(15,'(6(1x,e12.4))') x(iz-5),x(iz-4),x(iz-3), & x(jz-5),x(jz-4),x(jz-3) istop = 0 return else !i.e. St=0 write (15,201)ii,jj,dr,rtime,iwhere,iwhat(ii),iwhat(jj) c WRITE(6,201)ii,jj,dr,rtime,iwhere,iwhat(ii),iwhat(jj) c WRITE(6,*) 'drlt2.0ii,jj',ii,jj dx=dx*dset0/dr dy=dy*dset0/dr dz=dz*dset0/dr dr = dset0 goto 11 endif endif c particles in the pair reaching minimumal allowed separation if (dr+tiny.lt.dset) then if (drmin.gt.dr+tiny) then drmin = dr c iwhat(ii)=2 c iwhat(jj)=2 write (15,'(1x,f8.4,2(1x,i5),2(1x,e15.8),i2)') > rtime,ii,jj,drmin,dset,iwhere endif endif 11 continue drtest=min(drtest,dr) driv = 1.0d0/dr dis(1,kk) = dx*driv dis(2,kk) = dy*driv dis(3,kk) = dz*driv dis(4,kk) = dr if (drmin.gt.dr) drmin=dr 10 continue c pairs within the distance of real summation cut off radius jjrs=1 do 20 kk=1,nsq1 c write(*,*)dis(4,kk),rcutrs,jjrs,kk if (dis(4,kk).le.rcutrs) then if (jjrs.ne.kk) then id(1,jjrs)=id(1,kk) id(2,jjrs)=id(2,kk) dis(1,jjrs)=dis(1,kk) dis(2,jjrs)=dis(2,kk) dis(3,jjrs)=dis(3,kk) dis(4,jjrs)=dis(4,kk) c** NOT YET RESOLVED c if(id(1,jjrs).eq.125.or.id(2,jjrs).eq.125)then c write(*,*)jjrs,id(1,jjrs),id(2,jjrs),dis(4,jjrs) c endif endif jjrs=jjrs+1 endif 20 continue nsq1=jjrs-1 c pairs within the distance of near field cut off radius jjnf=1 do kk=1,nsq1 if (dis(4,kk).le.rcutnf) then idnf(1,jjnf)=id(1,kk) idnf(2,jjnf)=id(2,kk) disnf(1,jjnf)=dis(1,kk) disnf(2,jjnf)=dis(2,kk) disnf(3,jjnf)=dis(3,kk) disnf(4,jjnf)=dis(4,kk) jjnf=jjnf+1 if(jjnf .gt. nn2nf) then write(15,'(''nsq1nf'',i3,''>nn2nf'',i3)') jjnf,nn2nf istop = 0 return endif endif enddo nsq1nf=jjnf-1 c pairs within a distance smaller than near field cut off radius c supposingly facilitate the search of colliding pairs jjnfc=1 do kk=1,nsq1nf if (disnf(4,kk).le.rcutnfc) then idnfc(1,jjnfc)=idnf(1,kk) idnfc(2,jjnfc)=idnf(2,kk) disnfc(1,jjnfc)=disnf(1,kk) disnfc(2,jjnfc)=disnf(2,kk) disnfc(3,jjnfc)=disnf(3,kk) disnfc(4,jjnfc)=disnf(4,kk) jjnfc=jjnfc+1 if(jjnfc .gt. nn2nfc) then c write(70,'(''nsq1nfc'',i3,''>nn2nfc'',i3)') c > jjnfc,nn2nfc write(15,'(''nsq1nfc'',i3,''>nn2nfc'',i3)') > jjnfc,nn2nfc istop = 0 return endif endif enddo nsq1nfc=jjnfc-1 c set up near field pair list call setnflist(istop) return 201 format('olp ii/jj/dr/t/w',2(i5,2x),f15.10,2x,f7.3,i2) end c***************************************************************************c c***************************************************************************c