CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                  C
C     Make a randomly connected triangulated mesh                  C 
C                                                                  C
C     Modified to be a subroutine by Chris Martin 7/24/93          C
C                                                                  C
C     Modified by PDC 11/11/93 so as not to re-initialize the      C
C     mesh every time, but rather thermalize at the beginning      C
C     and then continue evolving the mesh.                         C
C                                                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE make_mesh

      CALL flip_many_links

      CALL produce_mesh
      
      RETURN
      END

C----------------------------------------------------------------------C

      SUBROUTINE init(ntherm)

      INCLUDE 'irreg.h'
      INCLUDE 'makemesh.h'

      INTEGER i,ntherm 

C
C --- Calculate all inverse numbers
C
      inv_n = 1.0/float(n)
      inv_3n = 1.0/float(3*n)
      inv_n6 = 1.0/float(n6)
      third_two_pi = two_pi/3.0

C
C --- Set up the initial mesh.
C
      CALL set_up_mesh

C
C --- Thermalize the mesh.
C
      DO i=1,ntherm
         CALL flip_many_links
      ENDDO

      RETURN
      END

C----------------------------------------------------------------------C


      SUBROUTINE set_up_mesh

      INCLUDE 'irreg.h'
      INCLUDE 'makemesh.h'

      INTEGER site, site1, s0, s2, m_start, m_start2
      integer nc,nc1,i,j
      
      DO site = 0 , nm
        m_start_s(site) = 6 * site
        n_connectivity_s(site) = 6
c$$$        owner_of_s(site) = (site/nprocs)
      ENDDO


      DO site = 0 , nm
        s_neighbor_m(6*site+0) = site + 1
        s_neighbor_m(6*site+1) = site + nx + 1
        s_neighbor_m(6*site+2) = site + nx
        s_neighbor_m(6*site+3) = site - 1
        s_neighbor_m(6*site+4) = site - nx - 1
        s_neighbor_m(6*site+5) = site - nx

        DO site1 = 0, 5
          m_next_m  (6*site+site1) = 6*site + site1 + 1
          m_former_m(6*site+site1) = 6*site + site1 - 1
          s_site_m  (6*site+site1) = site
        ENDDO
        m_next_m   ( 6*site + 5 ) = 6*site
        m_former_m ( 6*site     ) = 6*site + 5

        IF(site .LT. nx)THEN
          s_neighbor_m(6*site+5) = site + nx*(ny-1)
          s_neighbor_m(6*site+4) = site + nx*(ny-1) - 1
        ENDIF
        IF(nx*ny-1-site .lt. nx)then
          s_neighbor_m(6*site+2) = site - nx*(ny-1)
          s_neighbor_m(6*site+1) = site - nx*(ny-1) + 1
        ENDIF
        IF(mod(site,nx) .eq. nx-1)then
          s_neighbor_m(6*site+0) = site - nx + 1
          s_neighbor_m(6*site+1) = site + 1
        ENDIF
        IF(mod(site,nx) .eq. 0)then
          s_neighbor_m(6*site+3) = site + nx - 1
          s_neighbor_m(6*site+4) = site - 1
        ENDIF
      ENDDO

      s_neighbor_m(0+4) = nx * ny - 1
      s_neighbor_m(6*(nx*ny-1)+1) = 0

      DO s0=0, nm
        m_start = m_start_s(s0)
        nc = n_connectivity_s(s0)
        do i=1,nc
           s2 = s_neighbor_m(m_start)
           m_start2 = m_start_s(s2)
           nc1 = n_connectivity_s(s2)
           do j=1,nc1
              IF(s0.EQ.s_neighbor_m(m_start2))GOTO 52
              m_start2 = m_next_m(m_start2)
           enddo
 52        CONTINUE

CC    here m_start and m_start2 are the two opposite links.
           m_back_m(m_start ) = m_start2
           m_back_m(m_start2) = m_start
           m_start = m_next_m(m_start)
        enddo
      ENDDO

      RETURN
      END

C----------------------------------------------------------------------C

      SUBROUTINE produce_mesh

      INCLUDE 'irreg.h'
      INCLUDE 'makemesh.h'
      
      integer maxverts, nc
      parameter (maxverts = 40)
      integer tmp_s(maxverts)
      integer i,j,l
      integer m_start
      integer temp_owner,temp1
      integer mesh_n(0:nm)
      integer n_index(0:n6m)
      common /makesend/ mesh_n,n_index


CC    produce item 1.
c$$$      print *, 'number of node = ', n
c$$$      print *, 'number of triangular element = ', 2*n
c$$$
CC    produce item 2.
c$$$      print *,'node         its_coordinates'
c$$$      do i = 0, nm
c$$$         print *, i,'   ',(field(j,i),j=0,dm)
c$$$      enddo
c$$$
c$$$      do i =0, nm
c$$$         do j=0,dm
c$$$            coord_of_s(j,i)=field(j,i)
c$$$         enddo
c$$$      enddo
c$$$
c$$$      print *,'node         owner'
c$$$      do i = 0, nm
c$$$         print *, i,'   ',owner_of_s(i)
c$$$         info_of_s(0,i)=owner_of_s(i) ! owner of the node.
c$$$      enddo
c$$$
CC    produce item 3.
c$$$      print *,'node         its_adjacent_nodes'

      do i = 0,nm
         nc = n_connectivity_s(i)
         info_of_s(1,i)=n_connectivity_s(i)
         mesh_n(i)=n_connectivity_s(i)
         m_start = m_start_s(i)
         do j = 1,nc
            tmp_s(j) = s_neighbor_m(m_start)
	    temp_owner = (tmp_s(j)/nprocs)
            m_start = m_next_m(m_start)
            info_of_n_s(i,(j-1),0)=tmp_s(j)
	    info_of_n_s(i,(j-1),1)=temp_owner
         enddo
