      PROGRAM GRSOS
        IMPLICIT   NONE
        INTEGER    KDIM,N,M,NRUN,KRC,IN,I,J,IMK,KK,KK250
c--        PARAMETER  (NRUN=1,KDIM=800,M=1,N=10,KRC=KDIM/2)
        PARAMETER  (NRUN=1,KDIM=1000,M=1,N=5,KRC=KDIM/2)
        REAL*8     SHR(KRC,N)
        PARAMETER  (KK=KDIM**2,KK250=KDIM**2+250)
        INTEGER    JJ1,JJ2,II,JJ
        INTEGER    IJKL
        INTEGER    KP1,KM1,KP2,KM2
        INTEGER    K1,K2,LMAX,IND,L
        INTEGER    KP(KDIM+1),KM(KDIM)
        INTEGER    KPH1,KMH1,KPH2,KMH2
        REAL       RKRN
        REAL*8     DSEED
        INTEGER    IZ(KDIM,KDIM),INSEED(KK250),INDUM(KK250)
        INTEGER    IZI
        REAL       R1(KK),R2(KK),Y
        REAL*8     SMX,SMM(N),RHE,RKK,SNRR,SWW(N),SCX,SHH
        INTEGER    KCHECK(KDIM+1)
c
        common /area1/ IZ,INSEED,INDUM
        common /area2/ r1,r2
                                            
        EXTERNAL R250,RKRN,HTWD,CORR1   

        call initim()
c     open files
        OPEN(UNIT=11, FILE='d1000.dat',STATUS='unknown')
        OPEN(UNIT=12, FILE='c1000.dat',STATUS='unknown')

C       NORMALIZATION FACTORS
C
        SNRR=1.D0/DFLOAT(NRUN)
        RKK=1.D0/DFLOAT(KK)
C
C       RKK = INVERSE NUMBER OF CO-ORDINATES
C      
C====================================================================
C        INITIALIZE FIRST 250 SEEDS FOR R250 RANDOM                  
C        NUMBER GENERATOR                                            
C--------------------------------------------------------------------
C                              
                                                 
        DSEED = 128071.
C                                                                               
C     **INITIALIZE THE SEED ARRAY FOR R250 USING SUBROUTINE RKRN     
C     **LARGEST INTEGER IN CRAY XM-P IS 2**46 - 1, BUT IN CONVEX
C     ** 2**36 - 1
                                                                                
       DO 1 I=1,KK250
       INDUM(I)=INT(2147483647.D0 * RKRN(DSEED)) 
       IF(INDUM(I).LT.0.0) THEN
            WRITE(6,*)'ERROR IN RKRN: ',DSEED,INDUM(I)          
            STOP
       ENDIF 
1      CONTINUE                                                                 
C                             
C      SET UP TO SHUFFLE THE SEED ARRAY AND PUT THE RESULTS TO INSEED
C
       LMAX = KK250
       IND = 0

88     Y = RKRN(DSEED)
       L = Y*KK250 + 1
       IF ( L .GT. LMAX ) GOTO 88

       IND = IND + 1
       INSEED(IND) = INDUM(L)
       INDUM(L) = INDUM(LMAX)
       LMAX = LMAX - 1
       IF ( LMAX .EQ. 1 ) THEN
                          INSEED(KK250) = INDUM(1)
                          GOTO 888
       ENDIF
       GOTO 88

888     CONTINUE
C
C      *******************************************************
C
C  SET VECTORS FOR PERIODIC BOUNDARY CONDITIONS
C
        call shotim()
         DO 3 I=1,KDIM
         KP(I)=I+1
         KM(I)=I-1
         KCHECK(I)=I
3        CONTINUE
         KP(KDIM)=1
         KP(KDIM+1)=1
         KM(1)=KDIM    
C
         DO 7 I=1,N
         SMM(I)=0.D0
         SWW(I)=0.D0
7        CONTINUE                                                      
C
         DO 9876 IJKL=1,NRUN

C====================================================================
C   INITIALIZE THE LATTICE

         DO 5 JJ1=1,KDIM
         DO 5 JJ2=1,KDIM
         IZ(JJ1,JJ2)=0  
5        CONTINUE
C====================================================================
C
C====================================================================
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C          MONTE CARLO -LOOPS STARTS
C
C   SAMPLING LOOP; N=NUMBER OF SAMPLES:
C
          DO 70 IN=1,N     
C                      
C   SAMPLING INTERVAL LOOP; M=SAMPLING INTERVAL:
C
          DO 80 IMK=1,M
C
C   GENERATION OF RANDOM NUMBERS; ALTOGETHER KK NUMBERS
C
C     *******************************************
          CALL R250(INSEED,R1,KK) 
          CALL R250(INSEED,R2,KK) 
C     *******************************************
C
          DO 80 JJ2=1,KK 

C
C     LATTICE LOOP :
C
C--------------------------------------------------------------------
C   RANDOMLY CHOSEN CO-ORDINATES
C
          K1=KDIM*R1(JJ2)+1
          K2=KDIM*R2(JJ2)+1
          K1=KCHECK(K1)
          K2=KCHECK(K2)        

C--------------------------------------------------------------------
C   AND ITS NEIGHBORS
C
        KP1=KP(K1)
        KM1=KM(K1)
        KP2=KP(K2)
        KM2=KM(K2)

        IZI=IZ(K1,K2)+1

        KPH1=IABS(IZ(KP1,K2)-IZI)
        IF(KPH1.LE.1)THEN
          KMH1=IABS(IZ(KM1,K2)-IZI)
          IF(KMH1.LE.1)THEN
             KPH2=IABS(IZ(K1,KP2)-IZI)
             IF(KPH2.LE.1)THEN          
                KMH2=IABS(IZ(K1,KM2)-IZI)
                IF(KMH2.LE.1)THEN
                         IZ(K1,K2)=IZI
                        ENDIF
                      ENDIF
                    ENDIF               
                  ENDIF

C====================================================================
80      CONTINUE
   
        CALL HTWD(IZ,SMM,SWW,IN,RKK,KDIM)
        CALL CORR1(IZ,SHR,IN,KDIM,KRC)

        call shotim()
70      CONTINUE
        WRITE(6,*)IJKL
9876    CONTINUE
C
C====================================================================
C          MONTE CARLO -LOOP ENDS
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C       
C
C  * CALCULATE THE AVERAGES AND WRITE THE DATA INTO FILES
C
        
        DO 168 I=1,N
        SMM(I)=SMM(I)*SNRR
        SWW(I)=SWW(I)*SNRR
        WRITE(11,*)I*M,SWW(I),SMM(I)
 168    CONTINUE

        CLOSE(11)

        WRITE(6,*)'FINAL N'
        WRITE(6,*)N

        DO 169 J=1,N
        DO 169 I=1,KRC
        SHR(I,J)=SHR(I,J)*SNRR
        WRITE(12,*)I,SHR(I,J)
169     CONTINUE       
        CLOSE(12)

        STOP
        END
