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

      Program Ising_3D

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

      L=10
      H=0.0
      V=1.0
      TauMin=4.0
      TauMax=5.0
      Nruns=6
      Nsweeps=1000

      Call Initialize(Lat,L)

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

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

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

      Do 3 i=1,Nsweeps

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

    3 Continue

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

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

    4 Continue

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

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

      Subroutine Sweep(Lat,L)

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

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

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

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

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

      End If

1     Continue
      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,L)
      Common/Params/H,V,Tau,M,E,Iran

      Iran=99991

      Open(Unit=1,File='Ising_3D.out')
      Write(1,*)' L, H, V =',L,H,V
      Write(1,*)
      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
      Do 1 k=1,L

      Lat(i,j,k)=1

1     Continue

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

      End

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