      program               kmcbcc            ! Binary alloy (BCC mesh).

      implicit              none

      integer               l
      parameter            (l = 64)           ! Box-size: V = l^3.
      integer               ly_2
      parameter            (ly_2 = l / 2)     ! Alpha & Beta half-size.
      integer               isteps
      parameter            (isteps=20000000)  ! Number of MC steps.

      real*8                cb
      parameter            (cb = 0.40d0)      ! Fractional conc. of B.
      real*8                eps1
      parameter            (eps1 = 0.03d0)    ! Eps1 in units of "eV".
      real*8                eps1_T
      real*8                eps2
      parameter            (eps2 = -0.04d0)   ! Eps2 in units of "eV".
      real*8                eps2_T
      real*8                u1
      parameter            (u1 = eps2)        ! U1 in units of "eV".
      real*8                u1_T
      real*8                u2
      parameter            (u2 = 0.0d0)       ! U2 in units of "eV".
      real*8                u2_T
      real*8                k_b               ! Boltzman constant(eV/K).
      parameter            (k_b = 0.02585d0/300.0d0)
      real*8                temp              ! T in standard Kelvins(K)

      integer               loop              ! KMC loop index.
      integer               ageing            ! Ageing sequence loop.
      integer               n                 ! Total number of sites.
      integer               na                ! Number of A atoms.
      integer               nb                ! Number of B atoms.
      integer               mesh3d(l,l,l)     ! Simulation space.
      integer               alpha (l,ly_2,l)  ! Physical space (alpha).
      integer               beta  (l,ly_2,l)  ! Physical space (beta).
      integer               lmnvac(3)         ! Position of vacancy.

      integer               iseed             ! RNG seed (initialize -1)
      integer               unum              ! File unit number.

      real*8                eta, c_avg        ! Scalar order parameters.
      external              eta, c_avg
      external              init3d            ! Seeding routine.
      external              evolve            ! KMC step routine.
      external              lattic            ! Generates physical mesh.

c-----------------------------------------------------------------------
c                            Initialization: (seed lattice and ran3)
c-----------------------------------------------------------------------

      open (unit = 11, file = 'kmcbcc.out', status = 'new')

      call init3d(l,cb,n,na,nb,mesh3d,lmnvac,iseed)

      unum = 0

      call                   lattic(mesh3d,alpha,beta,l,ly_2,unum)

      write(11,*) 'eta = ',  eta(alpha,beta,l,ly_2,unum)

      write(11,*) 'cb  = ',c_avg(alpha,beta,l,ly_2,unum)

c-----------------------------------------------------------------------
c                            Warmup: (start system in B2 phase)
c-----------------------------------------------------------------------

      temp   = 800.0d0

      eps1_T = eps1 / (k_b * temp)           ! Scale energy parameters
      eps2_T = eps2 / (k_b * temp)           ! by temperature to produce
      u1_T   = u1   / (k_b * temp)           ! dimensionless inputs for
      u2_T   = u2   / (k_b * temp)           ! KMC evaluations...

      do 11 loop = 1, isteps
         call evolve(eps1_T,eps2_T,u1_T,u2_T,l,mesh3d,lmnvac,iseed)
11    continue

      unum = unum + 1

      call                   lattic(mesh3d,alpha,beta,l,ly_2,unum)

      write(11,*) 'eta = ',  eta(alpha,beta,l,ly_2,unum)

      write(11,*) 'cb  = ',c_avg(alpha,beta,l,ly_2,unum)

c-----------------------------------------------------------------------
c                            Warmup: (start system in B2 phase)
c-----------------------------------------------------------------------

      temp   = 750.0d0

      eps1_T = eps1 / (k_b * temp)           ! Scale energy parameters
      eps2_T = eps2 / (k_b * temp)           ! by temperature to produce
      u1_T   = u1   / (k_b * temp)           ! dimensionless inputs for
      u2_T   = u2   / (k_b * temp)           ! KMC evaluations...

      do 12 loop = 1, isteps
         call evolve(eps1_T,eps2_T,u1_T,u2_T,l,mesh3d,lmnvac,iseed)
12    continue

      unum = unum + 1

      call                   lattic(mesh3d,alpha,beta,l,ly_2,unum)

      write(11,*) 'eta = ',  eta(alpha,beta,l,ly_2,unum)

      write(11,*) 'cb  = ',c_avg(alpha,beta,l,ly_2,unum)

c-----------------------------------------------------------------------
c                            Warmup: (start system in B2 phase)
c-----------------------------------------------------------------------

      temp   = 700.0d0

      eps1_T = eps1 / (k_b * temp)           ! Scale energy parameters
      eps2_T = eps2 / (k_b * temp)           ! by temperature to produce
      u1_T   = u1   / (k_b * temp)           ! dimensionless inputs for
      u2_T   = u2   / (k_b * temp)           ! KMC evaluations...

      do 13 loop = 1, isteps
         call evolve(eps1_T,eps2_T,u1_T,u2_T,l,mesh3d,lmnvac,iseed)
13    continue

      unum = unum + 1

      call                   lattic(mesh3d,alpha,beta,l,ly_2,unum)

      write(11,*) 'eta = ',  eta(alpha,beta,l,ly_2,unum)

      write(11,*) 'cb  = ',c_avg(alpha,beta,l,ly_2,unum)

c-----------------------------------------------------------------------
c                            Warmup: (start system in B2 phase)
c-----------------------------------------------------------------------

      temp   = 650.0d0

      eps1_T = eps1 / (k_b * temp)           ! Scale energy parameters
      eps2_T = eps2 / (k_b * temp)           ! by temperature to produce
      u1_T   = u1   / (k_b * temp)           ! dimensionless inputs for
      u2_T   = u2   / (k_b * temp)           ! KMC evaluations...

      do 14 loop = 1, isteps
         call evolve(eps1_T,eps2_T,u1_T,u2_T,l,mesh3d,lmnvac,iseed)
