c
c***********************************************************************
c
c        program   : 2dnew
c
c        2d ising model
c
c        algorithm : swendsen and wang based on an idea by hu
c
c        version   : 0.1 as of september
c
c        dieter w. heermann
c
c***********************************************************************
c
      integer spin(0:30,0:30)
      integer lbond(1:30,1:30),ubond(1:30,1:30)
      integer stack(1:5000)
      integer nsd(1:1000)
      integer nsc(1:1000)
      integer nsa(1:1000)
      integer nsa1(1:1000)
c
      real    jkt
      real    mbond,sbond
c
c-----set and/or read the input parameter
c
      write(*,*) 'lattice size = '
      read(*,*) l
      write(*,*) 'temp   = '
      read(*,*) ttc
      write(*,*) 'mcsmax = '
      read(*,*) mcsmax
      write(*,*) 'iskip  = '
      read(*,*) iskip
      write(*,*) 'iseed  = '
      read(*,*) iseed
c
      write(*,*)
      write(*,*) 'program 2dnew2'
      write(*,*)
      write(9,*) 'simulation parameters'
      write(9,*) '---------------------'
      write(9,*)
      write(9,*) 'lattice dim        = ',l
      write(9,*) 'temperature        = ',ttc
      write(9,*) 'number of samples  = ',mcsmax
      write(9,*) 'iskip              = ',iskip
      write(9,*) 'seed for generator = ',iseed
      write(9,*)
c
c-----set the other parameters
c
      arg     = 1.0 + sqrt( 2.0 )
      jkt     = 0.5 * alog( arg ) / ttc
      pb      = 1.0 - exp( - 2.0 * jkt )
      write(9,*) 'pb is ',pb
      lsq     = l * l
      iseed   = 4 * iseed + 1
      rcum    = 0.0
      sus     = 0.0
c
      maxcla  = 0
      suscl   = 0.0
      prob    = 0.0
c
      mbond   = 0.0
      sbond   = 0.0
c
      stackm  = 5000
c
c-----set initial configuration ( random )
c
      do  10  j=1,l
       do  10  i=1,l
        r         = sunif(iseed) - 0.5
        spin(i,j) = sign(1.,r)
 10   continue
c
c-----set boundray bonds
c
      do  20  j=1,l
        lbond(j,1) = 0
        ubond(1,j) = 0
 20   continue
c
c        ***************************************************
c        *                                                 *
c        *           Monte Carlo Part                      *
c        *                                                 *
c        ***************************************************
c
      do  1000  mcs=1,mcsmax
c
c                                      generate the bonds
c
      cb = 0.0
      do  300  j=2,l
       do  300  i=1,l
        rb         = pb - sunif(iseed)
        ib         = sign(1.,rb)
        ib         = ( ib + 1 ) / 2
        ileft      = (spin(i,j-1) + spin(i,j)) / 2
        ileft      = abs(ileft)
        ib         = ib * ileft
        cb         = cb + ib
        lbond(i,j) = ib
 300  continue
      do  310  j=1,l
       do  310  i=2,l
        rb         = pb - sunif(iseed)
        ib         = sign(1.,rb)
        ib         = ( ib + 1 ) / 2
        iup        = (spin(i-1,j) + spin(i,j)) / 2
        iup        = abs(iup)
        ib         = ib * iup
        cb         = cb + ib
        ubond(i,j) = ib
 310  continue
      mbond = mbond + cb
      sbond = sbond + cb * cb
c
c                                      assign new random orientation
c                                      to the clusters
c
      icl = 0
      do  400  i=1,l
       do  400  j=1,l
c
        ileft = lbond(i,j)
        iup   = ubond(i,j)
        if ((ileft .eq. 1).and.(iup  .eq. 0)) spin(i,j) = spin(i,j-1)
        if ((iup   .eq. 1).and.(ileft.eq. 0)) spin(i,j) = spin(i-1,j)
        if ((ileft .eq. 1).and.(iup  .eq. 1)) then
