c MC Calculation of scattering in a beam emitted in a square pulse c from a rectangular heater c 4/12/99 move all floats to single precision c 4/22/99 Vector product formula for rsqm c 5/5/99 Fixed logic error: icmin might not be iearliest(i) when c passed to removecoll c 5/20/99 nmcoll, remove mirrors: mirrors provide no benefit c 5/26/99 Commented out ars[], we no longer use it and it c occasionally gives a divide by zero error. c 15 may 2001-- above comment not true c 5/26/99 Fixed problem writing arrays when inpp not divisible by 5 c 5/26/99 nmcoll2.for : change newcolls to calculate both c particles simultaneously c 5/28/99 nmcoll3.for : change v(i,1)=v1(i)... r0(i,1)=r1(i).... c 12 mar 2001 adding mpi calls to run on beowulf (sidd) c 24 may 2001 calls to ran() changed to rand() c 25 may 2001 REMINDER: i,j,inpp,iNCM etc must be made long integer c 4july2001: protect against simultaneous collisions?single precision ? c 2 sep 2001: attempt to eliminate messaging in newcolls c 21 May 2002: try to make program take any value of inpp c Although nproc must still be odd c 16 July 2002: include zd in messaging at end of step 1 c .. correcting bug that caused nonlocal zds to b zero c 23 Oct 2001:--------------------------------------------- c procedure to run jobs on osc beowulf cluster c we submit jobs to the cluster with a program named qsub c thus c qsub mp_33d99.job c 'Aha!' u say, 'but what is an mp_39d99.job ?' c -------- begin excerpt from beowulf session ------- c /b/osu2786/mp/mp_39d.work$ cat mp_39d99.job c #PBS -l cput=200:00 c #PBS -l walltime=5:00 c #PBS -l nodes=25:ppn=4 c #PBS -l mem=256mb c #PBS -N mp_39d99 c #PBS -S /bin/bash c cd $HOME/mp/mp_39d.work/ c mpiexec ./mp_39d99 > mp_39d99.out 2>mp_39d99.err c -------- end excerpt from beowulf session ------- c c first off u notice that the first line sez c /b/osu......ork$ at the beginning. All that bit is the c current directory. the next bit, 'cat mp....job' is c the command to print the job file to the monitor c which results in the following lines beginning with c #PBS -l cput ... c thru c mpiexec .... c c the first line of the jobfile specifies the total c CPU time over all processors, the second is the c CPU walltime, should be greater than the first # c divided by the # of cpus .. which should b specified c by the next line as the product of nodes and c processors per node (ppn). I say should, because altho c in this case the product is 100, in reality the numner c of cpus doin work is 99 (hence the 99 in the title). c we ask for 100 cpus bcoz there is no ez way to ask for c 99 (there are 32 nodes with 4 processors per node....) c c the next line specifies a memory limit (per processor c i believe) and the lines fokllowing specify a job name c and a shell, in this case /bin/bash. a shell is what u c type at and gives u a prompt (i typed the cat...job c command into the bash shell earlier) c c the penultimate line instructs the qsub program to c change directory to our working directory c in this case cd $HOME/mp/mp_39d.work/ ( $HOME c is the login directory, in this case /b/osu2786 ) c c the last line does the actual work, now that the c environment has hopefully been set up for it c mpiexec takes the executable file specified c in this case mp_39d99 (the ./ in the beginning is to c specify that mp_39d99 is in the current directory c and is not really neccessary), and redirect the output c to mp_39d99.out and eroors to mp_39d99.err c c so how do u make an executable like mp_39d99 from c the input fortran file: the magic words are c c -------begin excerpt --------------------------- c /b/osu2786/mp/mp_39d.work$ mpif77 -o mp_39d99 mp_39d99.for c /usr/local/pgi-3.2.3/linux86/lib/libpgftnrtl.a(cnfg.o): c In function `__fio_scratch_name': c cnfg.o(.text+0x33): the use of `tempnam' is dangerous, c better use `mkstemp' c /b/osu2786/mp/mp_39d.work$ c -------end excerpt --------------------------- c c we blandly ignore the errors. c c in any event.. this set of comments is actually bakward c coz u would in reality do the mpif77 command b4 the c qsub .. c c and when one does do the qsub, as i did on oct 19 c it does not run for many (4 days till today..) days c one investigates the status of ones jobs with c qstat c thus c ---- begin excerpt -------- c /b/osu2786/mp/mp_39d.work$ qstat | grep osu2786 c 87418.oscbw01 mp_39d99 osu2786 0 Q parallel c /b/osu2786/mp/mp_39d.work$ c ---end excerpt -------------- c c still dolefully queued. if one does not put in the | grep... c at the end one sees all jobs queued; the grep serves to c look for my jobs (i masquerade as osu2786 here) c and a c qstat -f c gives all gory details c c to look at remaining research units the command is c OSCusage c --- begin excerpt ------ c /b/osu2786/mp/mp_39d.work$ OSCusage c c Usage Statistics for project PAS0786 c for 10/24/2001 to 10/24/2001 c c RU Balance: 56.10416 c c Username Dates RUs Status c c osu2785 10/24/2001 to 10/24/2001 0.00000 ACTIVE c osu2786 10/24/2001 to 10/24/2001 0.00000 ACTIVE c osu1264 10/24/2001 to 10/24/2001 0.00000 RESTRICT c osu1263 10/24/2001 to 10/24/2001 0.00000 ACTIVE c ============= PAS0786 TOTAL 0.00000 c /b/osu2786/mp/mp_39d.work$ c ----------------------------end excerpt ------------ implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' integer err,rank,size,nproc,nloc,nloc1,nloc2 integer status(MPI_STATUS_SIZE) character*8 time_day, date_,time_now character*12 filename common/sigmao/rsg08p,cv2,cv2er0,a1,radsq common/ffo/gam,m4etc common/calph/alph common/cevap/Nevap,tau,aH,wH,lH,Rad0,Rad0sq,vbar,nu common/constants/NA,kb,hbar,pi,m4 common/pta/bm1,b0,b1,b2,b3,b4,b5,b6 c zd:=array(1..npp): # height of impact on dome c v:= array(1..npp,1..4): # particle velocity; v[i,4]:= vperp = vp c deadcoll:= array(1..ncm): # array of used collision numbers available c r0:= array(1..npp,1..3): # particle position at t=0 c tbirth:= array(1..npp): # birth times c iearliest:= array(1..npp): c number of earliest collision in particle list c vav:= array(1..6,[0 $ i=1..6] ): # 1..3=vav, 4..6=vavsq c inpp=83160=2*3^3*5*7*11 c inpp=332640=2^3*3^3*5*7*11 PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/ms/ierr common/cdatai/icollp(iNCM),jcollp(iNCM),iprev(iNCM), & inext(iNCM),jprev(iNCM),jnext(iNCM),iearliest(inpp), & nloc1,nloc2,rank,igposs common/cdataf/tcurr,sumsig,maxsig,tcoll(iNCM),vr2coll(iNCM) common/cdead/iposscoll,itotcoll,imcoll,idc0,ideadcoll(iNCM) common/pdataf/tdead(inpp),tdead0(inpp),tbirth(inpp),zd(inpp), & v1(inpp),v2(inpp),v3(inpp),v4(inpp),r1(inpp),r2(inpp),r3(inpp), & ars(inpp,3),vav(6),idist(iNd*2) c v4(i) = vperp common/cinit/imcoll0,iskip,rnear0,rfar0,crnear0,crfar0, & ratnear0,ratfar0,inear0,ifar0,ircol(200),izcol(200),ihard c iabin = number of bins in theta (angle on the sphere) c itbin = Number of bins in time (from 0 to 10 msec) PARAMETER(iabin=10,itbin=100) DIMENSION bin(itbin,iabin),bin2(itbin,iabin),isdist(2*iNd), & ibin(itbin,iabin) DIMENSION edist(iNd*2) c Initialise MPI library IF (nproc.EQ.1) rank=0 CALL MPI_INIT(ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr c Find rank, size and error code CALL MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr CALL MPI_COMM_SIZE(MPI_COMM_WORLD,size,ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr c Dismiss unnecessary processors IF (rank.GE.nproc) GOTO 999 nloc1=rank*nloc+1 nloc2=(rank+1)*nloc c last cpu has smaller number of local particles: IF (rank.EQ.(nproc-1)) nloc2=inpp ihard=0 iskip=0 ifile=1 filename='m332k1e3.x1' pi=4.0*atan(1.0) data hbar/1.0546D-27/,NA/6.022D+23/,kb/1.3806D-16/ m4=4.0028D0/NA call time(time_now) WRITE(6,*)'mp_n started on ',rank,' at: ',time_now c a0 = Scattering length (SAPT2) c r0 = Effective range. r02 = r0 / 2; c a0:=85.25: er0:=7.256: # in angstrom for SAPT2 (Janzen & Aziz) data a0/85.25/,er0/7.256D0/,ipot/0/ c data a0/-176.32592D0/,er0/7.9567759D0/,ipot/3/ sg0 = 8.0*pi*a0**2 c # Collision distance squared radsq = rad**2 lambda = radsq*pi/sg0*1D16 c converts vr from cm/sec to inverse angstrom cv = 1D-8*m4/2.0/hbar rsg08p = radsq/sg0*8.0*pi cv2 = cv**2 cv2er0 = cv2*er0/2.0 c a1 = 1/a needed in sigma function a1 = 1.0d0/a0 c tau:= 15D-6: # halflength of pulse in seconds tau = 15D-6 c tau = .25E-6 c halfwidth and halflength of heater in cm wH = .4D0 lH = .4D0 aH = wH*lH*4.0 c Radius of the Dome Rad0 = 3.5D0 c Rad0 = 8.2 Rad0sq = Rad0**2 nu = 1.0 tmax=5D-3*Rad0/5.0D0 c tmax=2D-3*Rad0/5.0D0 c Angstroms of film evaporated is Devap vv4 = 27.58D0/NA Nevap = inpp*lambda Devap=Nevap*vv4/aH*1D8 c R.L. Rusby JLTP 58, 203 (1985) Parameters for vapor pressure DATA bm1/-7.41816D0/, b0/5.42128D0/, b1/9.903203D0/ DATA b2/-9.617095D0/, b3/6.804602D0/, b4/-3.0154606D0/ DATA b5/.7461357D0/, b6/-.0791791D0/ c Bisection to find the beam evaporation temperature T1 = .25 T2 = 2.0 Tmm = (T1+T2)/2.0 fT1 = P_Nevap( T1 ) - PT( T1) fT2 = P_Nevap( T2 ) - PT( T2) DO 101 i=1 ,30 fTm = P_Nevap( Tmm ) - PT( Tmm) IF ((fT1*fTm) .LT. 0.0 ) THEN T2 = Tmm fT2 = fTm ELSE T1 = Tmm fT1 = fTm ENDIF Tmm = (T1+T2)/2.0 101 CONTINUE T = Tmm alph = m4/2/kB/T vbar = sqrt(8.0*kb*T/pi/m4) c Nominal heat into heater in ergs,calculated from Nevap L4 = 7.17D0 Q0 = Nevap*kb*(L4+2.0*T) WRITE(6,*)'T,Q0=',T,Q0 WRITE(6,*)'characteristic length in cm, tau*vbar=', tau*vbar WRITE(6,*)'ncm, max # of collisions per particle=',iNCM,iNCM/nloc c # of dead collision numbers in deadcoll idc0 = 0 c # = total no of collisions considered iposscoll = 0 c # = total no of actual collisions itotcoll = 0 c # number of last collision in list imcoll = 0 c # sum of crossections for collisions sumsig = 0 c # max of crossections for collisions maxsig = 0 imcoll0 = 0 inear0 = 0 ifar0 = 0 rfar0 = 0.0 rnear0 = 0.0 crfar0 = 0.0 crnear0 = 0.0 ratnear0 = 0.0 ratfar0 = 0.0 DO 139 i=1,200 ircol(i)=0 izcol(i)=0 139 CONTINUE call time(time_day) call date(date_) DO 108 i=1,itbin DO 109 j=1,iabin bin(i,j)=0.0 bin2(i,j)=0.0 ibin(i,j)=0 109 END DO 108 END DO DO 119 i=1,2*iNd edist(i) = 0.0 isdist(i) = 0 119 END DO isummcoll0=0 isumtotcol=0 iloops=1 DO 118 ii=1,iloops WRITE(6,*)'rank,iter loop ',rank,ii,'/',iloops dummy = Func(ii) IF (itotcoll .GT. 0) avrsq = sumsig/itotcoll WRITE(6,*)'end of iter loop' isummcoll0=isummcoll0+imcoll0 isumtotcol=isumtotcol+itotcoll WRITE(6,*)'Binning results rank:',rank DO 117 i=1,inpp cs = zd(i)/Rad0 ia = nint(iabin*zd(i)/Rad0+.5) IF (ia .GT. iabin) THEN WRITE(6,*)'rank,ERROR(main): ia > iabin : &ia,iabin,zd(i)',rank,ia,iabin,zd(i) GOTO 117 END IF it = int(itbin*tdead(i)/tmax+.5) cs = zd(i)/Rad0 id = nint(iNd*cs+.5) +iNd IF ((id.EQ.101).AND.(rank.EQ.0)) THEN WRITE(6,*)'zd=',zd(i),' vz=',v3(i),' tdead=', &tdead(i),' tbirth=',tbirth(i) END IF E = kb*L4+.5*m4*(v1(i)**2 + v2(i)**2 + v3(i)**2) El = E*lambda/iloops edist(id) = edist(id) + El isdist(id) = isdist(id) + 1 IF (ia .LT. 1) GOTO 117 c Particle hits the floor IF ((it .LT. 1) .OR. (it .GT. itbin) ) GOTO 117 c particle hits outside of the recording time bin(it,ia) = bin(it,ia) + El bin2(it,ia) = bin2(it,ia) + El**2 ibin(it,ia) = ibin(it,ia) + 1 117 END DO c DO 120 i=1,2*iNd c isdist(i)=isdist(i)+idist(i) c 120 END DO 118 END DO IF ((rank.EQ.0).AND.(ifile.EQ.1)) THEN open(file=filename,type='new',unit=3) c 'new'gives means: if outbeam;n exists, outbeam;n+1 is written WRITE(3,*)'# ',filename,' data run at ',time_day, date_ WRITE(3,*)'# nproc,ihard,filename',nproc,ihard,filename inp = inpp IF (inp .GT. 2000) THEN inp = 2000 ENDIF WRITE(3,901)inpp,inp,iNCM,iNd 901 FORMAT(1x,'npp:=',I6,': nnp:=',I6,': NCM:=',I8,': Nd:=',I4,':') WRITE(3,905)nu,wH,lH,Rad0sq,Rad0 905 FORMAT(1x,'nu:=',F14.10,': wH:=',F14.10,': lH:=',F14.10, & ': R0sq:=',F14.8,': R0:=',F14.8,':') WRITE(3,906)rad,radsq,tau,Q0 906 FORMAT(1x,'rad:=',F14.10,': radsq:=',F14.10,': tau:=', & F14.10,': Q0:=',F14.10,':') WRITE(3,907)Devap,Nevap,T,Q0 907 FORMAT(1x,'Devap:=',F12.6,': Nevap:=',E14.8,': T:=',F14.10, & ': Q0:=',F14.10,':') WRITE(3,908)iposscoll,itotcoll,imcoll 908 FORMAT(1x,'posscoll:=',I12,': totcoll:=',I10,': mcoll:=',I10, & ':') IF (itotcoll .GT. 0) avrsq = sumsig/itotcoll WRITE(3,909)tcurr,avrsq,maxsig 909 FORMAT(1x,'tcurr:=',F14.10,': avrsq:=',F16.10,': maxrsq:=', & F16.10,':') cpp = 2.*itotcoll/inpp WRITE(3,911)itotcoll,imcoll0,cpp 911 FORMAT(1x,'totcoll:=',I12,': mcoll0:=',I12,': cpp:=',F16.10,':') WRITE(3,*)'summcoll0:=',isummcoll0,': sumtotcoll:=', & isumtotcol,':' WRITE(3,922)sqrt(rnear0/inear0),sqrt(crnear0/inear0), & sqrt(rfar0/ifar0),sqrt(crfar0/ifar0) 922 FORMAT(1x,'rnear0:=',F14.10,': crnear0:=',F14.10,': rfar0:=', & F14.10,': crfar0:=',F14.10,':') WRITE(3,923)ratnear0/inear0,ratfar0/ifar0,inear0,ifar0 923 FORMAT(1x,'ratnear0:=',F14.10,': ratfar0:=',F14.10, & ': inear0:=',I12,': ifar0:=',I12,':') WRITE(3,910)iabin,itbin,tmax 910 FORMAT(1x,'abin:=',I6,': tbin:=',I6,': tmax:=',F14.10,':') WRITE(3,*)'iloops:=',iloops,': ipot:=', & ipot,':' WRITE(3,*)'ars:=array(1..nnp, 1..3, [' DO 102 i=1,inp-1 WRITE(3,916)ars(i,1),ars(i,2),ars(i,3) 102 END DO WRITE(3,917)ars(inp,1),ars(inp,2),ars(inp,3) WRITE(3,*)'v:=array(1..nnp, 1..4, [' DO 103 i=1,inp-1 WRITE(3,912)v1(i),v2(i),v3(i),v4(i) 103 END DO WRITE(3,913)v1(inp),v2(inp),v3(inp),v4(inp) WRITE(3,*)'zd:=array(1..nnp, [' DO 104 i=1,inp-5,5 WRITE(3,924)zd(i),zd(i+1),zd(i+2),zd(i+3),zd(i+4) 104 END DO IF (MOD(inp,5) .EQ. 0 ) THEN WRITE(3,925)zd(inp-4),zd(inp-3),zd(inp-2),zd(inp-1),zd(inp) ELSE imin=inp-MOD(inp,5)+1 IF (imin .LT. inp) THEN DO 121 i=imin,inp-1 WRITE(3,*)zd(i),',' 121 END DO END IF WRITE(3,*)zd(inp),' ] ):' END IF WRITE(3,*)'tdead1:=array(1..nnp, [' DO 105 i=1,inp-5,5 WRITE(3,914)tdead(i),tdead(i+1),tdead(i+2),tdead(i+3),tdead(i+4) 105 END DO IF (MOD(inp,5) .EQ. 0 ) THEN WRITE(3,915)tdead(inp-4),tdead(inp-3),tdead(inp-2), & tdead(inp-1),tdead(inp) ELSE imin=inp-MOD(inp,5)+1 IF (imin .LT. inp) THEN DO 122 i=imin,inp-1 WRITE(3,*)tdead(i),',' 122 END DO END IF WRITE(3,*)tdead(inp),' ] ):' END IF WRITE(3,*)'tdead0:=array(1..nnp, [' DO 106 i=1,inp-5,5 WRITE(3,914)tdead0(i),tdead0(i+1),tdead0(i+2),tdead0(i+3), & tdead0(i+4) 106 END DO IF (MOD(inp,5) .EQ. 0 ) THEN WRITE(3,915)tdead0(inp-4),tdead0(inp-3),tdead0(inp-2), & tdead0(inp-1),tdead0(inp) ELSE imin=inp-MOD(inp,5)+1 IF (imin .LT. inp) THEN DO 123 i=imin,inp-1 WRITE(3,*)tdead0(i),',' 123 END DO END IF WRITE(3,*)tdead0(inp),' ] ):' END IF WRITE(3,*)'vav:=array(1..6):' WRITE(3,902)vav(1),vav(2) WRITE(3,903)vav(3),vav(4) WRITE(3,904)vav(5),vav(6) 902 FORMAT(1x,'vav[1]:=',F16.7,': vav[2]:=',F16.7,':') 903 FORMAT(1x,'vav[3]:=',F16.7,': vav[4]:=',F16.7,':') 904 FORMAT(1x,'vav[5]:=',F16.7,': vav[6]:=',F16.7,':') WRITE(3,*)'dist:=array(-(Nd-1)..Nd, [' DO 107 i=1,iNd*2-5,5 WRITE(3,918)isdist(i),isdist(i+1),isdist(i+2),isdist(i+3), & isdist(i+4) 107 END DO WRITE(3,919)isdist(iNd*2-4),isdist(iNd*2-3),isdist(iNd*2-2), & isdist(iNd*2-1),isdist(iNd*2) WRITE(3,*)'Edist:=array(-(Nd-1)..Nd, [' DO 116 i=1,iNd*2-5,5 WRITE(3,932)edist(i),edist(i+1),edist(i+2),edist(i+3), & edist(i+4) 116 END DO WRITE(3,933)edist(iNd*2-4),edist(iNd*2-3),edist(iNd*2-2), & edist(iNd*2-1),edist(iNd*2) it4 = itbin-4 WRITE(3,*)'bin:=array(1..abin, 1..tbin, [ [' DO 110 j=1,iabin-1 DO 111 i=1,it4,5 IF (i .LT. it4) THEN WRITE(3,920)bin(i,j),bin(i+1,j),bin(i+2,j),bin(i+3,j), & bin(i+4,j) ELSE WRITE(3,921)bin(i,j),bin(i+1,j),bin(i+2,j),bin(i+3,j), & bin(i+4,j) END IF 111 END DO WRITE(3,*)'], [' 110 END DO DO 112 i=1,it4,5 IF (i .LT. it4) THEN WRITE(3,920)bin(i,iabin),bin(i+1,iabin),bin(i+2,iabin), & bin(i+3,iabin),bin(i+4,iabin) ELSE WRITE(3,921)bin(i,iabin),bin(i+1,iabin),bin(i+2,iabin), & bin(i+3,iabin),bin(i+4,iabin) END IF 112 END DO WRITE(3,*)'] ] ):' WRITE(3,*)'ibin:=array(1..abin, 1..tbin, [ [' DO 113 j=1,iabin-1 DO 114 i=1,it4,5 IF (i .LT. it4) THEN WRITE(3,930)ibin(i,j),ibin(i+1,j),ibin(i+2,j),ibin(i+3,j), & ibin(i+4,j) ELSE WRITE(3,931)ibin(i,j),ibin(i+1,j),ibin(i+2,j),ibin(i+3,j), & ibin(i+4,j) END IF 114 END DO WRITE(3,*)'], [' 113 END DO DO 115 i=1,it4,5 IF (i .LT. it4) THEN WRITE(3,930)ibin(i,iabin),ibin(i+1,iabin),ibin(i+2,iabin), & ibin(i+3,iabin),ibin(i+4,iabin) ELSE WRITE(3,931)ibin(i,iabin),ibin(i+1,iabin),ibin(i+2,iabin), & ibin(i+3,iabin),ibin(i+4,iabin) END IF 115 END DO WRITE(3,*)'] ] ):' WRITE(3,*)'rcol:=array(1..200, [' DO 138 i=1,191,5 WRITE(3,918)ircol(i),ircol(i+1),ircol(i+2),ircol(i+3), & ircol(i+4) 138 END DO WRITE(3,919)ircol(196),ircol(197),ircol(198), & ircol(199),ircol(200) WRITE(3,*)'zcol:=array(1..200, [' DO 137 i=1,191,5 WRITE(3,918)izcol(i),izcol(i+1),izcol(i+2),izcol(i+3), & izcol(i+4) 137 END DO WRITE(3,919)izcol(196),izcol(197),izcol(198), & izcol(199),izcol(200) 912 FORMAT(1x,'[',F14.5,',',F14.5,',',F14.5,',',F14.5,'],') 913 FORMAT(1x,'[',F14.5,',',F14.5,',',F14.5,',',F14.5,'] ] ):') 914 FORMAT(1x,F14.10,',',F14.10,',',F14.10,',',F14.10,',',F14.10,',') 924 FORMAT(1x,F14.6,',',F14.6,',',F14.6,',',F14.6,',',F14.6,',') 915 FORMAT(1x,F14.10,',',F14.10,',',F14.10,',',F14.10,',', & F14.10,' ] ):') 925 FORMAT(1x,F14.6,',',F14.6,',',F14.6,',',F14.6,',', & F14.6,' ] ):') 916 FORMAT(1x,'[',F14.10,',',F14.10,',',F14.10,'],') 917 FORMAT(1x,'[',F14.10,',',F14.10,',',F14.10,'] ] ):') 918 FORMAT(1x,I8,',',I8,',',I8,',',I8,',',I8,',') 919 FORMAT(1x,I8,',',I8,',',I8,',',I8,',',I8,' ] ):') 920 FORMAT(1x,E15.7,',',E15.7,',',E15.7,',',E15.7,',',E15.7,',') 921 FORMAT(1x,E15.7,',',E15.7,',',E15.7,',',E15.7,',',E15.7) 930 FORMAT(1x,I12,',',I12,',',I12,',',I12,',',I12,',') 931 FORMAT(1x,I12,',',I12,',',I12,',',I12,',',I12) 932 FORMAT(1x,E15.7,',',E15.7,',',E15.7,',',E15.7,',',E15.7,',') 933 FORMAT(1x,E15.7,',',E15.7,',',E15.7,',',E15.7,',',E15.7,' ] ):') END IF 999 call time(time_now) CALL MPI_FINALIZE(ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr WRITE(6,*)'mp_n finished on ',rank,' at: ',time_now END c a0 = Scattering length, er0 = effective range. c sig = cross-section in Angstrom^2; k in reciprocal angstrom: c k:='k': sig:= (i,k)-> 8*pi/(k^2 + ( er0[i]/2*k^2 - 1/A[i])^2 ): c A[1]:=87.5: er0[1]:=8.65339: # in angstrom for HFD-B2(HE) (Meyer thesis) c From Gutierrez et al PRB 29, 5211 (1984): c A[2]:=125.0810: er0[2]:=7.397569: # in angstrom for HFDHE2 c A[3]:=-176.32529: er0[3]:=7.9567759: # for Lennard-Jones (Meyer thesis) FUNCTION rsg0(vrsq) implicit real*8 (a-h,k-z), integer (i,j) common/sigmao/rsg08p,cv2,cv2er0,a1,radsq rsg0 = rsg08p/(cv2*vrsq + (vrsq*cv2er0-a1 )**2) RETURN END c Vapor pressure during pulse FUNCTION P_Nevap(T) implicit real*8 (a-h,k-z), integer (i,j) common/cevap/Nevap,tau,aH,wH,lH,Rad0,Rad0sq,vbar,nu common/constants/NA,kb,hbar,pi,m4 P_Nevap = kb*T*4.0*Nevap/(2.0*tau)/aH/sqrt(8.0*kb*T/pi/m4) RETURN END c Vapor pressure formula FUNCTION PT(T) implicit real*8 (a-h,k-z), integer (i,j) common/pta/bm1,b0,b1,b2,b3,b4,b5,b6 PT = 10.0*exp(bm1/T + b0 + T*(b1 + T*(b2 + T*(b3 +T*(b4 & + T*(b5 + T*b6) ) ) ) )) RETURN END FUNCTION vF(X) C Function MUST remain in DOUBLE PRECISION implicit real*8 (a-h,k-z), integer (i,j) common/calph/alph c Returns a velocity from a Maxwellian beam distribution c argument must be a DOUBLE PRECISION c random number between 0 and 1 c Binary search adapted from Maple file `Maxwell 7 Aug 02` c f(v) = 2*v^3*exp(-v^2): is beam distribution, c with v in units of (sqrt(2kBT/m) c Integral is F(v) = 1-(1+v^2)exp(-v^2) c Suppose F(v)=1-X where 0back0 x = 2.*RAND()*wH-wH y = 2.*RAND()*lH-lH vp = vv*sqrt(1.0-cs**2) c # velocity components v3(i) = vv*cs sendp(3,ib)=v3(i) v1(i) = vp*cos(ph) sendp(1,ib)=v1(i) v2(i) = vp*sin(ph) sendp(2,ib)=v2(i) c # vector r0 = position at t=0 plus vector v define trajectory r1(i) = x-v1(i)*ts sendp(4,ib)=r1(i) r2(i) = y-v2(i)*ts sendp(5,ib)=r2(i) r3(i) = -v3(i)*ts sendp(6,ib)=r3(i) c # tdead = time for particle to hit dome tdead(i) = twall(i) sendp(7,ib)=tdead(i) sendp(8,ib)=tbirth(i) sendp(9,ib)=zd(i) tdead0(i) = tdead(i) 201 CONTINUE CALL time(time_now) WRITE(6,*)'Processor ',rank,' finished local particle & list at ',time_now WRITE(6,*)'rank,particle nloc2 ',rank,r1(nloc2),v1(nloc2), &tbirth(nloc2),tdead(nloc2) DO 12 it=1,nproc-1 irk(it)=MOD(rank+it,nproc) CALL MPI_ISEND(sendp,9*nloc,MPI_DOUBLE_PRECISION &,irk(it),11,MPI_COMM_WORLD,ira(it),ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr 12 CONTINUE CALL time(time_now) WRITE(6,*)'processor ',rank,' broadcast data at:', &time_now DO 11 it=1,nproc-1 CALL MPI_RECV(recp,9*nloc,MPI_DOUBLE_PRECISION &,MPI_ANY_SOURCE,11,MPI_COMM_WORLD,status,ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr irank=status(MPI_SOURCE) iloc1=irank*nloc+1 iloc2=(irank+1)*nloc c Each cpu has nloc local particles except the last:. IF (irank.EQ.(nproc-1)) iloc2=inpp DO 11 i=iloc1,iloc2 ib=i-iloc1+1 v1(i)=recp(1,ib) v2(i)=recp(2,ib) v3(i)=recp(3,ib) r1(i)=recp(4,ib) r2(i)=recp(5,ib) r3(i)=recp(6,ib) tdead0(i)=recp(7,ib) tdead(i)=recp(7,ib) tbirth(i)=recp(8,ib) zd(i)=recp(9,ib) 11 CONTINUE CALL time(time_now) WRITE(6,*)'Finished global particle list: processor ',rank, &'at ',time_now WRITE(6,*)'rank,particle inpp',rank,r1(inpp),v1(inpp), &tbirth(inpp),tdead(inpp) WRITE(6,*)'rank,particle 1',rank,r1(1),v1(1), &tbirth(1),tdead(1) DO 16 it=1,nproc-1 CALL MPI_WAIT(ira(it),status,ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr 16 CONTINUE c STEP 2: Initial collision list:----------------------------------- c icoll is an index in the local collision list c imcoll is the max value of icoll c imcoll0 is imcoll at the end of initial collision list icoll = 0 idc0 = 0 tlocmin = 1.D10 c icmin is earliest local collision icmin = 0 imid=(nproc-1)/2 c nproc must be odd IF (rank.LT.imid) itarg=nloc2-1+imid*nloc IF (rank.EQ.imid) itarg=inpp-1 IF (rank.GT.imid) itarg=inpp-1+(rank-imid)*nloc C May 22 Old line was: itarg=nloc2-1+(inpp-nloc)/2 c itarg is used in DO 211 jdum=i,itarg, j=1+MOD(jdum,inpp)in c STEP 2 so that j runs from i+1 to 1+itarg if itarg iNCM : icoll,iNCM',icoll,iNCM END IF IF ( (r0z+tc*(vz+v3(j))+r3(j) .GT. dchar) .OR. & (tc .GT. t2) ) THEN rfar0 = rfar0 + rsqm/vrsq crfar0 = crfar0 + rsg ratfar0 = ratfar0 + rsqm/rsg/vrsq ifar0 = ifar0 + 1 ELSE rnear0 = rnear0 + rsqm/vrsq crnear0 = crnear0 + rsg ratnear0 = ratnear0 + rsqm/rsg/vrsq inear0 = inear0 + 1 END IF IF (tc .LT. tlocmin) THEN icmin = icoll tlocmin = tc END IF icollp(icoll) = i jcollp(icoll) = j c WRITE(6,*)'setting jcollp: j,jcollp,icoll=', c & j,jcollp(icoll),icoll tcoll(icoll) = tc vr2coll(icoll) = vrsq c iearliest(i,j) the numbers of the earliest local collisions c for local particle i and target j are set in prevnext dummy = prevnext(icoll,i,tc) c # sets prev and next in coll dummy = prevnext(icoll,j,tc) c od: od: # ends of loops in j,i 211 CONTINUE 210 CONTINUE WRITE(6,*)'rnear0,crnear0,rfar0,crfar0=',rnear0,crnear0, & rfar0,crfar0 WRITE(6,*)'inear0,ifar0=',inear0,ifar0 c # max number in list imcoll = icoll imcoll0 = icoll c iposscoll is the # of possible local collisions iposscoll = icoll WRITE(6,*)'rank,initial mcoll,tlocmin(musec),icmin', & rank,imcoll,1.D6*tlocmin,icmin CALL time(time_now) WRITE(6,*)'Finished collision list at', time_now c STEP 2A: c Fortran from maple file `message` of 9 Aug 01` c Calculate STEP 3 messaging instructions jmax=10 c # Max number of possible instructions per cpu c # inst(i,j) = jth instruction for cpu i # NOTE cpus numbered from 1 to npr c # Instruction - i means receive from i; +i means send to i; c $ 0 means do nothing c # inf(i)=1 means cpu i has information DO 57 i=1,nproc 57 inf(i)=1 isent=-1 c isent = -1 means no messages waiting c isent=sp means message waiting from cpu sp DO 58 j=1,jmax jflag=0 c jflag=0 means no info left except in cpu npr c initialize jth column of inst DO 59 i=1,nproc 59 inst(i,j)= 0 DO 61 i=nproc,1,-1 IF (inf(i).EQ.1) THEN c cpu i has information IF (i.NE.1) jflag=1 IF (isent.EQ.-1) THEN c no unreceived message isent=i c cpu i sends info to a later one inf(i)=0 c no info in cpu i after sending message ELSE inst(isent,j)=i c cpu isent sends info to i inst(i,j)=-isent c cpu i receives info from cpu isent isent=-1 ENDIF ENDIF 61 CONTINUE IF (isent.NE.-1) THEN c last message can't be received inf(isent)= 1 c # restore info to cpu isent isent=-1 ENDIF c instmax is the max # of instructions for any processor instmax=j-1 IF (jflag.NE.1) GOTO 65 58 CONTINUE c 58 is end of loop in j 65 CONTINUE DO 62 j=1,instmax 62 instr(j)=inst(rank+1,j) c instructions for cpu rank+1 in hunt for global tmin c instructions are reversed & inverted to broadcast the result DO 63 j=instmax,1,-1 IF (instr(j).NE.0) THEN instm=j c instm is # of instructions for local processor GOTO 64 ENDIF 63 CONTINUE 64 WRITE(6,*) 'rank, instm,instmax=', rank,instm,instmax IF (rank.EQ.0) THEN WRITE(6,*)((inst(i,j),i=1,nproc),j=1,instmax) END IF c STEP 3: Advance time & generate new collisions--------------- iseed=0 c initialize processors to same series of random #s CALL SRAND(iseed) WRITE(6,*)'begin step 3: rank iseed ',rank,iseed c itotcoll is the # of actual collisions itotcoll = 0 c Inirialize tcmin, the globally earliest collision time ... tcmin=1.D9 IF (iskip.EQ.0) THEN DO 220 WHILE (tcmin.LT.1.D10) c Initial values for sende: sende(1)= 1.0d0*icollp(icmin) sende(2)= 1.0d0*jcollp(icmin) sende(3)= tlocmin sende(4)= vr2coll(icmin) sende(5)= 1.D0*iposscoll c Start message instructions DO 67 ji=1,instm ins= instr(ji) IF (ins.GT.0) THEN CALL MPI_SEND(sende,5,MPI_DOUBLE_PRECISION &,ins-1,itotcoll+100,MPI_COMM_WORLD,ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr ELSEIF (ins.LT.0) THEN CALL MPI_RECV(rece,5,MPI_DOUBLE_PRECISION &,-ins-1,itotcoll+100,MPI_COMM_WORLD,status,ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr sende(5)=sende(5)+rece(5) IF (rece(3).LT.sende(3)) THEN DO 68 it=1,4 68 sende(it)=rece(it) END IF END IF c End of loop in ji: results are in sende 67 CONTINUE c # Broadcast result by reversing and inverting instructions DO 69 ji=instm,1,-1 ins= instr(ji) IF (ins.LT.0) THEN CALL MPI_SEND(sende,5,MPI_DOUBLE_PRECISION &,-ins-1,itotcoll+200,MPI_COMM_WORLD,ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr ELSEIF (ins.GT.0) THEN CALL MPI_RECV(sende,5,MPI_DOUBLE_PRECISION &,ins-1,itotcoll+200,MPI_COMM_WORLD,status,ierr) IF (ierr.NE.0) WRITE(6,*)'rank,ierr',rank,ierr END IF c End of inverse loop in ji: 69 CONTINUE c Now interpret result igposs=sende(5) tcmin=sende(3) imin=INT(sende(1)+.5D0) jmin=INT(sende(2)+.5D0) vrsq=sende(4) c no collisions left IF (tcmin.EQ.1.D10) GOTO 220 c # Advancing time: tcurr = tcmin c # total actual collisions itotcoll = itotcoll+1 c # i and j collide: c to compare results with different nproc, i is set less c than j so that reorientation of velocity is the same c for different nproc i = MIN(imin,jmin) j = MAX(imin,jmin) IF ((MOD(itotcoll,1000).EQ.1).OR.(itotcoll.LE.100).AND. &(rank.EQ.0)) THEN WRITE(6,*)'rank=',rank,' totcoll=',itotcoll,' pcoll=' &,' tcurr=',tcurr,' i=',i,' j=',j,' igposs=',igposs &,' iposscoll=',iposscoll END IF vr2 = .5*sqrt(vrsq) cs = rsg0(vrsq) sumsig = sumsig+cs IF (cs .GT. maxsig) maxsig = cs c sites of collisions for particles i,j xi=r1(i)+tcmin*v1(i) yi=r2(i)+tcmin*v2(i) zi=r3(i)+tcmin*v3(i) IF (zi .LT. 0.0) THEN WRITE(6,*)'ERR:rank,i,xi,yi,zi=',rank,i,zi,yi,zi WRITE(6,*)'rank,totcol,tcmin=',rank,itotcoll,tcmin END IF xj=r1(j)+tcmin*v1(j) yj=r2(j)+tcmin*v2(j) zj=r3(j)+tcmin*v3(j) IF (zj .LT. 0.0) THEN WRITE(6,*)'ERR:rank,i,xi,yi,zi=',rank,i,zi,yi,zi WRITE(6,*)'rank,totcol,tcmin=',rank,itotcoll,tcmin END IF IF (ihard.EQ.0) THEN c # CM velocity vcmx = .5*(v1(i)+v1(j)) vcmy = .5*(v2(i)+v2(j)) vcmz = .5*(v3(i)+v3(j)) c # random direction cs = RAND()*2.0-1.0 ph = pi*RAND()*2.0 - pi c # New trajectory for i: vrz = vr2*cs vp = vr2*sqrt(1.0-cs**2) vpc=vp*cos(ph) vps=vp*sin(ph) v1(i) = vcmx+vpc v2(i) = vcmy+vps v3(i) = vcmz+vrz v1(j) = vcmx-vpc v2(j) = vcmy-vps v3(j) = vcmz-vrz ELSE c 'hard sphere' calculation begins here c distance of closest approach rsq=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 c time of true hard sphere collision tsp=tcmin-sqrt((cs-rsq)/vrsq) vp=sqrt(vrsq*(cs-rsq))/cs vx=vp*(r1(i)-r1(j)+tsp*(v1(i)-v1(j))) vy=vp*(r2(i)-r2(j)+tsp*(v2(i)-v2(j))) vz=vp*(r3(i)-r3(j)+tsp*(v3(i)-v3(j))) v1(i)=v1(i)+vx v2(i)=v2(i)+vy v3(i)=v3(i)+vz v1(j)=v1(j)-vx v2(j)=v2(j)-vy v3(j)=v3(j)-vz END IF r1(i) = xi-v1(i)*tcmin r2(i) = yi-v2(i)*tcmin r3(i) = zi-v3(i)*tcmin r1(j) = xj-v1(j)*tcmin r2(j) = yj-v2(j)*tcmin r3(j) = zj-v3(j)*tcmin c CM site of the collision xcm=0.5*(xi+xj) ycm=0.5*(yi+yj) zcm=0.5*(zi+zj) c Save collision distance ir=int(sqrt(xcm**2 + ycm**2 + zcm**2)*200/Rad0) +1 iz=int(zcm*200/Rad0) +1 ircol(ir)=ircol(ir)+1 izcol(iz)=izcol(iz)+1 c # Revise deadtimes: tdead(i) = twall(i) tdead(j) = twall(j) c # Revise "birth" times tbirth(i) = tcmin tbirth(j) = tcmin c # Remove stale collisions with particle i: iei=iearliest(i) IF (iei .GT. 0) dummy = removecoll(iei,i) c # Remove stale collisions with particle j: iej=iearliest(j) IF (iej .GT. 0) dummy = removecoll(iej,j) c IF (rank.EQ.0) WRITE(6,*)'tcurr,i,j',tcurr,i,j c WRITE(6,*)'rank=',rank,' totcoll=',itotcoll,' pcoll=' c &,' tcurr=',tcurr,' i=',i,' j=',j,' igposs=',igposs c &,' iposscoll=',iposscoll c Find new colls for i and j, revise collision lists dummy = newcolls(i,j) c # Find new local earliest collision & tcmin: tlocmin =1.D10 DO 221 i=nloc1,nloc2 ie = iearliest(i) IF (ie .GT. 0) THEN c # particle i has a collision tc = tcoll(ie) IF (tc .LT. tlocmin) THEN tlocmin = tc icmin =ie END IF END IF 221 CONTINUE c If no particle has a collision left in earliest list,tlocmin=1.D11 c od: # end of while tcmin<1.D10 220 END DO END IF c ------------------END OF STEP 3--------------------------------- WRITE(6,*)'rank,ncm,mcoll,totcoll,pcoll,posscoll=',rank, &incm,imcoll,itotcoll,iposscoll WRITE(6,*)'rank,end of main procedure, tcurr,deadcoll=', &rank,1.E6*tcurr,idc0 WRITE(6,*)'rank, Find means and distributions',rank c for i to npp do cs = (zd(i)/R0): id = ceil(Nd*cs): DO 225 i=1,iNd*2 idist(i) = 0 225 END DO DO 230 i=1,inpp cs = zd(i)/Rad0 c id = int(iNd*cs) +1+iNd id = nint(iNd*cs+.5) +iNd IF ((id .LT. 1) .OR. (id .GT. 2*iNd)) THEN WRITE(6,*)'ERROR:rank, i,cs,id=',rank,i,cs,id ELSE idist(id) = idist(id)+1 END IF c # Distribution of apparent source, vperp : IF ( abs(v3(i)) .GT. 1E-6) THEN ts = -r3(i)/v3(i) ars(i,1) = abs(r1(i)+v1(i)*ts) ars(i,2) = abs(r2(i)+v2(i)*ts) ars(i,3) = ts ELSE ars(i,1) = 1D6 ars(i,2) = 1D6 ars(i,3) = -1D10 END IF c # vperp v4(i) = sqrt(v1(i)**2+v2(i)**2) 230 CONTINUE WRITE(6,*)'Mean velocities rank:',rank c # WARNING: These averages are weighted by speed: DO 235 i=1,6 vav(i) = 0 235 END DO DO 241 i=1,inpp vav(1) = vav(1)+v1(i) vav(4) = vav(4)+v1(i)**2 vav(2) = vav(2)+v2(i) vav(5) = vav(5)+v2(i)**2 vav(3) = vav(3)+v3(i) vav(6) = vav(6)+v3(i)**2 241 CONTINUE vav(1) = vav(1)/inpp vav(2) = vav(2)/inpp vav(3) = vav(3)/inpp vav(4) = vav(4)/inpp vav(5) = vav(5)/inpp vav(6) = vav(6)/inpp Func = 0.0 RETURN END c end: # end of main procedure c # Procedure for time to hit spherical dome: FUNCTION twall(i) implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/cevap/Nevap,tau,aH,wH,lH,Rad0,Rad0sq,vbar,nu common/pdataf/tdead(inpp),tdead0(inpp),tbirth(inpp),zd(inpp), & v1(inpp),v2(inpp),v3(inpp),v4(inpp),r1(inpp),r2(inpp),r3(inpp), & ars(inpp,3),vav(6),idist(iNd*2) c #scalar products: r0sq = r1(i)**2+r2(i)**2+r3(i)**2 vsq = v1(i)**2+v2(i)**2+v3(i)**2 r0v = r1(i)*v1(i)+r2(i)*v2(i)+r3(i)*v3(i) c # time of impact at sphere t = (-r0v+sqrt(r0v**2+vsq*(Rad0sq-r0sq) ) )/vsq c # height of impact z = r3(i)+v3(i)*t zd(i) = z c # For a hemisphere: IF (z .LT. 0.0) THEN t = -r3(i)/v3(i) c # leave zd as is END IF twall=t RETURN END c # end of twall c PREVNEXT:---------------------------------------------------- c procedure to get previous and next collisions FUNCTION prevnext(ic,i,tc) implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' integer err,rank,size,nproc,nloc,nloc1,nloc2 PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/ms/ierr common/cdatai/icollp(iNCM),jcollp(iNCM),iprev(iNCM), & inext(iNCM),jprev(iNCM),jnext(iNCM),iearliest(inpp), & nloc1,nloc2,rank,igposs common/cdataf/tcurr,sumsig,maxsig,tcoll(iNCM),vr2coll(iNCM) IF ( (i .NE. icollp(ic)) .AND. (i .NE. jcollp(ic)) ) THEN WRITE(6,*)'ERROR(prevnext): i,ic,icollp(ic),jcollp(ic)=', & i,ic,icollp(ic),jcollp(ic) END IF ie=iearliest(i) IF (ie .EQ. 0) THEN c # i did not have a collision before iearliest(i) = ic iprv = 0 inxt = 0 ELSE IF (tcoll(ie) .GT. tc) THEN c # new coll is earlier than current earliest c # insert new coll in sequence: iprv = 0 inxt = ie iearliest(i) = ic c # sets prev of collision ie to be the new one, ic dummy = setprev(ie,i,ic) ELSE c # new coll is later than current first collision for particle i c # set initial values for loop: iprv = ie inxt = inextcoll(ie,i) DO 300 WHILE ((inxt .GT. 0) .AND. (tcoll(inxt) .LT. tc)) iprv = inxt inxt = inextcoll(inxt,i) 300 END DO IF ( (ic .EQ. inxt) .OR. (ic .EQ. iprv) ) THEN WRITE(6,*)'ERROR(prevnext): ic,iprv,inxt,i,ie=', & ic,iprv,inxt,i,ie END IF c # insert new collision in sequence of lists: IF (inxt .GT. 0) THEN dummy = setprev(inxt,i,ic) END IF dummy = setnext(iprv,i,ic) END IF END IF IF ( (ic .EQ. inxt) .OR. (ic .EQ. iprv) ) THEN WRITE(6,*)'ERROR(prevnext):ic,iprv,inxt,i=',ic,iprv,inxt,i END IF IF (icollp(ic) .EQ. i) THEN iprev(ic) = iprv inext(ic) = inxt ELSE jprev(ic) = iprv jnext(ic) = inxt END IF prevnext=0.0 RETURN END c # end of prevnext c # set prev(i) in coll ic to be val FUNCTION setprev(ic,i,ival) implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' integer err,rank,size,nproc,nloc,nloc1,nloc2 PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/ms/ierr common/cdatai/icollp(iNCM),jcollp(iNCM),iprev(iNCM), & inext(iNCM),jprev(iNCM),jnext(iNCM),iearliest(inpp), & nloc1,nloc2,rank,igposs IF (i .LT. 1) THEN WRITE(6,*)'ERROR in setprev, i,ic=', i,ic setprev = 1.0 RETURN END IF IF (icollp(ic) .EQ. i) THEN iprev(ic) = ival ELSE IF (jcollp(ic) .NE. i) THEN WRITE(6,*)'ERROR(setprev): i,icollp,jcollp,ic=', & i,icollp(ic),jcollp(ic),ic END IF jprev(ic) = ival END IF setprev = 0.0 RETURN END c # end of setprev FUNCTION setnext(ic,i,ival) implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' integer err,rank,size,nproc,nloc,nloc1,nloc2 PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/ms/ierr common/cdatai/icollp(iNCM),jcollp(iNCM),iprev(iNCM), & inext(iNCM),jprev(iNCM),jnext(iNCM),iearliest(inpp), & nloc1,nloc2,rank,igposs setnext = 0.0 IF (i .LT. 1) THEN WRITE(6,*)'ERROR in setnext, i,ic=', i,ic RETURN END IF IF (icollp(ic) .EQ. i) THEN inext(ic) = ival ELSE IF (jcollp(ic) .NE. i) THEN WRITE(6,*)'ERROR(setnext): i,icollp,jcollp,ic=', & i,icollp(ic),jcollp(ic),ic END IF jnext(ic) = ival END IF RETURN END c # end of setnext FUNCTION iprevcoll(ic,i) implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' integer err,rank,size,nproc,nloc,nloc1,nloc2 PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/ms/ierr common/cdatai/icollp(iNCM),jcollp(iNCM),iprev(iNCM), & inext(iNCM),jprev(iNCM),jnext(iNCM),iearliest(inpp), & nloc1,nloc2,rank,igposs IF (i .LT. 1) THEN WRITE(6,*)'ERROR in iprevcoll, i,ic=', i,ic iprevcoll=0 RETURN END IF IF (icollp(ic) .EQ. i) THEN iprevcoll=iprev(ic) ELSE IF (jcollp(ic) .NE. i) THEN WRITE(6,*)'ERROR(iprevcoll): i,icollp,jcollp,ic=', & i,icollp(ic),jcollp(ic),ic END IF iprevcoll=jprev(ic) END IF IF ( (iprevcoll .NE. 0 ) .AND. & (i .NE. icollp(iprevcoll)) .AND. & (i .NE. jcollp(iprevcoll)) ) THEN WRITE(6,*)'ERROR(iprevcoll): prev. coll. does not contain i!!' WRITE(6,*)':: i,ic,icollp(iprv),jcollp(iprv),iprv=', & i,ic,icollp(iprevcoll),jcollp(iprevcoll),iprevcoll END IF RETURN END c end:# end of prevcoll FUNCTION inextcoll(ic,i) implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' integer err,rank,size,nproc,nloc,nloc1,nloc2 PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/ms/ierr common/cdatai/icollp(iNCM),jcollp(iNCM),iprev(iNCM), & inext(iNCM),jprev(iNCM),jnext(iNCM),iearliest(inpp), & nloc1,nloc2,rank,igposs IF (i .LT. 1) THEN WRITE(6,*)'ERROR in inextcoll, i,ic=', i,ic inextcoll=0 RETURN END IF IF (icollp(ic) .EQ. i) THEN inextcoll=inext(ic) ELSE IF (jcollp(ic) .NE. i) THEN WRITE(6,*)'ERROR(inextcoll): i,icollp,jcollp,ic=', & i,icollp(ic),jcollp(ic),ic END IF inextcoll=jnext(ic) END IF IF ( (inextcoll .NE. 0 ) .AND. & (i .NE. icollp(inextcoll)) .AND. & (i .NE. jcollp(inextcoll)) ) THEN WRITE(6,*)'ERROR(inextcoll): next coll. does not contain i!!' WRITE(6,*)':: i,ic,icollp(inxt),jcollp(inxt),inxt=', & i,ic,icollp(inextcoll),jcollp(inextcoll),inextcoll END IF RETURN END c # end of next coll c REMOVECOLL -------------------------------------------------------- c # After collision ic, involving particle i, later collisions with c # i are removed from list ( if i>0) FUNCTION removecoll(ic,i) implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' integer err,rank,size,nproc,nloc,nloc1,nloc2 PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/ms/ierr common/cdatai/icollp(iNCM),jcollp(iNCM),iprev(iNCM), & inext(iNCM),jprev(iNCM),jnext(iNCM),iearliest(inpp), & nloc1,nloc2,rank,igposs common/cdead/iposscoll,itotcoll,imcoll,idc0,ideadcoll(iNCM) removecoll = 0.0 IF ( (i .NE. icollp(ic)) .AND. (i .NE. jcollp(ic)) ) THEN WRITE(6,*)'ERROR in removecoll, i,ic,icollp,jcollp= & rank,itotcoll',i,ic,icollp(ic),jcollp(ic),rank,itotcoll END IF DO 400 WHILE (ic .GT. 0) c # ic is now dead idc0 = idc0 + 1 ideadcoll(idc0) = ic j = icollp(ic) IF (j .EQ. i) j = jcollp(ic) iprv = iprevcoll(ic,j) inxt = inextcoll(ic,j) IF (iprv .GT. 0) THEN dummy = setnext(iprv,j,inxt) ELSE iearliest(j) = inxt END IF IF (inxt .GT. 0) THEN dummy = setprev(inxt,j,iprv) END IF ic = inextcoll(ic,i) 400 END DO iearliest(i) = 0 RETURN END c # end of removecoll c NEWCOLLS:----------------------------------------------------- c Procedure to revise collision list for particle ii excluding its c previous collidee jj and vice versa c 2 sep 2001 -- attempt to remove messaging FUNCTION newcolls(ii,jj) implicit real*8 (a-h,k-z), integer (i,j) include 'mpif.h' integer status(MPI_STATUS_SIZE),request integer err,rank,size,nproc,nloc,nloc1,nloc2 PARAMETER(rad=.1d0,nproc=25,inpp=332640,nloc=inpp/nproc+1, &iNCM=(rad*inpp)**2/40/nproc,iNd=100) common/ms/ierr common/cdatai/icollp(iNCM),jcollp(iNCM),iprev(iNCM), & inext(iNCM),jprev(iNCM),jnext(iNCM),iearliest(inpp), & nloc1,nloc2,rank,igposs common/cdataf/tcurr,sumsig,maxsig,tcoll(iNCM),vr2coll(iNCM) common/cdead/iposscoll,itotcoll,imcoll,idc0,ideadcoll(iNCM) common/pdataf/tdead(inpp),tdead0(inpp),tbirth(inpp),zd(inpp), & v1(inpp),v2(inpp),v3(inpp),v4(inpp),r1(inpp),r2(inpp),r3(inpp), & ars(inpp,3),vav(6),idist(iNd*2) common/cevap/Nevap,tau,aH,wH,lH,Rad0,Rad0sq,vbar,nu common/constants/NA,kb,hbar,pi,m4 common/sigmao/rsg08p,cv2,cv2er0,a1,radsq newcolls = 0.0 tdedi = tdead(ii) tbiri = tbirth(ii) tdedj = tdead(jj) tbirj = tbirth(jj) r0xi = r1(ii) r0yi = r2(ii) r0zi = r3(ii) vxi = v1(ii) vyi = v2(ii) vzi = v3(ii) dxji = r1(jj)-r0xi dyji = r2(jj)-r0yi dzji = r3(jj)-r0zi dvxji = v1(jj)-vxi dvyji = v2(jj)-vyi dvzji = v3(jj)-vzi DO 500 i=nloc1,nloc2 IF ( (tdead(i) .LE. tcurr) .OR. (i .EQ. ii) .OR. & (i .EQ. jj) ) GOTO 500 dx0i = r1(i)-r0xi dy0i = r2(i)-r0yi dz0i = r3(i)-r0zi vrxi = v1(i)-vxi vryi = v2(i)-vyi vrzi = v3(i)-vzi dx0j = dx0i-dxji dy0j = dy0i-dyji dz0j = dz0i-dzji vrxj = vrxi-dvxji vryj = vryi-dvyji vrzj = vrzi-dvzji vrsqi = vrxi**2+vryi**2+vrzi**2 vrsqj = vrxj**2+vryj**2+vrzj**2 c # distance of closest approach squared * vrsq rsqmi = (dx0i*vryi-dy0i*vrxi)**2 + (dx0i*vrzi-dz0i*vrxi)**2 & + (dy0i*vrzi-dz0i*vryi)**2 rsqmj = (dx0j*vryj-dy0j*vrxj)**2 + (dx0j*vrzj-dz0j*vrxj)**2 & + (dy0j*vrzj-dz0j*vryj)**2 c Calculate possible collision with particle ii IF (rsqmi .GT. radsq*vrsqi) GOTO 510 dr0vr = dx0i*vrxi+dy0i*vryi+dz0i*vrzi c conjectured collision time tc = -dr0vr/vrsqi IF ( (tc .GE. tdedi) .OR. (tc .LE. tbiri) .OR. &(tc.GE.tdead(i)).OR.(tc.LE.tbirth(i)) ) GOTO 510 IF (rsqmi .GT. rsg0(vrsqi)*vrsqi) GOTO 510 c # there is a collision at tc IF (idc0 .EQ. 0) THEN imcoll = imcoll+1 icoll = imcoll IF (icoll .GT. iNCM) THEN WRITE(6,*)'Error(newcolls): icoll > iNCM &: icoll,iNCM',icoll,iNCM END IF ELSE c # use a number from deadcoll list: icoll = ideadcoll(idc0) idc0 = idc0 - 1 END IF iposscoll = iposscoll + 1 icollp(icoll) = i jcollp(icoll) = ii tcoll(icoll) = tc vr2coll(icoll) = vrsqi dummy = prevnext(icoll,i,tc) dummy = prevnext(icoll,ii,tc) c Check for collision with particle jj 510 CONTINUE IF (rsqmj .GT. radsq*vrsqj) GOTO 500 c # conjectured collision time dr0vr = dx0j*vrxj+dy0j*vryj+dz0j*vrzj tc = -dr0vr/vrsqj IF ( (tc .LE. tbirth(i)) .OR. (tc .LE. tbirj) .OR. &(tc.GE.tdead(i)).OR.(tc.GE.tdedj)) GOTO 500 IF (rsqmj .GT. rsg0(vrsqj)*vrsqj) GOTO 500 c # there is a collision at tc IF (idc0 .EQ. 0) THEN imcoll = imcoll+1 icoll = imcoll IF (icoll .GT. iNCM) THEN WRITE(6,*)'Err(newcolls):icoll>iNCM:icoll,iNCM',icoll,iNCM END IF ELSE c # use a number from deadcoll list: icoll = ideadcoll(idc0) idc0 = idc0 - 1 END IF iposscoll = iposscoll + 1 icollp(icoll) = i jcollp(icoll) = jj tcoll(icoll) = tc vr2coll(icoll) = vrsqj dummy = prevnext(icoll,i,tc) dummy = prevnext(icoll,jj,tc) c # end of loop in i 500 CONTINUE END c end: # end of newcolls