c    ===================================================================
        program mc
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Monte Carlo of 2d Ising magnet, constant-T ensemble            |
c    ===================================================================

c    ** Subprograms referenced **
        external start, runner, finish

c    -------------------------------------------------------------------

        call start
	call runner
	call finish

	stop
	end
	
		
	
c    ===================================================================
        subroutine start
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Reads in run parameters and initial configuration              |
c    |  Calculates and prints out initial values                       |
c    ===================================================================

c    ** Local variables **
        integer         i, j, ll

c    ** Subprograms referenced **
        external getintg, getreal, engall
		
c    -------------------------------------------------------------------

c    ** Output program identification **	
        print *, '-----------------------------------------------------'	
	print *, 'Monte Carlo simulation of 2d Ising magnet'
	print *, 'Constant-T ensemble'
	print *, 'Square periodic boundary conditions'

c    ** Input run parameters **	
        print *, '-----------------------------------------------------'	
	call getintg ( 'number of blocks     ', nblock )
	call getintg ( 'number of steps/block', nstep  )
	call getintg ( 'random number seed   ', seed   )
        call getreal ( 'temperature          ', temp   )

c    ** Input starting configuration **	
        print *, '-----------------------------------------------------'	
        print *, 'Reading initial configuration from file mc.old'
	open(unit=3,file='mc.old',status='old')
	read(3,*) ll
	if ( ll .ne. l ) stop 'Error in number of spins'
        do j = 1, l
	   read(3,*) (spin(i,j),i=1,l)
	enddo
	close(unit=3)

c    ** Output initial values **
	call engall 
        print *, '-----------------------------------------------------'	
	print *, 'Initial values'
	write(*,'(1x,''Energy/spin        = '',f10.4)') 
     :     real(eng)/real(n)
	write(*,'(1x,''Magnetization/spin = '',f10.4)') 
     :     real(mag)/real(n)
	
	return
	end
	
	
	
c    ===================================================================
        subroutine finish
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Calculates and prints out final values                         |
c    |  Writes out final configuration                                 |
c    ===================================================================

c    ** Local variables **
	integer i, j

c    ** Subprograms referenced **
        external engall
	
c    -------------------------------------------------------------------

c    ** Output final values **
	call engall 
        print *, '-----------------------------------------------------'	
	print *, 'Final values'
	write(*,'(1x,''Energy/spin        = '',f10.4)') 
     :    real(eng)/real(n)
	write(*,'(1x,''Magnetization/spin = '',f10.4)') 
     :    real(mag)/real(n)

c    ** Output final configuration **
        print

 *, '-----------------------------------------------------'	
        print *, 'Writing final configuration to file mc.new'
	open(unit=3,file='mc.new',status='unknown')
	write(3,'(1x,i5)') l
        do j = 1, l
	   write(3,'(1x,20i3)') (spin(i,j),i=1,l)
	enddo
	close(unit=3)
	
        print *, '-----------------------------------------------------'	
	print *, 'Program ends'
        print *, '-----------------------------------------------------'	

        return
	end



c    ===================================================================
        subroutine runner
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Conducts Monte Carlo run in nblock blocks of nstep steps       |
c    |  Each step consists of N attempted flips                        |
c    ===================================================================

C    ** Local variables **
        integer block, step, i, j, turn, nacc
        integer im, ip, jm, jp, spinn, spinij, delta, magnew

c    ** Subprograms referenced **
        real     ran
	external ran, engall
	external blkzer, blkacc, blkout
	external runzer, runacc, runout
	
c    -------------------------------------------------------------------

        write(*,'(/1x,50(''=''))')
        write(*,'(1x,''          Avg.Period  Mov.Rate'',
     :               ''    Energy  Magnetzn'')')
        write(*,'(1x,50(''=''))')
	
c    ** Zero run accumulators **
        call runzer

c    ** Perform simulation block-by-block **
        do block = 1, nblock

c       ** Recalculate total energy **
           call engall 
	   
c       ** Zero block accumulators **
           call blkzer

c       ** Perform required steps in block **
           do step = 1, nstep

              nacc = 0
	      
c          ** Attempt N flips  **
              do turn = 1, n

c             ** Pick spin at random **
                 i = 1 + ran(seed)*real(l)
                 if ( i .lt. 1 ) i = 1
                 if ( i .gt. l ) i = l
                 j = 1 + ran(seed)*real(l)
                 if ( j .lt. 1 ) j = 1
                 if ( j .gt. l ) j = l
                 spinij = spin(i,j)

c             ** Locate neighbours & sum spins**
                 im = i-1
                 if ( im .lt. 1 ) im = l
                 ip = i+1
                 if ( ip .gt. l ) ip = 1
                 jm = j-1
                 if ( jm .lt. 1 ) jm = l
                 jp = j+1
                 if ( jp .gt. l ) jp = 1
                 spinn = spin(i,jp) + spin(i,jm) 
     :                 + spin(ip,j) + spin(im,j)

c             ** Compute energy change and new magnetization **
                 delta = 2 * ( spinij * spinn )
                 magnew = mag - 2*spinij

c             ** Check for acceptance **
		 if ( exp(-real(delta)/temp) .gt. ran(seed) ) then
c                ** ... accept **
                    spin(i,j) = - spinij
                    eng  = eng + delta
		    mag  = magnew
                    nacc = nacc + 1
                 endif
             
              enddo
c          ** End of loop over spins **

c          ** Values for this step **
              mratio = real(nacc)/real(n)
	      
c          ** Accumulate block averages **
              call blkacc

           enddo
