implicit none
      integer i,j,d,N,L,step,dU,MCS,ii,f,licznik,a
      parameter(L=40,MCS=600000)
      integer S(L,L),ni(L),pi(L)
      real r,ran1,k,T,P,w,m,suma,ms,ms2,ms4,x,E,Es,C,Es2,y,U
      real ex(5)

      d=-1
      open(1,file='L=40B.txt')
      open(2,file='data.txt')

      do i=1,L
       ni(i)=i+1
       pi(i)=i-1
      enddo
      ni(L)=1
      pi(1)=L

      do i=1,L
        do j=1,L
          k=ran1(d)
          if(k.lt.0.5) then
          S(i,j)=-1
          elseif(k.ge.0.5) then
          S(i,j)=1
          endif
        enddo
      enddo

      !T=1.68
      DO a=1,51
      !x=-5.95+(a-1)*0.16
      !T=2.269*(1-exp(x)/L)
      T=2.0+(a-1)*0.007
      licznik=0
      ms=0
      ms2=0
      ms4=0
      Es=0
      Es2=0

      do i=1,5
        ex(i)=exp(-(-8+(i-1)*4)/T)
      enddo

      do ii=1,MCS
       suma=0
       E=0

       do i=1,L
        do j=1,L
          dU=2*S(i,j)*(S(ni(i),j)+S(pi(i),j)+S(i,ni(j))+S(i,pi(j)))
          f=(dU+8)/4+1
          w=min(1.0,ex(f))
          P=ran1(d)
            if(P.le.w) then
              S(i,j)=-S(i,j)
            endif
           suma=suma+S(i,j)
           E=E+0.5*S(i,j)*(S(ni(i),j)+S(pi(i),j)+S(i,ni(j))+S(i,pi(j)))
        enddo
       enddo

       if(ii.gt.30000) then
         if(mod(ii,100).eq.0) then
          licznik=licznik+1
          m=suma/(L*L)
          ms=ms+abs(m)
          ms2=ms2+m**2
          ms4=ms4+m**4
          E=E/(L*L)
          Es=Es+E
          Es2=Es2+E**2
          !write(1,*) ii, m
          !write(*,*) ii, m
         endif
       endif

      !if(ii.lt.30) then         !termalizacja
         !m=suma/(L*L)
          !write(1,*) ii, m
          !write(*,*) ii, m
       !endif

      enddo

      ms=ms/float(licznik)
      ms2=ms2/float(licznik)
      ms4=ms4/float(licznik)
      x=L*L*(ms2-ms**2)/T
      Es=-Es/float(licznik)
      Es2=Es2/float(licznik)
      C=(Es2-Es**2)/(T*T*L*L)
      !x=log(abs(1-T/2.269)*L)
      !y=log(ms*L**0.125)
      U=1-ms4/(3*(ms2**2))

      write(1,*) T, U
      write(*,*) T, U

      ENDDO

        !do j=1,L
          !write(2,1000)(S(i,j),i=1,200)
        !enddo
1000        !format(1x,200I4)

      close(1)
      close(2)

      read(*,*)
      end

      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM

      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END