c 220605 c if(encmin.eq.0) then nerby gravity +impacts ignored c better to do something with impacts then! c fft-array definations moved to ssbox.inc 31/05/2003 c contains individual eromax 6/4/2003 c gravinew567: c gmet=5 old c gmet=6 da/dt c gmet=7 gmet=6 + impact list separately c gmet=8 only far-gravity c uses also glim: delta > glim*ssum c gravinew2 -> evaluate the mutual gravity of c potentially colliding particles at dimpacts c -> excluded here c gravinew -> evaluate also dA/dt for nearby gravity c remove all except gmet=5 c include non-colliding tracer-particles c in this case use encmin=eromax + non-zero softening c ~hsalo/POTSDAM/SSBOX/OCT2000/ssgravi_d2.f c correct bug in defination of dyp with gmet=5 c add individual version-number to gidl and ggidl-files c remove velocity,xgrid from gidl output c add ss-output to gidl output c update ggidl-output (gmet=5) c _d2 -> include gravity from replicas (to gmet=5) c determined by mxgravi,mygravi c fz on particles outside grid: linear interpolation c 17.10.2000 c _d3 -> include possible gravity softening c esoft=max(encsoft,ssum) c*********************************************************** SUBROUTINE GRAVIGRID(ival) c*********************************************************** c calculate separately grid-gravity and c that due to nearby particles: HS 11/08/98 c every step: c 1) calculate direct forces c from particles with DIST < ENCMIN c 2) add contribution of gravity c from particles with ENCMIN < DIST < ENCLIM c using FFT on xyz-grid: nxg*nyg*nyz cells c inclined grid c 3) add vertical gravity due to infinite sheet c ival determines test-output c ival=0 - initialization c ival=1 - force profiles c ival=2 - force profiles + individual forces c IVAL=2 -> output forces to DGIDL-file, UNIT 30 c R,AGRAV_near c V,AGRAV_tot c xgrid,ygrid,zgrid, fxgrid,fygrid,fzgrid c -> calculation of gravity can be checked c IVAL=1,2 -> output grid-gravity cuts to DGGIDL, UNIT 31 c for 20 locations along y,z=0 c for 10 location along x,z=0 c for 60 locations along x,y=0 c cpugrav1 = CPU spent in sorting c cpugrav2 = CPU spent in tabulating distant gravity c cpugrav3 = CPU spent in calculation of nearby gravity c cpugrav4 = CPU spent in interpolation of distant gravity C*********************************************************** INCLUDE 'ssbox.inc' c for near gravity real xvar(nstar2),yvar(nstar2) dimension indx(nstar2) c for tests real xtest(160),ytest(160),ztest(160) real fxapu31(160),fyapu31(160),fzapu31(160) c dimension agrav_sheet(nstar) c for FFT c 3D: maximum true grid-size 512*256*8 c -> 176 mega c parameter (nxmax=512,nymax=256,nzmax=8) c parameter (nxmax=128,nymax=128,nzmax=32) c parameter (nzmax2=2*nzmax,ndatamax=4*nxmax*nymax*nzmax) c parameter declaration moved to include-file 310503 c real rhogrid3(nxmax,nymax,nzmax) c real fxgrid3(nxmax,nymax,nzmax) c real fygrid3(nxmax,nymax,nzmax) c real fzgrid3(nxmax,nymax,nzmax) c real frrho3(nxmax,nymax,nzmax2) c real firho3(nxmax,nymax,nzmax2) c real distx3(nxmax,nymax,nzmax2) c real disty3(nxmax,nymax,nzmax2) c real distz3(nxmax,nymax,nzmax2) c real frdist3(nxmax,nymax,nzmax2) c real fidist3(nxmax,nymax,nzmax2) real data(ndatamax) dimension nnn(3) VVAK=(KAPPA2-4.D0*OMEGA2)/OMEGA*XL CALL CPU(TIME1) encmin2=encmin**2 tgrid=0.d0 svak=bvak/2.d0/xl c interval for f_checks (in steps) i200=100 c*********************************************** c*********************************************** c TABULATE DISTANT GRAVITY c forces from xp,yp,zp density with FFT c*********************************************** c*********************************************** c encmin = limit for nearby gravity c enclim = limit for near+far gravity c encmin < enclim -> use fft for the zone from encmin to enclim svak=bvak/2.d0/xl n128x=nxg n128y=nyg n128z=2*nzg n64x=n128x/2 n64y=n128y/2 n64z=n128z/2 dxg=2.d0*xl/n128x dyg=2.d0*yl/n128y c----------------------------------------------- c vertical extent defined by fzrange c fzrange > 0 use zmax=fzrange*stdev(zz) c fzrange < 0 use zmax=abs(fzrange) apu=0.d0 apu2=0.d0 do 400 i=1,lkm apu2=apu2+r(i,3)**2 apu=apu+r(i,3) 400 continue apu=apu/lkm apu2=sqrt(apu2/lkm) apu_zdisp=apu2 apu_zmean=apu if(fzrange.gt.0) then zmax=fzrange*apu2 endif if(fzrange.lt.0) zmax=-fzrange zmax_stamp=zmax zmin=-zmax dzg=(zmax-zmin)/nzg c use this variable for storing zmax (output to dgidl and dggidl) tgrid=zmax c-------------------------------------------------------------- c moved here -> above test-quantities are always calculated if(encmin.lt.enclim) then CALL CPU(TTIME1) c----------------------------------------------- c clear density arrays if(time.lt.period/20.) then CALL CPU(TTTIME1) endif do 401 i=1,n128x do 402 j=1,n128y do 403 k=1,n128z rhogrid3(i,j,k)=0.d0 distx3(i,j,k)=0. disty3(i,j,k)=0. distz3(i,j,k)=0. 403 continue 402 continue 401 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-GRID',time/period, + tttime2-tttime1 write(6,*) rhogrid3(1,1,nzg/2),rhogrid3(1,1,nzg/2+1) CALL CPU(TTTIME1) endif c------------------------------------------------------- c calculate density rhogrid3 c cic=0 ngp c cic=2 cic cic=2 c ngp-assignment if(cic.lt.1) then do 410 i=1,lkm xp=r(i,1) yp=r(i,2)-svak*xp zp=r(i,3) if(yp.gt.yl) yp=yp-2.d0*yl if(yp.lt.-yl) yp=yp+2.d0*yl if(yp.gt.yl) yp=yp-2.d0*yl if(yp.lt.-yl) yp=yp+2.d0*yl if(xp.gt.xl) xp=xp-2.d0*xl if(xp.lt.-xl) xp=xp+2.d0*xl ixp=(xp+xl)/dxg+1 iyp=(yp+yl)/dyg+1 izp=(zp-zmin)/dzg+1 if(izp.ge.1.and.izp.le.n64z) + rhogrid3(ixp,iyp,izp)= + rhogrid3(ixp,iyp,izp)+gg*xmass(i) 410 continue endif c cic-assignment if(cic.gt.1) then do 411 i=1,lkm xp=r(i,1) yp=r(i,2)-svak*xp zp=r(i,3) if(yp.gt.yl) yp=yp-2.d0*yl if(yp.lt.-yl) yp=yp+2.d0*yl if(yp.gt.yl) yp=yp-2.d0*yl if(yp.lt.-yl) yp=yp+2.d0*yl if(xp.gt.xl) xp=xp-2.d0*xl if(xp.lt.-xl) xp=xp+2.d0*xl px=(xp+xl)/2./xl*nxg+n128x-0.5 py=(yp+yl)/2./yl*nyg+n128y-0.5 pz=(zp-zmin)/(zmax-zmin)*nzg+n128z-0.5 iix=int(px)-n128x iiy=int(py)-n128y iiz=int(pz)-n128z px=px-n128x py=py-n128y pz=pz-n128z dx=px-iix dy=py-iiy dz=pz-iiz tx=1.-dx ty=1.-dy tz=1.-dz ix1=iix+1 ix2=iix+2 iy1=iiy+1 iy2=iiy+2 iz1=iiz+1 iz2=iiz+2 if(ix1.gt.n128x) ix1=ix1-n128x if(ix2.gt.n128x) ix2=ix2-n128x if(iy1.gt.n128y) iy1=iy1-n128y if(iy2.gt.n128y) iy2=iy2-n128y if(ix1.lt.1) ix1=ix1+n128x if(ix2.lt.1) ix2=ix2+n128x if(iy1.lt.1) iy1=iy1+n128y if(iy2.lt.1) iy2=iy2+n128y if(ix2.gt.n128x.or.ix1.lt.1) then write(6,*) 'x-outgrid1',time/period,i, + r(i,1),ix1,ix2,n128x goto 199 endif if(iy2.gt.n128y.or.iy1.lt.1) then write(6,*) 'y-outgrid1',time/period,i, + r(i,2),iy1,iy2,n128y goto 199 endif c mass falling outside the grid if(iz2.gt.nzg) tz=0. if(iz2.lt.1) tz=0. if(iz1.gt.nzg) dz=0. if(iz1.lt.1) dz=0. c these just for convenience if(iz1.lt.1) iz1=1 if(iz2.gt.nzg) iz2=nzg if(iz2.lt.1) iz2=1 if(iz1.gt.nzg) iz1=nzg rhogrid3(ix2,iy2,iz2)= + rhogrid3(ix2,iy2,iz2)+gg*xmass(i)*dx*dy*dz rhogrid3(ix2,iy2,iz1)= + rhogrid3(ix2,iy2,iz1)+gg*xmass(i)*dx*dy*tz rhogrid3(ix2,iy1,iz2)= + rhogrid3(ix2,iy1,iz2)+gg*xmass(i)*dx*ty*dz rhogrid3(ix2,iy1,iz1)= + rhogrid3(ix2,iy1,iz1)+gg*xmass(i)*dx*ty*tz rhogrid3(ix1,iy2,iz2)= + rhogrid3(ix1,iy2,iz2)+gg*xmass(i)*tx*dy*dz rhogrid3(ix1,iy2,iz1)= + rhogrid3(ix1,iy2,iz1)+gg*xmass(i)*tx*dy*tz rhogrid3(ix1,iy1,iz2)= + rhogrid3(ix1,iy1,iz2)+gg*xmass(i)*tx*ty*dz rhogrid3(ix1,iy1,iz1)= + rhogrid3(ix1,iy1,iz1)+gg*xmass(i)*tx*ty*tz 199 continue 411 continue endif if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-RHO',tttime2-tttime1,cic write(6,*) rhogrid3(1,1,nzg/2),rhogrid3(1,1,nzg/2+1) CALL CPU(TTTIME1) endif c------------------------------------------------------- c construct distx3,disty3,distz3 c this makes sure that forces are symmetric if only one box is used i_yla=n64x-1 j_yla=n64y-1 mxgravi=0 mygravi=0 if(mxgravi.gt.0) i_yla=n64x if(mygravi.gt.0) j_yla=n64y do 420 i=-n64x+1,i_yla do 421 j=-n64y+1,j_yla dx=i*dxg dy=j*dyg+svak*dx apu2=dx*dx+dy*dy do 422 k=-n64z+1,n64z-1 dz=k*dzg if(apu2.gt.encmin2.and.apu2.le.enclim2) then apu=dx*dx+dy*dy+dz*dz ii=i+1 if(ii.lt.1) ii=ii+n128x jj=j+1 if(jj.lt.1) jj=jj+n128y kk=k+1 if(kk.lt.1) kk=kk+n128z apu=apu+encsoft2 apu3=1./sqrt(apu)/apu apux=dx*apu3 apuy=dy*apu3 apuz=dz*apu3 distx3(ii,jj,kk)=apux disty3(ii,jj,kk)=apuy distz3(ii,jj,kk)=apuz endif c _d2: if(mxgravi.gt.0.or.mygravi.gt.0) then do 1420 imx=-mxgravi,mxgravi do 1421 imy=-mygravi,mygravi if(imx.ne.0.or.imy.ne.0) then dx=i*dxg+2.*xl*imx dy=j*dyg+svak*dx+2.*yl*imy apu=dx*dx+dy*dy+dz*dz apu2=dx*dx+dy*dy if(apu2.gt.encmin2. + and.apu.le.enclim2) then apu=apu+encsoft2 apu3=1./sqrt(apu)/apu apux=dx*apu3 apuy=dy*apu3 apuz=dz*apu3 distx3(ii,jj,kk)= + distx3(ii,jj,kk)+apux disty3(ii,jj,kk)= + disty3(ii,jj,kk)+apuy distz3(ii,jj,kk)= + distz3(ii,jj,kk)+apuz endif endif 1421 continue 1420 continue endif 422 continue 421 continue 420 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-DIST',tttime2-tttime1 CALL CPU(TTTIME1) endif c------------------------------------------------------- c FFT of rhogrid3 ndim=3 nnn(1)=n128x nnn(2)=n128y nnn(3)=n128z lask=0 do 430 k=1,n128z do 431 j=1,n128y do 432 i=1,n128x lask=lask+1 data(lask)=rhogrid3(i,j,k) lask=lask+1 data(lask)=0. 432 continue 431 continue 430 continue ISIGN=+1 CALL FOURN(DATA,NNN,NDIM,ISIGN) lask=0 do 433 k=1,n128z do 434 j=1,n128y do 435 i=1,n128x lask=lask+1 frrho3(i,j,k)=data(lask) lask=lask+1 firho3(i,j,k)=data(lask) 435 continue 434 continue 433 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-FFT(RHO)',tttime2-tttime1 CALL CPU(TTTIME1) endif c---------------------------------- c only vertical forces? c 09/10/2000 <- 24/1/2000: not done if encz=1 if(encz.eq.1) then do 793 k=1,n64z do 794 j=1,n128y do 795 i=1,n128x fxgrid3(i,j,k)=0. fygrid3(i,j,k)=0. 795 continue 794 continue 793 continue endif if(encz.ne.1) then c------------------------------------------------------- c CALCULATION OF FXGRID3 c------------------------------------------------------- c FFT of distx3 lask=0 do 440 k=1,n128z do 441 j=1,n128y do 442 i=1,n128x lask=lask+1 data(lask)=distx3(i,j,k) lask=lask+1 data(lask)=0. 442 continue 441 continue 440 continue ISIGN=+1 CALL FOURN(DATA,NNN,NDIM,ISIGN) lask=0 do 443 k=1,n128z do 444 j=1,n128y do 445 i=1,n128x lask=lask+1 frdist3(i,j,k)=data(lask) lask=lask+1 fidist3(i,j,k)=data(lask) 445 continue 444 continue 443 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-FFT(distx3)',tttime2-tttime1 CALL CPU(TTTIME1) endif c------------------------------------------------------- c multiply fft(rhogrid3) by fft(dist3) c and insert into data lask=0 do 450 k=1,n128z do 451 j=1,n128y do 452 i=1,n128x lask=lask+1 data(lask)=frdist3(i,j,k)*frrho3(i,j,k)- + fidist3(i,j,k)*firho3(i,j,k) lask=lask+1 data(lask)=frdist3(i,j,k)*firho3(i,j,k)+ + fidist3(i,j,k)*frrho3(i,j,k) 452 continue 451 continue 450 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-MULTIPLY',tttime2-tttime1 CALL CPU(TTTIME1) endif c------------------------------------------------------- c obtain FX by inverse transform c just the real part isign=-1 CALL FOURN(DATA,NNN,NDIM,ISIGN) lask=0 do 453 k=1,n64z do 454 j=1,n128y do 455 i=1,n128x lask=lask+1 fxgrid3(i,j,k)=-data(lask)/n128x/n128y/n128z lask=lask+1 455 continue 454 continue 453 continue if(time.lt.period/20.) then CALL CPU(TTTIME2) write(6,*) 'GMET=5 CPU-FXGRID',tttime2-tttime1 CALL CPU(TTTIME1) endif c------------------------------------------------------- c CALCULATION OF FYGRID3 c------------------------------------------------------- c FFT of disty3 lask=0 do 460 k=1,n128z do 461 j=1,n128y do 462 i=1,n128x lask=lask+1 data(lask)=disty3(i,j,k) lask=lask+1 data(lask)=0. 462 continue 461 continue 460 continue ISIGN=+1 CALL FOURN(DATA,NNN,NDIM,ISIGN) lask=0 do 463 k=1,n128z do 464 j=1,n128y do 465 i=1,n128x lask=lask+1 frdist3(i,j,k)=data(lask) lask=lask+1 fidist3(i,j,k)=data(lask) 465 continue 464 continue 463 continue c------------------------------------------------------- c multiply fft(rhogrid3) by fft(dist3) c and insert into data lask=0 do 470 k=1,n128z do 471 j=1,n128y do 472 i=1,n128x lask=lask+1 data(lask)=frdist3(i,j,k)*frrho3(i,j,k)- + fidist3(i,j,k)*firho3(i,j,k) lask=lask+1 data(lask)=frdist3(i,j,k)*firho3(i,j,k)+ + fidist3(i,j,k)*frrho3(i,j,k) 472 continue 471 continue 470 continue c------------------------------------------------------- c obtain FY by inverse transform c just the real part isign=-1 CALL FOURN(DATA,NNN,NDIM,ISIGN) lask=0 do 473 k=1,n64z do 474 j=1,n128y do 475 i=1,n128x lask=lask+1 fygrid3(i,j,k)=-data(lask)/n128x/n128y/n128z lask=lask+1 475 continue 474 continue 473 continue endif c encz ne 1? c------------------------------------------------------- c CALCULATION OF FZGRID3 c------------------------------------------------------- c FFT of distz3 lask=0 do 480 k=1,n128z do 481 j=1,n128y do 482 i=1,n128x lask=lask+1 data(lask)=distz3(i,j,k) lask=lask+1 data(lask)=0. 482 continue 481 continue 480 continue ISIGN=+1 CALL FOURN(DATA,NNN,NDIM,ISIGN) lask=0 do 483 k=1,n128z do 484 j=1,n128y do 485 i=1,n128x lask=lask+1 frdist3(i,j,k)=data(lask) lask=lask+1 fidist3(i,j,k)=data(lask) 485 continue 484 continue 483 continue c------------------------------------------------------- c multiply fft(rhogrid3) by fft(dist3) c and insert into data lask=0 do 490 k=1,n128z do 491 j=1,n128y do 492 i=1,n128x lask=lask+1 data(lask)=frdist3(i,j,k)*frrho3(i,j,k)- + fidist3(i,j,k)*firho3(i,j,k) lask=lask+1 data(lask)=frdist3(i,j,k)*firho3(i,j,k)+ + fidist3(i,j,k)*frrho3(i,j,k) 492 continue 491 continue 490 continue c------------------------------------------------------- c obtain FZ by inverse transform c just the real part isign=-1 CALL FOURN(DATA,NNN,NDIM,ISIGN) lask=0 do 493 k=1,n64z do 494 j=1,n128y do 495 i=1,n128x lask=lask+1 fzgrid3(i,j,k)=-data(lask)/n128x/n128y/n128z lask=lask+1 495 continue 494 continue 493 continue CALL CPU(TTIME2) cpugrav2=cpugrav2+ttime2-ttime1 c encmin = limit for nearby gravity c enclim = limit for near+far gravity c encmin < enclim -> use fft for the zone from encmin to enclim c encmin exclude gravity of particles in impact list c**************************************************** c**************************************************** raja1_sum=0. raja1_lkm=0. raja1_max=0. ix=0 if(encmin.gt.1d-6) then c if(simustep.eq.1) then c write(97,*) 'ii,jj,nntemp,mmtemp,ssum,ero,dero,raja1_use' c write(97,*) 'dx,dy,dz,dvx,dvy,dvz' c endif c sort particles according to radius with padding CALL CPU(TTIME1) do 5 i=1,lkm xvar(i)=r(i,1) xvar(lkm+i)=r(i,1)+2.*xl yvar(i)=r(i,2) yvar(lkm+i)=r(i,2)+bvak 5 continue lkm2=2*lkm call indexr(lkm2,xvar,indx) CALL CPU(TTIME2) cpugrav1=cpugrav1+ttime2-ttime1 npotsum0=0 npotsum1=0 npotsum2=0 npotsum3=0 npotsum4=0 c calculate gravity + find nearby pairs CALL CPU(TTIME1) ix=0 rajatmp=encmin do 20 i=1,lkm ii=indx(i) xtemp1=xvar(ii) xrajatmp=xtemp1+rajatmp do 30 j=i+1,lkm2 jj=indx(j) xtemp2=xvar(jj) npotsum0=npotsum0+1 if(xtemp2.gt.xrajatmp) goto 20 mmtemp=99 dy=yvar(ii)-yvar(jj) if(abs(dy).lt.rajatmp) mmtemp=0 if(abs(dy+2.*yl).lt.rajatmp) mmtemp=1 if(abs(dy-2.*yl).lt.rajatmp) mmtemp=-1 if(mmtemp.eq.99) goto 30 nntemp=0 if(jj.gt.lkm) nntemp=-1 jjt=jj if(jj.gt.lkm) jjt=jj-lkm dx=r(ii,1)-r(jjt,1)+nntemp*avak dy=r(ii,2)-r(jjt,2)+nntemp*bvak+mmtemp*cvak dz=r(ii,3)-r(jjt,3) npotsum1=npotsum1+1 ddist2=dx*dx+dy*dy ddist=ddist2+dz*dz ssum=s(ii)+s(jjt) dist=sqrt(ddist) ero=dist-ssum c******************************************* c nearby gravitational forces if(ddist2.le.encmin2) then npotsum2=npotsum2+1 dist1=dist if(dist1.le.ssum*glim) dist1=ssum*glim if(encsoft.gt.0) dist1=sqrt(dist1**2+encsoft2) dvx=rp(ii,1)-rp(jjt,1) dvy=rp(ii,2)-rp(jjt,2)+nntemp*vvak dvz=rp(ii,3)-rp(jjt,3) c------------------------------------------------------------------ c gravinew567 c gmet=8: no nearby gravity c this is calculated in dimpacts if(gmet.eq.8) goto 9876 c-------------------------------- c individual eromax c used if particles have experienced enough impacts c (together more than coll_eropart impacts since last update) c (coll_eropart negative -> do not use) c OTHERWISE MAKE THE SAME WITH EROPART_OLD,EROCOLL_OLD C LAST REMEDY IS TO USE RAJA1 raja1_use=raja1 raja1_use_ori=0. coll_eropart=-10 if(coll_eropart.ge.0) then C normal case c----------------------------------------------------- if(collpart(ii)+collpart(jjt).ge. + coll_eropart) then raja1_use=eropart(ii)*SAFE_EROPART raja1_use2=eropart(JJT)*SAFE_EROPART if(RAJA1_USE2.gt.raja1_use) + raja1_use=RAJA1_USE2 raja1_use_ori=raja1_use if(raja1_use.lt.raja1min) + raja1_use=raja1min goto 1234 endif C from previous step c this helps to avoid spurious spikes in ix? c----------------------------------------------------- if(collpart_old(ii)+collpart_old(jjt).ge. + coll_eropart) then raja1_use=eropart_old(ii)*SAFE_EROPART raja1_use2=eropart_old(JJT)*SAFE_EROPART if(RAJA1_USE2.gt.raja1_use) + raja1_use=RAJA1_USE2 raja1_use_ori=raja1_use if(raja1_use.lt.raja1min) + raja1_use=raja1min endif endif 1234 continue raja1_sum=raja1_sum+raja1_use raja1_lkm=raja1_lkm+1. if(raja1_use_ori.gt.raja1_max) + raja1_max=raja1_use_ori c gmet=7: exclude pairs which are inserted to collision list if(gmet.eq.7) then if(ero.lt.raja1_use) then dero=EROVAK1*sqrt(dvx**2+dvy**2+dvz**2)*dt+ + EROVAK2*raja1_use if(ero.lt.dero) then c write(6,*) '*gmet7:excl',time,ii,jj,dist goto 9876 endif endif endif c------------------------------------------------------------------ if(simustep.lt.1) then if(ii.eq.10.or.jjt.eq.10) then if(ii.lt.jjt) write(97,*) 'PARI: ',ii,jjt if(jjt.lt.ii) write(97,*) 'PARI: ',jjt,ii 1010 format(4i8,4f8.3) 1011 format(6f8.3) write(97,*) 'ii,jjt,nn,mm,dx,dy,dz,dist,Ri,Rj' write(97,1010) + ii,jjt,nntemp,mmtemp,dx,dy,dz,dist1 write(97,1011) + r(ii,1),r(ii,2),r(ii,3), + r(jjt,1),r(jjt,2),r(jjt,3) endif endif FAPU=GG/dist1/dist1/dist1 FAPU1=-FAPU*XMASS(JJT) FAPU2= FAPU*XMASS(II) agrav(ii,1)=agrav(ii,1)+fapu1*dx agrav(ii,2)=agrav(ii,2)+fapu1*dy agrav(ii,3)=agrav(ii,3)+fapu1*dz agrav(jjt,1)=agrav(jjt,1)+fapu2*dx agrav(jjt,2)=agrav(jjt,2)+fapu2*dy agrav(jjt,3)=agrav(jjt,3)+fapu2*dz c energy due to over-border gravity c s=-shear c but nntemp=-1 -> negative sign if(nntemp.ne.0) then enes_grav=enes_grav-xmass(ii)*2.*xl*shear* + fapu1*dy*dt e_obg_size(is(ii))=e_obg_size(is(ii))- + xmass(ii)*2.*xl*shear*fapu1*dy*dt endif c gmet=5: do not calculate dA/dt if(gmet.eq.5) goto 9876 c------------------------------------------------------------------ c dA/dt only for gmet=6 and gmet=7 cv=dvx*dx+dvy*dy+dvz*dz apux=dvx-3.d0*cv*dx/dist1/dist1 apuy=dvy-3.d0*cv*dy/dist1/dist1 apuz=dvz-3.d0*cv*dz/dist1/dist1 aagrav(ii,1)=aagrav(ii,1)+fapu1*apux aagrav(ii,2)=aagrav(ii,2)+fapu1*apuy aagrav(ii,3)=aagrav(ii,3)+fapu1*apuz aagrav(jjt,1)=aagrav(jjt,1)+fapu2*apux aagrav(jjt,2)=aagrav(jjt,2)+fapu2*apuy aagrav(jjt,3)=aagrav(jjt,3)+fapu2*apuz c no gravity or dA/dt for this pair 9876 continue c--------------------------------------------------- c dist use fft for the zone from encmin to enclim if(encmin.lt.enclim) then CALL CPU(TTIME1) c**************************************************** c first normal particles c**************************************************** do 60 i=1,lkm fxapu=0. fyapu=0. fzapu=0. x0p=r(i,1) y0p=r(i,2)-svak*x0p z0p=r(i,3) if(y0p.gt.yl) y0p=y0p-2.*yl if(y0p.lt.-yl) y0p=y0p+2.*yl px=(x0p+xl)/2./xl*nxg+n128x-0.5 py=(y0p+yl)/2./yl*nyg+n128y-0.5 pz=(z0p-zmin)/(zmax-zmin)*nzg+n128z-0.5 iix=int(px)-n128x iiy=int(py)-n128y iiz=int(pz)-n128z px=px-n128x py=py-n128y pz=pz-n128z dx=px-iix dy=py-iiy dz=pz-iiz tx=1.-dx ty=1.-dy tz=1.-dz ix1=iix+1 ix2=iix+2 iy1=iiy+1 iy2=iiy+2 iz1=iiz+1 iz2=iiz+2 if(ix1.gt.n128x) ix1=ix1-n128x if(ix2.gt.n128x) ix2=ix2-n128x if(iy1.gt.n128y) iy1=iy1-n128y if(iy2.gt.n128y) iy2=iy2-n128y if(ix1.lt.1) ix1=ix1+n128x if(ix2.lt.1) ix2=ix2+n128x if(iy1.lt.1) iy1=iy1+n128y if(iy2.lt.1) iy2=iy2+n128y if(ix2.gt.n128x.or.ix1.lt.1) then write(6,*) 'x-outgrid',time/period,i, + r(i,1),ix1,ix2,n128x goto 99 endif if(iy2.gt.n128y.or.iy1.lt.1) then write(6,*) 'y-outgrid',time/period,i, + r(i,2),iy1,iy2,n128y goto 99 endif c mass falling outside the grid c 9.5.2002 if(iz2.gt.nzg) tz=0. if(iz2.lt.1) tz=0. if(iz1.gt.nzg) dz=0. if(iz1.lt.1) dz=0. if(iz1.lt.1) iz1=1 if(iz2.lt.1) iz2=1 if(iz1.gt.nzg) iz1=nzg if(iz2.gt.nzg) iz2=nzg fxapu= + dx*dy*dz*fxgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fxgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fxgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fxgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fxgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fxgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fxgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fxgrid3(ix1,iy1,iz1) fyapu= + dx*dy*dz*fygrid3(ix2,iy2,iz2)+ + dx*dy*tz*fygrid3(ix2,iy2,iz1)+ + dx*ty*dz*fygrid3(ix2,iy1,iz2)+ + dx*ty*tz*fygrid3(ix2,iy1,iz1)+ + tx*dy*dz*fygrid3(ix1,iy2,iz2)+ + tx*dy*tz*fygrid3(ix1,iy2,iz1)+ + tx*ty*dz*fygrid3(ix1,iy1,iz2)+ + tx*ty*tz*fygrid3(ix1,iy1,iz1) fzapu= + dx*dy*dz*fzgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fzgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fzgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fzgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fzgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fzgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fzgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fzgrid3(ix1,iy1,iz1) c set grid-force to zero if outside grid c add later that due to a constant-density sheet if(abs(z0p).gt.zmax) fzapu=0. 99 continue agrav(i,1)=agrav(i,1)+fxapu agrav(i,2)=agrav(i,2)+fyapu agrav(i,3)=agrav(i,3)+fzapu if(encz.eq.1) agrav(i,1)=0. if(encz.eq.1) agrav(i,2)=0. 60 continue if(simustep.le.5) then write(*,*) 'PP +FFT Forces on particle #10 and #100' write(*,*) 'A(10) ',agrav(10,1),agrav(10,2),agrav(10,3) write(*,*) 'A(100)',agrav(100,1),agrav(100,2),agrav(100,3) write(*,*) '------------------------------------------' endif c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c add checks: sum of vertical gravity should be zero c contribution from near+far particles: fz_check2 if(mod(simustep,i200).eq.0) then fx_check2=0. fy_check2=0. fz_check2=0. axmin=1d8 aymin=1d8 azmin=1d8 axmax=-1d8 aymax=-1d8 azmax=-1d8 n_zover=0 do 102 i=1,lkm fx_check2=fx_check2+xmass(i)*agrav(i,1) fy_check2=fy_check2+xmass(i)*agrav(i,2) fz_check2=fz_check2+xmass(i)*agrav(i,3) if(agrav(i,1).lt.axmin) axmin=agrav(i,1) if(agrav(i,2).lt.aymin) aymin=agrav(i,2) if(agrav(i,3).lt.azmin) azmin=agrav(i,3) if(agrav(i,1).gt.axmax) axmax=agrav(i,1) if(agrav(i,2).gt.aymax) aymax=agrav(i,2) if(agrav(i,3).gt.azmax) azmax=agrav(i,3) if(abs(r(i,3)).gt.zmax) n_zover=n_zover 102 continue fx_check2=fx_check2/tmass/omega2 fy_check2=fy_check2/tmass/omega2 fz_check2=fz_check2/tmass/omega2 axmin2=axmin/omega2 aymin2=aymin/omega2 azmin2=azmin/omega2 axmax2=axmax/omega2 aymax2=aymax/omega2 azmax2=azmax/omega2 write(6,1104) 'f_check2:',orb,fx_check2,fy_check2,fz_check2 write(6,1105) 'f_check2: zmean, zdisp, zmax, nov', + apu_zmean,apu_zdisp,zmax,n_zover write(6,1104) 'f_amin2: ',orb,axmin2,aymin2,azmin2 write(6,1104) 'f_amax2: ',orb,axmax2,aymax2,azmax2 endif c**************************************************** c tracer-particles c**************************************************** c**************************************************** c**************************************************** c output slices of force to unit 31 c ival=1 or ival=2 c**************************************************** c**************************************************** if(ival.eq.1.or.ival.eq.2) then if(mod(simustep,5).eq.0) then n30=90 write(31) time,omega,tgrid,igrav write(31) nxg,nyg,n30 write(31) xl,yl,enclim,encmin do 61 i=1,20 xtest(i)=-xl+0.05*(i-0.5)*2.*xl ytest(i)=0. ztest(i)=0. 61 continue do 62 i=1,10 xtest(20+i)=0. ytest(20+i)=-yl+0.10*(i-0.5)*2.*yl ztest(20+i)=0. 62 continue do 621 i=1,60 xtest(30+i)=0. ytest(30+i)=0. ztest(30+i)=zmin*2+(i-0.5)*(zmax-zmin)*2/30. 621 continue do 63 i=1,n30 x0p=xtest(i) y0p=ytest(i)-svak*x0p z0p=ztest(i) if(y0p.gt.yl) y0p=y0p-2.*yl if(y0p.lt.-yl) y0p=y0p+2.*yl px=(x0p+xl)/2./xl*nxg+n128x-0.5 py=(y0p+yl)/2./yl*nyg+n128y-0.5 pz=(z0p-zmin)/(zmax-zmin)*nzg+n128z-0.5 iix=int(px)-n128x iiy=int(py)-n128y iiz=int(pz)-n128z px=px-n128x py=py-n128y pz=pz-n128z dx=px-iix dy=py-iiy dz=pz-iiz tx=1.-dx ty=1.-dy tz=1.-dz ix1=iix+1 ix2=iix+2 iy1=iiy+1 iy2=iiy+2 iz1=iiz+1 iz2=iiz+2 if(ix1.gt.n128x) ix1=ix1-n128x if(ix2.gt.n128x) ix2=ix2-n128x if(iy1.gt.n128y) iy1=iy1-n128y if(iy2.gt.n128y) iy2=iy2-n128y if(ix1.lt.1) ix1=ix1+n128x if(ix2.lt.1) ix2=ix2+n128x if(iy1.lt.1) iy1=iy1+n128y if(iy2.lt.1) iy2=iy2+n128y c mass falling outside the grid if(iz2.gt.nzg) tz=0. if(iz2.lt.1) tz=0. if(iz1.gt.nzg) dz=0. if(iz1.lt.1) dz=0. if(iz1.lt.1) iz1=1 if(iz2.lt.1) iz2=1 if(iz1.gt.nzg) iz1=nzg if(iz2.gt.nzg) iz2=nzg fxapu= + dx*dy*dz*fxgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fxgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fxgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fxgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fxgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fxgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fxgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fxgrid3(ix1,iy1,iz1) fyapu= + dx*dy*dz*fygrid3(ix2,iy2,iz2)+ + dx*dy*tz*fygrid3(ix2,iy2,iz1)+ + dx*ty*dz*fygrid3(ix2,iy1,iz2)+ + dx*ty*tz*fygrid3(ix2,iy1,iz1)+ + tx*dy*dz*fygrid3(ix1,iy2,iz2)+ + tx*dy*tz*fygrid3(ix1,iy2,iz1)+ + tx*ty*dz*fygrid3(ix1,iy1,iz2)+ + tx*ty*tz*fygrid3(ix1,iy1,iz1) fzapu= + dx*dy*dz*fzgrid3(ix2,iy2,iz2)+ + dx*dy*tz*fzgrid3(ix2,iy2,iz1)+ + dx*ty*dz*fzgrid3(ix2,iy1,iz2)+ + dx*ty*tz*fzgrid3(ix2,iy1,iz1)+ + tx*dy*dz*fzgrid3(ix1,iy2,iz2)+ + tx*dy*tz*fzgrid3(ix1,iy2,iz1)+ + tx*ty*dz*fzgrid3(ix1,iy1,iz2)+ + tx*ty*tz*fzgrid3(ix1,iy1,iz1) c _d2 c if(abs(z0p).gt.zmax) fzapu=fzapu*abs(z0p)/zmax c _tracer c extrapolation: include forces beyond encmin c encmin < enclim if(abs(z0p).gt.zmax) then finf=2.*pii*gg*tmass/4./xl/yl fac=(1.+(encmin/z0p)**2)**(-.5) fzapu=-finf*fac*z0p/abs(z0p) endif if(abs(z0p).lt.zmax) then finf=2.*pii*gg*tmass/4./xl/yl fac=(1.+(enclim/z0p)**2)**(-.5) fzapu=fzapu-finf*fac*z0p/abs(z0p) endif fxapu31(i)=fxapu/omega2 fyapu31(i)=fyapu/omega2 fzapu31(i)=fzapu/omega2 63 continue write(31) (xtest(i),i=1,n30) write(31) (ytest(i),i=1,n30) write(31) (ztest(i),i=1,n30) write(31) (fxapu31(i),i=1,n30) write(31) (fyapu31(i),i=1,n30) write(31) (fzapu31(i),i=1,n30) endif endif CALL CPU(TTIME2) cpugrav4=cpugrav4+ttime2-ttime1 c encmin = limit for nearby gravity c enclim = limit for near+far gravity c encmin < enclim -> use fft for the zone from encmin to enclim c encmin no gravity due to infinite sheet c gsheet=1 -> add gravity due to infinite sheet c***************************************************************** c***************************************************************** apuapu=0.d0 apuapu2=0.d0 if(gsheet.ge.1) then c normal particles finf=2.d0*pii*gg*tmass/4.d0/xl/yl do 80 i=1,lkm agrav_sheet(i)=0.d0 if(r(i,3).ne.0) then rajainf=enclim if(abs(r(i,3)).gt.zmax) rajainf=encmin fac=1.d0/sqrt((1.d0+(rajainf/r(i,3))**2)) agrav(i,3)=agrav(i,3)-finf*fac*r(i,3)/dabs(r(i,3)) agrav_sheet(i)=-finf*fac*r(i,3)/dabs(r(i,3)) apuapu=apuapu-finf*fac*r(i,3)/dabs(r(i,3))*xmass(i) apuapu2=apuapu2+agrav(i,3)*xmass(i) endif 80 continue c tracer particles c for particles beyond grid: add forces from > encmin c for other particles: add forces from > enclim c add checks: sum of vertical gravity should be zero c contribution from near+far particles+distgrav: fz_check3 if(mod(simustep,i200).eq.0) then fz_check3=0. azmin=1d8 azmax=-1d8 do 103 i=1,lkm fz_check3=fz_check3+xmass(i)*agrav(i,3) if(agrav(i,3).lt.azmin) azmin=agrav(i,3) if(agrav(i,3).gt.azmax) azmax=agrav(i,3) 103 continue fz_check3=fz_check3/tmass/omega2 azmin3=azmin/omega2 azmax3=azmax/omega2 write(6,1104) 'f_check3:',orb,fz_check3 write(6,1104) 'f_amin3: ',orb,azmin3 write(6,1104) 'f_amax3: ',orb,azmax3 write(6,*) 'apuapu',apuapu/omega2/tmass,apuapu2/tmass endif c gsheet=1 ? endif c***************************************************************** c***************************************************************** c CORRECT MEAN Az to zero? c addition of distgrav can give a slightly non-zero mean acceleration, c which can gradually accumulate to give large non-zero mean z! c correct this c unless it is caused by the moonlet c if(n_moon.eq.0) then c do 81 i=1,lkm c agrav(i,3)=agrav(i,3)-apuapu2/tmass c 81 continue c endif c the correction is now made optional c gcorr=0 -> no correction c gcorr=1 -> correction if(gcorr.ge.1) then do 81 i=1,lkm agrav(i,3)=agrav(i,3)-apuapu2/tmass 81 continue if(mod(simustep,i200).eq.0) then fz_check4=0. azmin=1d8 azmax=-1d8 do 104 i=1,lkm fz_check4=fz_check4+xmass(i)*agrav(i,3) if(agrav(i,3).lt.azmin) azmin=agrav(i,3) if(agrav(i,3).gt.azmax) azmax=agrav(i,3) 104 continue fz_check4=fz_check4/tmass/omega2 azmin4=azmin/omega2 azmax4=azmax/omega2 write(6,1104) 'f_check4:',orb,fz_check4 write(6,1104) 'f_amin4: ',orb,azmin4 write(6,1104) 'f_amax4: ',orb,azmax4 write(6,*) '-------------------------------------------' endif endif 1104 format(a10,f12.3,3e20.6) 1105 format(a35,g14.3,2f12.3,i6) c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c**************************************************** c**************************************************** c OUTPUT TO GIDL: 30: dgravi1_read.pro c includes near+grid+sheet c add separately that due to sheet c**************************************************** c**************************************************** if(ival.eq.2) then write(30) (agrav(i,1)/omega2,i=1,lkm) write(30) (agrav(i,2)/omega2,i=1,lkm) write(30) (agrav(i,3)/omega2,i=1,lkm) c added with version 3 write(30) (agrav_sheet(i)/omega2,i=1,lkm) write(30) tgrid,igrav nnxg=nxg nnyg=nyg nnzg=nzg write(30) + (((fxgrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg) write(30) + (((fygrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg) write(30) + (((fzgrid3(i,j,k)/omega2,k=1,nnzg),j=1,nnyg),i=1,nnxg) write(30) cpugrav1,cpugrav2,cpugrav3,cpugrav4 endif c ************************************************** c ************************************************** c data collection (moved from lim1) c ************************************************** c ************************************************** C ******* TABULATE IXMIN ETC IF(IX.LT.IXMIN) IXMIN=IX IF(IX.GT.IXMAX) IXMAX=IX IXL=IXL+1 IXS=IXS+IX C ******* WARNING MESSAGE IF TOO MANY PAIRS IF(IX.GT..9*NIX) WRITE(IDEV6,*)'WARN: IX > ',0.9*NIX,IX CALL CPU(TIME2) CPUenc=CPUenc+TIME2-TIME1 c cpugrav1 = CPU spent in sorting c cpugrav2 = CPU spent in tabulating distant gravity c cpugrav3 = CPU spent in calculation of nearby gravity c cpugrav4 = CPU spent in interpolation of distant gravity if(mod(simustep,100).eq.0.or.simustep.le.5) then WRITE(*,*) 'GRAV: orb= ',orb c write(6,9049) ' GRAV: sort,tabu,nearby,far ', write(6,9049) ' GRAV: 1dtab,pp, grid, interp ', + cpugrav1,cpugrav3,cpugrav2,cpugrav4 c write(49,4949) orb,cpugrav1,cpugrav2,cpugrav3,cpugrav4, c + cpupota,cpuint,cputot WRITE(*,*) 'GRAV: npot0,1,2,3,4 = ', + npotsum0,npotsum1,npotsum2,npotsum3,npotsum4 endif 4949 format(f8.3,f6.1,2f10.1,f6.1,3f10.1) 9049 format(A30, 4g10.3) c---------------------------------------------------------- c GRAVITY FELT BY THE MOONLET IS SET TO ZERO c added 15/01/03 n_moon=0 if(n_moon.ne.0) then agrav(1,1)=0. agrav(1,2)=0. agrav(1,3)=0. aagrav(1,1)=0. aagrav(1,2)=0. aagrav(1,3)=0. endif RETURN END c--------------------------------------------------------------------- SUBROUTINE FOURN(DATA,NN,NDIM,ISIGN) double precision WR,WI,WPR,WPI,WTEMP,THETA dimension NN(NDIM),DATA(*) NTOT=1 DO 11 IDIM=1,NDIM NTOT=NTOT*NN(IDIM) 11 CONTINUE NPREV=1 DO 18 IDIM=1,NDIM N=NN(IDIM) NREM=NTOT/(N*NPREV) IP1=2*NPREV IP2=IP1*N IP3=IP2*NREM I2REV=1 DO 14 I2=1,IP2,IP1 IF(I2.LT.I2REV)THEN DO 13 I1=I2,I2+IP1-2,2 DO 12 I3=I1,IP3,IP2 I3REV=I2REV+I3-I2 TEMPR=DATA(I3) TEMPI=DATA(I3+1) DATA(I3)=DATA(I3REV) DATA(I3+1)=DATA(I3REV+1) DATA(I3REV)=TEMPR DATA(I3REV+1)=TEMPI 12 CONTINUE 13 CONTINUE ENDIF IBIT=IP2/2 1 IF ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) THEN I2REV=I2REV-IBIT IBIT=IBIT/2 GO TO 1 ENDIF I2REV=I2REV+IBIT 14 CONTINUE IFP1=IP1 2 IF(IFP1.LT.IP2)THEN IFP2=2*IFP1 THETA=ISIGN*6.28318530717959D0/(IFP2/IP1) WPR=-2.D0*DSIN(0.5D0*THETA)**2 WPI=DSIN(THETA) WR=1.D0 WI=0.D0 DO 17 I3=1,IFP1,IP1 DO 16 I1=I3,I3+IP1-2,2 DO 15 I2=I1,IP3,IFP2 K1=I2 K2=K1+IFP1 TEMPR=SNGL(WR)*DATA(K2)-SNGL(WI)*DATA(K2+1) TEMPI=SNGL(WR)*DATA(K2+1)+SNGL(WI)*DATA(K2) DATA(K2)=DATA(K1)-TEMPR DATA(K2+1)=DATA(K1+1)-TEMPI DATA(K1)=DATA(K1)+TEMPR DATA(K1+1)=DATA(K1+1)+TEMPI 15 CONTINUE 16 CONTINUE WTEMP=WR WR=WR*WPR-WI*WPI+WR WI=WI*WPR+WTEMP*WPI+WI 17 CONTINUE IFP1=IFP2 GO TO 2 ENDIF NPREV=N*NPREV 18 CONTINUE RETURN END