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

c    ===================================================================
c    |  Initial configuration of 2d Ising magnet                       |
c    |  Allows user to prescribe magnetization                         |
c    ===================================================================

c    ** Local variables **
        integer         i, j, flip, nflip
	real            mags, amags
        real            tol
        parameter     ( tol = 1.0e-6 )

c    ** External subprograms **
        real     ran
        external ran

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

        print *, 'System dimension = ', l
        print *, 'Number of spins  = ', n

        print *, 'Enter desired magnetization per spin'
        print *, 'Allowed range is -1..1'
        read  *, mags
        amags = abs(mags)
        if ( amags .gt. 1.0 ) stop 'Value outside range'
        nflip = nint(real(n)*(1.0-amags)/2.0)
        amags = 1.0 - 2.0*real(nflip)/real(n)
        if ( abs(amags-abs(mags)) .gt. tol ) then
           print *, 'Closest we can get is ', sign(amags,mags)
        else
           print *, 'OK'
        endif
        print *, 'Enter random seed'
        read  *, seed

        do i = 1, l
           do j = 1, l
              spin(i,j) = 1
           enddo
        enddo

        do flip = 1, nflip
10         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
           if ( spin(i,j) .lt. 0 ) goto 10
           spin(i,j) = -1
        enddo

        if ( mags .lt. 0.0 ) then
           do i = 1, l
              do j = 1, l
                 spin(i,j) = -spin(i,j)
              enddo
           enddo
        endif

        print *, 'Writing configuration to file mc.old'
	open(unit=3,file='mc.old',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)

        stop
        end