implicit none
      integer i, j, k, n, ngolfer, ntourney, length, burn, ni(605),l
      integer ntot,id(6000)
      double precision tourney(6000), score(6000)
      double precision asig, bsig, atau, btau
      double precision mmu, s2mu, theta(605), sumy(605)
      double precision mu, sig2, tau2, sumtheta, ssq2, ssq1
      double precision astar, bstar, mustar, sigstar
      real zero, one, gennor, gengam
      external gennor, gengam, setall

      n = 5700
      ngolfer = 604
      open(1,FILE="fortrangolf.dat",STATUS="OLD")
      open(7,FILE="fortranout.txt",STATUS="UNKNOWN")

      call setall(12434,523)
      do i=1,n
         read(1,*) id(i), score(i), tourney(i)
      end do

      close(1)
      zero = 0.0
      one = 1.0
      length = 50000
      burn = 10000

c     priors

      asig = 83.0d0
      bsig = 0.0014d0
      atau = 18.0d0
      btau = 0.015d0
      mmu = 72.0d0
      s2mu = 2.0d0

c     starting values
c      sumy=0.0d0
      ssq1=0.0d0
      ssq2=0.0d0
      
      do i=1,ngolfer
         theta(i) = 72.0d0
      end do
      mu = 72.0
      sig2 = 9.0d0
      tau2 = 6.0d0
      
      do i=1,ngolfer
         sumy(i)=0.0d0
         ni(i)=0
        
      end do
      do i=1,n
         sumy(id(i))=sumy(id(i))+score(i)
         ni(id(i))=ni(id(i))+1
         
      end do
      ntot=length+burn
      do i=1,ntot
         sumtheta=0.0d0
         do j=1,ngolfer
            mustar=(sumy(j)*tau2+mu*sig2)/(ni(j)*tau2+sig2)
            sigstar=(tau2*sig2)/(ni(j)*tau2+sig2)
            theta(j)=mustar+dsqrt(sigstar)*dble(gennor(zero,one))
            sumtheta=sumtheta+theta(j)
         end do
         mustar=(sumtheta*s2mu+mmu*tau2)/(ngolfer*s2mu+tau2)
         sigstar=s2mu*tau2/(ngolfer*s2mu+tau2)
         mu=mustar+dsqrt(sigstar)*dble(gennor(zero,one)) 
         ssq1=0.0d0
         ssq2=0.0d0
         do j=1,n
            ssq1=ssq1+(score(j)-theta(id(j)))**2
         end do
         do j=1,ngolfer
            ssq2=ssq2+(theta(j)-mu)**2
         enddo
         
         astar=real(asig+n/2.0d0)
         bstar=real(1.0d0/bsig+ssq1/2.0d0)
         sig2=1.0d0/dble(gengam(bstar,astar))
         
         astar=real(atau+ngolfer/2.0d0)
         bstar=real(1.0d0/btau+ssq2/2.0d0)
         tau2=1.0d0/dble(gengam(bstar,astar))
         
         if(i.gt.burn) write(7,*), (theta(l), l=1,ngolfer), mu, tau2, 
     c        sig2

      end do
      stop
      end