c
c	INIT  :  Generates the initial conditions for the matrix
c	--------------------------------------------------------------
c	this subr. reads in a microstructure from a previous run from a file
c	called "micro.input"; the information can be copied in from a "ph" file
c
c	IMPORTANT!  the dimensions of the structure must match those
c	of the previous simulation.  Therefore a check is performed and
c	a fatal error generated if a mismatch is detected.
c
c	Grains with spin less than or equal to q/2 are unrecrystallized
c	grains.
c
	program uniquespin
	include 'common.f'
	integer site, mfile, nfile, lenfile, col, row
	real  linefile, nucfile, rbarfile, tfile, tempfile, radius, bdry
	character*25  run_name,fname
c
	write(*,*) 'This program takes an input and re-writes it with'
	write(*,*) ' one spin number per grain'
	write(*,*) 'This version set up for M x N = ',m,n
	write(*,*) 'If you need a different structure, change common.f'
	write(*,*)
	call prep
c
10	write(*,*) 'Enter the file name of the file to convert, in quotes'
	read(*,*) fname
c
	open(17,file=fname,status='old',err=10)
	read (17,*) mfile, nfile
c
	if((mfile.ne.m).or.(nfile.ne.n))    stop 'mismatched lattices'
	read (17,*)  run_name,tfile,rbarfile, tempfile
	read (17,*)  lenfile,linefile,nucfile
c	
	read (17,*) (spins(site), site=1,mn)
	close (17)
c
	call infiniq(icount,r2bar)
c
	write(*,*) 'check in INIT:  '
	write(*,*) 'No. grains = Q = ',icount
	write(*,*) 'Av. radius = ',rbar
	write(*,*) 'Av. Area = ',abar
c
20	write(*,*) 'Enter the file name of the file to write'
	read(*,*) fname
c
	open(17,file=fname,status='new',err=20)
	write (17,*) mfile, nfile, icount
	write (17,*)  "'",run_name,"'",tfile,rbarfile,tempfile
	write (17,*)  lenfile,linefile,nucfile
c	
	write (17,*) (spins(site), site=1,mn)
	close (17)
c
	call exit
	end
c
c  -------------------
c
	subroutine infiniq(icount,r2bar)
c  
c  converts a structure to infinite Q
c  calculates sizes also
c
	include 'common.f'
c
	integer site,sorted(mn),spinfix,nn,nlist(40000),listend
	real asum,rsum,r2bar
	real isize(40000)
	real pi
	data pi /3.1415926536/
c
	icount=1
	jcount=0
	abar=0.
	rbar=0.
	a2bar=0.
	r2bar=0.
c  ICOUNT is the spin number
c  JCOUNT is the grain area
c  ISIZE has the size of each grain
c
	do 5, i=1,mn
5	sorted(i)=spins(i)
c
	nn=0
10	call find(site,nindex,spinfix,sorted)
	write(*,*) 'in infiniq: success= ',nindex,'  site= ',
     &	site, ' spin= ',spinfix
	if(nindex.eq.0) then
		goto 200
c  could not find a new spin to burn, so complete
		else
		sorted(site)=-1*icount
		jcount=1
	endif
c
	write(*,*) 'in infiniq: icount= ',icount
	call findun(site,spinfix,nlist,listend,sorted)
c
c  LISTEND is the pointer to the last entry in the list of sites unburnt
c  NLIST is a list of site numbers to be burnt
c
20	if(listend.gt.0) then
		site=nlist(listend)
c  point to the last site on the list
		if(sorted(site).gt.0) then
			sorted(site)=-1*icount
c  if not already burnt, then do so
			jcount=jcount+1
c			write(*,*) 'in INFINIQ, burnt site ',site,'   grain no.',icount
		endif
		listend=listend-1
		call findun(site,spinfix,nlist,listend,sorted)
		goto 20
	else
		isize(icount)=jcount
c  record the size of the grain
		icount=icount+1
		jcount=0
		listend=0
		goto 10
	endif