14    continue

      unum = unum + 1

      call                   lattic(mesh3d,alpha,beta,l,ly_2,unum)

      write(11,*) 'eta = ',  eta(alpha,beta,l,ly_2,unum)

      write(11,*) 'cb  = ',c_avg(alpha,beta,l,ly_2,unum)

c-----------------------------------------------------------------------
c                           Main KMC Loop:
c-----------------------------------------------------------------------

      temp   = 348.0d0

      eps1_T = eps1 / (k_b * temp)           ! Scale energy parameters
      eps2_T = eps2 / (k_b * temp)           ! by temperature to produce
      u1_T   = u1   / (k_b * temp)           ! dimensionless inputs for
      u2_T   = u2   / (k_b * temp)           ! KMC evaluations...

      do 16 ageing = 1, 10
         do 15 loop = 1, isteps
            call evolve(eps1_T,eps2_T,u1_T,u2_T,l,mesh3d,lmnvac,iseed)
15       continue
         unum = unum + 1
         call                    lattic(mesh3d,alpha,beta,l,ly_2,unum)
         write(11,*) 'eta = ',   eta(alpha,beta,l,ly_2,unum)
         write(11,*) 'cb  = ', c_avg(alpha,beta,l,ly_2,unum)
16    continue

      temp   = 200.0d0

      eps1_T = eps1 / (k_b * temp)           ! Scale energy parameters
      eps2_T = eps2 / (k_b * temp)           ! by temperature to produce
      u1_T   = u1   / (k_b * temp)           ! dimensionless inputs for
      u2_T   = u2   / (k_b * temp)           ! KMC evaluations...

      do 18 ageing = 1, 10
         do 17 loop = 1, isteps
            call evolve(eps1_T,eps2_T,u1_T,u2_T,l,mesh3d,lmnvac,iseed)
17       continue
         unum = unum + 1
         call                    lattic(mesh3d,alpha,beta,l,ly_2,unum)
         write(11,*) 'eta = ',   eta(alpha,beta,l,ly_2,unum)
         write(11,*) 'cb  = ', c_avg(alpha,beta,l,ly_2,unum)
18    continue

      end file (unit = 11)

      close (unit = 11)

      stop
      end

c=======================================================================
c     subroutine            lattic
c
c    (returns the "physical" meshes of interest)
c
c     'mesh3d' was the original simulation cell, unphysically cubic:
c     _________________           ________________
c     \  rhombohedral  \         |     simple     |
c      \   simulation   \  ==>   |     cubic      |
c       \      cell      \       |    remapping   |    an l X l X l box.
c        \----------------\      |----------------| , 
c
c     'alpha' is the array of lattice points comprising half of the box.
c     'beta'  is the array of lattice points comprising the other half.
c     'l' is the lenth of these arrays in x and z directions.
c     'ly_2' is the half-as-long length along the y direction.
c     'unum' is the unit number of the formated file output.
c
c  !!!ALPHA AND BETA ARE HALF AS BIG AS MESH3D ALONG THE Y-DIRECTION!!!
c
c     To generate the complete simulation cell, 'alpha' and 'beta'
c     lattices must be added together as follows:
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c     In real "physical" space, the 'beta' lattice is offset from the
c     alpha lattice by a vector, (1/2,1/2,1/2).
c
c=======================================================================
      subroutine            lattic(mesh3d,alpha,beta,l,ly_2,unum)
      implicit              none
      integer               unum
      integer               i, j, k
      integer               l, ly_2
      integer               mesh3d(l,l,l)
      integer               alpha (l,ly_2,l)
      integer               beta  (l,ly_2,l)
      integer               temp_alpha (l,ly_2,l)
      integer               temp_beta  (l,ly_2,l)
      integer               pbc
      external              pbc

      do 13 k = 1, l
         do 12 j = 1, ly_2
            do 11 i = 1, l
               alpha(i,j,k) = 
     &                mesh3d(pbc(i+j-1,l),pbc(2*j-1,l),pbc(k+j-1,l))
               beta (i,j,k) = 
     &                mesh3d(pbc( i+j ,l),pbc( 2*j ,l),pbc( k+j ,l))
               if (alpha(i,j,k).GT.0) THEN
                  temp_alpha(i,j,k) = alpha(i,j,k)
               else
                  temp_alpha(i,j,k) = 2
               endif
               if (beta(i,j,k).GT.0) THEN
                  temp_beta(i,j,k) = beta(i,j,k)
               else
                  temp_beta(i,j,k) = 2
               endif
11          continue
12       continue
13    continue

      write (5000+unum*10+1,*) "Time = 1000000"
      write (5000+unum*10+1,*) "MATRIX SIZE (x,y,z)"
      write (5000+unum*10+1,*) l
      write (5000+unum*10+1,*) ly_2
      write (5000+unum*10+1,*) l
      write (5000+unum*10+1,*) "SPECIES"
      write (5000+unum*10+1,*) "Vac        0"
      write (5000+unum*10+1,*) "A          1"
      write (5000+unum*10+1,*) "B          2"

      write (5000+unum*10+2,*) "Time = 1000000"
      write (5000+unum*10+2,*) "MATRIX SIZE (x,y,z)"
      write (5000+unum*10+2,*) l
      write (5000+unum*10+2,*) ly_2
      write (5000+unum*10+2,*) l
      write (5000+unum*10+2,*) "SPECIES"
      write (5000+unum*10+2,*) "Vac        0"
      write (5000+unum*10+2,*) "A          1"
      write (5000+unum*10+2,*) "B          2"

      do k = 1,l
         do j = 1,ly_2
               write (5000+unum*10+1,*) (temp_alpha(i,j,k), i = 1,l)
         enddo
      enddo

      do k = 1,l
         do j = 1,ly_2
               write (5000+unum*10+2,*) (temp_beta(i,j,k), i = 1,l)
         enddo
      enddo
       
      return
      end

