c----------------------------------------- c simple030707.f = cleaned version c most of input paramaters removed c initial particle positions/velocities read from inipos c most of output subroutines removed c----------------------------------------- c*********************************************************************** PROGRAM simple c*********************************************************************** c Local Simulation Code (HS 18.02.99 ->) C*********************************************************************** INCLUDE 'ssbox.inc' real xani(nstar),yani(nstar),zani(nstar) character*80 dummyline double precision dummytab(500) CALL CPU(TALKU) codeid='simple_030707.f' c----------------------------------------------------------------------- c READ INPUT VARIABLES c------------------------------------------------------ c RUNID = run identification write(*,*) 'RUNID' read(*,*) runid do 1001 i=1,80 if(runid(i:i).eq.'') goto 1002 1001 continue 1002 nlen=i-1 write(*,*) RUNID,nlen c------------------------------------------------------ c Particle data from 'iniposnew'-file (fixed name) c This has the same format as the finposrot files c of the old version + some new data c time = starting time c lkm = number of particles c r(lkm,3) = particle position vectors c rp(lkm,3) = particle velocity vectors c myy(lkm) = inertial mass c s(lkm) = particle radius c xmass(lkm) = gravity mass c rota(lkm,3)= spin velocity (in finposrot) c new additions: c is(lkm) = size group index c eta = distance from saturn c kapparat,wzrat = kappa/omega, omegaz/omega c xl,yl = radial, tangential semiwidth c All values are given in SI units open(unit=50,file='iniposnew',status='old') read(50,*) time,lkm write(*,*) 'READ INITIAL POSITIONS FROM iniposnew' write(*,*) 'time,lkm',time,lkm read(50,*) (r(i,1),r(i,2),r(i,3),i=1,lkm) read(50,*) (rp(i,1),rp(i,2),rp(i,3),i=1,lkm) read(50,*) (myy(i),i=1,lkm) read(50,*) (s(i),i=1,lkm) read(50,*) (xmass(i),i=1,lkm) read(50,*) (ROTA(i,1),ROTA(i,2),ROTA(i,3),i=1,lkm) read(50,*) (is(i),i=1,lkm) read(50,*) eta read(50,*) kappafac,wzfac read(50,*) xl,yl write(6,*) 'inipos-READ SUCCESFULL' close(50) c------------------------------------------------------ c Internal density roo c Already set by xmass and s --> recalculate here c This allows easier testing of the influence of roo c TODO: same with eta ? read(*,*) roo c----------------------------------------------------- c calculate some constants C ******* CONSTANTS RELATED TO IMAGE-PARTICLES C X -> X + NN*AVAK C Y -> Y + NN*BVAK +MM*CVAK C VX -> VX C VY -> VY + NN*VVAK C ******* EQUATIONS OF MOTION C XPP - 2*OMEGA * YP + FACTEMP * X = FX C YPP + 2*OMEGA * XP = FY C ZPP + WZ2 * Z = FZ W0=DSQRT(GM/ETA/ETA/ETA) OMEGA=W0 W2=2.D0*OMEGA KAPPA=KAPPAFAC*OMEGA OMEGA2=OMEGA*OMEGA KAPPA2=KAPPA*KAPPA WZ=W0*WZFAC WZ2=WZ*WZ PERIOD=2.D0*PII/OMEGA FACTEMP=KAPPA2-4.D0*OMEGA2 SHEAR=-FACTEMP/2.d0/OMEGA AVAK=2.D0*XL CVAK=2.D0*YL VVAK=FACTEMP/OMEGA*XL BVAK=VVAK*TIME time0=time ORB=TIME/PERIOD orbori=orb do 1020 i=1,lkm xmass(i)=4.*pii/3.*roo*s(i)**3 myy(i)=xmass(i) if(roo.le.0) + myy(i)=4.*pii/3.*900.*s(i)**3 1020 continue tot_area=0.d0 tot_mass=0.d0 do 1030 i=1,lkm tot_mass=tot_mass+4.*pii/3.*roo*s(i)**3 tot_area=tot_area+pii*s(i)**2 1030 continue tautot=tot_area/(4.*xl*yl) sigmatot=tot_mass/(4.*xl*yl) c previous variable names tmass=tot_mass taumin=tautot tau0=taumin write(*,*) 'lkm=',lkm write(*,*) 'xl,yl=',xl,yl write(*,*) 'smin,smax',s(1),s(lkm) write(*,*) 'roo=',roo write(*,*) 'tau =',tautot write(*,*) 'Sigma=',sigmatot c---------------------------------------------------------------- c VARIABLES RELATED TO INTEGRATION c dt_orb = timestep for impact lists/gravity update c (in orbital periods) c orbmax = length of integration read(*,*) dt_orb, orbmax dt=dt_orb*period c---------------------------------------------------------------- c VARIABLES RELATED TO IMPACTS read(*,*) dummyline c FORCE-MODEL PARAMETERS C alfstep: divide each timestep into alfstep integ steps C alfome: frequency of harmonic force (in OMEGA'S) read(*,*) alfstep,alfome alfome=alfome*omega C ELASTICITY MODELS c IMODEL determines elasticity law EPS_N(v) c LAMDA and V0 parameters of EPS_N c BETA = 1-EPS_T tangential coefficient (actually limits friction) c use beta = 2 -> can not have effect c VCRITA set eps_n to 1 below this V_normal (WT88-trick) read(*,*) IMODEL,LAMDA,V0,BETA,VCRITA C FRICTION c lfric= coefficient of friction (negative) c fricmodel = model for friction 23.1.03 c 0 -> as before mu_fric=lfric, with beta limitation c 1 -> as above, but zigzag prevented c 2 -> mu_fric chosen to yield beta, ft=mu_fric*fn c 3 -> kmu_fric chosen to yield beta, ft=kmu_fric*gt*fn read(*,*) lfric, fricmodel c ADHESION read(5,*) gamma_adh c TESTING: IGNORE SOME IMPACTS? c nob_imp = 1 -> ignore impacts over all boundaries c 2 -> ignore impacts over radial boundaries c 3 -> ignore impacts over tangential boundaries c 4 -> ignore all impacts c 0 -> calculate all impacts (normal case) read(*,*) nob_imp c-------------------------------------------------------------- c VARIABLES RELATED TO THE SEARCH OF NEARBY PAIRS c RAJA1 = limiting distance for including into close pairs c initial value, updated automatically during the run c RAJA1MIN, RAJA1MAX: limits for RAJA1 c Additional adjustments: pair is included to impact list if c ero < EROVAK1*sqrt(dvx**2+dvy**2+dvz**2)*dt+EROVAK2*raja1 c where ero= dist-ssum (distance between surfaces) read(*,*) dummyline read(*,*) raja1, raja1min, raja1max read(*,*) erovak1, erovak2 c if erovak1=0 and erovak2=1 -> original condition ero old method using 1-dim ordering read(*,*) mx_mc,my_mc c---------------------------------------------------------------- c VARIABLES RELATED TO SELFGRAVITY c Method of gravity calculation: routines gravi, nongravi c (same routines make the calculation impact lists) c gmet=2 --> no self-gravity c gmet=5,6 --> with self-gravity c encmin = distance inside which particle-particle method c is used (accurate gravitational encounters) c enclim = limiting radii for gravity calculation c between encmin - enclim, FFT is used c enclim must be le min(xl,yl) c enclim < 0 -> given relative to min(xl,yl) c nxg,nyg,nzg = gravity grid c fzrange = vertically grid extends fzrange*stdev(zz) c encsoft = explicit gravity softening (typically zero) c test the effect of vertical gravity: c encz=1 -> only vertical gravity c encz=0 -> all components (normally choice) c avoid errors due to missed impacts: c glim>0 -> limit gravity to that with delta=ssum*glim c glim=0 -> no limit c gravity of infinite sheet outside enclim? c calculation made optional c gsheet = 1 -> vertical gravity due infinite sheet included c 0 -> not included c gcorr = 1 -> mean vertical gravity set to zero c 0 -> no setting read(*,*) dummyline read(*,*) gmet read(*,*) encmin,enclim read(*,*) nxg,nyg,nzg read(*,*) fzrange read(*,*) encsoft read(*,*) encz read(*,*) glim read(*,*) gsheet,gcorr if(enclim.lt.0) then enclim=-enclim*xl else if(enclim.gt.xl) enclim=xl if(enclim.gt.yl) enclim=yl endif enclim2=enclim*enclim encsoft2=encsoft*encsoft c---------------------------------------------------------------- c VARIABLES RELATED TO OUTPUT c (output intervals given in orbital periods) read(*,*) dummyline c VELOCITY DISPERSION (nidl-files) c new output c dt_nidl = collection interval c out_nidl = output interval (use this later) read(*,*) dt_nidl, out_nidl open(unit=35,file=runid(1:nlen)//'.nidl',status='new') call nidl_data(0) c POSITION SNAPSHOTS (pidl-files) c dt_pidl = output interval for pidl-files read(*,*) dt_pidl open(unit=25,file=runid(1:nlen)//'.pidl',status='new') call pidl_data(0) c TABULATION OF RADIAL FOURIER MODES (ufidl-files) c i60 number of tabulated modes c dt_ufidl = interval for output c dummytab needed to be compaible with old format do 1100 i=1,500 dummytab(i)=0.d0 1100 continue read(*,*) i60,dt_ufidl open(unit=71,file=runid(1:nlen)//'.ufidl',status='new', + form='unformatted') versio=1. write(71) versio write(71) i60 write(71) (dummytab(i),i=1,i60) write(71) (dummytab(i),i=1,i60) write(71) xl,taumin,lkm,omega c ANIMATION OUTPUT (aniidl-files) c collect output for making fast animations c animin: start at orbit c animax: stop at orbit c dt_ani: output interval in orbital periods read(*,*) animin,animax,dt_ani if(dt_ani.gt.0) then open(unit=90,file=runid(1:nlen)//'.aniidl', + status='new',form='unformatted') aniversio=0. write(90) aniversio write(90) animin,animax,dt_ani endif c CHECKING OF GRAVITY CALCULATIONS c dt_gidl = output interval for gidl-files read(*,*) dt_gidl open(unit=30,file=runid(1:nlen)//'.dgidl',status='new', + form='unformatted') open(unit=31,file=runid(1:nlen)//'.dggidl',status='new', + form='unformatted') versio=1 write(30) versio write(31) versio c FINPOS-FILE c dt_finpos = output interval for finposnew-files c (can be used for continuation) read(*,*) dt_finpos c--------------------------------------------------------------- c TIMESTAMP - output run info after every step c--------------------------------------------------------------- call timestamp(0) C************************************************************ C INITIALIZE COLLISION LIST COUNTERS IXL=0 IXS=0 IXMIN=100000 IXMAX=0 vero1max=0. ERO1S=0. ERO1L=0. msis=0 coll=0 C************************************************************ C INITIALIZE CPU-COUNTERS C CPU_RET = RETURN PARTICLES TO BOX C CPU_LIM = CONSTRUCTION OF LIST OF NEARBY PARTICLES c in dimpacts: C CPU_SAV = bookkeeping of ongoing impacts C CPU_POTA = construction of forces for overlapping particles C CPU_INT = integration of orbits C CPU_POTB = overhead in dimpacts C CPU_SG = CALCULATION OF GRAVITATIONAL FORCES C CPU_TOT = TOTAL CPU C 0 REFERS TO ITECHECK-OUTPUT PERIOD CPU_RET=0. CPU_LIM=0. CPU_SAV=0. CPU_POTA=0. CPU_INT=0. CPU_POTB=0. CPU_SG=0. CPU_TOT=0. CPU_RET0=0. CPU_LIM0=0. CPU_SAV0=0. CPU_POTA0=0. CPU_INT0=0. CPU_POTB0=0. CPU_SG0=0. CPU_TOT0=0. C*********************************************** C INITIALIZE OUTPUT-COUNTERS c previous animation storage orbold_ani=orb animui=1 c ufidl-tabulation orbold_ufidl=orb c finposnew orbold_finpos=orb c itecheck orbold_raja1=orb c velocity dispersion c new nidl output orbold_nidl=orb C ******************************************** C SET MISCELLENOUS VARIABLES ALFSUM=0. XL1=-XL XL2=XL XL0=0.0D0 YL0=0.0D0 c radial boundary crossing c nha = inward, nhy = outward nha=0 nhy=0 c tangential boundary crossing npa=0 npy=0 c-------------------------------------------------------------------- c INITIAL OUTPUT TO FILES c-------------------------------------------------------------------- c CONTROL FILE c write some input-file parameters to control-file c these will be checked on each time step c --> offers possibility to terminate the run open(unit=95,file=runid(1:nlen)//'.control',status='new') write(95,*) orbmax write(95,*) raja1min close(95) c------------------------------------------------------------------- c store position file call pidl_data(1) c radial fourier amplitudes call ufidl_data orbold_ufidl=orb c velocity dispersion call nidl_data(1) c---------------------------------------------------------- call timestamp(1) c---------------------------------------------------------- C*********************************************************************** C MAIN SIMULATION-LOOP C*********************************************************************** simustep=-1 1 CONTINUE ORB=TIME/PERIOD simustep=simustep+1 c------------------------------------------------- C CALCULATE GRAVITATIONAL FORCES C LIST OF POTENTIAL IMPACTORS c------------------------------------------------- c return particles to the calculation region CALL RETBOX c clear force arrays CALL ENCZERO c gmet=2 -> no sg, just collision lists c gmet=5,6 -> sg + collision lists c 5 -> previous method A=constant c 6 -> A=A0+dA/dt*dt if(gmet.eq.2) CALL NONGRAVI C GRAVITY CHECK? c ival=1 -profiles ival=2 profiles, individual forces IVAL_G=1 IF((ORB-ORBOLD_gidl).GE.0.999*dt_gidl) THEN ival_g=2 ORBOLD_gidl=ORB ENDIF if(gmet.ge.5) then if(mx_mc*my_mc.eq.0) CALL GRAVIGRID(ival_g) if(mx_mc*my_mc.gt.0) CALL GRAVIGRID_mc(ival_g) endif c------------------------------------------ c IMPACTS + integration c------------------------------------------ call dimpacts_rkhf4 c------------------------------------------ c UPDATE RAJA1 c output various counters c------------------------------------------ IF((orb-orbold_raja1).ge.0.99999*dt_raja1) then orbold_raja1=orb call itecheck ENDIF c------------------------------------------ c OUTPUT c------------------------------------------ c output current state of the system (unit 77) call timestamp(2) c ADD: velocity dispersion IF((orb-orbold_nidl).ge.0.99999*dt_nidl) then orbold_nidl=orb call nidl_data(1) ENDIF c projection files IF((orb-orbold_pidl).ge.0.99999*dt_pidl) then orbold_pidl=orb call retbox call pidl_data(1) ENDIF c radial fourier components if((orb-orbold_ufidl).ge.0.99999*dt_ufidl) then orbold_ufidl=orb call ufidl_data endif c animation files if(orb.ge.animin.and.orb.le.animax) then if((orb-orbold_ani).ge.0.99999*dt_ani) then call retbox do 9090 ii=1,lkm xani(ii)=r(ii,1) yani(ii)=r(ii,2) zani(ii)=r(ii,3) 9090 continue write(90) orb,omega,xl,yl,lkm write(90) (xani(ii),ii=1,lkm) write(90) (yani(ii),ii=1,lkm) write(90) (zani(ii),ii=1,lkm) orbold_ani=orb animui=2 endif endif if(animui.eq.2.and.orb.ge.animax) then animui=3 close(90) endif c--------------------------------------------------- c WRITE FINPOSNEW-FILE (for continuation)? c--------------------------------------------------- if((orb-orbold_finpos).ge.0.999*dt_finpos) then orbold_finpos=orb call retbox open(unit=50,file='finposnew',status='unknown') write(6,*) 'OUTPUT POSITIONS to finposnew' write(6,*) 'time,lkm',time,lkm write(50,*) time ,lkm write(50,*) (r(i,1),r(i,2),r(i,3),i=1,lkm) write(50,*) (rp(i,1),rp(i,2),rp(i,3),i=1,lkm) write(50,*) (myy(i),i=1,lkm) write(50,*) (s(i),i=1,lkm) write(50,*) (xmass(i),i=1,lkm) write(50,*) (rota(i,1),rota(i,2),rota(i,3),i=1,lkm) write(50,*) (is(i),i=1,lkm) write(50,*) eta write(50,*) kappafac,wzfac write(50,*) xl,yl close(50) endif c--------------------------------------------------------- c CHECK THE CONTROL FILE c (edit orbmax --> force run to terminate) c--------------------------------------------------------- open(unit=95,file=runid(1:nlen)//'.control') read(95,*) orbmax_new read(95,*) raja1min_new if(orbmax_new.ne.orbmax) then write(*,*) 'CONTROL-FILE: new orbmax',orbmax_new orbmax=orbmax_new endif if(raja1min_new.ne.raja1min) then write(*,*) 'CONTROL-FILE: new raja1min',raja1min_new raja1min=raja1min_new endif close(95) C*************************************************************** C TERMINATE RUN AFTER ORBMAX c (also if there are too many impact pairs) IF(ORB.GE.ORBMAX) GOTO 2 IF(IX.GT..9*NIX) then write(*,*) 'too many IX pairs',IX goto 2 endif GOTO 1 2 CONTINUE CALL TIMESTAMP(3) close(71) close(25) close(35) close(90) END C************************************************************** SUBROUTINE ENCZERO C************************************************************** INCLUDE 'ssbox.inc' c ***** set agrav to zero do 30 i=1,lkm agrav(i,1)=0.d0 agrav(i,2)=0.d0 agrav(i,3)=0.d0 c in connection with gravinew: aagrav(i,1)=0.d0 aagrav(i,2)=0.d0 aagrav(i,3)=0.d0 30 continue return end C*********************************************************************** C SUBROUTINE pidl_data(IDUM) c STORES PARTICLE POSITIONS TO pidl-file c idum=0 --> initial output c idum=1 --> snapshots C*********************************************************************** INCLUDE 'ssbox.inc' dimension rtemppi(3) idev=25 c initial output to pidl-file if(idum.eq.0) then IVERSIO=3 WRITE(IDEV,*) IVERSIO WRITE(IDEV,*) ENcsoft,ENClim,ENCZ etatkm=eta/1d6 nind=0 mjak=0 qq=0 nsig=0 ceks=0 cink=0 qzinit=0 idtmax=0 tsmin=0 tsmax=0 WRITE(IDEV,*) LKM,ETATKM,S(1),S(LKM),MJAK,QQ,NSIG WRITE(IDEV,*) IMODEL,LAMDA,V0,BETA,VCRITA WRITE(IDEV,*) CEKS,CINK,QZINIT,YL,XL,TAU0,TAUMIN WRITE(IDEV,*) TIME0,DT,IDTMAX,TSMIN,TSMAX WRITE(IDEV,*) RAJA3,RAJA1,xraja,RAJA1MIN,RAJA3MIN WRITE(IDEV,*) COLLMAX,COLLS,COLLL WRITE(IDEV,*) ORBMAX,ORBD,ORBE WRITE(IDEV,*) SIGTAR,SIGTER,WZFAC,KAPPAFAC WRITE(IDEV,*) DZCUM,NZCUM,ROO WRITE(IDEV,*) ORBZ1,ORBZ2,ORBSC WRITE(IDEV,*) NIX,NJX,NIND,NSTAR return endif c snapshot output to pidl-file ANGLE=W0*TIME SINA=dSIN(ANGLE) COSA=dCOS(ANGLE) WRITE(IDEV,3535) COLL,TIME,LKM 3535 FORMAT(I12,F20.8,I12) DO 333 I=1,LKM rtemppi(1)=r(i,1) rtemppi(2)=r(i,2) rtemppi(3)=r(i,3) if(rtemppi(1).gt.9999.) rtemppi(1)=9999. if(rtemppi(2).gt.9999.) rtemppi(2)=9999. if(rtemppi(3).gt.9999.) rtemppi(3)=9999. if(rtemppi(1).lt.-9999.) rtemppi(1)=-9999. if(rtemppi(2).lt.-9999.) rtemppi(2)=-9999. if(rtemppi(3).lt.-9999.) rtemppi(3)=-9999. WRITE(IDEV,3536) IS(I),S(I),(Rtemppi(j),j=1,3), + (RP(I,J),J=1,3) c spins are always expressed in rotating system if(beta.ne.0) then rosig1=rota(i,1)/s(i)/w0 rosig2=rota(i,2)/s(i)/w0 rosig3=rota(i,3)/s(i)/w0 WRITE(IDEV,3537) rosig1,rosig2,rosig3 endif 333 CONTINUE 3536 FORMAT(I3,f10.5,3f10.3,3e12.4) 3537 FORMAT(3E12.4) RETURN END C*********************************************************************** C SUBROUTINE nidl_data(IDUM) c store velocity dispersion etc c idum=0 --> initial output c idum=1 --> snapshots C*********************************************************************** INCLUDE 'ssbox.inc' idev=35 c initial output to nidl-file if(idum.eq.0) then IVERSIO=1 WRITE(IDEV,*) IVERSIO return endif dispvx=0. dispvy=0. dispvz=0. dispvxy=0. do 2001 I=1,lkm vxt=rp(i,1) vyt=rp(i,2)+shear*r(i,1) vzt=rp(i,3) dispvx=dispvx+vxt**2 dispvy=dispvy+vyt**2 dispvz=dispvz+vzt**2 dispvxy=dispvxy+vxt*vyt 2001 continue c dispersion, normalized by omega*1m dispvx=sqrt(dispvx/lkm)/omega dispvy=sqrt(dispvy/lkm)/omega dispvz=sqrt(dispvz/lkm)/omega dispvxy=sqrt(dispvxy/lkm)/omega cx_stamp=dispvx cy_stamp=dispvy cz_stamp=dispvz cxy_stamp=dispvxy orb=time/period write(idev,*) orb,coll write(idev,*) dispvx,dispvy,dispvz,dispvxy c-------------------------------------------------- c constants of motion FACTEMP=KAPPA2-4.D0*OMEGA2 SHEAR=-FACTEMP/2.d0/OMEGA DO 11 I=1,LKM UVAK=UVAK+RP(I,1)*MYY(I) WVAK=WVAK+(RP(I,2)+SHEAR*R(I,1))*MYY(I) ZVAK=ZVAK+R(I,3)*MYY(I) ZPVAK=ZPVAK+RP(I,3)*MYY(I) MSUM=MSUM+MYY(I) 11 CONTINUE UVAK=UVAK/MSUM/XL/W0 WVAK=UVAK/MSUM/YL/W0 ZVAK=ZVAK/MSUM ZPVAK=ZPVAK/MSUM/XL/W0 uvak_all=uvak wvak_all=wvak zvak_all=zvak zpvak_all=zpvak RETURN END c---------------------------------------------------------------------- c list of output files opened c------------------------------------------- c ssbox_em.f: open(unit=7 ,file=runid(1:nlen)//'.logi',status='new') c ssbox_em.f: open(unit=17 ,file='TIMESTAMP',status='unknown') c ssbox_em.f: open(unit=20,file=runid(1:nlen)//'.cidl',status='new') c ssbox_em.f: open(unit=25,file=runid(1:nlen)//'.pidl',status='new') c ssbox_em.f: open(unit=30,file=runid(1:nlen)//'.dgidl',status='new', c + form='unformatted') c ssbox_em.f: open(unit=31,file=runid(1:nlen)//'.dggidl',status='new', c + form='unformatted') c ssbox_em.f: open(unit=71,file=runid(1:nlen)//'.ufidl',status='new', c + form='unformatted') c ssbox_em.f: open(unit=90,file=runid(1:nlen)//'.aniidl', c + form='unformatted') c these are closed right after use: c ssbox_em.f: open(unit=50,file='inipos',status='old') c ssbox_em.f: open(unit=50,file='finpos',status='unknown')