c          search for the up label
           isrch = spin(i-1,j)
 401       inew  = stack(isrch)
           if ( inew .gt. 0 ) then
              isrch = inew
              goto 401
             else
              iu = isrch
            end if
c          search for the left label
           isrch = spin(i,j-1)
 402       inew  = stack(isrch)
           if ( inew .gt. 0 ) then
              isrch = inew
              goto 402
             else
              il = isrch
            end if
           mini        = min0(iu,il)
           maxi        = max0(iu,il)
           if ( mini .ne. maxi )  then
              stack(maxi) = mini
              spin(i-1,j) = iu
              spin(i,j-1) = il
             end if
           spin(i,j)   = mini
          end if
        if ((ileft.eq.0).and.(iup.eq.0)) then
           icl        = icl + 1
           if ( icl .gt. stackm ) stop 99
           r          = sunif(iseed) - 0.5
           inew       = sign(1.,r)
           stack(icl) = -2 + inew
           spin(i,j)  = icl
          end if
c
 400   continue
c
       rm = 0.0
       do  410  i=1,l
        do  410  j=1,l
           isrch = spin(i,j)
 420       inew  = stack(isrch)
           if ( inew .gt. 0 ) then
              isrch = inew
              goto 420
             else
              im         = inew + 2
              spin(i,j)  = im
              rm         = rm + im
              nsd(isrch) = nsd(isrch) + 1
            end if
 410   continue
c
       itest = (mcs/iskip) * iskip
       if ( itest .eq. mcs ) then
c                            compute the magnetic properties
       icount = icount + 1
       rm   = rm / lsq
       sus  = sus  + rm * rm
       rcum = rcum + abs(rm)
c      write(*,*) 'mc step ',mcs,' m = ',rm
c
c                            compute the cluster properties
c
       maxcl = 0
       do  500  i=1,icl
         ns     = nsd(i)
         nsd(i) = 0
         if ( ns .gt. 0 ) then
             nsc(ns) = nsc(ns) + 1
             nsa(ns) = nsa(ns) + 1
             if ( ns .gt. maxcl )  maxcl = ns
           end if
 500   continue
       prob = prob + maxcl
       if ( maxcl .gt. maxcla ) maxcla = maxcl
c
       cm = 0.0
       do  510  i=1,maxcl-1
         ns     = nsc(i)
         nsa1(i) = nsa1(i) + ns
         nsc(i) = 0
         if ( ns .ne. 0 ) then
c            write(9,*) i,ns
           end if
         cm = cm + i * i * ns
 510   continue
       cm = cm / lsq
       suscl = suscl + cm
c
c      write(9,*) 'mean cluster size (sus.) = ',cm
c
       end if
c
 1000  continue
c
c        ***************************************************
c        *                                                 *
c        *           A n a l y s i s                       *
c        *                                                 *
c        ***************************************************
c
       sus = sus - rcum * rcum /icount
       sus = sus / icount  
       write(9,*)
       write(9,*) 'results: program 2dnew2'
       write(9,*) '-----------------------'
       write(9,*)
       write(9,*) 'sample size              = ',icount
       write(9,*) 'average magnetization    = ',rcum/icount
       write(9,*) 'susceptibility           = ',sus * lsq
c
       prob = (prob / icount)/ lsq
       write(9,*) 'percolation probability  = ',prob
c
       suscl = suscl / icount
       write(9,*) 'average sus from cluster = ',suscl
c
       mbond = (mbond/lsq) / icount
       sbond = (sbond/lsq) / icount 
       write(9,*) 'mean number of bonds     = ',mbond
       write(9,*) 'fluctuation of bonds     = ',sbond
c
       write(9,*) 'average cluster size distribution'
       write(9,*)
       cm1 = 0.0
       cm2 = 0.0
       do  1510  i=1,maxcla-1
         ns     = nsa(i)
         if ( ns .ne. 0 ) then
             write(9,*) i,ns
             cm2 = cm2 + i * i * ns
           end if
 1510  continue
       cm2 = (cm2 / lsq) / icount
       write(9,*) 'mean cluster size (sus.) = ',cm2
       do  1520  i=1,maxcla-1
         ns     = nsa1(i)
         if ( ns .ne. 0 ) then
             write(9,*) i,ns
             cm2 = cm2 + i * i * ns
           end if
 1520  continue
       cm2 = (cm2 / lsq) / icount
       write(9,*) 'mean cluster size (sus.) = ',cm2