c=======================================================================
c     function              fcctyp
c
c    (tells whether or not you are on the fcc sublattice of a simple
c     cubic lattice which happens to include the point 1,1,1 at the
c     origin. Does that make sense? Think about it, and it will... )
c
c    (...or maybe it doesn't.)
c
c=======================================================================
      logical function      fcctyp(i,j,k)
      implicit              none
      integer               i, j, k
      if (mod(j,2) .eq. 0) then
         if (mod(i,2) .eq. 0) then
            if (mod(k,2) .eq. 0) then
               fcctyp = .false.
               return
            else
               fcctyp = .true.
               return
            endif
         else
            if (mod(k,2) .eq. 0) then
               fcctyp = .true.
               return
            else
               fcctyp = .false.
               return
            endif
         endif
      else
         if (mod(i,2) .eq. 0) then
            if (mod(k,2) .eq. 0) then
               fcctyp = .true.
               return
            else
               fcctyp = .false.
               return
            endif
         else
            if (mod(k,2) .eq. 0) then
               fcctyp = .false.
               return
            else
               fcctyp = .true.
               return
            endif
         endif
      endif
      end

c=======================================================================
c     function              n_k
c
c    (returns an atomic occupancy)
c    ('typmsh' = -1 iff we're talking about the alpha mesh.)
c    ('typmsh' = +1 iff we're talking about the beta  mesh.)
c
c=======================================================================
      integer function      n_k(i,j,k,typmsh,alpha,beta,l,ly_2)
      implicit              none
      integer               typmsh
      integer               i, j, k, l, ly_2
      integer               alpha(l,ly_2,l), beta(l,ly_2,l)
      if (typmsh .eq. -1) then
         if (alpha(i,j,k) .eq. 1) then
            n_k = 1
            return
         else
            n_k = 0
            return
         endif
      else
         if (beta(i,j,k) .eq. 1) then
            n_k = 1
            return
         else
            n_k = 0
            return
         endif
      endif
      end

c=======================================================================
c     function              s_k
c
c    (returns another atomic occupancy)
c    ('typmsh' = -1 iff we're talking about the alpha mesh.)
c    ('typmsh' = +1 iff we're talking about the beta  mesh.)
c
c=======================================================================
      integer function      s_k(i,j,k,typmsh,alpha,beta,l,ly_2)
      implicit              none
      integer               typmsh
      integer               i, j, k, l, ly_2
      integer               alpha(l,ly_2,l), beta(l,ly_2,l)
      if (typmsh .eq. -1) then
         if (alpha(i,j,k) .eq. -1) then
            s_k = 1
            return
         else
            s_k = -1
            return
         endif
      else
         if (beta(i,j,k) .eq. 1) then
            s_k = 1
            return
         else
            s_k = -1
            return
         endif
      endif
      end

c=======================================================================
c     subroutine            omega
c
c    (tells where all the nearest neighbors and next-nearest-neighbors
c     are located and which mesh they are on.)
c    ('typmsh' = -1 iff we're talking about the alpha mesh.)
c    ('typmsh' = +1 iff we're talking about the beta  mesh.)
c
c=======================================================================
      subroutine            omega(i,j,k,typmsh,omegai,l,ly_2)
      implicit              none
      integer               typmsh
      integer               i, j, k, l, ly_2
      integer               omegai(14,4)
      integer               pbc
      external              pbc

      omegai(1,1) = pbc(i,l)
      omegai(1,2) = pbc(j,ly_2)
      omegai(1,3) = pbc(k,l)
      omegai(1,4) = -1 * typmsh

      if (typmsh .eq. -1) then
         omegai(2,1) = pbc(i-1,l)
         omegai(2,2) = pbc(j,ly_2)
         omegai(2,3) = pbc(k,l)
         omegai(2,4) = 1
         omegai(3,1) = pbc(i,l)
         omegai(3,2) = pbc(j-1,ly_2)
         omegai(3,3) = pbc(k,l)
         omegai(3,4) = 1
         omegai(4,1) = pbc(i,l)
         omegai(4,2) = pbc(j,ly_2)
         omegai(4,3) = pbc(k-1,l)
         omegai(4,4) = 1
         omegai(5,1) = pbc(i-1,l)
         omegai(5,2) = pbc(j-1,ly_2)
         omegai(5,3) = pbc(k,l)
         omegai(5,4) = 1
         omegai(6,1) = pbc(i,l)
         omegai(6,2) = pbc(j-1,ly_2)
         omegai(6,3) = pbc(k-1,l)
         omegai(6,4) = 1
         omegai(7,1) = pbc(i-1,l)
         omegai(7,2) = pbc(j,ly_2)
         omegai(7,3) = pbc(k-1,l)
         omegai(7,4) = 1
         omegai(8,1) = pbc(i-1,l)
         omegai(8,2) = pbc(j-1,ly_2)
         omegai(8,3) = pbc(k-1,l)
         omegai(8,4) = 1
      else
         omegai(2,1) = pbc(i+1,l)
         omegai(2,2) = pbc(j,ly_2)
         omegai(2,3) = pbc(k,l)
         omegai(2,4) = -1
         omegai(3,1) = pbc(i,l)
         omegai(3,2) = pbc(j+1,ly_2)
         omegai(3,3) = pbc(k,l)
         omegai(3,4) = -1
         omegai(4,1) = pbc(i,l)
         omegai(4,2) = pbc(j,ly_2)
         omegai(4,3) = pbc(k+1,l)
         omegai(4,4) = -1
         omegai(5,1) = pbc(i+1,l)
         omegai(5,2) = pbc(j+1,ly_2)
         omegai(5,3) = pbc(k,l)
         omegai(5,4) = -1
         omegai(6,1) = pbc(i,l)
         omegai(6,2) = pbc(j+1,ly_2)
         omegai(6,3) = pbc(k+1,l)
         omegai(6,4) = -1
         omegai(7,1) = pbc(i+1,l)
         omegai(7,2) = pbc(j,ly_2)
         omegai(7,3) = pbc(k+1,l)
         omegai(7,4) = -1
         omegai(8,1) = pbc(i+1,l)
         omegai(8,2) = pbc(j+1,ly_2)
         omegai(8,3) = pbc(k+1,l)
         omegai(8,4) = -1
      endif

      omegai(9,1)  = pbc(i+1,l)
      omegai(9,2)  = pbc(j,ly_2)
      omegai(9,3)  = pbc(k,l)
      omegai(9,4)  = typmsh
      omegai(10,1) = pbc(i-1,l)
      omegai(10,2) = pbc(j,ly_2)
      omegai(10,3) = pbc(k,l)
      omegai(10,4) = typmsh
      omegai(11,1) = pbc(i,l)
      omegai(11,2) = pbc(j+1,ly_2)
      omegai(11,3) = pbc(k,l)
      omegai(11,4) = typmsh
      omegai(12,1) = pbc(i,l)
      omegai(12,2) = pbc(j-1,ly_2)
      omegai(12,3) = pbc(k,l)
      omegai(12,4) = typmsh
      omegai(13,1) = pbc(i,l)
      omegai(13,2) = pbc(j,ly_2)
      omegai(13,3) = pbc(k+1,l)
      omegai(13,4) = typmsh
      omegai(14,1) = pbc(i,l)
      omegai(14,2) = pbc(j,ly_2)
      omegai(14,3) = pbc(k-1,l)
      omegai(14,4) = typmsh

      return
      end

c=======================================================================
c     function              eta
c
c    (returns the B2 order parameter, averaged over all sites.)
c
c     'alpha' is the array of lattice points comprising half of the box.
c     'beta'  is the array of lattice points comprising the other half.
c     'l' is the lenth of these arrays in x and z directions.
c     'ly_2' is the half-as-long length along the y direction.
c     'unum' is the unit number of the formated file output.
c
c     To generate the complete simulation cell, 'alpha' and 'beta'
c     lattices must be added together as follows:
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c     In real "physical" space, the 'beta' lattice is offset from the
c     alpha lattice by a vector, (1/2,1/2,1/2).
c
c=======================================================================
      real*8 function       eta(alpha,beta,l,ly_2,unum)
      implicit              none
      integer               unum
      integer               total
      integer               i, j, k, l, ly_2
      integer               alpha(l,ly_2,l)
      integer               beta (l,ly_2,l)
      real*8                temper
      real*8                temper_1
      real*8                eta_i
      external              eta_i

      eta = 0.0d0
      temper = 0.0d0
      temper_1 = 0.0d0

      total = 0

      do 13 k = 1, l
         do 12 j = 1, ly_2
            do 11 i = 1, l
               total = total + 1
               temper = eta_i(i,j,k,-1,alpha,beta,l,ly_2)
               eta = eta + temper
               temper_1 = (temper*temper)
               if ((temper_1.gt.0.5).and.(j.eq.8)) then
                  write (5000+unum*10+3,*) i , k
               endif
11          continue
12       continue
13    continue

      do 16 k = 1, l
         do 15 j = 1, ly_2
            do 14 i = 1, l
               total = total + 1
               eta = eta + eta_i(i,j,k, 1,alpha,beta,l,ly_2)
14          continue
15       continue
16    continue

      eta = eta / dfloat(total)

      return
      end

c=======================================================================
c     function              eta_i (Equation 10.)
c
c    (returns the B2 order field as a function of position)
c
c    ('position' is specified by 4 coordinates: i, j, k, typmsh)
c
c     'i' is the index of position along the x-direction of a lattice.
c     'j' is the index of position along the y-direction of a lattice.
c     'k' is the index of position along the z-direction of a lattice.
c     'typmsh' is the mesh coordinate. (-1 = A mesh, +1 = B mesh)
c
c     'alpha' is the array of lattice points comprising half of the box.
c     'beta'  is the array of lattice points comprising the other half.
c     'l' is the lenth of these arrays in x and z directions.
c     'ly_2' is the half-as-long length along the y direction.
c
c     To generate the complete simulation cell, 'alpha' and 'beta'
c     lattices must be added together as follows:
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c     In real "physical" space, the 'beta' lattice is offset from the
c     alpha lattice by a vector, (1/2,1/2,1/2).
c
c     uses: s_k, omega, pbc
c
c=======================================================================
      real*8 function       eta_i(i,j,k,typmsh,alpha,beta,l,ly_2)
      implicit              none
      integer               i, j, k, l, ly_2, typmsh
      integer               n, total
      integer               omegai(14,4)
      integer               alpha(l,ly_2,l), beta(l,ly_2,l)
      integer               s_k
      external              s_k, omega

      call                  omega(i,j,k,typmsh,omegai,l,ly_2)

      eta_i = 0.0d0

      total  = 0

      do 11 n = 1, 14
         total = total + 1
         eta_i = eta_i + dfloat(s_k(
     &                  omegai(n,1),omegai(n,2),omegai(n,3),omegai(n,4),
     &                  alpha,beta,l,ly_2))
11    continue

      eta_i = eta_i / dfloat(total)

      return
      end

c=======================================================================
c     function              c_avg
c
c    (returns the B2 concentration parameter, averaged over all sites.)
c
c     'alpha' is the array of lattice points comprising half of the box.
c     'beta'  is the array of lattice points comprising the other half.
c     'l' is the lenth of these arrays in x and z directions.
c     'ly_2' is the half-as-long length along the y direction.
c     'unum' is the unit number of the formated file output.
c
c     To generate the complete simulation cell, 'alpha' and 'beta'
c     lattices must be added together as follows:
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c     In real "physical" space, the 'beta' lattice is offset from the
c     alpha lattice by a vector, (1/2,1/2,1/2).
c
c=======================================================================
      real*8 function       c_avg(alpha,beta,l,ly_2,unum)
      implicit              none
      integer               unum
      integer               total
      integer               i, j, k, l, ly_2
      integer               alpha(l,ly_2,l)
      integer               beta (l,ly_2,l)
      real*8                temper
      real*8                c_i
      external              c_i

      c_avg = 0.0d0
      temper = 0.0d0

      total = 0

      do 13 k = 1, l
         do 12 j = 1, ly_2
            do 11 i = 1, l
               total = total + 1
               temper = c_i(i,j,k,-1,alpha,beta,l,ly_2)
               c_avg = c_avg + temper
11         continue
12       continue
13    continue

      do 16 k = 1, l
         do 15 j = 1, ly_2
            do 14 i = 1, l
               total = total + 1
               c_avg = c_avg + c_i(i,j,k, 1,alpha,beta,l,ly_2)
14          continue
15       continue
16    continue

      c_avg = c_avg / dfloat(total)

      return
      end

c=======================================================================
c     function              c_i (Equation 10.)
c
c    (returns the B2 concentration field as a function of position)
c
c    ('position' is specified by 4 coordinates: i, j, k, typmsh)
c
c     'i' is the index of position along the x-direction of a lattice.
c     'j' is the index of position along the y-direction of a lattice.
c     'k' is the index of position along the z-direction of a lattice.
c     'typmsh' is the mesh coordinate. (-1 = A mesh, +1 = B mesh)
c
c     'alpha' is the array of lattice points comprising half of the box.
c     'beta'  is the array of lattice points comprising the other half.
c     'l' is the lenth of these arrays in x and z directions.
c     'ly_2' is the half-as-long length along the y direction.
c
c     To generate the complete simulation cell, 'alpha' and 'beta'
c     lattices must be added together as follows:
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c         B       B       B       B       B       B       B       B
c
c     A       A       A       A       A       A       A       A
c
c     In real "physical" space, the 'beta' lattice is offset from the
c     alpha lattice by a vector, (1/2,1/2,1/2).
c
c     uses: n_k, omega, pbc
c
c=======================================================================
      real*8 function       c_i(i,j,k,typmsh,alpha,beta,l,ly_2)
      implicit              none
      integer               i, j, k, l, ly_2, typmsh
      integer               n, total
      integer               omegai(14,4)
      integer               alpha(l,ly_2,l), beta(l,ly_2,l)
      integer               n_k
      external              n_k, omega

      call                  omega(i,j,k,typmsh,omegai,l,ly_2)

      c_i = 0.0d0

      total  = 0

      do 11 n = 1, 14
         total = total + 1
         c_i = c_i + dfloat(n_k(
     &               omegai(n,1),omegai(n,2),omegai(n,3),omegai(n,4),
     &               alpha,beta,l,ly_2))
11    continue

      c_i = c_i / dfloat(total)

      return
      end

c=======================================================================
c     subroutine            init3d
c
c    (randomly populate the 3d mesh with A atoms, B atoms and a vacancy)
c    (calculates the initial position of the vacancy, "lmnvac")
c    (initializes the time and the random-number-generator seed.)
c
c=======================================================================
      subroutine            init3d(l,cb,n,na,nb,mesh3d,lmnvac,iseed)
      implicit              none
      integer               i, j, k, l
      integer               i_ran, j_ran, k_ran, ncount
      integer               n
      integer               na
      integer               nb
      integer               mesh3d(l,l,l)
      integer               lmnvac(3)
      integer               iseed
      real*8                ca, cb
      real*8                ran3
      external              ran3
      iseed = -137
c     rseed = 98693.0d0
      n = l**3
      ca = 1.0d0 - cb
      na = max0(min0(nint(ca*(n-1)),n-1),0)
      nb = n - na - 1
      if (na .ge. nb) then
         do 13 k = 1, l
            do 12 j = 1, l
               do 11 i = 1, l
                  mesh3d(i,j,k) = -1
11             continue
12          continue
13       continue
         if (nb .eq. 0) goto 3
      else
         do 16 k = 1, l
            do 15 j = 1, l
               do 14 i = 1, l
                  mesh3d(i,j,k) = 1
14             continue
15          continue
16       continue
         if (na .eq. 0) goto 3
      endif

      ncount = 0

      if (na .ge. nb) then
1        continue
         i_ran = 1 + nint((l-1)*ran3(iseed))
         j_ran = 1 + nint((l-1)*ran3(iseed))
         k_ran = 1 + nint((l-1)*ran3(iseed))
         if (mesh3d(i_ran,j_ran,k_ran) .eq. -1) then
            mesh3d(i_ran,j_ran,k_ran) = 1
            ncount = ncount + 1
         endif
         if (ncount .eq. nb) goto 3
         goto 1
      else
2        continue
         i_ran = 1 + nint((l-1)*ran3(iseed))
         j_ran = 1 + nint((l-1)*ran3(iseed))
         k_ran = 1 + nint((l-1)*ran3(iseed))
         if (mesh3d(i_ran,j_ran,k_ran) .eq. 1) then
            mesh3d(i_ran,j_ran,k_ran) = -1
            ncount = ncount + 1
         endif
         if (ncount .eq. na) goto 3
         goto 2
      endif

3     continue

      if (na .ge. nb) then
4        continue
         i_ran = 1 + nint((l-1)*ran3(iseed))
         j_ran = 1 + nint((l-1)*ran3(iseed))
         k_ran = 1 + nint((l-1)*ran3(iseed))
         if (mesh3d(i_ran,j_ran,k_ran) .eq. -1) then
            mesh3d(i_ran,j_ran,k_ran) = 0
            lmnvac(1) = i_ran
            lmnvac(2) = j_ran
            lmnvac(3) = k_ran
            goto 6
         endif
         goto 4
      else
5        continue
         i_ran = 1 + nint((l-1)*ran3(iseed))
         j_ran = 1 + nint((l-1)*ran3(iseed))
         k_ran = 1 + nint((l-1)*ran3(iseed))
         if (mesh3d(i_ran,j_ran,k_ran) .eq. 1) then
            mesh3d(i_ran,j_ran,k_ran) = 0
            lmnvac(1) = i_ran
            lmnvac(2) = j_ran
            lmnvac(3) = k_ran
            goto 6
         endif
         goto 5
      endif

6     continue

      return
      end

c=======================================================================
c     subroutine            evolve
c
c    (performs a single Kinetic Monte Carlo step)
c
c=======================================================================
      subroutine            evolve(eps1,eps2,u1,u2,l,
     &                             mesh3d,lmnvac,iseed)
      implicit              none
      integer               j, l
      integer               mesh3d(l,l,l)
      integer               lmnvac(3)
      integer               swap
      integer               pbc
      integer               iseed
      real*8                eps1, eps2, u1, u2
      real*8                ttmp, rantmp, sumgam
      real*8                gam_nn(8)
      real*8                tau, ran3
      external              tau, ran3, pbc, gamslv

      call gamslv(gam_nn,eps1,eps2,u1,u2,mesh3d,l,
     &            lmnvac(1),lmnvac(2),lmnvac(3))

      ttmp = tau(gam_nn)

      rantmp = ran3(iseed) / ttmp

      sumgam = 0.0d0

      j = 1

      if (rantmp .lt. gam_nn(1)) goto 1

      do 11 j = 2, 8
         sumgam = sumgam + gam_nn(j-1)
         if ((rantmp .ge. sumgam) .and. 
     &       (rantmp .lt. sumgam+gam_nn(j))) goto 1
11    continue

      j = 8

1     continue

      if (j .eq. 1) then
       swap = 
     & mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2)+1,l),pbc(lmnvac(3),l))
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2)+1,l),pbc(lmnvac(3),l))
     & = 0
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2),l),pbc(lmnvac(3),l))=swap
       lmnvac(1) = pbc(lmnvac(1),l)
       lmnvac(2) = pbc(lmnvac(2)+1,l)
       lmnvac(3) = pbc(lmnvac(3),l)
      else if (j .eq. 2) then
       swap = 
     & mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2)+1,l),pbc(lmnvac(3)+1,l))
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2)+1,l),pbc(lmnvac(3)+1,l))
     & = 0
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2),l),pbc(lmnvac(3),l))=swap
       lmnvac(1) = pbc(lmnvac(1),l)
       lmnvac(2) = pbc(lmnvac(2)+1,l)
       lmnvac(3) = pbc(lmnvac(3)+1,l)
      else if (j .eq. 3) then
       swap = 
     & mesh3d(pbc(lmnvac(1)+1,l),pbc(lmnvac(2)+1,l),pbc(lmnvac(3)+1,l))
       mesh3d(pbc(lmnvac(1)+1,l),pbc(lmnvac(2)+1,l),pbc(lmnvac(3)+1,l))
     & = 0
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2),l),pbc(lmnvac(3),l))=swap
       lmnvac(1) = pbc(lmnvac(1)+1,l)
       lmnvac(2) = pbc(lmnvac(2)+1,l)
       lmnvac(3) = pbc(lmnvac(3)+1,l)
      else if (j .eq. 4) then
       swap = 
     & mesh3d(pbc(lmnvac(1)+1,l),pbc(lmnvac(2)+1,l),pbc(lmnvac(3),l))
       mesh3d(pbc(lmnvac(1)+1,l),pbc(lmnvac(2)+1,l),pbc(lmnvac(3),l))
     & = 0
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2),l),pbc(lmnvac(3),l))=swap
       lmnvac(1) = pbc(lmnvac(1)+1,l)
       lmnvac(2) = pbc(lmnvac(2)+1,l)
       lmnvac(3) = pbc(lmnvac(3),l)
      else if (j .eq. 5) then
       swap = 
     & mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2)-1,l),pbc(lmnvac(3),l))
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2)-1,l),pbc(lmnvac(3),l))
     & = 0
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2),l),pbc(lmnvac(3),l))=swap
       lmnvac(1) = pbc(lmnvac(1),l)
       lmnvac(2) = pbc(lmnvac(2)-1,l)
       lmnvac(3) = pbc(lmnvac(3),l)
      else if (j .eq. 6) then
       swap = 
     & mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2)-1,l),pbc(lmnvac(3)-1,l))
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2)-1,l),pbc(lmnvac(3)-1,l))
     & = 0
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2),l),pbc(lmnvac(3),l))=swap
       lmnvac(1) = pbc(lmnvac(1),l)
       lmnvac(2) = pbc(lmnvac(2)-1,l)
       lmnvac(3) = pbc(lmnvac(3)-1,l)
      else if (j .eq. 7) then
       swap = 
     & mesh3d(pbc(lmnvac(1)-1,l),pbc(lmnvac(2)-1,l),pbc(lmnvac(3)-1,l))
       mesh3d(pbc(lmnvac(1)-1,l),pbc(lmnvac(2)-1,l),pbc(lmnvac(3)-1,l))
     & = 0
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2),l),pbc(lmnvac(3),l))=swap
       lmnvac(1) = pbc(lmnvac(1)-1,l)
       lmnvac(2) = pbc(lmnvac(2)-1,l)
       lmnvac(3) = pbc(lmnvac(3)-1,l)
      else if (j .eq. 8) then
       swap = 
     & mesh3d(pbc(lmnvac(1)-1,l),pbc(lmnvac(2)-1,l),pbc(lmnvac(3),l))
       mesh3d(pbc(lmnvac(1)-1,l),pbc(lmnvac(2)-1,l),pbc(lmnvac(3),l))
     & = 0
       mesh3d(pbc(lmnvac(1),l),pbc(lmnvac(2),l),pbc(lmnvac(3),l))=swap
       lmnvac(1) = pbc(lmnvac(1)-1,l)
       lmnvac(2) = pbc(lmnvac(2)-1,l)
       lmnvac(3) = pbc(lmnvac(3),l)
      endif

      return
      end

