c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c this program generates a random distribution of points and then grows circles
c  or spheres to fill a domain (set idimen).
c To view grains launch gnuplot and type: load 'fort.13'
c
      parameter(igrains=40000)  !total number of grains
c
      integer igroup(igrains), icontact(10), imat, index(igrains)
      integer islice(3),ibox(3),idimen,imovie
      double precision xyzr(4,igrains),x,y,xdelta,ydelta,zdelta
      double precision box(3),zslice,volumei,pi, grains, rmag
      double precision dt,vmag, davg, vfrac, L, rdist(igrains),layer
      logical grow
      common /time/ imovie,itime
c
      open(3, file='meso.cth', status='REPLACE')
c
      idimen  = 2        ! this is the dimension either 2 or 3
      davg    = 0.020d0  ! cm  0.0032 cm = 32 um 
      box(1)  = 0.250d0  ! cm  - this is the longitudinal direction
      box(2)  = 0.250d0  ! cm
      box(3)  = 0.250d0  ! cm  - this value is ignored if idimen=2
      vfrac   = 0.5962264151d0  ! dry=0.596d0  ! volume fraction filled
      nummat  = 10 !number of different materials
      imat    = 1 !initate material assignment
      imoist  = 1  ! this turns on the moistue layer, mesomoist.cth
      layer   = 0.03d0*davg  !layer thickness
c
      pi = dacos(-1.d0)
      itime = 0
      imovie=101
c
c the size of the grains will grow slightly over davg to exactly achieve
c  the vfrac specified above.
c
      if (idimen . eq. 2) then
       grains = vfrac*box(1)*box(2)/(pi*(davg/2.d0)**2) ! 2D
       igrain = grains
c       grains = dsqrt(vfrac*box(1)*box(2)/(dfloat(igrain)*pi)) ! 2D
       davg   = ( vfrac*box(1)*box(2)/
     &  (dfloat(igrain)*pi ) )**(1.d0/2.d0)
       davg   = 2.d0 * davg ! davg above is actually a radius
c
      elseif (idimen .eq. 3) then
       grains = vfrac*box(1)*box(2)*box(3)/(4.d0*pi*(davg/2.d0)**3/3.d0)
       igrain = grains
       davg = ( vfrac*box(1)*box(2)*box(3)/
     &  (dfloat(igrain)*4.d0*pi/3.d0) )**(1.d0/3.d0)
       davg = 2.d0 * davg
       print *,'estimate igrain= ',
     &  vfrac*box(1)*box(2)*box(3)*6.d0/dacos(-1.d0)/davg**3.d0
c
      else
       print *,'Error: The variable idimen must be either 2 or 3'
       pause
      endif
       igrain = grains
c      if (grains .gt. igrain) igrain = igrain +1
c      box(1) = float(igrain)* pi*(davg/2.d0)**2/vfrac/box(2)
      print *,'This is a ',idimen,'-D calculation.'
      print *,'Grains= ',igrain
      print *,'davg  = ',davg
      pause
c
      igroup(1)      = imat
c
c  place grains in domain
         do i = 1,igrain
c          igrain = igrain + 1
c          print *,'igrain',iend,igrain
          xyzr(1,i) = rand(0)*box(1)
          xyzr(2,i) = rand(0)*box(2)
c      zslice = zdelta/2.d0 -davg/2.d0 + dfloat(i-1)*davg/10.d0
         if (idimen .eq. 2) then
          xyzr(3,i) = box(3)/2.d0  ! 2D: fix the z direction to be zero
         else 
          xyzr(3,i) = rand(0)*box(3) ! 3D
         endif
          xyzr(4,i) = 0.d0  !davg/2.d0
         enddo

      do i=1,igrain
       igroup(i) = 3  !make everyone blue
      enddo
c
      print *,'Done particle placement'
c=======================================================================
c grow one grain at a time
      grow = .true.
      igrow = 0
      grownum = 0
      volume = 0.d0
c
      do while ( grow )
      itime = itime + 1
      igrow = igrow + 1 ! this is the grain to grow for this time step
      grownum = grownum + 1
      imovie  = imovie  + 1
      ierr    = 0  ! there are no overlap errors
      if (igrow .gt. igrain) igrow = 1
      xyzr(4,igrow) = davg/2.0d0
c
c write to gnuplot
      islice(1) = 1     !islice(1) = xplane
      islice(2) = 2     !islice(2) = yplane
      islice(3) = 3     !islice(3) = cut plane
      zslice    = zdelta/2.d0 + dfloat(ibox(islice(3))/2-1)*zdelta
      zslice    = xyzr(3,1)
      call gnuwrite(imovie,box,islice,zslice,igrain,xyzr,igroup,
     &     layer,0)
c end gnuplot
        iplot=0  !do not plot gnumovie files as you remove overlap
        iter = 0  !number of overlap iterations, 0 means remove all overlap
        rmag = 0.5d0 + (vfrac-volume)/(2.d0*vfrac) 
        if(imovie .eq. 151) rmag = 0.5d0
        call overlap(igrow,igrain,box,xyzr,iplot,igroup,iter,rmag,ierr)
        call perturb2(igrow,igrain,ibox,box,xyzr,davg,vmag,iseed,idimen,
     &   igroup,imovie)
        if (ierr .ne. 0) then
         print *,'Error: error in overlap...'
         grow = .false.
         goto 99
        endif
c
c calculate volume fraction
 99     iigrow = 0
      volume = 0.d0
      do i=1,igrain
      if (xyzr(4,i) .lt. davg/2.d0) iigrow = 1
c       if (xyzr(1,i) .lt. box(1) .and. xyzr(1,i) .gt. 0.d0 .and.
c     &     xyzr(2,i) .lt. box(2) .and. xyzr(2,i) .gt. 0.d0 .and.
c     &     xyzr(3,i) .lt. box(3) .and. xyzr(3,i) .gt. 0.d0 ) then
       if (idimen .eq. 2) then
        volume=volume + pi*xyzr(4,i)**2  ! 2D
       else
        volume=volume + 4.d0*pi*xyzr(4,i)**3/3.d0  ! 3D
       endif ! idimen
c       endif
      enddo
c
      if (idimen .eq. 2) then
       volumne=volume/(box(1)*box(2))   !*box(3)),idumb
       print *,'Current vfrac= ',volume,'  Grownum =',igrow !num
      else
       volume = volume/(box(1)*box(2)*box(3))
       print *,'Current vfrac= ',volume
       print *,'Average Radius= ',xyzr(4,1)   
      endif