c
c
200	icount=icount-1
c
	write(*,*) 'Total number of grains= ',icount
	do 210, i=1,icount
		write(*,*) i,isize(i)
210	continue
c
	do 300, i=1,icount
	rsize=isize(i)
	abar=abar+rsize
	rbar=rbar+sqrt(rsize/pi)
	a2bar=a2bar+rsize**2
	r2bar=r2bar+rsize/pi
300	continue
	rtemp=float(icount)
	abar=abar/rtemp
	rbar=rbar/rtemp
	a2bar=a2bar/rtemp
	r2bar=r2bar/rtemp
c
	write(*,*) 'INFINIQ is COMPLETE'
	write(*,*) 'No. grains = Q = ',icount
	write(*,*) 'Av. radius = ',rbar
	write(*,*) 'Av. Area = ',abar
c
c	if(icount.gt.q) stop 'too many grains, > Q !!!!'
c
	do 400, i=1,mn
		ijk=iabs(sorted(i))
c		if(ijk.eq.1) then
c			spins(i)=2
c		else
c			spins(i)=q2+ijk
c			spins(i)=q2
			spins(i)=ijk
c		endif
400	continue
c
c	write(*,*) 'hit return when ready!'
c	read (*,*) dummy
	return
	end
c
c  -------------------
c
	subroutine findun(zsite,spinfix,nlist,listend,sorted)
c  to add unburnt sites neighboring isite to the list to be burned
c
	include 'common.f'
	logical qtest
	integer zsite,nlist(40000),nn,spinfix,sorted(mn),listend
c
c	write(*,*) 'in findun: site,spinfix,listend= ',zsite,spinfix,listend
c	if(nbors.eq.6) then
		limit=nbors
c		endif
c	if(nbors.eq.8) then
c		limit=4
c		endif
	if(nbors.eq.18) then
		limit=6
		endif
c	if(nbors.eq.4) then
c		limit=4
c		endif
	if(nbors.eq.8) limit=4
c	nn=0
	do 100, i=1,limit
		if(sorted(neighs(zsite,i)).eq.spinfix) then
			qtest=.true.
c  if QTEST is true, then the neighbor site is not already in NLIST
			do 50, j=1,listend
			   if(neighs(zsite,i).eq.nlist(j)) qtest=.false.
50			continue
			if(qtest) then
				listend=listend+1
				nlist(listend)=neighs(zsite,i)
			endif
		endif
100	continue
c	write(*,*) 'end of findun: no. of like sites= ',listend
	return
	end
c
c  -------------------
c
	subroutine find(site,isuccess,spinfix,sorted)
	include 'common.f'
	integer site,spinfix,isuccess,sorted(mn)
	isuccess=0
	site=0
	do 100, i=1,mn
	if(sorted(i).gt.0) then
		site=i
		spinfix=spins(i)
		isuccess=1
		return
	endif
100	continue
	return
	end
c
c
c	PREP :  Prepares the array and parameters for the run
c	-----------------------------------------------------------------
c	copyright 1992 Elizabeth A. Holm
c
c	-sets up a table of the neighbors of sites
c	-initializes the grain structure
c	-places second phase particles in the matrix
c	-places recrystallized grain nuclei in the matrix
c	-determines time increment for the classical Monte Carlo routine
c
c	To change lattice type, uncomment the appropriate neighbor
c	tabulator below
c
	subroutine prep
	include 'common.f'
	integer i,site,kjp1,kjm1,kip1,kim1
	character*10 lattype
c
       PI=3.14159265
