************************************************************************

      Program Ising_1

      Integer*2 Lat(10,10)
      Common/Params/H,V,Tau,M,E,Iran
      Real*8 SumM, SumE, SumM2, SumE2

      L=10
      H=0.1
      V=0.5
      TauMin=1.0
      TauMax=10.0
      Nruns=10
      Nmoves=4000

      Call Initialize(Lat,L)

      DelTau=(TauMax-TauMin)/(Nruns-1)
      Do 4 Tau=TauMin, TauMax, DelTau

      Do 2 i=1,100*L**2
    2 Call MakeMove(Lat,L)

      SumM=0
      SumM2=0
      SumE=0
      SumE2=0

      Do 3 i=1,Nmoves

      Call MakeMove(Lat,L)
      SumM=SumM+Abs(M)
      SumM2=SumM2+M**2
      SumE=SumE+E
      SumE2=SumE2+E**2

    3 Continue

      AveM=SumM/Nmoves
      AveM2=SumM2/Nmoves
      AveE=SumE/Nmoves
      AveE2=SumE2/Nmoves
      DelSqM=AveM2-(AveM)**2
      DelSqE=AveE2-(AveE)**2
      AveM=AveM/L**2
      AveE=AveE/L**2
      Chi=DelSqM/(Tau**2*L**2)
      SpecHeat=DelSqE/(Tau**2*L**2)

      Write(1,100)Tau, AveM, AveE, Chi, SpecHeat

    4 Continue

      Close(1)
100   Format(F6.2,4F18.6)
      End

************************************************************************

      Subroutine MakeMove(Lat,L)

      Integer*2 Lat(L,L)
      Common/Params/H,V,Tau,M,E,Iran

      i=1+L*Ran(Iran)
      j=1+L*Ran(Iran)

      DelE=2*Lat(i,j)*(H+V*(Lat(i,nm1(j,L))+Lat(i,np1(j,L))
     #+Lat(nm1(i,L),j)+Lat(np1(i,L),j)))

      If ((DelE.LT.0).OR.(Ran(Iran).LT.Exp(-DelE/Tau))) then

        Lat(i,j)=-Lat(i,j)
        M=M+2*Lat(i,j)
        E=E+DelE

      End If

      End

************************************************************************

      Function np1(n,L)
      If(n.EQ.L) then
        np1=1
      Else
        np1=n+1
      End If
      Return
      End

************************************************************************

      Function nm1(n,L)
      If(n.EQ.1) then
        nm1=L
      Else
        nm1=n-1
      End If
      Return
      End

************************************************************************

      Subroutine Initialize(Lat,L)

      Integer*2 Lat(L,L)
      Common/Params/H,V,Tau,M,E,Iran

      Iran=99991

      Open(Unit=1,File='Ising_1.out')
      Write(1,*)'L, H, V =',L,H,V
      Write(1,50)
   50 Format(//'   Tau',14X,'AveM',14X,'AveE',15X,'Chi',10X,'SpecHeat')
      Write(1,*)

      Do 1 i=1,L
      Do 1 j=1,L

      Lat(i,j)=1

1     Continue

      M=L**2
      E=-(H+2*V)*L**2

      End

************************************************************************