c
      if (iigrow .eq. 0) grow = .false.
      enddo ! do while ( grow )
c      pause
c=======================================================================
c write to gnuplot
      islice(1) = 1     !islice(1) = xplane
      islice(2) = 2     !islice(2) = yplane
      islice(3) = 3     !islice(3) = cut plane
      zslice    = zdelta/2.d0 + dfloat(ibox(islice(3))/2-1)*zdelta
      zslice    = xyzr(3,1)
      call gnuwrite(13,box,islice,zslice,igrain,xyzr,igroup,
     &     layer,0)

c=======================================================================
      iter = 0 !number of overlap iterations, 0 means remove all overlap
      rmag = 0.5d0 ! this means grains overlap will be completely removed
c      call overlap(igrain,box,xyzr,iplot,igroup,iter,rmag,ierr)
c
c      call perturb(igrain,ibox,box,xyzr,davg,vmag,iseed,idimen,
c     &   igroup,imovie)
c
c=======================================================================
c write to gnuplot
      islice(1) = 1     !islice(1) = xplane
      islice(2) = 2     !islice(2) = yplane
      islice(3) = 3     !islice(3) = cut plane
      zslice    = zdelta/2.d0 + dfloat(ibox(islice(3))/2-1)*zdelta
      zslice    = xyzr(3,i)
      call gnuwrite(14,box,islice,zslice,igrain,xyzr,igroup,
     &     layer,0)
      iter=0
c      call overlap(igrain,box,xyzr,1,igroup,iter,rmag,ierr)
c
c      islice(1) = 1     !islice(1) = xplane
c      islice(2) = 2     !islice(2) = yplane
c      islice(3) = 3     !islice(3) = cut plane
c      zslice    = box(3)/2.d0+zdelta
c      print *,'zslice',zslice,zdelta,ibox(islice(3)/2)
c      call gnuwrite(14,box,islice,zslice,igrain,xyzr,igroup,
c     &     layer,0)
cc
c      islice(1) = 1     !islice(1) = xplane
c      islice(2) = 3     !islice(2) = yplane
c      islice(3) = 2     !islice(3) = cut plane
c      zslice    = ydelta/2.d0
c      print *,'zslice',zslice,zdelta,ibox(islice(3)/2)
c      call gnuwrite(15,box,islice,zslice,igrain,xyzr,igroup,
c     &     layer,0)
c      
      call materialassign(imat,nummat,xyzr,igroup,igrain)
c      
c estimate final pack volume
      volume = 0.d0
      do i=1,igrain
       if (xyzr(1,i) .lt. box(1) .and. xyzr(1,i) .gt. 0.d0 .and.
     &     xyzr(2,i) .lt. box(2) .and. xyzr(2,i) .gt. 0.d0 .and.
     &     xyzr(3,i) .lt. box(3) .and. xyzr(3,i) .gt. 0.d0 ) then
        if (idimen .eq. 2) then
         volume=volume + pi*xyzr(4,i)**2
        else
         volume=volume + 4.d0*pi*xyzr(4,i)**3/3.d0
        endif
       endif
       enddo
c
      if (idimen .eq. 2) then
       volume = volume/(box(1)*box(2))
      else 
       volume = volume/(box(1)*box(2)*box(3))
      endif
      print *,'Final vfrac= ',volume
c generate several gnuplot files     
      do i=2,10
      zslice = box(3)/2.d0 -davg/2.d0 + dfloat(i-1)*davg/10.d0
c      print *,'zslice= ',zslice
      call gnuwrite(15+i,box,islice,zslice,igrain,xyzr,igroup,
     &     layer,0)
      enddo
c
      call cthwrite(3,idimen,igrain,nummat,box,xyzr,igroup)
      if (imoist .eq. 1) then
      print *,'Writing moisture layer...'
c      pause
      call moist(layer,idimen,igrain,nummat,box,xyzr,igroup)
      endif
c      
      end
c
c======================================================================
      subroutine overlap(i,igrain,box,xyzr,iplot,igroup,iter,rmag,ierr)
c
      parameter (igrains=40000)
c
      double precision box(3), xyzr(4,igrain), 
     .                 davg, rr(igrains), vv(3,igrains), f(3,igrains)
      double precision   pi, k, L,
     .                 theta, phi,t,dt,tmax, vmag, r, xfac,yfac,zfac,
     .                 fmag, rdist, roverlap, disp, vfrac,
     .                 v, vmax, xvmax, yvmax, mag,totalover,
     &                 xshift, yshift,zshift,zslice,rmag
      INTEGER i,j,igrain,iplot, ioverlap(igrains), icount, ivmax
     &        istart, igroup(igrains),ipert, ifixed(igrains),fix,
     &        imovie,islice(3),itan
c
      LOGICAL loverlap
      common /time/ imovie,itime
c
      ignuprint = 0
      ierr = 0
      pi = dacos(-1.d0)
c
c clear ifixed
      fix = igrain+1
c      do i=1,igrain
c       ifixed(i)=0
c      enddo      
      ifixnum   = 1
      igroup(fix)=2
      ifixed(1) = fix  !make the grain that grew fixed in space and
c                        force those around it to grow      
      loverlap = .TRUE.
      icount = 0
      do while ( loverlap )
       icount    = icount + 1
       icontact  = 0
       totalover = 0.d0  ! clear total amount of overlap
c
c       WRITE(*,'(70("="),/, 20X, i6, /,70("="))') icount
c       do i=1,igrain
c        ioverlap(i) = 0
c       enddo
c
c   see which circles overlap with others and by how much
c       do i=1,igrain
        ii0 = 0  !-1  !this was done in an effort to make the search smarter
        ii1 = 0  ! 1 ! and therefore faster
        jj0 = 0  !-1 
        jj1 = 0  ! 1
        kk0 = 0  !-1
        kk1 = 0  ! 1
        if (xyzr(1,i) .lt. 0.d0   + 3.d0*xyzr(4,i) ) ii1 =   1
        if (xyzr(1,i) .gt. box(1) - 3.d0*xyzr(4,i) ) ii0 =  -1
        if (xyzr(2,i) .lt. 0.d0   + 3.d0*xyzr(4,i) ) jj1 =   1
        if (xyzr(2,i) .gt. box(2) - 3.d0*xyzr(4,i) ) jj0 =  -1
        if (xyzr(3,i) .lt. 0.d0   + 3.d0*xyzr(4,i) ) kk1 =   1
        if (xyzr(3,i) .gt. box(3) - 3.d0*xyzr(4,i) ) kk0 =  -1
