c	SOD - Seed Orchard Design Program - Mk 2.2

c	J.K. Vanclay,		October 1984

	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	logical ok
	character*80 title
	data rate/5./

	vacant=-999
	write(6,101)
101	format(/' #SOD  Mk 2.2')
	ind=0 
	nrams=1
	call input(title)
	maxram=ntot/nclone+1
	write(6,*)
	write(6,*) ntot,' ramets required'
10	num=0
	diff=0.
	call doit(ok,diff)
	if(ok) goto 50

c	this attempt failed - alter parameters and try again
	call clear
c
c	adjust weights
	do 40 i=1,nclone
	no=mclone(2,i)-idat(i,i)
	nexp=mclone(2,i)*float(num)/float(ntot)+.5
	weight(i)=weight(i)*(float(nexp)+rate)/(float(no)+rate)
40	continue
c
c	adjust horizon
	if(ind.eq.0.or.ind.eq.4) then
		horiz=horiz-1
		write(6,402) horiz
402		format(/' Horizon reduced to',i3)
		do 45 i=1,3
		if(horiz.gt.ns(i)) goto 45
		ind=i
		goto 10
45		continue
c
	elseif(ind.ge.1.and.ind.le.3) then
		if(.not.adjsep) then
			write(6,403) nrams
403   format(/' Ramets available for each clone increased by',i4)
			do 48 i=1,nclone
			if(mclone(2,i).gt.0) mclone(2,i)=mclone(2,i)+nrams
48			idat(i,i)=mclone(2,i)
			nrams=nrams*2
			if(nrams.gt.maxram) adjsep=.true.
		else
			if(ns(ind).ge.horiz) ns(ind)=ns(ind)-1
			ind=ind+1
			write(6,401) ns
401   format(/' Separations reduced to',3i3)
			if(horiz.eq.1.and.ind.eq.4) adjsep=.false.
		endif
	if(horiz.eq.1.and.ind.eq.4.and.nrams.gt.maxram) adjsep=.true.
	endif 
	goto 10

50	diff=diff*100.
	write(6,501) diff
501	format(/' Degree of Difficulty',f6.1,'%')
	call output(title)
	write(6,*) char(12)
	stop
	end



 
	subroutine doit (ok,diff)

c	starting from the centre and spiralling outwards,
c	take each position in turn and try to find a suitable ramet

	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	integer radius,circum
	logical ok

	ok=.true.
	nundo=0
	nsav=0
c	start from centre
	i=(is+1)/2
	j=(js+1)/2
	radius=0
	if(iplan(i,j).ne.vacant) goto 10
	num=num+1 
	call fill(i,j,nch,iclone)
	if(debug) write(6,101) num,i,j,nch,iclone

c	and spiralling outwards
10	radius=radius+1
	circum=0
20	circum=circum+1
	if(circum.gt.radius*8) goto 10

c	get rectangular coordinates
	call convrt(radius,circum,i,j)
	if(i.lt.1.or.i.gt.is) goto 20
	if(j.lt.1.or.j.gt.js) goto 20
	if(iplan(i,j).ne.vacant) goto 20
	num=num+1 
	call fill(i,j,nch,iclone)
	if(debug) write(6,101) num,i,j,nch,iclone
101	format(5i5)
	if(nch.gt.0) then
c		if orchard is full, we're finished
		if(num.ge.ntot) return
c		save any alternatives
		if(nch.gt.1) call savit(radius,circum,nsav)
		goto 20
	elseif(nsav.gt.0.and.nundo.lt.num/2) then
c		can't fill this position,
c		take a previous position and use an alternative
		diff=max(diff,nundo*2./num)
		call undo(radius,circum,nsav,ok)
		nundo=nundo+1
		if(ok) goto 20
	endif 

c	either we've run out of alternatives,
c	or we've undone it many times and it's time to start afresh
	num=num-1 
	write(6,*) 'Failed after',num,'ramets.'
	ok=.false.
	return
	end



 
	subroutine convrt(radius,circum,i,j)

c	given polar coordinates (radius, circum),
c	convert to rectangular coordinates (i,j)
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	integer radius,circum,side
	jstart=(js+1)/2
	istart=(is+1)/2
	if(circum.le.radius*2) then
		j=jstart+radius
		i=istart-radius+circum
	elseif(circum.le.radius*4) then
		i=istart+radius
		side=circum-radius*2
		j=jstart+radius-side
	elseif(circum.le.radius*6) then
		j=jstart-radius
		side=circum-radius*4
		i=istart+radius-side
	elseif(circum.le.radius*8) then
		i=istart-radius
		side=circum-radius*6
		j=jstart-radius+side
	else
		call abxrt('BAD CONV    ')
	endif 
	return
	end



 
	subroutine input(title)
	parameter (nokeys=17)

