PROGRAM MCMC
C  This is a driver program for weighted Markov chain Monte Carlo 
C  simulation of the posterior distribution.
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (LPAR=250)
      COMMON/AMCMC/WT,PRI,PDAT,PAR(LPAR),ITER,NPAR,LRECRD,LWRITE,LERR
      COMMON/BMCMC/NITER,NSKIP
      LOGICAL LRECRD,LWRITE,LERR
      COMMON/FILES/HFILE(4)
      CHARACTER*40 HFILE
      CHARACTER*75 HDESC
      REAL*4 CPSEC
      CALL NOWPR(6)
      READ(5,*,END=1090,ERR=1100)NITER,NSKIP,NWRITE
      CALL PINTI('Number of iterations',NITER,1,2147483647)
      CALL PINTI('Skips between iteerations',NSKIP,0,NITER)
      LWRITE=.TRUE.
      IF(NWRITE.EQ.0)LWRITE=.FALSE.
      NS=NSKIP+1
      MITER=NITER*NS
      CALL INIT
      CALL PD0
      WRITE(6,5)NPAR,NITER,NSKIP
    5 FORMAT(/,'  Markov chain Monte Carlo posterior simulation',/,
     1       4X,'Parameters written in each record:',I10,/,
     2       27X,'Iterations:',I10,/,
     3       13X,'Skips between iterations:',I10)
      IF(LWRITE)CALL FOPEN('Posterior simulation matrix',3,.TRUE.)
      READ(5,10,END=50)HDESC
   10 FORMAT(A75)
      WRITE(6,20)HDESC
   20 FORMAT(/,'  Description:',/,2X,A75)
   30 READ(5,10,END=50)HDESC
      WRITE(6,40)HDESC
   40 FORMAT(2X,A75)
      GO TO 30
   50 IF(.NOT.LWRITE)GO TO 60
      CALL PINTI('Number of parameters',NPAR,1,LPAR)
      WRITE(3,55)NITER,NPAR
   55 FORMAT(2I15)
   60 TIME=CPSEC()
      ITER=0
      DO 120 LITER=1,NITER
      DO 110 ISKIP=1,NS
      ITER=ITER+1
      LRECRD=.FALSE.
      IF(ISKIP.EQ.NS)LRECRD=.TRUE.
      CALL CPD
  110 IF(LERR)GO TO 1010
  120 IF(LWRITE)WRITE(3,130)ITER,WT,PRI,PDAT,(PAR(IPAR),IPAR=1,NPAR)
  130 FORMAT(I16,3(1PE16.7),/,(5(1PE16.7)))
  140 TIME=CPSEC()
      WRITE(6,150)ITER,TIME,TIME/ITER
  150 FORMAT(/,'   Iterations completed:',I12,/,
     1         '               CPU time:',F12.4,/,
     2         '  CPU time per iteraton:',F12.4)
      CALL PD1
      CALL NOWPR(6)
      STOP
 1010 WRITE(6,1020)ITER,(PAR(IPAR),IPAR=1,NPAR)
 1020 FORMAT(/,' Fatal error from CPD on iteration',I10,/,
     1       ' Parameter vector:',/,(1X,5(1PE15.6)))
      GO TO 140
 1090 CALL PEND('Unexpected end of file, line 1 of control file')
 1100 CALL PEND('Read error, line 1 of control file')
      END