c       ** End of loop over steps in block

c       ** Normalize and print block averages **
           call blkout ( block )

c       ** Accumulate run averages **
           call runacc

        enddo
c    ** End of loop over blocks

c    ** Normalize and print run averages **
        write(*,'(1x,50(''=''))')
        call runout
        write(*,'(1x,50(''='')/)')
	
	return
	end


c    ===================================================================
        subroutine engall 
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Calculates energy of entire configuration                      |
c    ===================================================================

c    ** Local variables **
	integer i, j, ip, jp, spinij

c    -------------------------------------------------------------------
	
c    ** Zero energy **
        eng = 0
	mag = 0

c    ** Loop over spins **
        do j = 1, l
           do i = 1, l
              spinij = spin(i,j)
              mag = mag + spinij

c          ** Get neighbour spins and add in **
              ip = i+1
              if ( ip .gt. l ) ip = 1
              jp = j+1
              if ( jp .gt. l ) jp = 1
              eng = eng - spinij*(spin(i,jp)+spin(ip,j))
	      
	   enddo
	enddo
c    ** End of loop over spins **

	return
	end
	

c    ===================================================================
        subroutine blkzer
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Sets to zero all block accumulators                            |
c    ===================================================================

c    -------------------------------------------------------------------

        blkmr  = 0.0d0
        blkeng = 0.0d0
        blkmag = 0.0d0
	
	return
	end
	

c    ===================================================================
        subroutine runzer
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Sets to zero all run accumulators                              |
c    ===================================================================

c    -------------------------------------------------------------------

        runmr  = 0.0d0
        runeng = 0.0d0
        runmag = 0.0d0
	
        errmr  = 0.0d0
        erreng = 0.0d0
        errmag = 0.0d0
	
	return
	end

	

c    ===================================================================
        subroutine blkacc
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Increments all block accumulators                              |
c    ===================================================================

c    -------------------------------------------------------------------

        blkmr  = blkmr  + mratio
        blkeng = blkeng + real(eng)/real(n)
        blkmag = blkmag + real(mag)/real(n)
	
	return
	end

	

c    ===================================================================
        subroutine runacc
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Increments all run accumulators (including error estimators)   |
c    ===================================================================

c    -------------------------------------------------------------------

        runmr  = runmr  + blkmr 
        runeng = runeng + blkeng
        runmag = runmag + blkmag

        errmr  = errmr  + blkmr **2
        erreng = erreng + blkeng**2
        errmag = errmag + blkmag**2
	
	return
	end


	

c    ===================================================================
        subroutine blkout ( block )
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Normalizes and prints out block averages                       |
c    |  Supplied argument: block  .............. current block number  |
c    ===================================================================

c    ** Arguments **
        integer block
	
c    -------------------------------------------------------------------

        blkmr  = blkmr  / dble(nstep)
        blkeng = blkeng / dble(nstep)
        blkmag = blkmag / dble(nstep)
	
	write(*,'(1x,''Block  '',i3,i10,3f10.4)') 
     :     block, nstep, blkmr , blkeng, blkmag
     
	return
	end

	

c    ===================================================================
        subroutine runout
	implicit none
	include 'mc.inc'
c    ===================================================================

c    ===================================================================
c    |  Normalizes and prints out run averages with error estimates    |
c    ===================================================================

c    -------------------------------------------------------------------

        runmr  = runmr  / dble(nblock)
        runeng = runeng / dble(nblock)
        runmag = runmag / dble(nblock)

        errmr  = errmr  / dble(nblock)
        erreng = erreng / dble(nblock)
        errmag = errmag / dble(nblock)
        errmr  = ( errmr  - runmr **2 ) / dble(nblock)
        erreng = ( erreng - runeng**2 ) / dble(nblock)
        errmag = ( errmag - runmag**2 ) / dble(nblock)
	if (errmr .gt.0.0d0) errmr =sqrt(errmr )
	if (erreng.gt.0.0d0) erreng=sqrt(erreng)
	if (errmag.gt.0.0d0) errmag=sqrt(errmag)
	
	write(*,'(1x,''Run avg.  '',i10,3f10.4)') 
     :     nblock*nstep, 
     :     runmr , runeng, runmag
	write(*,'(1x,''Error est.'',i10,3f10.4)') 
     :     nblock*nstep,
     :     errmr , erreng, errmag

	return
	end



c    ===================================================================
        subroutine getintg ( string, value )
	implicit none
c    ===================================================================

c    ===================================================================
c    |  Prompts for an integer value with a specified string           |
c    |  Echoes the supplied value in the output file                   |
c    |  Supplied arguments: string        .... string used in prompt   |
c    |  Returned values:    value         .... the input value         |
c    ===================================================================

c    ** Arguments **
        character*(*) string
	integer       value
	
c    -------------------------------------------------------------------

        print *, 'Enter ',string
	read  *, value
	print *, 'Value of ',string,' = ', value
	return
	end
	


c    ===================================================================
        subroutine getreal ( string, value )
	implicit none
c    ===================================================================

c    ===================================================================
c    |  Prompts for an real value with a specified string              |
c    |  Echoes the supplied value in the output file                   |
c    |  Supplied arguments: string        .... string used in prompt   |
c    |  Returned values:    value         .... the input value         |
c    ===================================================================

c    ** Arguments **
        character*(*) string
	real          value
	
c    -------------------------------------------------------------------

        print *, 'Enter ',string
	read  *, value
	print *, 'Value of ',string,' = ', value
	return
	end