c=======================================================================
c     subroutine            gamslv
c
c    (solves for the 8 nearest-neighbor transition-amplitudes)
c
c=======================================================================
      subroutine            gamslv(gam_nn,eps1,eps2,u1,u2,
     &                             mesh3d,l,lx,mx,nx)
      implicit              none
      integer               l, lx, mx, nx
      integer               mesh3d(l,l,l)
      integer               lmnvac(3)
      real*8                eps1, eps2, u1, u2
      real*8                gam_nn(8)
      integer               neigh1, neigh2, pbc
      real*8                actbcc
      external              pbc, actbcc, neigh1, neigh2

      gam_nn(1)=exp(
     &        -actbcc(mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx,l)),
     &                eps1,eps2,u1,u2,
     &                neigh1(lx,mx+1,nx,
     &                mesh3d,l,
     &                mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx,l))),
     &                neigh2(lx,mx+1,nx,
     &                mesh3d,l,
     &                mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx,l))),
     &                neigh1(lx,mx+1,nx,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx,l))),
     &                neigh2(lx,mx+1,nx,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx,l))))
     &                                                                 )

      gam_nn(2)=exp(
     &        -actbcc(mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx+1,l)),
     &                eps1,eps2,u1,u2,
     &                neigh1(lx,mx+1,nx+1,
     &                mesh3d,l,
     &                mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx+1,l))),
     &                neigh2(lx,mx+1,nx+1,
     &                mesh3d,l,
     &                mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx+1,l))),
     &                neigh1(lx,mx+1,nx+1,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx+1,l))),
     &                neigh2(lx,mx+1,nx+1,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx+1,l))))
     &                                                                 )

      gam_nn(3)=exp(
     &        -actbcc(mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx+1,l)),
     &                eps1,eps2,u1,u2,
     &                neigh1(lx+1,mx+1,nx+1,
     &                mesh3d,l,
     &                mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx+1,l))),
     &                neigh2(lx+1,mx+1,nx+1,
     &                mesh3d,l,
     &                mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx+1,l))),
     &                neigh1(lx+1,mx+1,nx+1,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx+1,l))),
     &                neigh2(lx+1,mx+1,nx+1,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx+1,l))))
     &                                                                 )

      gam_nn(4)=exp(
     &        -actbcc(mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx,l)),
     &                eps1,eps2,u1,u2,
     &                neigh1(lx+1,mx+1,nx,
     &                mesh3d,l,
     &                mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx,l))),
     &                neigh2(lx+1,mx+1,nx,
     &                mesh3d,l,
     &                mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx,l))),
     &                neigh1(lx+1,mx+1,nx,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx,l))),
     &                neigh2(lx+1,mx+1,nx,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx,l))))
     &                                                                 )

      gam_nn(5)=exp(
     &        -actbcc(mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx,l)),
     &                eps1,eps2,u1,u2,
     &                neigh1(lx,mx-1,nx,
     &                mesh3d,l,
     &                mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx,l))),
     &                neigh2(lx,mx-1,nx,
     &                mesh3d,l,
     &                mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx,l))),
     &                neigh1(lx,mx-1,nx,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx,l))),
     &                neigh2(lx,mx-1,nx,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx,l))))
     &                                                                 )

      gam_nn(6)=exp(
     &        -actbcc(mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx-1,l)),
     &                eps1,eps2,u1,u2,
     &                neigh1(lx,mx-1,nx-1,
     &                mesh3d,l,
     &                mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx-1,l))),
     &                neigh2(lx,mx-1,nx-1,
     &                mesh3d,l,
     &                mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx-1,l))),
     &                neigh1(lx,mx-1,nx-1,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx-1,l))),
     &                neigh2(lx,mx-1,nx-1,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx-1,l))))
     &                                                                 )

      gam_nn(7)=exp(
     &        -actbcc(mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx-1,l)),
     &                eps1,eps2,u1,u2,
     &                neigh1(lx-1,mx-1,nx-1,
     &                mesh3d,l,
     &                mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx-1,l))),
     &                neigh2(lx-1,mx-1,nx-1,
     &                mesh3d,l,
     &                mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx-1,l))),
     &                neigh1(lx-1,mx-1,nx-1,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx-1,l))),
     &                neigh2(lx-1,mx-1,nx-1,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx-1,l))))
     &                                                                 )

      gam_nn(8)=exp(
     &        -actbcc(mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx,l)),
     &                eps1,eps2,u1,u2,
     &                neigh1(lx-1,mx-1,nx,
     &                mesh3d,l,
     &                mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx,l))),
     &                neigh2(lx-1,mx-1,nx,
     &                mesh3d,l,
     &                mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx,l))),
     &                neigh1(lx-1,mx-1,nx,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx,l))),
     &                neigh2(lx-1,mx-1,nx,
     &                mesh3d,l,
     &               -mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx,l))))
     &                                                                 )

      return
      end