c	solicit the necessary input
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	character title*80,fname1*40,fname3*40
	integer ramets
	character status*8,form*12
	character*4 keynam,keys(nokeys)
	logical octal
	data lu/15/,keys/'ROT1','ROT2','ROT3','REFL','OCTA'
&	,'WEIG','SIMP','DETE','INCE','DEBU','RAME','SEPA'
&	,'CIRC','SQUA','OFFR','OFFT','    '/
	data status/'old     '/,form/'formatted   '/

	read(5,100,end=91) title
100	format(a80)
	write(6,103) title
103	format(///1x,a80)

	ibuf(81:132)=''
	read(5,101,end=92) ibuf
101	format(a80)
	i=1
	call getnum(ibuf,i,80,nclone)
	if(i.lt.80) call getnum(ibuf,i,80,is)
	if(i.lt.80) call getnum(ibuf,i,80,js)
	do 123 j=1,3
	if(i.lt.80) call getnum(ibuf,i,80,ns(j))
123	continue
	if(i.lt.80) call getnum(ibuf,i,80,horiz)
	if(i.lt.80) call getnum(ibuf,i,80,iseed)
	if(i.gt.80) goto 92
	write(6,106) nclone,is,js
106	format(5x,'Number of clones available:',i5/
&	,5x,'Size of orchard:',i5,' rows *',i5,' trees.')
	if(nclone.le.0) stop
	if(nclone.gt.mc) stop 'Too many clones'
	msep=sqrt(float(nclone))
	if(is.le.0.or.is.gt.mpi) stop 'Too many rows'
	if(js.le.0.or.js.gt.mpj) stop 'Too many cols'
c
	if(horiz.lt.0) horiz=4
	if(iseed.lt.0) iseed=0
	if(horiz.gt.msep) horiz=msep
	do 5 i=1,3
	if(ns(i).lt.0) ns(i)=2
	if(ns(i).ge.msep) ns(i)=msep-1
	if(ns(i).gt.horiz) horiz=ns(i)
5	continue
	write(6,107) ns
107	format(5x,'Separation to be maintained between half-sibs:',i2/
&	,41x,'full-sibs:',i2/,41x,'ramets of same clone:',i2)
	write(6,*) '    Horizon used in design:',horiz

	read(5,101,end=93) ibuf
	fname1=''
	fname2=''
	fname3=''
	do 51 i=1,80
	if(ibuf(i:i).gt.' ') goto 511
51	continue
	goto 93
511	do 512 j=i,80
	if(ibuf(i:i).eq.' '.or.ibuf(j:j).eq.',') goto 513
512	continue
	j=81
513	if(j-i+1.gt.40) j=i+40-1
	fname1(1:j-i+1)=ibuf(i:j)
	do 52 i=j+1,80
	if(ibuf(i:i).gt.' ') goto 521
52	continue
	goto 54
521	do 522 j=i,80
	if(ibuf(i:i).eq.' '.or.ibuf(j:j).eq.',') goto 523
522	continue
	j=81
523	if(j-i+1.gt.40) j=i+40-1
	fname2(1:j-i+1)=ibuf(i:j)
	do 53 i=j+1,80
	if(ibuf(i:i).gt.' ') goto 531
53	continue
	goto 54
531	do 532 j=i,80
	if(ibuf(i:i).eq.' '.or.ibuf(j:j).eq.',') goto 533
532	continue
	j=81
533	if(j-i+1.gt.40) j=i+40-1
	fname3(1:j-i+1)=ibuf(i:j)

54	write(6,105) fname1
105	format(5x,'Data files: ',a40)
	if(fname2.ne.'') write(6,104) fname2
	if(fname3.ne.'') write(6,104) fname3
104	format(17x,a40)

	nrot=0
	refl=.false.
	octal=.false.
	debug=.false.
	adjsep=.false.
	stratg=keys(8)
	circle=.true.
	ibuf(77:80)=''
10	read(5,102,end=90) keynam,ibuf
102	format(a4,a76)
	do 11 i=1,4
	if(keynam(i:i).ge.'a')
&		keynam(i:i)=char(ichar(keynam(i:i))-ichar('a')+ichar('A'))
11	continue
c**OLD	call intuch(4,keynam)
	write(6,108) keynam
108	format(5x,a4)
	do 19 i=1,nokeys
	if(keynam.eq.keys(i)) goto 192
19	continue
	i=0
c**OLD	i=iscan(keynam,nokeys,keys)
192	goto (20,20,20,30,40,50,50,50,50,60,70,75,80,85,86,87,90) i
	write(6,109)
109	format(5x,'----> Invalid Option')
	goto 10
20	nrot=i
	goto 10
30	refl=.true.
	goto 10
40	octal=.true.
	goto 10
50	stratg=keynam
	goto 10
60	debug=.true.
	goto 10
70	adjsep=.false.
	goto 10
75	adjsep=.true.
	goto 10
80	circle=.true.
	goto 10
85	circle=.false.
	goto 10
86	i=1
	call getnum(ibuf,i,80,ioff)
	goto 10
87	i=1
	call getnum(ibuf,i,80,joff)
	goto 10

90	open(lu,file=fname1,status=status,form=form)
	rewind lu
	call getcln(lu,ramets)
	call getpln(lu,octal)
	close (lu)
	if(ntot.gt.ramets) stop 'Too few ramets'
	call calink

	if(fname3.eq.'') goto 95
	open(lu,file=fname3,status=status,form=form)
	rewind lu
	call exist(lu,octal)
	close(lu)

95	if(stratg.eq.keys(8).or.stratg.eq.keys(9)) return
	write(6,*) '    Random selection initialized with',iseed
c*****	home version - rand(3f)
	call srand(iseed)
c*****	portable version
c*****	x=rand(iseed)
	return
91	stop 'No data'
92	stop 'Insufficient data'
93	stop 'No filename'
	end





	subroutine exist(lu,octal)

c	collect info on existing section of orchard

	include 'sod.incl'
	dimension list(mpj)
	logical octal

	jso=js
	if(octal) jso=(js+2)/3

	do 20 i=1,is
	read(lu,101)(list(j),j=1,jso)
101	format(80i1)
	if(octal) then
		do 11 j=jso,1,-1
		do 10 k=3,1,-1
		n=3*(j-1)+k
		if(n.gt.js) goto 10
		list(n)=mod(int(list(j)/2**(3-k)),2)
10		continue
11		continue
	endif
	do 15 j=1,js
	if(list(j).ne.0.and.iplan(i,j).ne.0) stop 'design conflict'
	if(list(j).ne.0) iplan(i,j)=1
15	continue
20	continue
	do 50 i=1,is
	no=0
	do 30 j=1,js
	if(iplan(i,j).eq.1) no=no+1
30	continue
	if(no.eq.0) goto 50
	read(lu,*)(list(j),j=1,no)
	no=0
	do 40 j=1,js
	if(iplan(i,j).eq.1) then
		no=no+1
		if(list(no).le.0.or.list(no).gt.nclone)
&			stop 'bad existing data'
		iplan(i,j)=-list(no)
	endif
40	continue
50	continue
	return
	end





	subroutine calink
c	compute possible links in orchard, number of cells not related,
c	and thus derive expected number of links per cell.
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2
	integer relatd

	n1=0
	n2=0
	do 25 ii=1,is
	do 20 jj=1,js
	if(iplan(ii,jj).eq.0) goto 20
	do 15 i=0,1
	ix=ii+i
	if(ix.lt.1.or.ix.gt.is) goto 15
	do 10 j=-1,1
	if(i.eq.0.and.j.le.0) goto 10
	jx=jj+j
	if(jx.lt.1.or.jx.gt.js) goto 10
	if(iplan(ix,jx).eq.0) goto 10
	n1=n1+1
10	continue
15	continue
20	continue
25	continue

	do 40 i=1,nclone-1
	do 30 j=i+1,nclone
	if(mclone(2,i).eq.0.or.mclone(2,j).eq.0) goto 30
	if(relatd(i,j).gt.0) goto 30
	n2=n2+1
30	continue
40	continue

	if(n2.gt.0) then
		nlpc=n1/n2
	else
		nlpc=0
	endif
	return
	end



 
	subroutine getcln(lu,ramets)

c	get clone names and numbers of ramets from file
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	character*4 iclone(2,mc)
	equivalence(mclone,iclone)
	integer ramets

	write(6,*)
	ramets=0
	do 5 i=1,nclone
5	weight(i)=2.
	do 10 i=1,nclone
	read(lu,101,end=90) clonam(i),iclone(1,i),mclone(2,i)
101	format(a8,a4,i4)
	write(6,102) clonam(i),iclone(1,i),mclone(2,i)
102	format(5x,a8,2x,a4,i4)
	ramets=ramets+mclone(2,i)
10	idat(i,i)=mclone(2,i)

c	find the related clones
	do 25 i=2,nclone
	do 20 j=1,i-1
	do 12 k=1,3,2
	l=4-k
	if(iclone(1,i)(k:k+1).eq.'  ') goto 12
	if(iclone(1,i)(k:k+1).eq.iclone(1,j)(k:k+1).or.
&		iclone(1,i)(k:k+1).eq.iclone(1,j)(l:l+1))goto 15
12	continue
	goto 20
15	idat(i,j)=1
	if((iclone(1,i).eq.iclone(1,j).or.
&		iclone(1,i)(1:2).eq.iclone(1,j)(3:4).and.
&		iclone(1,i)(3:4).eq.iclone(1,j)(1:2)).and.
&		iclone(1,i)(1:2).ne.'  '.and.iclone(1,i)(3:4).ne.'  ')
&		idat(i,j)=2
	weight(i)=weight(i)+1
	weight(j)=weight(j)+1
20	continue
25	continue
	return
90	stop 'Insufficient clone names'
	end



 
	subroutine getpln(lu,octal)

c	get the orchard plan from file
c	if no plan given, assume all positions to be grafted
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	logical octal

	ntot=0
	jso=js
	if(octal) jso=(js+2)/3
	do 25 i=1,is
	read(lu,101,end=80)(iplan(i,j),j=1,jso)
101	format(80i1)
	if(octal) then
		do 11 j=jso,1,-1
		do 10 k=3,1,-1
		n=3*(j-1)+k
		if(n.gt.js) goto 10
		iplan(i,n)=mod(int(iplan(i,j)/2**(3-k)),2)
10		continue
11		continue
	endif
	do 20 j=1,js
	if(iplan(i,j).eq.0) goto 20
	ntot=ntot+1
	iplan(i,j)=vacant
20	continue
25	continue
	return
80	do 95 j=i,is
	do 90 k=1,js
	iplan(j,k)=vacant
90	ntot=ntot+1
95	continue
	return
	end



 
	subroutine clear

c	houskeeping routine to clear arrays
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	do 15 i=1,is
	do 10 j=1,js
	if(iplan(i,j).gt.0) iplan(i,j)=vacant
10	continue
15	continue
	idat(1,1)=mclone(2,1)
	do 25 i=2,nclone
	idat(i,i)=mclone(2,i)
	do 20 j=1,i-1
20	idat(j,i)=0
25	continue
	do 35 i=1,ms2
	do 30 j=1,ms
30	isave(j,i)=0
35	continue
	return
	end



 
	subroutine fill (ii,jj,ich,iclone)

c	find a ramet suitable for position (ii,jj)
c	using the nominated strategy
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	dimension score(mc)

	iclone=0
c	calculate scores for each clone
	do 10 i=1,nclone
10	score(i)=idat(i,i)*weight(i)
	call look(ii,jj,score)


c	find score total and number of alternatives
	ich=0
	stot=0.0
	do 15 i=1,nclone
	if(score(i).gt.0.0) ich=ich+1
15	stot=stot+score(i)
	if(ich.eq.0) return

	if(stratg.eq.'DETE') then
c		choose the best
		best=0.0
		iclone=0
		do 20 j=1,nclone
		if(score(j).gt.best) then
			iclone=j 
			best=score(j)
		endif
20		continue
		if(iclone.ne.0) goto 45
	elseif(stratg.eq.'WEIG') then
c		use weighted random selection
		z=rand()*stot
c*****portable	z=rand(0)*stot
		zz=0.0
		do 25 i=1,nclone
		zz=zz+score(i)
		if(z.le.zz.and.score(i).gt.0.0) goto 40
25		continue
	elseif(stratg.eq.'SIMP') then
c		use simple random selection
		n=rand()*ich +1
c*****portable	n=rand(0)*ich +1
		nn=0 
		do 30 i=1,nclone
		if(score(i).gt.0.0) nn=nn+1
		if(nn.eq.n) goto 40
30		continue
	elseif(stratg.eq.'INCE') then
c		choose the worst which meets minimum criteria
		best=stot
		iclone=0
		do 35 j=1,nclone
		if(score(j).gt.0.0.and.score(j).le.best) then
			iclone=j 
			best=score(j)
		endif
35		continue
		if(iclone.ne.0) goto 45
	endif 
	call abxrt('ERR IN FILL ')
40	iclone=i
45	iplan(ii,jj)=iclone
	idat(iclone,iclone)=idat(iclone,iclone)-1
	call put(ii,jj,iclone,1)
	return
	end



 
	subroutine undo(i,j,nsav,ok)

c	we've reached an impasse at i,j
c	take a previous position, and try one of the alternatives
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	dimension temp(ms3)
	logical ok

	ok=.true.
10	ii=isave(nsav,1)
	jj=isave(nsav,2)
	do 25 ll=i,ii,-1
	do 20 mm=ll*8,1,-1
	if(ll.eq.i.and.mm.gt.j) goto 20
	if(ll.eq.ii.and.mm.lt.jj) goto 20
	call convrt(ll,mm,l,m)
	if(l.lt.1.or.l.gt.is) goto 20
	if(m.lt.1.or.m.gt.js) goto 20
	if(iplan(l,m).le.0) goto 20
	n=iplan(l,m)
	iplan(l,m)=vacant
	idat(n,n)=idat(n,n)+1
	call put(l,m,n,-1)
	if(debug)write(6,*)'Remove ',num,l,m,n
	num=num-1
	if(num.le.0) call abxrt('BAD NUM     ')
20	continue
25	continue
	do 30 kk=1,ms3
	k=kk+2
	if(isave(nsav,k).eq.0) isave(nsav,k)=n
	isv=isave(nsav,k)
	temp(kk)=weight(isv)
	weight(isv)=0.
	if(isv.eq.n) goto 31
30	continue
	kk=ms3
31	call convrt(ii,jj,l,m)
	call fill(l,m,nch,iclone)
	if(debug) write(6,*) 'Undo to ',num,l,m,nch,iclone,kk
	do 40 k=1,kk
40	weight(isave(nsav,k+2))=temp(k)
	if(kk.ge.ms3.or.nch.eq.0) then
		isave(nsav,1)=0
		nsav=nsav-1
		if(nsav.eq.0) nsav=ms1
		if(isave(nsav,1).eq.0) nsav=0
	endif 
	if(nch.eq.0.and.nsav.ne.0) goto 10
	if(nch.eq.0) ok=.false.
	i=ii
	j=jj
	return

	entry savit(ix,jx,nsv)

c	save position which has alternatives
	nsv=nsv+1
	if(nsv.gt.ms1) nsv=1
	isave(nsv,1)=ix
	isave(nsv,2)=jx
	do 50 it=3,ms2
50	isave(nsv,it)=0
	return
	end



 
	subroutine put(ii,jj,k,inc)

c	do housekeeping associated with putting clone k in position i,j
c		(or with removing it when called from undo)
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	do 15 i=-1,1
	l=ii+i
	if(l.lt.1.or.l.gt.mpi) goto 15
	do 10 j=-1,1
	m=jj+j
	if(m.lt.1.or.m.gt.mpj) goto 10
	if(i.eq.0.and.j.eq.0) goto 10
	if(iplan(l,m).eq.vacant.or.iplan(l,m).eq.0) goto 10
	k1=abs(iplan(l,m))
	if(k.eq.k1) goto 10
	n1=min(k,k1)
	n2=max(k,k1)
	idat(n1,n2)=idat(n1,n2)+inc
10	continue
15	continue
	return
	end



 
	subroutine output(title)

c	arrange output of the successful design
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	character*80 title
	dimension nrels(mc,5),nrtot(5)

	do 10 j=1,js,44
	jj=j+43
	if(jj.gt.js) jj=js
	call oplan(j,jj,title)
10	continue

	call comrel(nrels,nrtot)
	call pramet(title,nrels,nrtot,score)

	n1=0
	n2=0
	do 30 j=2,nclone,23
	jj=j+22
	if(jj.gt.nclone) jj=nclone
	call plinks(j,jj,title,n1,n2)
30	continue

	write(6,*) '  For current section of orchard only:'
	write(6,*) '    Attainable links per cell =',nlpc
	if(n1.gt.0) then
		t=100*(1.-float(n2)/float(n1))
		write(6,301) t
301		format(5x,'Panmixia =',f5.1,' %')
		t=t*score/100
		write(6,302) t
302		format(//5x,'Composite Score =',f5.1,' %')
	endif

	call grafti
	do 40 i=1,nclone,7
	j=i+6
	if(j.gt.nclone) j=nclone
	call pclone(i,j,title)
40	continue

	do 80 i=1,2
	do 70 j=1,js,20
	jj=j+19
	if(jj.gt.js) jj=js
	call pplan(j,jj,title,i)
70	continue
80	continue
	return
	end




	subroutine comrel(nrels,nrtot)
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	dimension nrels(mc,5),nrtot(5)
	integer relatd

	do 6 j=1,5
	nrtot(j)=0
	do 5 i=1,mc
5	nrels(i,j)=0
6	continue

	do 35 i=1,is
	do 30 j=1,js
	k=iplan(i,j)
	ka=abs(k)
	if(k.eq.0) goto 30

	do 20 k1=0,5
	l1=i+k1
	if(l1.lt.1.or.l1.gt.is) goto 20

	do 10 k2=-5,5
	if(k1.eq.0.and.k2.le.0) goto 10
	l2=j+k2
	if(l2.lt.1.or.l2.gt.js) goto 10
	if(circle.and.(k1*k1+k2*k2).ge.36) goto 10
	l=iplan(l1,l2)
	if(l.eq.0) goto 10
	if(k.lt.0.and.l.lt.0) goto 10
	la=abs(l)

	if(circle) then
		m=int(sqrt(float(k1*k1+k2*k2)))
		if(m.gt.5) goto 10
	else
		m=max(iabs(k1),iabs(k2))
	endif
	nrtot(m)=nrtot(m)+2
	if(relatd(ka,la).eq.0) goto 10
	nrels(ka,m)=nrels(ka,m)+1
	nrels(la,m)=nrels(la,m)+1
10	continue
20	continue
30	continue
35	continue
	return
	end



 
	integer function relatd(i,j)
c	if i and j are unrelated, relatd=0
c			halfsibs, relatd=1
c			fullsibs, relatd=2
c			identical relatd=3
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	if(i-j) 10,20,30

10	relatd=idat(j,i)
	return
20	relatd=3
	return
30	relatd=idat(i,j)
	return
	end



 
	subroutine grafti

c	computations necessary to produce grafing instructions
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2
	logical dump
	character status*8,form*12
	data status/'unknown '/,form/'formatted   '/
	data lu/16/

	dump=.false.
	if(fname2.ne.'') dump=.true.
	if(dump) then
		open(lu,file=fname2,status=status,form=form)
		rewind lu
	endif

	do 6 i=1,mc
	do 5 j=1,mc
5	idat(i,j)=0
6	continue
	do 10 i=1,nclone
10	mclone(2,i)=0
	do 25 i=1,is
	do 20 j=1,js
	if(iplan(i,j).le.0) goto 20
	call adjust(i,j,ii,jj)
	k=iplan(i,j)

	if(dump) write(lu,101) ii,jj,k
101	format(3i4)

	mclone(2,k)=mclone(2,k)+1
	n=mclone(2,k)
	idat(n,k)=ii*1000+jj
20	continue
25	continue

	if(dump) close(lu)

	if(.not.refl.and.nrot.eq.0) return
	do 30 i=1,nclone
	call qkrsi(idat(1,i),mclone(2,i))
30	continue
	return
	end



 
	subroutine adjust(i,j,ii,jj)

c	calculate the necessary transformation corresponding to the
c	specified reflection and rotation of the basic design
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	ii=i
	jj=j
	if(refl) jj=js-jj+1
	goto (1,2,3) nrot
	goto 10

1	ii=js-jj+1
	jj=i
	goto 10

2	ii=is-ii+1
	jj=js-jj+1
	goto 10

3	ii=jj
	jj=is-i+1

10	ii=ii+ioff
	jj=jj+joff
	return


	entry adjcol(j1,j2)
	j2=j1
	if(refl) j2=js-j1+1
	if(nrot.eq.1.or.nrot.eq.2) j2=js-j2+1
	if(nrot.eq.0.or.nrot.eq.2) then
		j2=j2+joff
	else
		j2=j2+ioff
	endif
	return

	entry adjrow(i1,i2)
	i2=i1
	if(nrot.gt.1) i2=is-i1+1
	if(nrot.eq.0.or.nrot.eq.2) then
		i2=i2+ioff
	else
		i2=i2+joff
	endif
	return
	end



 
	subroutine pclone(i,j,title)

c	print grafting instructions
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	character title*80
	data nl/60/

	iblank(1)='    '
	ibuf='    '
	m=0
	do 5 k=i,j
5	m=max(mclone(2,k),m)
	if(m.eq.0) return

	nl=nl+m+6
	if(nl.gt.60) then
		nl=m+11
		write(6,*) char(12)
		write(6,10004) title
10004	format(//40x,a80//40x,'Positions of Clones for Grafting')
	endif
	write(6,10001)(k,k=i,j)
10001	format(///7(3x,'Clone',i5:,6x))
	write(6,10002) (clonam(k),k=i,j)
10002	format(7(5x,a8:,6x))
	write(6,10005)(iblank(1),k=i,j)
10005	format(7(a4,'# Row Tree':,5x))

	do 20 n=1,m
	write(ibuf,10003)(n,int(idat(n,k)/1000),mod(idat(n,k),1000),k=i,j)
10003	format(7(1x,3i4:,6x))
	do 15 k=i,j
	l=k-i
15	if (n.gt.mclone(2,k)) ibuf(l*19+2:l*19+13)='            '
20	write(6,10006) ibuf
10006	format(a132)
	return
	end



 
	subroutine oplan(ii,jj,title)

c	print basic plan
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2
	character*80 title

	call reduce(ii,i,1)
	call reduce(jj,j,-1)
	if(j-i.ge.44.or.j-i.lt.0) call abxrt('OPLAN ERROR ')
	ibuf='    '

	do 20 kk=1,is,50
	write(6,*) char(12)
	write(6,10001)title
10001	format(//40x,a80//40x,'Basic Plan'/)

	do 10 k=kk,kk+49
	if(k.le.0.or.k.gt.is) goto 10
	write(ibuf,10004) (iplan(k,l),l=i,j)
10004	format(44i3)
	call zot(132,2,ibuf)
	write(6,10005) ibuf
10005	format(a132)
10	continue
20	continue
	return
	end



 
	subroutine pramet(title,rels,nrtot,score)

c	print list of clones with ramets
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2
	character*80 title
	integer relatd,rels
	dimension list(mc),rels(mc,5),nrtot(5)

	write(6,*) char(12)
	write(6,101) title,(i,i=1,5)
101	format(//40x,a80//40x,'List of Clones with Numbers of Ramets '
&	,'and Relationships'//' Clone and Pedigree    Ramets',5x
&	,'Related Ramets at Distance      Related Clones'/32x,5i5/)

	score=0.
	do 40 i=1,nclone
	mclone(2,i)=mclone(2,i)-idat(i,i)
	sc=1.
	do 20 j=1,5
	sc=sc/2.
20	score=score+sc*rels(i,j)
	n=0
	do 30 j=1,nclone
	if(i.eq.j) goto 30
	if(relatd(i,j).eq.0) goto 30
	n=n+1
	list(n)=j
30	continue
40	write(6,103)i,clonam(i),(mclone(j,i),j=1,2),(rels(i,j),j=1,5)
&	,(list(k),k=1,n)
103	format('  #',i2,2x,a8,2x,a4,i6,5x,5i5,6x,16i4,(/63x,16i4))

	stot=0.
	sc=1.
	do 50 j=1,5
	sc=sc/2.
50	stot=stot+sc*nrtot(j)
	score=100*(1.-score/stot)
	write(6,104) score
104	format(/10x,'Separation =',f5.1,' %')
	return
	end



 
	subroutine plinks(i,j,title,n1,n2)

c	print table showing panmixia
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2
	character title*80,kk*4
	character*1 char(4)
	data char/' ','H','F','.'/

	write(6,*) char(12)
	write(6,101) title
101	format(//40x,a80//40x,'Adjacencies between Clones ',
&	 '(H = Half sibs, F = Full sibs)'/)
	write(6,105) (k,k=i,j)
105	format(1x,25(' #',i2:,1x))
	write(6,102) (clonam(k)(1:4),k=i,j)
	write(6,102) (clonam(k)(5:8),k=i,j)
102	format(1x,25(1x,a4))

	do 10 k=1,j-1
	ibuf=' '
	n=i
	if(k.ge.i) n=k+1
	do 5 l=n,j
	m=(l-i)*5+4
	ibuf(m:m)=char(1+idat(l,k))
	m=m+1 
	if(idat(k,l).gt.0) then
		write(char(4),103) idat(k,l)
103		format(i1)
		ibuf(m:m)=char(4)
	else
		ibuf(m:m)='.'
	endif 
	if(idat(l,k).eq.0.and.mclone(2,l).ne.0.and.mclone(2,k).ne.0)then
		n1=n1+max(1,nlpc)
		if(idat(k,l).lt.nlpc) n2=n2+(nlpc-idat(k,l))
	endif 
5	continue
	write(kk,106) k
106	format(' #',i2)
	ibuf(5*(j-i)+7:5*(j-i)+10)=kk
	ibuf(5*(j-i)+12:5*(j-i)+19)=clonam(k)
10	write(6,104) ibuf
104	format(a132)
	return
	end



 
	subroutine pplan(ii,jj,title,ind)

c	print detailed design
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2
	character*8 clnam
	dimension iclone(2,mcx),clnam(mcx)
	equivalence (iblank,iclone),(bcl,clnam)

	character title*80
	dimension list(mpj)
	character*4 lable(2)
	character*12 lab(2)
	dimension linmax(2)
	data linmax/18,14/
	data lab/'Pedigree    ','Clone Names '/
	data lable/'Tree','Row '/

	clnam(1)='        '
	call reduce(ii,ir,1)
	call reduce(jj,jr,-1)
	if(jr-ir.ge.20.or.jr-ir.lt.0) call abxrt('PPLAN ERROR ')
	ibuf='    '
	write(6,*) char(12)
	write(6,10001)title,lab(ind)
10001	format(//40x,a80//40x,'Design showing',1x,a12/)
	l=mod(nrot,2)+1
	write(6,10008) lable(3-l),(lable(l),k=ir,jr),lable(3-l)
10008	format(/1x,a4,21(1x,a4,1x))
	do 5 k=ir,jr
	call adjcol(k,ka)
5	list(k)=ka
	write(6,10002) (list(k),k=ir,jr)
10002	format(3x,20i6)
	line=1

	do 10 k=1,is
	call adjrow(k,ka)
	line=line+1
	if(line.le.linmax(ind)) goto 15
	line=0
	write(6,*) char(12)
	write(6,10010) lable(3-l),(lable(l),m=ir,jr),lable(3-l)
10010	format(//a4,21(1x,a4,1x))
	write(6,10002) (list(m),m=ir,jr)
15	write(ibuf,10003) ka,(abs(iplan(k,m)),m=ir,jr),ka
10003	format(i4,21(i5,1x))
	call zot(132,1,ibuf)
	write(6,10004) ibuf
10004	format(/a132)
	if(ind.eq.1) then
		write(6,10006) (iclone(1,abs(iplan(k,m))+1),m=ir,jr)
	else
		write(6,10006) (clnam(abs(iplan(k,m))+1)(1:4),m=ir,jr)
		write(6,10006) (clnam(abs(iplan(k,m))+1)(5:8),m=ir,jr)
	endif
10006	format(4x,20(2x,a4))
10	continue
	return
	end



 
	subroutine reduce(jj,j,inc)

c	eliminate redundant entries from output
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2

	j=jj
10	do 20 i=1,is
	if(iplan(i,j).eq.0) goto 20
	return
20	continue
	j=j+inc
	goto 10
	end



 
	subroutine zot(j,i,buf)

c	remove surplus zeroes from print line
	character*132 buf
	do 10 k=i,j-1
	if(buf(k:k+1).eq.' 0') buf(k:k+1)='  '
10	continue
	return
	end



	subroutine abxrt(msg)
	character*12 msg
	write(6,*) msg
	write(6,*)
	write(6,*)'You have exposed a previously undetected error',
&	'in the program.'
	write(6,*)'Please report the error with all details you can',
&	'supply to:'
	write(6,*)'   J.K. Vanclay,'
	write(6,*)'   c/- Department of Forestry,'
	write(6,*)'   Box 944, G.P.O.,'
	write(6,*)'   Brisbane, Qld, 4000.'
	stop
	end





	subroutine look(ii,jj,score)
c	examine all positions in the neighbourhood of position (ii,jj),
c	and compute a score for each clone
	parameter (mc=100, mpi=50, mpj=50, ms1=5, ms2=12)
	parameter (ms3=ms2-2, mcx=mc+1)
	logical refl,debug,adjsep,circle
	character stratg*4, bcl*8,clonam*8,ibuf*132,iblank*4,fname2*40
	integer horiz,vacant
	common/zero/nclone,is,js,ns(3),horiz,num,ntot,nlpc,nrot
     2,refl,debug,adjsep,circle,stratg,ioff,joff,vacant
	common/one/iblank(2),mclone(2,mc)
	common/two/iplan(mpi,mpj)
	common/three/idat(mc,mc),weight(mc)
	common/four/isave(ms1,ms2)
	common/five/bcl,clonam(mc)
	common/edcom/ibuf,fname2
	dimension score(mc),sep(3)
	integer relatd,sep

	msep2=(horiz+1)**2
	do 50 k=1,nclone
	if(score(k).le.0) goto 50
	do 5 m=1,3
5	sep(m)=msep2

	do 30 i=-horiz,horiz
	ix=ii+i
	if(ix.lt.1.or.ix.gt.is) goto 30

	do 20 j=-horiz,horiz
	if(i.eq.0.and.j.eq.0) goto 20
	jx=jj+j
	if(jx.lt.1.or.jx.gt.js) goto 20

	if(circle.and.(i*i+j*j).ge.msep2) goto 20
	l=iplan(ix,jx)
	if(l.eq.0.or.l.eq.vacant) goto 20
	l=abs(l)
	m=relatd(l,k)
	if(m.gt.0) goto 10

c	ramets not related - if adjacent, adjust scores
	if(i*i+j*j.gt.2) goto 20
	if(k-l) 7,20,8
7	d=idat(k,l)
	goto 9
8	d=idat(l,k)
9	score(k)=score(k)*2.*(d+5.)/(d+1.)
	goto 20

c	related ramets
10	if(circle) then
		sep(m)=min(i*i+j*j,sep(m))
	else
		sep(m)=min(max(i*i,j*j),sep(m))
	endif
20	continue
30	continue
	do 40 m=1,3
	if(sep(m).lt.(ns(m)+1)**2) then
		score(k)=0.
	else
		score(k)=score(k)*sep(m)
	endif
40	continue
50	continue
	return
	end




	subroutine getnum(arg,i,j,k)
	character arg*80,fmt*4

c	read a number from string of digits in arg,
c	from between the positions i and j (inclusive)
c	and return the value in k (-1 indicates invalid)
c	Update i, so that it points to space following the number

	k=-1

c	find first digit
	do 10 m=i,j
	if(arg(m:m).ge.'0'.and.arg(m:m).le.'9') goto 20
10	continue
	i=j+1
	return

20	do 25 l=m,j
	if(arg(l:l).lt.'0'.or.arg(l:l).gt.'9') goto 30
25	continue
	l=j+1
30	i=l
	l=l-1
	write(fmt,101) l-m+1
101	format('(i',i1,')')
	read(arg(m:l),fmt) k
	return
	end



	subroutine qkrsi(m,n)
	dimension m(n)
	do 50 i=1,n
	icnt=0
	do 10 j=2,n
	if(m(j-1).gt.m(j)) then
		k=m(j)
		m(j)=m(j-1)
		m(j-1)=k
		icnt=icnt+1
	endif
10	continue
	if(icnt.le.0) return
50	continue
	stop 'error in sort'
	end