what Sidd does sometimes

is flog a beowulf

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