c        print *,'ii0,ii1,jj0,jj1,kk0,kk1',ii0,ii1,jj0,jj1,kk0,kk1  
        do ii = ii0 , ii1
        do jj = jj0 , jj1
        do kk = kk0 , kk1
         xshift = dfloat(ii)*box(1)
         yshift = dfloat(jj)*box(2)
         zshift = dfloat(kk)*box(3)
c         print *,'ii,jj,kk ',i,ii,jj,kk
         do j=1, igrain
         if (j .ne. i) then
          r = DSQRT(( (xyzr(1,i)+xshift) - xyzr(1,j) )**2 
     &            + ( (xyzr(2,i)+yshift) - xyzr(2,j) )**2
     &            + ( (xyzr(3,i)+zshift) - xyzr(3,j) )**2)
          if (xyzr(4,i)+xyzr(4,j) - r .gt. 
     &                    1.0d-4*(xyzr(4,i)+xyzr(4,j))/2.d0) then
c           icontact  = icontact  + 1
           ignuprint = ignuprint + 1
           totalover = totalover + xyzr(4,i)+xyzr(4,j)-r
c           
           disp      = xyzr(4,i)+xyzr(4,j)-r
           xfac      = (xyzr(1,i)+xshift - xyzr(1,j) )/r
           yfac      = (xyzr(2,i)+yshift - xyzr(2,j) )/r 
           zfac      = (xyzr(3,i)+zshift - xyzr(3,j) )/r           
c
           if (i .eq. j) then
                   print *,'error: i = j'
                   pause
           endif
c
           imove = 1
           jmove = 0
c
           if    (imove .eq. 0 .and. jmove .eq. 1) then
             xyzr(1,j) =  xyzr(1,j) - 1.d0*disp*xfac
             xyzr(2,j) =  xyzr(2,j) - 1.d0*disp*yfac 
             xyzr(3,j) =  xyzr(3,j) - 1.d0*disp*zfac
c keep particles in domain
         if (xyzr(1,j) .lt. 0.00d0) xyzr(1,j) = xyzr(1,j) + box(1)
         if (xyzr(1,j) .gt. box(1)) xyzr(1,j) = xyzr(1,j) - box(1)
         if (xyzr(2,j) .lt. 0.00d0) xyzr(2,j) = xyzr(2,j) + box(2)
         if (xyzr(2,j) .gt. box(2)) xyzr(2,j) = xyzr(2,j) - box(2)
         if (xyzr(3,j) .lt. 0.00d0) xyzr(3,j) = xyzr(3,j) + box(3)
         if (xyzr(3,j) .gt. box(3)) xyzr(3,j) = xyzr(3,j) - box(3)
c         print *,'Going to scooch...'
c      call scooch(igrain,box,xyzr,i,j,xfac,yfac,zfac,disp)
c            pause
c            ifixnum = ifixnum + 1
c            ifixed(ifixnum) = j
            elseif(imove .eq. 1 .and. jmove .eq. 0) then
             xyzr(1,i) =  xyzr(1,i) + 1.d0*disp*xfac
             xyzr(2,i) =  xyzr(2,i) + 1.d0*disp*yfac
             xyzr(3,i) =  xyzr(3,i) + 1.d0*disp*zfac
c            pause
c            ifixnum = ifixnum + 1
c            ifixed(ifixnum) = i
           elseif (imove .eq. 1 .and. jmove .eq. 1) then
c            disp=rmag*disp !+ 0.5d0*disp/dfloat(icount)
c            if (icontact .lt. 10) then
c            disp = 0.5d0*disp + 0.5d0*disp/dfloat(icontact)
c            else
            disp = 0.5d0*disp
c            endif
            xyzr(1,i) = xyzr(1,i) + disp*xfac
            xyzr(2,i) = xyzr(2,i) + disp*yfac
            xyzr(3,i) = xyzr(3,i) + disp*zfac
            xyzr(1,j) = xyzr(1,j) - disp*xfac
            xyzr(2,j) = xyzr(2,j) - disp*yfac 
            xyzr(3,j) = xyzr(3,j) - disp*zfac
c            print *,i,j,disp,xyzr(1,i),xyzr(2,i),xfac,yfac
c            print *,j,i,disp,xyzr(1,j),xyzr(2,j)
c            ifixnum = ifixnum + 1
c            ifixed(ifixnum) = i 
c            ifixnum = ifixnum + 1
c            ifixed(ifixnum) = j 
           else
            print *,'Error: no grains to move' 
            pause
           endif
c keep particles in domain
         if (xyzr(1,i) .lt. 0.00d0) xyzr(1,i) = xyzr(1,i) + box(1)
         if (xyzr(1,i) .gt. box(1)) xyzr(1,i) = xyzr(1,i) - box(1)
         if (xyzr(2,i) .lt. 0.00d0) xyzr(2,i) = xyzr(2,i) + box(2)
         if (xyzr(2,i) .gt. box(2)) xyzr(2,i) = xyzr(2,i) - box(2)
         if (xyzr(3,i) .lt. 0.00d0) xyzr(3,i) = xyzr(3,i) + box(3)
         if (xyzr(3,i) .gt. box(3)) xyzr(3,i) = xyzr(3,i) - box(3)
c         if (xyzr(1,j) .lt. 0.00d0) xyzr(1,j) = xyzr(1,j) + box(1)
c         if (xyzr(1,j) .gt. box(1)) xyzr(1,j) = xyzr(1,j) - box(1)
c         if (xyzr(2,j) .lt. 0.00d0) xyzr(2,j) = xyzr(2,j) + box(2)
c         if (xyzr(2,j) .gt. box(2)) xyzr(2,j) = xyzr(2,j) - box(2)
c         if (xyzr(3,j) .lt. 0.00d0) xyzr(3,j) = xyzr(3,j) + box(3)
c         if (xyzr(3,j) .gt. box(3)) xyzr(3,j) = xyzr(3,j) - box(3)
c
c write to gnuplot
        if (iplot       .eq. 1 .and. 
     &      ignuprint/1 .eq. float(ignuprint)/1. .and.
     &      1 .eq. 1 ) then
         islice(1) = 1     !islice(1) = xplane
         islice(2) = 2     !islice(2) = yplane
         islice(3) = 3     !islice(3) = cut plane
