           program parisi
           integer spin(4,4),sb(2,2)
           integer s1,s2,s3
           real*8 k1,k2,k3,kp1,kp2,kp3
           real*8 b1,b2,b3,b4,b
           real*8 factor(2,2,-1:1)
            

           write(*,*) 'give the couplings k1,k2,k3'
           read*, k1,k2,k3

C     initalize the boltzmann factors for the blocks:
         b1=0.0D0
         b2=0.0D0
         b3=0.0D0
         b4=0.0D0
C     loop over all configurations
         do  j1=-1,1,2
                        spin(1,1)=j1
         do  j2=-1,1,2
                        spin(1,2)=j2
         do  j3=-1,1,2
                        spin(1,3)=j3
         do  j4=-1,1,2
                        spin(1,4)=j4
         do  j5=-1,1,2
                        spin(2,1)=j5
         do  j6=-1,1,2
                        spin(2,2)=j6
         do  j7=-1,1,2
                        spin(2,3)=j7
         do  j8=-1,1,2
                        spin(2,4)=j8
         do  j9=-1,1,2
                        spin(3,1)=j9
         do  j10=-1,1,2
                        spin(3,2)=j10
         do  j11=-1,1,2
                        spin(3,3)=j11
         do  j12=-1,1,2
                        spin(3,4)=j12
         do  j13=-1,1,2
                        spin(4,1)=j13
         do  j14=-1,1,2
                        spin(4,2)=j14
         do  j15=-1,1,2
                        spin(4,3)=j15
         do  j16=-1,1,2
                        spin(4,4)=j16

           s1=0
           s2=0
           s3=0

C       summation over the lattice
           do i1=1,4
             ir1=i1+1
C          periodic boundary conditions:
             if(ir1.eq.5) ir1=1
             il1=i1-1
             if(il1.eq.0) il1=4
             do i2=1,4
               ir2=i2+1
C          periodic boundary conditions:
               if(ir2.eq.5) ir2=1

               s1=s1+spin(i1,i2)*(spin(ir1,i2)+spin(i1,ir2))
               s2=s2+spin(i1,i2)*(spin(ir1,ir2)+spin(il1,ir2))
               s3=s3+spin(i1,i2)*spin(ir1,i2)*spin(i1,ir2)*spin(ir1,ir2)

             enddo
           enddo

C           the Boltzmann:

               b=exp(k1*s1+k2*s2+k3*s3)

C           compute the blocks:  
               sb(1,1)=spin(1,1)+spin(1,2)+spin(2,1)+spin(2,2)          
               sb(1,2)=spin(1,3)+spin(1,4)+spin(2,3)+spin(2,4)          
               sb(2,1)=spin(3,1)+spin(3,2)+spin(4,1)+spin(4,2)          
               sb(2,2)=spin(3,3)+spin(3,4)+spin(4,3)+spin(4,4)     

               do ii1=1,2
               do ii2=1,2
               if(sb(ii1,ii2).lt.0) then
                 factor(ii1,ii2,-1)=1.0D0
                 factor(ii1,ii2,1)=0.0D0  
               else if(sb(ii1,ii2).eq.0) then
                 factor(ii1,ii2,-1)=0.5D0
                 factor(ii1,ii2,1)= 0.5D0
               else
                 factor(ii1,ii2,-1)=0.0D0
                 factor(ii1,ii2,1)= 1.0D0
               endif
               enddo
               enddo   


C           compute the Boltzmann factors for the block configs:
C           +  +           
C           +  +
          b1=b1+factor(1,1,1)*factor(1,2,1)
     &         *factor(2,1,1)*factor(2,2,1)*b
C           +  +           
C           +  -
          b2=b2+factor(1,1,1)*factor(1,2,1)
     &         *factor(2,1,1)*factor(2,2,-1)*b
C           +  -           
C           +  -
          b3=b3+factor(1,1,1)*factor(1,2,-1)
     &         *factor(2,1,1)*factor(2,2,-1)*b
C           -  +           
C           +  -
          b4=b4+factor(1,1,-1)*factor(1,2,1)
     &         *factor(2,1,1)*factor(2,2,-1)*b


         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo

C        compute the  new couplings form the block boltzmann factors:

C  eq(1)      b1=const*exp(+kp1*8+kp2*8+kp3*4)
C  eq(2)      b2=const*exp(+kp1*0+kp2*0-kp3*4)
C  eq(3)      b3=const*exp(+kp1*0-kp2*8+kp3*4)
C  eq(4)      b4=const*exp(-kp1*8+kp2*8+kp3*4)

C      take the logarithm and solve with respect to ln(const),kp1,kp2,kp3

         b1=log(b1)
         b2=log(b2)
         b3=log(b3)
         b4=log(b4)

        kp1=(b1-b4)/16.0D0
        kp2=(b4+8.0D0*kp1-b3)/16.0D0
        kp3=( b1-8.0D0*(kp1+kp2)-b2)/8.0D0
        const=b2+4.0D0*kp3



         write(*,*) 'kp1,kp2,kp3'
         write(*,*) kp1,kp2,kp3
         k1=kp1
	 k2=kp2
	 k3=kp3

        end