c=======================================================================
c     function              tau
c
c    (calculates the average residence time for the current config)
c
c=======================================================================
      real*8 function       tau(gam_nn)
      implicit              none
      integer               i
      real*8                gam_nn(8)
      tau = 0.0d0
      do 11 i = 1, 8
         tau = tau + gam_nn(i)
11    continue
      tau = 1.0d0 / max(abs(tau),1.0d-36)
      return
      end

c=======================================================================
c     function              actbcc
c
c    (calculates the activation energy)
c
c=======================================================================
      real*8 function       actbcc(typx,eps1,eps2,u1,u2,n1x,n2x,n1y,n2y)
      implicit              none
      integer               typx
      integer               n1x, n2x, n1y, n2y
      real*8                eps1, eps2
      real*8                u1, u2
      actbcc = 0.5d0 * (eps1*n1y + eps2*n2y + typx*(u1*n1x + u2*n2x))
      return
      end

c=======================================================================
c     function              neigh1
c
c    (calculates the number of first nearest neighbors of type x.)
c
c=======================================================================
      integer function      neigh1(lx,mx,nx,mesh3d,l,typx)
      implicit              none
      integer               typx
      integer               l
      integer               mesh3d(l,l,l)
      integer               lx, mx, nx
      integer               pbc
      external              pbc
      neigh1 = 
     &         mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx,l)) + 
     &            
     &         mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx,l)) + 
     &            
     &         mesh3d(pbc(lx,l),pbc(mx+1,l),pbc(nx+1,l)) + 
     &            
     &         mesh3d(pbc(lx+1,l),pbc(mx+1,l),pbc(nx+1,l)) + 
     &            
     &         mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx,l)) + 
     &            
     &         mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx,l)) + 
     &            
     &         mesh3d(pbc(lx,l),pbc(mx-1,l),pbc(nx-1,l)) + 
     &            
     &         mesh3d(pbc(lx-1,l),pbc(mx-1,l),pbc(nx-1,l))

      neigh1 = (8 + neigh1*typx - mod(neigh1,2)) / 2

      return
      end