c         zslice    = zdelta/2.d0 + dfloat(ibox(islice(3))/2-1)*zdelta
         zslice    = xyzr(3,1)
         imovie = imovie + 1
         call gnuwrite(imovie,box,islice,zslice,igrain,xyzr,igroup,
     &       0.d0,0)
         print *,'----------------- overlap movie write --------------'
         print *,'imovie =',imovie,' icontact = ',icontact
        endif
c
           endif  ! overlap
          endif 
          enddo   ! j loop
          enddo   !kk loop
         enddo    !jj loop
         enddo    !ii loop
c
c        enddo     ! i loop = 1,igrain
c
c keep particles in domain
c        do i=1,igrain
c         if (xyzr(1,i) .lt. 0.00d0) xyzr(1,i) = xyzr(1,i) + box(1)
c         if (xyzr(1,i) .gt. box(1)) xyzr(1,i) = xyzr(1,i) - box(1)
c         if (xyzr(2,i) .lt. 0.00d0) xyzr(2,i) = xyzr(2,i) + box(2)
c         if (xyzr(2,i) .gt. box(2)) xyzr(2,i) = xyzr(2,i) - box(2)
c         if (xyzr(3,i) .lt. 0.00d0) xyzr(3,i) = xyzr(3,i) + box(3)
c         if (xyzr(3,i) .gt. box(3)) xyzr(3,i) = xyzr(3,i) - box(3)
c        enddo ! ii loop to keep particles in domain
c         
      print *,'Iter: ',icount,' Gains:',icontact
     &   ,totalover
     & ,' Avg overlap ',totalover/dfloat(icontact)
c
         icontact = 0
      if(icontact .eq. 0) loverlap = .FALSE.
        if(icount .eq. iter .and. itan .ne. 0) then 
         loverlap = .FALSE.
         write(*,*) 'WARNING:  Terminating with overlap...'
c         pause
        endif
      enddo !overlap = true then continue to iterate
c
c      do ii = 1, igrain
c         igroup(ii) = 3
c      enddo
 99   return
      end
c======================================================================
      subroutine materialassign(imat,nummat,xyzr,igroup,igrain)
c
      parameter(igrains=40000)
c
      integer ii,jj,isave,igrain,igroup(igrains),index(igrains)
      integer icontact(igrains)
      double precision xyzr(4,igrains),rdist(igrains)
c
c  randomize the grain material numbering scheme in order to insure good mixture
c 
      do i = 1,igrain
       index(i)  = i !index will be the randimized array of igrain
       igroup(i) = 0
       enddo   ! igrain
c
      do i=1,igrain*4
       ii = rand(0)*igrain
       jj = rand(0)*igrain
       if (ii .le. 0      ) ii = 1
       if (ii .ge. igrain ) ii = igrain
       if (jj .le. 0      ) jj = 1
       if (jj .ge. igrain ) jj = igrain
       isave     = index(ii) !randomize by switching index
       index(ii) = index(jj)
       index(jj) = isave
      enddo
c 
c  insure that circles touching aren't in the same group
      do i=1,igrain
       jj = index(i)  ! this is the grain # we are assigning a group to
        do ii=1,nummat
           icontact(ii) = 0
        enddo
        do ii=1, igrain-1
           rdist(ii) = dsqrt( (xyzr(1,ii)-xyzr(1,jj))**2
     &                      + (xyzr(2,ii)-xyzr(2,jj))**2
     &                      + (xyzr(3,ii)-xyzr(3,jj))**2 )
           if (rdist(ii) .le. 2.d0 * (xyzr(4,ii)+xyzr(4,jj)) ) then
              icontact(igroup(ii)) = 1
           endif
        enddo
        do k=1,nummat
           if(icontact(imat) .eq. 1) then
              imat=imat+1
              if(imat .eq. nummat+1) imat=1
           endif
        enddo
c
        igroup(jj) = imat
c        igroup(jj) = 3 ! this will make all circles blue
        imat = imat + 1
        if(imat .eq. nummat+1) imat=1
      enddo  ! i-loop
c
      return
      end
c======================================================================
      subroutine periodic(iperiodic,igrain,box,xyzr,igroup)
c
      parameter(igrains=40000)
c
      integer ii,jj,kk,igrain,islice(3),igroup(igrains),iflag
      double precision box(3),xyzr(4,igrains)
      double precision xshift,yshift,zshift
c
c add periodic particles
c
      if (1 .eq. 1 ) then
      iflag = 0
      do i = 1,igrain
      if (xyzr(1,i) .gt. box(1) .or. xyzr(1,i) .lt. 0.d0) then
       print *,'xyzr(x,i) ',i,xyzr(1,i),xyzr(2,i),xyzr(3,i)
       iflag=iflag+1
      endif
      if (xyzr(2,i) .gt. box(2) .or. xyzr(2,i) .lt. 0.d0) then
       iflag=iflag+1
       print *,'xyzr(x,i) ',i,xyzr(1,i),xyzr(2,i),xyzr(3,i)
      endif
      if (xyzr(3,i) .gt. box(3) .or. xyzr(3,i) .lt. 0.d0) then
       iflag=iflag+1
       print *,'xyzr(x,i) ',i,xyzr(1,i),xyzr(2,i),xyzr(3,i)
      endif
      enddo
      if (iflag .ne. 0) then
       print *,'Error: ',iflag,' particles out of domain'
       print *,i,xyzr(1,i),xyzr(2,i),xyzr(3,i),box(1),box(2),box(3)
c       pause
      endif
