	program ising
c
c ***********************************************************************
c       Simple Ising model simulation on square lattice, no tricks.
c ***********************************************************************
c
	parameter (maxn=101,maxt=100001)
	implicit real*8 (a-h,o-z)
	common ispin(0:maxn,0:maxn)
	dimension order(maxt)
	character*12 of1,of2,ans,yes
	real etime,time(2)
	yes='y'
10	print *,'Ising model on square lattice:'
	print *,'nedge, Ti, Tf, #Tdivisions, init[1=ordered], seed ?'
	read(5,*) nedge,t1,t2,nt,init,isd
	if(nedge.gt.(maxn-1).or.nt.lt.0) go to 10
	print *,'ms0 (transient), ms (interval), ns (how many) ?'
	read(5,*) ms0,ms,ns
	print *,'periodic (1) or free (2) boundary ?'
	read(5,*) ib
	print *,'output file for magnetization and energy ?'
	read(5,15) of1
15	format(a12)
	print *,'record magnetization time series ?'
	read(5,16) ans
16	format(a1)
	if(ans.eq.yes) then
		print *,'output file for magnetization time series ?'
		read(5,15) of2
		if((ms0+(ns-1)*ms+1).gt.maxt) go to 10
	endif
	time0=etime(time)
	call initialize(t1,nedge,isd,init,ms0,ib,order)
	open(1,file=of1)
	rewind(1)
	if(nt.gt.0) then
		dt=(t2-t1)/nt
	else
		dt=0.d0
	endif
	write(6,35)
35	format(/2x,'Temperature',1x,'Magnetization',3x,'Energy',2x,
     c         'Susceptibility',1x,'Heat Capacity'/)
	m1=ms0
	do 20 i=1,nt+1
		t=t1+(i-1)*dt
		call calc(t,ord,en,chi,ch,nedge,m1,ms,ns,isd,ib,order)
		write(1,25) t,ord,en,chi,ch
		write(6,25) t,ord,en,chi,ch
25		format(1x,1p,5(e12.5,1x))
20	continue
	close(1)
	if(ans.eq.yes) then
		open(2,file=of2)
		rewind(2)
		dn=1.d0/dfloat(nedge**2)
		do 50 m=1,m1
			write(2,55) m,order(m)*dn
50		continue
55		format(1x,i6,1x,1p,e12.5)
		close(2)
	endif
	time1=etime(time)
	write(6,45) time1-time0
45	format(/1x,'CPU time = ',1p,e12.5,' sec')
	stop
	end
c
c
c
	subroutine initialize(t,nedge,isd,init,ms0,ib,order)
	parameter (maxn=101,maxt=100001)
	implicit real*8 (a-h,o-z)
	real ranf
	common ispin(0:maxn,0:maxn)
	dimension order(1)
	do 70 j=0,nedge+1
	ispin(0,j)=0
70	ispin(nedge+1,j)=0
	do 80 i=1,nedge
	ispin(i,0)=0
80	ispin(i,nedge+1)=0
	if(init.eq.1) then
		do 10 i=1,nedge
			do 15 j=1,nedge
				ispin(i,j)=1
15			continue
10		continue
		order(1)=nedge**2
	else
		order(1)=0.d0
		do 20 i=1,nedge
			do 25 j=1,nedge
				ispin(i,j)=1
				if(ranf(isd).gt.0.5) ispin(i,j)=-1
				order(1)=order(1)+ispin(i,j)
25			continue
20		continue
	endif
	do 30 m=2,ms0+1
		m1=m
		if(m1.gt.maxt) then
                   order(1)=order(m1-1)
                   m1=2
                endif
		order(m1)=order(m1-1)
		do 40 i=1,nedge
			i1=i-1
			if(ib.eq.1.and.i.eq.1) i1=nedge
			i2=i+1
			if(ib.eq.1.and.i.eq.nedge) i2=1
			do 50 j=1,nedge
				j1=j-1
				if(ib.eq.1.and.j.eq.1) j1=nedge
				j2=j+1
				if(ib.eq.1.and.j.eq.nedge) j2=1
				isum=ispin(i1,j)+ispin(i2,j)
     c					+ispin(i,j1)+ispin(i,j2)			
				de=dfloat(2*ispin(i,j)*isum)
				if(de.lt.0) then
					ispin(i,j)=-ispin(i,j)
					order(m1)=order(m1)+2*ispin(i,j)
				elseif(ranf(isd).le.dexp(-de/t)) then
					ispin(i,j)=-ispin(i,j)
					order(m1)=order(m1)+2*ispin(i,j)
				endif
50			continue
40		continue
30	continue
	return
	end
c
c
c
	subroutine calc(t,ord,en,chi,ch,nedge,m1,ms,ns,isd,ib,order)
	parameter (maxn=101,maxt=100001)
	implicit real*8 (a-h,o-z)
	real ranf
	common ispin(0:maxn,0:maxn)
	dimension order(1)
c
c initialize sum,sum2,sume,sume2 for the quantities to compute
c
	sum=0.d0
	sum2=0.d0
	sume=0.d0
	sume2=0.d0
	do 10 k=1,ns
	xmag=0.d0
	xe=0.d0
	if(k.eq.1) then
		m2=1
	else
		m2=ms+1
	endif
	do 30 m=1,m2
		if(m.lt.m2) then
			m1=m1+1
		        if(m1.gt.maxt) then
                           order(1)=order(m1-1)
                           m1=2
                        endif
			order(m1)=order(m1-1)
		endif
		do 40 i=1,nedge
			i1=i-1
			if(ib.eq.1.and.i.eq.1) i1=nedge
			i2=i+1
			if(ib.eq.1.and.i.eq.nedge) i2=1
			do 50 j=1,nedge
				j1=j-1
				if(ib.eq.1.and.j.eq.1) j1=nedge
				j2=j+1
				if(ib.eq.1.and.j.eq.nedge) j2=1
				isum=ispin(i1,j)+ispin(i2,j)
     c					+ispin(i,j1)+ispin(i,j2)			
				if(m.lt.m2) then
				   de=dfloat(2*ispin(i,j)*isum)
				   if(de.lt.0) then
					ispin(i,j)=-ispin(i,j)
					order(m1)=order(m1)+2*ispin(i,j)
				   elseif(ranf(isd).le.dexp(-de/t)) then
					ispin(i,j)=-ispin(i,j)
					order(m1)=order(m1)+2*ispin(i,j)
				   endif
				else
				   xmag=xmag+dfloat(ispin(i,j))
				   xe=xe-0.5d0*dfloat(ispin(i,j)*isum)
				endif
50			continue
40		continue
30	continue
c
c compute mag,energy add these to the sum,sum2
c
	sum=sum+xmag
	sum2=sum2+xmag**2
	sume=sume+xe
	sume2=sume2+xe**2
10	continue
c
c average and compute mag,energy,chi,ch and put them in the common
c
	ord=sum/dfloat(ns)
	chi=(sum2/dfloat(ns)-ord**2)/(t*dfloat(nedge**2))
	ord=ord/dfloat(nedge**2)
	en=sume/dfloat(ns)
	ch=(sume2/dfloat(ns)-en**2)/(t*nedge)**2
	en=en/dfloat(nedge**2)
	return
	end
c
c ***********************************************************************
c       Random number generator
c ***********************************************************************
c
      real function ranf(iseed)
      iseed=iseed*1566083941+1
      ranf=float(iseed)*2.328306e-10+0.5
      return
      end