c
       stop
       end
c        mz_soft : random number generators
c
C**********************************************************************CSUN  600
C                                                                      CSUN  610
C                                                                      CSUN  620
C       0 , 1   - U N I F O R M  DISTRIBUTION                          CSUN  630
C                                                                      CSUN  640
C                                                                      CSUN  650
C**********************************************************************CSUN  660
C**********************************************************************CSUN  670
C                                                                      CSUN  680
C     FOR DETAILS SEE:                                                 CSUN  690
C                                                                      CSUN  700
C               AHRENS, J.H., DIETER, U. AND GRUBE, A.                 CSUN  710
C               PSEUDO-RANDOM NUMBERS:  A NEW PROPOSAL                 CSUN  720
C                     FOR THE CHOICE OF MULTIPLICATORS                 CSUN  730
C               COMPUTING, 6 (1970), 121 - 138                         CSUN  740
C                                                                      CSUN  750
C**********************************************************************CSUN  760
C                                                                       SUN  770
      REAL FUNCTION SUNIF(IR)                                           SUN  780
      DOUBLE PRECISION R,FACTOR,TWO28                                   SUN  790
C                                                                       SUN  800
C     FACTOR - INTEGER OF THE FORM 8*K+5 AS CLOSE AS POSSIBLE           SUN  810
C              TO  2**26 * (SQRT(5)-1)/2     (GOLDEN SECTION)           SUN  820
C     TWO28  = 2**28  (I.E. 28 SIGNIFICANT BITS FOR DEVIATES)           SUN  830
C                                                                       SUN  840
      DATA FACTOR /41475557.0D0/, TWO28 /268435456.0D0/                 SUN  850
C                                                                       SUN  860
C     RETURNS SAMPLE U FROM THE  0,1 -UNIFORM DISTRIBUTION              SUN  870
C     BY A MULTIPLICATIVE CONGRUENTIAL GENERATOR OF THE FORM            SUN  880
C        R := R * FACTOR (MOD 1) .                                      SUN  890
C     IN THE FIRST CALL R IS INITIALIZED TO                             SUN  900
C        R := IR / 2**28 ,                                              SUN  910
C     WHERE IR MUST BE OF THE FORM  IR = 4*K+1.                         SUN  920
C     THEN R ASSUMES ALL VALUES  0 < (4*K+1)/2**28 < 1 DURING           SUN  930
C     A FULL PERIOD 2**26 OF SUNIF.                                     SUN  940
C     THE PARAMETER IR IS USED ONLY IN THE FIRST CALL FOR               SUN  950
C     INITIALIZATION OF SUNIF. THEREAFTER (WHEN NEGATIVE)               SUN  960
C     IR BECOMES A DUMMY VARIABLE.                                      SUN  970
C                                                                       SUN  980
      IF (IR .GE. 0) GO TO 1                                            SUN  990
C                                                                       SUN 1000
C     STANDARD CASE:  SAMPLING                                          SUN 1010
C                                                                       SUN 1020
      R=DMOD(R*FACTOR,1.0D0)                                            SUN 1030
      SUNIF=SNGL(R)                                                     SUN 1040
      RETURN                                                            SUN 1050
C                                                                       SUN 1060
C     FIRST CALL: INITIALIZATION                                        SUN 1070
C                                                                       SUN 1080
1     R=DBLE(FLOAT(IR))/TWO28                                           SUN 1090
      R=DMOD(R*FACTOR,1.0D0)                                            SUN 1100
      SUNIF=SNGL(R)                                                     SUN 1110
      IR=-1                                                             SUN 1120
      RETURN                                                            SUN 1130
      END                                                               SUN 1140
C                                                                       SEX   10
C**********************************************************************CSEX   20