c
      iperiodic = 0
      do i=1,igrain
      do ii = -1,1
      do jj = -1,1
      do kk = -1,1
      if (ii .eq. 0 .and. jj .eq. 0 .and. kk .eq. 0) then
      else
       xshift = dfloat(ii)*box(1)
       yshift = dfloat(jj)*box(2)
       zshift = dfloat(kk)*box(3)
       if( xyzr(1,i) + xshift - dfloat(ii)*xyzr(4,i) .gt. 0.0d0  .and.
     &     xyzr(1,i) + xshift - dfloat(ii)*xyzr(4,i) .lt. box(1) .and.
     &     xyzr(2,i) + yshift - dfloat(jj)*xyzr(4,i) .gt. 0.0d0  .and.
     &     xyzr(2,i) + yshift - dfloat(jj)*xyzr(4,i) .lt. box(2) .and.
     &     xyzr(3,i) + zshift - dfloat(kk)*xyzr(4,i) .gt. 0.0d0  .and.
     &     xyzr(3,i) + zshift - dfloat(kk)*xyzr(4,i) .lt. box(3)
     &   ) then
        iperiodic      = iperiodic + 1
        xyzr(1,igrain+iperiodic) = xyzr(1,i) + xshift
        xyzr(2,igrain+iperiodic) = xyzr(2,i) + yshift
        xyzr(3,igrain+iperiodic) = xyzr(3,i) + zshift
        xyzr(4,igrain+iperiodic) = xyzr(4,i)
        igroup(igrain+iperiodic) = igroup(i) ! this should be igroup(i), 1 = red
       endif ! is shifted particle in the box
      endif
      enddo  ! kk-loop
      enddo  ! jj-loop
      enddo  ! ii-loop
      enddo  !  i-loop
      endif
c      
      return
      end
c
c======================================================================
      subroutine gnuwrite(ifile,box,islice,zslice,
     &       igrain,xyzr,igroup,layer,iflag)
c
      parameter(igrains=40000)
c
      integer in,jn,kn,igrain,islice(3),igroup(igrains),iflag
      double precision box(3),xyzr(4,igrains), radius, radius2, zslice
      double precision layer
      character*48 name
c
      if (ifile .le. 9) then
      write(name,'(a9,I1)') 'gnumovie.',ifile
      elseif (ifile .le. 99) then
      write(name,'(a8,I2)') 'gnuplot.',ifile
      elseif (ifile .le. 999) then
      write(name,'(a9,I3)') 'gnumovie.',ifile
      elseif (ifile .le. 9999) then
      write(name,'(a9,I4)') 'gnumovie.',ifile
      endif
c
      open(10, file=name, status='REPLACE')
c write output file
c
      print *,'ifile= ',ifile
c      write(10,'(a)')         'set term X11'   !aqua'
      write(10,'(a)')         'set parametric'
      write(10,'(a)')         'set key off'
      write(10,'(a)')         'pi = acos(-1)'
      write(10,'(a,f10.5)')   'set size ratio ',
     &              box(islice(2))/box(islice(1))
      write(10,'(a,f10.5,a)') 'set xrange [ 0.00:',box(islice(1)),']'
      write(10,'(a,f10.5,a)') 'set yrange [ 0.00:',box(islice(2)),']'
      write(10,'(a)')         'set trange [0:2*pi]'
      write(10,'(a)')         'plot \\'

c
      isave = 0
      call periodic(iperiodic,igrain,box,xyzr,igroup) ! this adds the periodic grains
c
       do i=1,igrain+iperiodic
c
        if (xyzr(islice(3),i)+xyzr(4,i) .gt. zslice .and.
     &      xyzr(islice(3),i)-xyzr(4,i) .lt. zslice       ) then
         radius = dsqrt( xyzr(4,i)**2 -
     &   (dabs(xyzr(islice(3),i)-zslice))**2 )
         radius2= dsqrt((xyzr(4,i)+layer)**2 -
     &   (dabs(xyzr(islice(3),i)-zslice))**2 )

         isave  = i
c         if (radius .gt. 0.0032d0/2.d0) then
c                 print *,'big radius',xyzr(4,i),radius
c                 pause
c         endif
        write(10,901) xyzr(islice(1),i),'+',radius,'*cos(t),',
     &                   xyzr(islice(2),i),'+',radius,'*sin(t)
     &                  ls ',igroup(i) ,',\\' !igroup(i)
        if (iflag .eq. 1) then  !this shows moisture layer
        write(10,901) xyzr(islice(1),i),'+',radius2,'*cos(t),',
     &                   xyzr(islice(2),i),'+',radius2,'*sin(t)
     &                  ls ',igroup(i) ,',\\' !igroup(i)
        endif
c
        if (xyzr(islice(1),i) .gt. box(islice(1)) .or. 
     &      xyzr(islice(1),i) .lt. 0.00d00        .or.
     &      xyzr(islice(2),i) .gt. box(islice(2)) .or. 
     &      xyzr(islice(2),i) .lt. 0.00d00       ) then
c         print *,'BoxOut',xyzr(islice(1),i),xyzr(islice(2),i),
c     &                    xyzr(islice(3),i)
        endif
       endif !zslice check
       enddo ! grain loop
c
       if ( isave .eq. 0 ) then
c re-write last one just to get the gnuplot end-of-file right
        print *,'Warning: No particles written to gnuplot',ifile
        print *,islice(1),islice(2),islice(3)
        print *,box(1),box(2),box(3)
        pause
       else
        write(10,901) xyzr(islice(1),isave),'+',radius,'*cos(t),',   !re-writes last particle to end gnufile correctly
     &                   xyzr(islice(2),isave),'+',radius,'*sin(t)
     &                   ls ', igroup(isave)
       endif
c
  901 format(f8.6,a,f8.6,a,3x,f8.6,a,f8.6,a,i2,a,2x)
c
      close(10)
c
      return
      end
c=========================================================================
      subroutine cthwrite(ifile,idimen,igrain,nummat,box,xyzr,igroup)
c
      parameter(igrains=40000)
c
      integer igrain,igroup(igrains),iperiodic,idimen
      double precision box(3),xyzr(4,igrains)
c
      iperiodic = 0
      call periodic(iperiodic,igrain,box,xyzr,igroup) ! this adds the periodic grains
c
c write output file
c
      open(3, file='meso.cth', status='REPLACE')
c
      do imat = 1,nummat
       print *,'cthwrite',imat
       write(3,'(/, 6X, A,I1)') 'package wc0', imat
       write(3,'(9X, A, I3)') 'material', imat
       write(3,'(9X, A)') 'nsub 100'