c=======================================================================
c     function              neigh2
c
c    (calculates the number of second nearest neighbors of type x.)
c
c=======================================================================
      integer function      neigh2(lx,mx,nx,mesh3d,l,typx)
      implicit              none
      integer               typx
      integer               l
      integer               mesh3d(l,l,l)
      integer               lx, mx, nx
      integer               pbc
      external              pbc
      neigh2 = 
     &         mesh3d(pbc(lx,l),pbc(mx,l),pbc(nx+1,l)) + 
     &            
     &         mesh3d(pbc(lx,l),pbc(mx,l),pbc(nx-1,l)) + 
     &            
     &         mesh3d(pbc(lx+1,l),pbc(mx,l),pbc(nx,l)) + 
     &            
     &         mesh3d(pbc(lx-1,l),pbc(mx,l),pbc(nx,l)) + 
     &            
     &         mesh3d(pbc(lx+1,l),pbc(mx+2,l),pbc(nx+1,l)) + 
     &            
     &         mesh3d(pbc(lx-1,l),pbc(mx-2,l),pbc(nx-1,l))

      neigh2 = (6 + neigh2*typx - mod(neigh2,2)) / 2

      return
      end

c=======================================================================
c     function              pbc
c
c    (takes care of the periodic boundary conditions)
c
c=======================================================================
      integer function      pbc(nx,l)
      implicit              none
      integer               nx, l
      pbc = mod(nx+l-1,l)+1
      return
      end