c$$$         print *, i,'   ', (tmp_s(j),j=1,nc)
       enddo

C-------------Open the files and write the data------------
C----Important ****** Change the format Statement when using the high
c***
c***         open(unit=45,file='n_info',form='formatted',
c***     &      access='sequential',status='replace')
c***
c***         open(unit=46,file='n_count',form='formatted',
c***     &      access='sequential',status='replace')
c$$$
c$$$         open(unit=47,file='n_coord',form='formatted',
c$$$     &      access='sequential',status='new')
c$$$
c$$$         open(unit=48,file='np_coord',form='formatted',
c$$$     &      access='sequential',status='new')

         l=0
      do i =0, nm
c$$$         print *,'---------------------------------------------'
c$$$         print *,'I am node',i
c$$$         print *,'My owner is ',info_of_s(0,i)
c$$$         print *,'My neighbours are',info_of_s(1,i)
c***         write(unit=45,*) info_of_s(1,i)
c***c 67      format('',I3)
         nc = info_of_s(1,i)
c$$$         print *,'---------------------------------------------'
c$$$         print *,'neighbour        its owner   coordinates'
         do j=1,nc
           temp1=info_of_n_s(i,(j-1),0)
c$$$           temp2=info_of_n_s(i,(j-1),1)
c$$$         print *,temp1,'  ',temp2,' ',(coord_of_s(k,temp1),k=0,2)
c$$$         write(unit=48,fmt=66) coord_of_s(0,temp1),coord_of_s(1,temp1),
c$$$     &        coord_of_s(2,temp1)
c$$$ 66      format(F9.6,3X,F9.6,3X,F9.6)
c***         write(unit=46,*) temp1
         n_index(l)=temp1
         l=l+1
         enddo
      enddo

      RETURN
      END

C----------------------------------------------------------------------C

      BLOCK DATA
      REAL u(97), c, cd, cm
      INTEGER i97, j97, ivec
      LOGICAL test
      COMMON /raset1/ u, c, cd, cm, i97, j97, ivec, test
      DATA test /.false./
      END

C----------------------------------------------------------------------C

      SUBROUTINE flip_many_links

      INCLUDE 'irreg.h'
      INCLUDE 'makemesh.h'

      INTEGER link, s1, s2, s3, s4, i
      INTEGER c1, c2, c3, c4
      INTEGER i_flip
      INTEGER m_start, site, ma ,mb 
      INTEGER m_f1, m_f2, m_f3, m_f4
      INTEGER m_b1, m_b2, m_b3, m_b4
      
      REAL ranmar

      DO i_flip = 1, 3 * n

        link = n6 * (1.0-ranmar())

        IF(link.EQ.n6)THEN
          link = n6m
        ENDIF

        s1 = s_site_m(link)
        c1 = n_connectivity_s(s1)
        IF(c1 .LT. 4)GOTO 10
        s2 = s_neighbor_m(link)
        c2 = n_connectivity_s(s2)
        IF(c2 .LT. 4)GOTO 10

        m_f1 = m_former_m(link)
        m_b4 = m_next_m(link) 
        m_b1 = m_back_m(m_f1)
        m_f2 = m_former_m(m_b1)
        m_f4 = m_back_m(m_b4)
        m_b3 = m_next_m(m_f4)
        m_f3 = m_back_m(m_b3)
        m_b2 = m_back_m(m_f2)

        s3 = s_neighbor_m(m_f1)
        s4 = s_neighbor_m(m_b4)

        c3 = n_connectivity_s(s3)
        c4 = n_connectivity_s(s4)

        m_start = m_start_s(s3)
        do i=1,c3
           site = s_neighbor_m(m_start)
           IF(site .EQ. s4)GOTO 10
           m_start = m_next_m(m_start)
        enddo

        ma = link
        mb = m_back_m(link)

        m_former_m(m_b4) = m_f1
        m_next_m(m_f1) = m_b4
        if(link .eq. m_start_s(s1))m_start_s(s1) = m_f1

        m_former_m(m_b2) = m_f3
        m_next_m(m_f3) = m_b2
        if(mb .eq. m_start_s(s2))m_start_s(s2) = m_f3

        m_former_m(link) = m_f2
        m_next_m(link) = m_b1
        m_next_m(m_f2) = link
        m_former_m(m_b1) = link

        m_former_m(mb) = m_f4
        m_next_m(mb) = m_b3
        m_next_m(m_f4) = mb
        m_former_m(m_b3) = mb

        n_connectivity_s(s1) =  c1 - 1
        n_connectivity_s(s2) =  c2 - 1
        n_connectivity_s(s3) =  c3 + 1
        n_connectivity_s(s4) =  c4 + 1

        s_site_m(ma) = s3
        s_site_m(mb) = s4
        s_neighbor_m(ma) = s4
        s_neighbor_m(mb) = s3
        
 10     CONTINUE
        
      ENDDO
      
      RETURN
      END
      