c
       do i=1,igrain+iperiodic
        if (   igroup(i) .eq. imat ) then
         if (   idimen    .eq.    2 ) then
          write(3,902) 'center = ',xyzr(1,i),xyzr(2,i),
     &                 'radius = ',xyzr(4,i)
         elseif(idimen    .eq.    3 ) then
          write(3,903) 'center = ',xyzr(1,i),xyzr(2,i),xyzr(3,i),
     &                 'radius = ',xyzr(4,i)
         endif
        endif  ! igroup = imat
       enddo ! grain loop
       write(3,'(6X,A)') 'endpackage'
      enddo ! imat loop
c
  902    format(9X, 'insert circle', /, 12X, A, F12.8,',',F12.8,
     &   /,12X, A, F12.8, /, 9X, 'endinsert')
  903    format(9X, 'insert sphere', /, 12X, A, F12.8,',',F12.8, ',', 
     &    F12.8, /,12X, A, F12.8, /, 9X, 'endinsert')
c
      return
      end
c======================================================================
      subroutine moist(layer,idimen,igrain,nummat,box,xyzr,igroup)
c           call moist(layer,idimen,igrain,nummat,box,xyzr,igroup)
c
      parameter(igrains=40000)
c
      integer igrain,igroup(igrains),iperiodic,idimen,islice(3)
      double precision layer,box(3),xyzr(4,igrains),zslice
c
      iperiodic = 0
      call periodic(iperiodic,igrain,box,xyzr,igroup) ! this adds the periodic grains
c
c write output file
c
      print *,'layer',layer
c      pause
      open(3, file='mesomoist.cth', status='REPLACE')
c
      do imat = 1,nummat
       print *,'cthwrite',imat
       write(3,'(/, 6X, A,I2)') 'package wc11', 11
       write(3,'(9X, A, I2)') 'material ', 14  !imat
       write(3,'(9X, A)') 'nsub 100'
c
       do i=1,igrain+iperiodic
        if (   igroup(i) .eq. imat ) then
         if (   idimen    .eq.    2 ) then
          write(3,902) 'center = ',xyzr(1,i),xyzr(2,i),
     &                 'radius = ',xyzr(4,i)+layer
         elseif(idimen    .eq.    3 ) then
          write(3,903) 'center = ',xyzr(1,i),xyzr(2,i),xyzr(3,i),
     &                 'radius = ',xyzr(4,i)+layer
         endif
        endif  ! igroup = imat
       enddo ! grain loop
       write(3,'(6X,A)') 'endpackage'
      enddo ! imat loop
c
c write to gnuplot
      islice(1) = 1     !islice(1) = xplane
      islice(2) = 2     !islice(2) = yplane
      islice(3) = 3     !islice(3) = cut plane
c      zslice    = zdelta/2.d0 + dfloat(ibox(islice(3))/2-1)*zdelta
      zslice    = xyzr(3,1)   !box(3)/2.d0
      call gnuwrite(99,box,islice,zslice,igrain,xyzr,igroup,
     &     layer,1)
c end gnuplot
  902    format(9X, 'insert circle', /, 12X, A, F12.8,',',F12.8,
     &   /,12X, A, F12.8, /, 9X, 'endinsert')
  903    format(9X, 'insert sphere', /, 12X, A, F12.8,',',F12.8, ',', 
     &    F12.8, /,12X, A, F12.8, /, 9X, 'endinsert')
c
      return
      end
c======================================================================
      subroutine perturb2(i,igrain,ibox,box,xyzr,davg,vmag,iseed,idimen,
     &   igroup,imovie)
c
c move particle towards the void centroid
c
      parameter (igrains=40000)
c
      double precision box(3), xyzr(4,igrain),  
     .                 davg, rr(igrains), vv(3,igrains), f(3,igrains)
      double precision   pi, k, L, rmax,rmin,rrmin(100,100,100),
     .                 theta, phi,t,dt,tmax, vmag, r, xfac,yfac,zfac,
     .                 fmag, rdist, roverlap, disp, vfrac,
     .                 v, vmax, xvmax, yvmax, mag,
     &                 xshift, yshift,zshift,zslice
      INTEGER i,j,igrain, iseed(2), ioverlap(igrains), icount, ivmax
     &        istart, igroup(igrains),ipert,islice(3),ibox(3)
c
      LOGICAL overlap
c
      print *,'in pert',i
         vv(1,i) = 0.d0  
         vv(2,i) = 0.d0
         vv(3,i) = 0.d0
          f(1,i) = 0.d0
          f(2,i) = 0.d0
          f(3,i) = 0.d0
c     
      do ii=1,igrain
       rmin = 9.d30
       do jj=1,igrain
        if (ii .ne. jj) then
         r = dsqrt( (xyzr(1,ii)-xyzr(1,jj) )**2
     &           +  (xyzr(2,ii)-xyzr(2,jj) )**2
     &           +  (xyzr(3,ii)-xyzr(3,jj) )**2 )
c        
         if (r .lt. rmin ) rmin = r
        endif
       enddo ! jj-loop
       rsum = rsum + rmin
      enddo ! ii-loop
      L = rsum/dfloat(igrain-1)
      print *,'L here',L
c
      istep  = 100
      mfp    = 2 
      istart = 1
      vmag  = (L/2.d0-davg)*dfloat(mfp)/dfloat(isteps) !this is the velocity in order to complete the transits
      vmag = davg/50.d0
      dt    = 1.d0  !vmag*dfloat(mfg)/dfloat(isteps)  !time step
      tmax  = dt*dfloat(isteps)
      k = 0.025d0/(dt*dt)
c
      imove = 0
      overlap = .true.
      do while ( overlap )
      icontact = 0
       imove=imove+1
       if (imove .eq. 50) then
        print *,'Finding void center...'
        do ii = 1,100
         do jj = 1,100
          do kk = 1,100
           rrmin(ii,jj,kk) = 9.d9
           do iii = 1,igrain
             r = sqrt( (xyzr(1,iii)-box(1)/dfloat(ii))**2
     &               + (xyzr(2,iii)-box(2)/dfloat(jj))**2
     &               + (xyzr(3,iii)-box(3)/dfloat(kk))**2 )
            if (r .lt. rrmin(ii,jj,kk)) then
             rrmin(ii,jj,kk) = r
            endif
           enddo  ! iii
          enddo
         enddo
        enddo
        rmax = 0.d0
        do ii = 1,100
         do jj = 1,100
          do kk = 1,100
           if (rrmin(ii,jj,kk) .gt. rmax) then
            rmax = rrmin(ii,jj,kk)
            ivoid = ii
            jvoid = jj
            kvoid = kk
           endif
           enddo  ! iii
          enddo
         enddo
        print *,'Max void',box(1)/dfloat(ivoid),
     &       box(2)/dfloat(jvoid),box(3)/dfloat(kvoid)
        pause
        endif