c=======================================================================
c     function              rantim
c
c    (general function which calls one of the two following rng funcs.)
c
c=======================================================================
      real*8 function       rantim(iseed,rseed)
      implicit              none
      logical               switch
      parameter            (switch = .true.)
      integer               iseed
      real*8                rseed
      real*8                ran3, rand
      external              ran3, rand
      if (switch) then
         rantim = ran3(iseed)
         return
      else
         rantim = rand(rseed)
         return
      endif
      end

c=======================================================================
c     This is a multiplicative congruential random number generator
c     with modulus 2**31-1
c     Taken from the book "Monte Carlo: Concepts, Algorithms, and 
c     Applications" by George S. Fishman  P. 604
c
c     N = # of random numbers to be generated, <=NDIM
c     SEED = initial positive seed on entering, final seed on exit
c     DUNIF = vector of uniform (0,1) random numbers on exit
c     NDIM = dimension of DUNIF
c=======================================================================
      REAL*8 FUNCTION RAND(SEED)

      REAL*8 SEED
      REAL*8 DZ,DOVER,DZ1,DZ2,DOVER1,DOVER2
      REAL*8 DTWO31,DMDLS,DA1,DA2
      REAL*8 Q,T,XP,NUM,DEN

c      PRINT *, SEED,RAND
c      OPEN (UNIT=20,FILE='RNG_OUT',STATUS='UNKNOWN',ACCESS='APPEND')
      DATA DTWO31/2147483648D+0/,DMDLS/2147483647D+0/
      DATA DA1/41160.0D+0/,DA2/950665216.0D+0/

      DZ=SEED
         DZ=IDINT(DZ)
         DZ1=DZ*DA1
         DZ2=DZ*DA2
         DOVER1=IDINT(DZ1/DTWO31)
         DOVER2=IDINT(DZ2/DTWO31)
         DZ1=DZ1-DOVER1*DTWO31
         DZ2=DZ2-DOVER2*DTWO31
         DZ=DZ1+DZ2+DOVER1+DOVER2
         DOVER=IDINT(DZ/DMDLS)
         DZ=DZ-DOVER*DMDLS
         RAND=DZ/DMDLS
         SEED=DZ
