what Sidd does sometimes
here is what he be floggin it with currently
and for u skeptics, it actually runs
and here is a picture with colors
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 00
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
back