c  check for contact between cylinders
       ii0 = 0 !-1  !this was done in an effort to make the search smarter
       ii1 = 0 ! 1 ! and therefore faster
       jj0 = 0 !-1 
       jj1 = 0 ! 1
       kk0 = 0 !-1
       kk1 = 0 ! 1
       if (xyzr(1,i) .lt. 0.d0   + 3.d0*xyzr(4,i) ) ii1 =   1
       if (xyzr(1,i) .gt. box(1) - 3.d0*xyzr(4,i) ) ii0 =  -1
       if (xyzr(2,i) .lt. 0.d0   + 3.d0*xyzr(4,i) ) jj1 =   1
       if (xyzr(2,i) .gt. box(2) - 3.d0*xyzr(4,i) ) jj0 =  -1
       if (xyzr(3,i) .lt. 0.d0   + 3.d0*xyzr(4,i) ) kk1 =   1
       if (xyzr(3,i) .gt. box(3) - 3.d0*xyzr(4,i) ) kk0 =  -1
c       print *,'ii0,ii1,jj0,jj1,kk0,kk1',ii0,ii1,jj0,jj1,kk0,kk1  
        do ii = ii0 , ii1
        do jj = jj0 , jj1
        do kk = kk0 , kk1   !this translates i particle for periodic contact check
         xshift = dfloat(ii)*box(1)
         yshift = dfloat(jj)*box(2)
         zshift = dfloat(kk)*box(3)
         do j=1,igrain !start at i+1 so particles don't get double loaded.
          if(i .ne. j) then
           r = dsqrt( (xyzr(1,i)+xshift-xyzr(1,j) )**2 
     &             +  (xyzr(2,i)+yshift-xyzr(2,j) )**2
     &             +  (xyzr(3,i)+zshift-xyzr(3,j) )**2 )
          if (r .lt. xyzr(4,i)+xyzr(4,j) .and. 
     &          xyzr(4,j) .gt. 0.d0) then
           icontact = icontact + 1
c            print *,'chub',i,j,r
c            print *,xyzr(1,i),xyzr(2,i),xyzr(3,i),xyzr(4,i)
c            print *,xyzr(1,j),xyzr(2,j),xyzr(3,j),xyzr(4,j)
c            pause
c actual contact
           if (imove .le. 50) then
            fmag = 0.5*k*(dabs(r - xyzr(4,i) - xyzr(4,j)) + 5.d-2)
            xfac = (xyzr(1,i)+xshift-xyzr(1,j) )/r
            yfac = (xyzr(2,i)+yshift-xyzr(2,j) )/r
            zfac = (xyzr(3,i)+zshift-xyzr(3,j) )/r
            f(1,i) = f(1,i) + fmag*xfac
            f(2,i) = f(2,i) + fmag*yfac
            if (idimen .eq. 2) then
             f(3,i) = 0.d0
            elseif (idimen .eq. 3) then
             f(3,i) = f(3,i) + fmag*zfac
            endif
            igroup(j) = 4
            print *,'fac',igrain,i,j,xfac,yfac,zfac
            print *,'xyz',xyzr(1,j),xyzr(2,j),xyzr(3,j),xyzr(4,j)
           endif
          endif
c push towards void space                   
c          if (imove .gt. 50) then
c           fmag = 0.5*k*1.d-3  !dabs((r - xyzr(4,i) - xyzr(4,j)))
c           xfac = (xyzr(1,i)+xshift-xyzr(1,j) )/r
c           yfac = (xyzr(2,i)+yshift-xyzr(2,j) )/r
c           zfac = (xyzr(3,i)+zshift-xyzr(3,j) )/r
c           f(1,i) = f(1,i) + fmag*xfac
c           f(2,i) = f(2,i) + fmag*yfac
c           if (idimen .eq. 2) then
c            f(3,i) = 0.d0
c           elseif (idimen .eq. 3) then
c            f(3,i) = f(3,i) + fmag*zfac
c           endif
c          endif

            endif ! i .ne. j
c           endif ! r < xyzr(4,i)+xyzr(4,j)
          ENDDO !  j loop
         enddo  ! kk loop
         enddo  ! jj loop
         enddo  ! ii loop
c push towards void space
         if (imove .gt. 50) then
            r = sqrt( (xyzr(1,i)-box(1)/dfloat(ivoid))**2
     &              + (xyzr(2,i)-box(2)/dfloat(jvoid))**2
     &              + (xyzr(3,i)-box(3)/dfloat(kvoid))**2 )
            fmag = 0.5*k*dabs(r)
            xfac = (xyzr(1,i)-box(1)/dfloat(ivoid) )/r
            yfac = (xyzr(2,i)-box(2)/dfloat(jvoid) )/r
            zfac = (xyzr(3,i)-box(3)/dfloat(kvoid) )/r
            f(1,i) = f(1,i) - fmag*xfac
            f(2,i) = f(2,i) - fmag*yfac
            if (idimen .eq. 2) then
             f(3,i) = 0.d0
            elseif (idimen .eq. 3) then
             f(3,i) = f(3,i) - fmag*zfac
            endif
         endif
c
         if (imove .gt. 200) icontact = 0
         if (icontact .eq. 0) then
          overlap = .false.
          goto 99
         endif
c
c update velocities
        vmax = 0.d0
         vv(1,i) = vv(1,i) + f(1,i)*dt
         vv(2,i) = vv(2,i) + f(2,i)*dt
         if (idimen .eq. 2) then
          vv(3,i) = 0.d0
         elseif (idimen .eq. 3) then
          vv(3,i) = vv(3,i) + f(3,i)*dt
         else
          print *,'Error: check idimen'
          pause
         endif

c   calculate position explicitly for new time
         xyzr(1,i) = xyzr(1,i) + vv(1,i)*dt
         xyzr(2,i) = xyzr(2,i) + vv(2,i)*dt
         xyzr(3,i) = xyzr(3,i) + vv(3,i)*dt