c      WRITE(20,*) SEED,RAND
c      CLOSE(UNIT=20)
      RETURN
      END

c=======================================================================
c     Knuth random-number-generator from Numerical Recipes, section 7.1
c=======================================================================
      real*8 function       ran3(idum)
      implicit              none
      integer               idum
      integer               mbig, mseed, mz
      real*8                fac
      parameter            (mbig=1000000000,mseed=161803398,mz=0,
     &                      fac=1.0d0/mbig)
      integer               i, iff, ii, inext, inextp, k
      integer               mj, mk, ma(55)
      save                  iff, inext, inextp, ma
      data iff /0/ 
      if(idum.lt.0.or.iff.eq.0)then 
        iff=1 
        mj=mseed-iabs(idum) 
        mj=mod(mj,mbig) 
        ma(55)=mj 
        mk=1 
        do 11 i=1,54 
          ii=mod(21*i,55) 
          ma(ii)=mk 
          mk=mj-mk 
          if(mk.lt.mz)mk=mk+mbig 
          mj=ma(ii) 
11      continue 
        do 13 k=1,4 
          do 12 i=1,55 
            ma(i)=ma(i)-ma(1+mod(i+30,55)) 
            if(ma(i).lt.mz)ma(i)=ma(i)+mbig 
12        continue 
13      continue 
        inext=0 
        inextp=31 
        idum=1 
      endif 
      inext=inext+1 
      if(inext.eq.56)inext=1 
      inextp=inextp+1 
      if(inextp.eq.56)inextp=1 
      mj=ma(inext)-ma(inextp) 
      if(mj.lt.mz)mj=mj+mbig 
      ma(inext)=mj 
      ran3=mj*fac 
      return 
      end