c
c
c	for each matrix site:
c	List the neighbors of 'site' in a table  with entries of the 
c	form 'neighs(site,neighbor)'.  Periodic boundary conditions 
c	are used.
c
c	in the tri(1) lattice, 'neighbor' is given by
c
c			3	2
c		4      site	1
c		5	6
c
c	lattype='tri(1)'
c	do 20 j=0,mn-m,m
c	kjp1=mod(j+m,mn)
c	kjm1=mod(j-m+mn,mn)
c		do 10 i=1,n,1
c		kip1=mod(i,n)
c		kim1=mod(i-2+n,n)
c		site=i+j
c		neighs(site,1)=j+kip1+1
c		neighs(site,2)=kjm1+kip1+1
c		neighs(site,3)=kjm1+i
c		neighs(site,4)=j+kim1+1
c		neighs(site,5)=kjp1+kim1+1
c		neighs(site,6)=kjp1+i
c10		continue
c20	continue
c
c	in the sq(1,2) lattice, 'neighbor' is given by
c
c		6	2	5
c		3      site	1
c		7	4	8
c
c	lattype='sq(1,2)'
c	do 20 j=0,mn-m,m
c		do 10 i=1,m
c			site=j+i
c			neighs(site,1)=j+mod(i,m)+1
c			neighs(site,2)=mod(j-m+mn,mn)+i
c			neighs(site,3)=j+mod(i-2+m,m)+1
c			neighs(site,4)=mod(j+m,mn)+i
c			neighs(site,5)=mod(j-m+mn,mn)+mod(i,m)+1
c			neighs(site,6)=mod(j-m+mn,mn)+mod(i-2+m,m)+1
c			neighs(site,7)=mod(j+m,mn)+mod(i-2+m,m)+1
c			neighs(site,8)=mod(j+m,mn)+mod(i,m)+1
c10		continue
c20	continue
c
c	in the tri(1,2) lattice, 'neighbor' is given by
c
c                       11      10      9
c               12      3       2       8
c	13	4      site     1       7
c	14	5       6       18
c	15	16      17
c
	noffset=(m/2)*n			! offset for kip below
	nmoffset=mn-noffset		! can use this for subtraction
c
c  the new lines that change the offsets at the left
c and right borders of the lattice ensure that there is
c  no "skew" in the geometry of the domain;  adr iv 01
	lattype='tri(1,2)'
	do 20 j=0,mn-m,m
	kjp1=mod(j+m,mn)
	kjm1=mod(j-m+mn,mn)
		do 10 i=1,n,1
		kip1=mod(i,n)
		kim1=mod(i-2+n,n)

		site=i+j
		neighs(site,1)=j+kip1+1
		neighs(site,2)=kjm1+kip1+1
		neighs(site,3)=kjm1+i
		neighs(site,4)=j+kim1+1
		neighs(site,5)=kjp1+kim1+1
		neighs(site,6)=kjp1+i
10		continue
20	continue
	do 25 site=1,mn
		n1=neighs(site,1)
		n2=neighs(site,2)
		if(i.eq.n) then
			neighs(site,1)=mod(j+noffset,mn)+kip1+1
c			write(*,*) 'i,j,site,neighs(site,1) ',i,j,site,neighs(site,1)
			neighs(site,2)=mod(kjm1+noffset,mn)+kip1+1
c			write(*,*) 'i,j,site,neighs(site,2) ',i,j,site,neighs(site,2)
		endif
		n3=neighs(site,3)
		n4=neighs(site,4)
		n5=neighs(site,5)
		if(i.eq.1) then
			neighs(site,4)=mod(j+nmoffset,mn)+kim1+1
c			write(*,*) 'i,j,site,neighs(site,4) ',i,j,site,neighs(site,4)
			neighs(site,5)=mod(kjp1+nmoffset,mn)+kim1+1
c			write(*,*) 'i,j,site,neighs(site,5) ',i,j,site,neighs(site,5)
		endif
		n6=neighs(site,6)
		neighs(site,7)=neighs(n1,1)
		neighs(site,8)=neighs(n2,1)
		neighs(site,9)=neighs(n2,2)
		neighs(site,10)=neighs(n3,2)
		neighs(site,11)=neighs(n3,3)
		neighs(site,12)=neighs(n4,3)
		neighs(site,13)=neighs(n4,4)
		neighs(site,14)=neighs(n5,4)
		neighs(site,15)=neighs(n5,5)
		neighs(site,16)=neighs(n6,5)
		neighs(site,17)=neighs(n6,6)
		neighs(site,18)=neighs(n1,6)
25	continue 
c
	return
	end