c keep particles in box domain
         if (xyzr(1,i) .lt. 0.00d0) xyzr(1,i) = xyzr(1,i) + box(1)
         if (xyzr(1,i) .gt. box(1)) xyzr(1,i) = xyzr(1,i) - box(1)
         if (xyzr(2,i) .lt. 0.00d0) xyzr(2,i) = xyzr(2,i) + box(2)
         if (xyzr(2,i) .gt. box(2)) xyzr(2,i) = xyzr(2,i) - box(2)
         if (xyzr(3,i) .lt. 0.00d0) xyzr(3,i) = xyzr(3,i) + box(3)
         if (xyzr(3,i) .gt. box(3)) xyzr(3,i) = xyzr(3,i) - box(3)
        if (icontact .gt. 0) then
        print *,'Brownian motion ',i,imove,icontact
c        print *,f(1,i),f(2,i),f(3,i),vv(1,i),vv(2,i),vv(3,i)
        endif
c
            f(1,i) = 0.d0
            f(2,i) = 0.d0
            f(3,i) = 0.d0
           vv(1,i) = 0.d0
           vv(2,i) = 0.d0
           vv(3,i) = 0.d0
c write to gnuplot
c      if (int(t)/3 .eq. t/3.d0) then
      islice(1) = 1     !islice(1) = xplane
      islice(2) = 2     !islice(2) = yplane
      islice(3) = 3     !islice(3) = cut plane
      zslice    = zdelta/2.d0 + dfloat(ibox(islice(3))/2-1)*zdelta
      zslice    = xyzr(3,1)
      imovie = imovie + 1
      igroup(i) = 2
      call gnuwrite(imovie,box,islice,zslice,igrain,xyzr,igroup,
     &     0.d0,0)
      igroup(i)=3
c      endif
c end gnuplot
c
      do ii = 1,igrain
       igroup(ii) = 3
      enddo
      if (imove .eq. 50) pause
 99   enddo !do while
      if (imove .eq. 3) pause
      return
      end
c===============================================================================
      subroutine scooch(igrain,box,xyzr,i,j,xfac,yfac,zfac,disp)
c this subroutine schooches all partiles over 
c
      parameter (igrains=40000)
c
      double precision box(3), xyzr(4,igrain),xfac,yfac,zfac,disp
      INTEGER i,j,igrain,iscooch(igrains)
c
      LOGICAL overlap
c
      do ic=1,igrain
       iscooch(ic) = 0  ! clear out the scooch array, this array identifies all grains to be scooched
      enddo
      iscooch(i)  = -1  !this is connected but should not be scooched
      iscooch(j)  =  1  !this is the first grain to check
c
      overlap = .true.
      do while ( overlap )
       do ii=1,igrain   ! identify grain to be be investigated
        if (iscooch(ii) .eq. 1 ) then
         ic = ii
         iscooch(ii) = 2  ! mark it as having been checked for additional connectivity
         goto 5
        endif
       enddo
       overlap = .false.
       goto 99
c
 5      ii0 = 0  !-1  !this was done in an effort to make the search smarter
        ii1 = 0  ! 1 ! and therefore faster
        jj0 = 0  !-1 
        jj1 = 0  ! 1
        kk0 = 0  !-1
        kk1 = 0  ! 1
        if ( xyzr(1,ic) .lt. 0.d0   + 3.d0*xyzr(4,ic) ) ii1 =   1
        if ( xyzr(1,ic) .gt. box(1) - 3.d0*xyzr(4,ic) ) ii0 =  -1
        if ( xyzr(2,ic) .lt. 0.d0   + 3.d0*xyzr(4,ic) ) jj1 =   1
        if ( xyzr(2,ic) .gt. box(2) - 3.d0*xyzr(4,ic) ) jj0 =  -1
        if ( xyzr(3,ic) .lt. 0.d0   + 3.d0*xyzr(4,ic) ) kk1 =   1
        if ( xyzr(3,ic) .gt. box(3) - 3.d0*xyzr(4,ic) ) kk0 =  -1
c        print *,'ii0,ii1,jj0,jj1,kk0,kk1',ii0,ii1,jj0,jj1,kk0,kk1  
        do ii = ii0 , ii1
        do jj = jj0 , jj1
        do kk = kk0 , kk1
         xshift = dfloat(ii)*box(1)
         yshift = dfloat(jj)*box(2)
         zshift = dfloat(kk)*box(3)
c         print *,'ii,jj,kk ',i,ii,jj,kk
         do jc=1, igrain
          r = DSQRT( ( (xyzr(1,ic)+xshift) - xyzr(1,jc) )**2 
     &            +  ( (xyzr(2,ic)+yshift) - xyzr(2,jc) )**2
     &            +  ( (xyzr(3,ic)+zshift) - xyzr(3,jc) )**2)
          if (xyzr(4,ic)+xyzr(4,jc) - r .ge. 0.d0 .and.
     &        jc .ne. ic .and. jc .ne. i .and. 
     &        iscooch(jc) .eq. 0) then
          icontact  = icontact  + 1
           iscooch(jc) = 1  ! thus jc is connected to ic and must itself be checked
c          print *,'scooch',i,j,ic,jc
          endif
         enddo  ! jc
         enddo  ! kk
         enddo  ! jj
         enddo  ! ii
        enddo   ! do while
c
 99   do ic = 1,igrain
       if (iscooch(ic) .eq. 2 .and. ic .ne. j) then
        xyzr(1,ic) =  xyzr(1,ic) - 1.d0*disp*xfac
        xyzr(2,ic) =  xyzr(2,ic) - 1.d0*disp*yfac 
        xyzr(3,ic) =  xyzr(3,ic) - 1.d0*disp*zfac
      endif
c keep particles in domain
         if (xyzr(1,ic) .lt. 0.00d0) xyzr(1,ic) = xyzr(1,ic) + box(1)
         if (xyzr(1,ic) .gt. box(1)) xyzr(1,ic) = xyzr(1,ic) - box(1)
         if (xyzr(2,ic) .lt. 0.00d0) xyzr(2,ic) = xyzr(2,ic) + box(2)
         if (xyzr(2,ic) .gt. box(2)) xyzr(2,ic) = xyzr(2,ic) - box(2)
         if (xyzr(3,ic) .lt. 0.00d0) xyzr(3,ic) = xyzr(3,ic) + box(3)
         if (xyzr(3,ic) .gt. box(3)) xyzr(3,ic) = xyzr(3,ic) - box(3)
      enddo
c      
      return
      end