!###############################BEGIN GMCMC############################# 
! General Markov Chain Monte Carlo routines
!  author: Guthrie Miller guthriemiller@gmail.com
!  reference: http://www.lulu.com/content/paperback-book/probabilistic-interpretation-of-data/13620353
      module GMCMC_MOD
      character*10, save :: tapefile='tape.txt'
! real*8 esumu or esumv is written on tape file in hex
! real*8 thetu or thetv is written on tape file in hex
!  one extra character to provide spacing and readability
      character*16, save :: cesum
      character*17, save, dimension(:), allocatable :: cthet
      character*39, save :: info ! date, time, size info about file
      character*39, save :: infoid ! file info about datafile
      character*80, save :: dataid  ! name of datafile
      logical, save:: accept,ranwalk,slow0,slow=.true.
      logical, save, dimension(:), allocatable :: bad
      integer, save :: ialg,icalcvx,icycle,i
     1 ,idorun,ifocus,ihypu,ihypv
     1 ,incr(2),incracpt(2),incrhyp(2,2),iprbmax,irecmax(2),irpt
     1 ,iseed,iv,ivmin,nburn,ndebug,ngrps,nhyp=1,nprm
     1 ,nprmsgrp,nrpts,nsumlike,ntape,nv,nvgood
      integer, save ,dimension(:), allocatable::iv1good,nprmgrp,iwrap
      integer, save, dimension(:,:), allocatable::iiprm
      integer*8, save :: irun,irunacpt,irundbug
     1 ,nrun,nrun0,nthin
      real*8, save :: burnfrac,c2(2,2),c2f,esumu,eminv
     1 ,eminv0,c2ndref=1,esumref,esumrefp,fracrw,phyp(2)
     1 ,proddthet,prodsum,pswitch=0.1d0,sumlike,xsec,xsec0
     1 ,xrand,One=1
      real*8 dranr250
      external dranr250
      real*8, save, dimension(:), allocatable :: esumv
     1 ,pfocus,pfocusn,pfocusc,prbnd,prbcd,thetu,dthet
      real*8, save, dimension(:,:), allocatable :: thetv
      end module GMCMC_MOD
!-----------------------------------------------------------------------
!  the user calls subroutine GMCMC below and needs to supply following five routines:
!                    calcvx        calculate extra v (candidate) variables xv,
!                                   which  are the natural variables of the problem,
!                                   from the fundamental theta variables
!                                   and calculate energy = chi**2/2
!                    movex         move xv variables to xu variables
!                                   when candidate is accepted and chain moved
!                    mvbackvx(iv1) restore xv variables from xu for candidate iv1
!                    increm        increment tallies
!                    zero          zero tallies
!  theta variables written to tapefile (TAPE.TXT),
!   which contains NTAPE records uniformly spaced throughout entire run
!   Other output also written to file GMCMC.TXT as well as to terminal
!   random numbers should be generated using e.g. xrand=dranr250(0)
!-----------------------------------------------------------------------
      subroutine GMCMC(nrunIN,ntapeIN,nrptsIN,ndebugIN
     1 ,mdataIN,iseedIN,icycleIN,nvIN,ialgIN
     1 ,slowIN,burnfracIN,fracrwIN
     1 ,nprmIN,ngrpsIN,nprmsgrpIN,nprmgrpIN,iiprmIN
     1 ,pfocusIN,iwrapIN,thetIN,dthetIN
     1 ,datafileIN)
! this program validates input parameters, initializes variables
!  and initiates run by calling runchain
!  arguments:
!   NRUN  number of chain iterations
!   NTAPE  number of chain iterations written on file tapefile
!   NRPTS  number of reports written out during the run
!   NDEBUG  for debugging, number of detailed print outs
!   MDATA  number of data points (for calculation of ChiSq/mDATA)
!   ISEED  random number seed
!   ICYCLE = 1 for parameter movement by regular cycling through groups
!          = 2 for randomly chosen next group
!          = 3 for randomly chosen with probability of group PFOCUS(IGRP)
!   NV  number of candidates
!   IALG = 1 for MRT algorithm
!        = 2 for B algorithm
!        = 3 for multiple candidate (MC) algorithm
!   SLOW  logical variable that if true uses slow but sure calculations
!    without tricky speedups
!   BURNFRAC  the fraction of total iterations that are ignored in tallies
!    because of possible initialization bias
!   FRACRW  fraction of time random walk is used rather than choosing
!    candidate from entire space
!   NPRM  number of parameters
!   NGRPS  number of parameter groups
!   NPRMSGRP maximum value of NPRMGRP(IGRP) for array allocation
!   NPRMGRP(IGRP), IGRP=1,NGRPS  number of parameters withing group IGRP
!   IIPRM(IPRMGRP,IGRP), IPRMGRP=1,NPRMGRP(IGRP) parameters withing group IGRP
!   PFOCUS(IGRP) if NGRPS>1, probability of focusing on group IGRP, when ICYCLE=3
!   IWRAP(IPRM)  wraping of parameter IPRM, 1 no wrapping, 2 cyclical wrapping
!    when parameter goes outside of theta interval (0,1)
!   THET(IPRM)  intial value of theta parameter IPRM
!   DTHET(IPRM)  random walk half width for parameter IPRM
!   DATAFILE  name of datafile, used with checksum to see if data have changed
!
!   These parameters for NHYP=2 can be changed by referencing them in GMCMC_MOD
!   NHYP = 1 for no alternate (reference) hypotheses (default)
!        = 2 for reference hypothesis
!   if NHYP = 2, these parameters used:
!     C2NDREF  equivalent Chisq/nData to give constant energy when in the reference hypothesis
!     PSWITCH  probability to have proposal switch hypothesis from normal to reference

      use GMCMC_MOD
      implicit none

      integer ialgIN,icycleIN,iseedIN,mdataIN,ndebugIN,ngrpsIN
     1 ,nprmIN,nprmsgrpIN,nrptsIN,ntapeIN,nvIN,ireturn
      integer iiprmIN(nprmsgrpIN,ngrpsIN),iwrapIN(nprmIN)
     1 ,nprmgrpIN(ngrpsIN)
      integer*8 nrunIN
      logical slowIN
      real*8 burnfracIN,fracrwIN
      real*8 dthetIN(nprmIN),pfocusIN(ngrpsIN),thetIN(nprmIN)

      character*(*) datafileIN

      character*80 ermess1,ermess2
      logical ex
      integer ih,igrp,iprm,iprmgrp,iv1

      open(1,file='gmcmc.out')
      write(*,*)'-------------------begin GMCMC---------------------'
      write(1,*)'-------------------begin GMCMC---------------------'

      if(iseedIN.eq.0) then
         write(ermess1,*)'iseed=',iseedIN
         ermess2='random number seed must be positive integer'
         go to 999
      endif
      iseed=iseedIN

      if(ndebugIN.gt.0) then
         write(*,*)'GMCMC:  ndebug=',ndebugIN
     1    ,' number of debug printouts'
         write(1,*)'GMCMC:  ndebug=',ndebugIN
     1    ,' number of debug printouts'
      endif
      ndebug=ndebugIN

      if(slowIN) then
         write(*,*)'GMCMC: REMINDER slow=',slowIN
     1    ,' means always use slow but sure calculations' 
         write(1,*)'GMCMC: REMINDER slow=',slowIN
     1    ,' means always use slow but sure calculations' 
      endif
      slow=slowIN
      slow0=slow
 
      if(burnfracIN.eq.0.or.burnfracIN.gt.1.) then
         write(ermess1,*)'burnfrac=',burnfracIN
         ermess2='burn fraction should be small fraction, i.e. 0.05'
         go to 999
      endif
      burnfrac=burnfracIN

      if(fracrwIN.lt.0.or.fracrwIN.gt.1.) then
         write(ermess1,*)'fracrw=',fracrwIN
         ermess2='fraction of time to random walk, between 0 and 1'
         go to 999
      endif
      fracrw=fracrwIN

      if(nrunIN.eq.0) then
         write(ermess1,*)'nrun=',nrunIN
         ermess2='run length, number of mcmc iterations, not set'
         go to 999
      endif
      nrun=nrunIN

      if(ntapeIN.eq.0.or.ntapeIN.gt.nrunIN) then
         write(ermess1,*)'ntape=',ntapeIN
         ermess2='number of iterations to write on tape bad'
         go to 999
      endif
      ntape=ntapeIN

      if(nrptsIN.eq.0.or.nrptsIN.gt.ntapeIN) then
         write(*,*)'GMCMC: WARNING  nrpts bad, nrpts=',nrptsIN
         write(1,*)'GMCMC: WARNING  nrpts bad, nrpts=',nrptsIN
      endif
      nrpts=nrptsIN

      if(mdataIN.eq.0) then
         write(ermess1,*)'mdata=',mdataIN
         ermess2='number of data points'
         go to 999
      endif

!  c2f is coefficient that converts energy (sum of squared residuals)
!   to chisq per data point
      c2f=2.d0/mdataIN
 
      if(nhyp.eq.2) then
         write(*,*)'GMCMC:  nhyp=',nhyp
     1    ,' zero-dimensional reference hypothesis'
         write(1,*)'GMCMC:  nhyp=',nhyp
     1    ,' zero-dimensional reference hypothesis'
         if(pswitch.lt.0.or.pswitch.gt.1) then
            write(ermess1,*)'ih,pswitch=',ih,pswitch
            ermess2='switch probabiltiy bad'
     1       //' should be >= 0 and <= 1'
            go to 999
         endif
         if(c2ndref.lt.0) then
            write(ermess1,*)'c2ndref=',c2ndref
            ermess2='need positive c2ndref=c2f*esumref'
            go to 999
         endif
         esumref=c2ndref/c2f
         esumrefp=esumref+log(pswitch) ! esumrefp=-log(P0/pswitch)
         write(*,*)'GMCMC:  c2ndref (Equivalent Chisq/nData) ='
     1    ,c2ndref,' esumref=',esumref
         write(1,*)'GMCMC:  c2ndref (Equivalent Chisq/nData) ='
     1    ,c2ndref,' esumref=',esumref
         write(*,*)' esumrefp=esumref+log(pswitch)=',esumrefp
     1    ,' Equivalent Chisq/nData=',c2f*esumrefp
         write(1,*)' esumrefp=esumref+log(pswitch)=',esumrefp
     1    ,' Equivalent Chisq/nData=',c2f*esumrefp
      elseif(nhyp.eq.1) then
      else
         write(ermess1,*)'nhyp bad, nhyp=',nhyp
         ermess2='nhyp = 1 or 2'
         go to 999
      endif

c     write(*,*)'datafile=',trim(datafile)
      if(trim(datafileIN).ne.'') then
         inquire(file=datafileIN,exist=ex)
         if(.not.ex) then
            write(ermess1,*)'datafile=',trim(datafileIN)
            ermess2="datafile doesn't exist"
            go to 999
         endif
      else
         write(*,*)'GMCMC: WARNING datafile='
     1    ,trim(datafileIN),' blank'
         write(1,*)'GMCMC: WARNING datafile='
     1    ,trim(datafileIN),' blank'
      endif
      dataid=datafileIN

c     write(*,*)'ngrps=',ngrps
      if(nprmIN.le.0) then
         write(ermess1,*)'nprm=',nprmIN
         ermess2='nprm bad'
         go to 999
      endif
      nprm=nprmIN

      if(nvIN.le.0) then
         write(ermess1,*)'nv=',nvIN
         ermess2='number of candidates bad'
         go to 999
      endif
      nv=nvIN

      if(ngrpsIN.le.0) then
         write(ermess1,*)'ngrps=',ngrpsIN
         ermess2='ngrps bad'
         go to 999
      endif
      ngrps=ngrpsIN
c      nogrps=.true.

      if(nprmsgrpIN.le.0) then
         write(ermess1,*)'nprmsgrp=',nprmsgrpIN
         ermess2='nprmsgrp bad'
         go to 999
      endif
      nprmsgrp=nprmsgrpIN

      if(allocated(iiprm))deallocate(iiprm)
      if(allocated(nprmgrp))deallocate(nprmgrp)
      if(allocated(thetv))deallocate(thetv)
      if(allocated(thetu))deallocate(thetu)
      if(allocated(cthet))deallocate(cthet)
      if(allocated(bad))deallocate(bad)
      if(allocated(dthet))deallocate(dthet)
      if(allocated(iwrap))deallocate(iwrap)
      if(allocated(iv1good))deallocate(iv1good)
      if(allocated(esumv))deallocate(esumv)
      if(allocated(prbnd))deallocate(prbnd)
      if(allocated(prbcd))deallocate(prbcd)

      allocate(iiprm(nprmsgrp,ngrps),nprmgrp(ngrps))
      allocate(thetv(nprm,nv),thetu(nprm),cthet(nprm),bad(nv))
      allocate(dthet(nprm),iwrap(nprm))
      allocate(iv1good(nv),esumv(nv))
      allocate(prbnd(nv+1),prbcd(nv+1))

c     write(*,*)'GMCMC: after allocate'

      if(ngrps.gt.1.and.
     1   (icycleIN.le.0.or.icycleIN.gt.3)) then
         write(ermess1,*)'ngrps,icycle='
     1    ,ngrps,icycleIN
         ermess2='icycle=1 (cycle), 2 (random)'//
     1    ', 3 (prob specified by pfocus)'
         go to 999
      endif
      icycle=icycleIN

      igrpLoop: do igrp=1,ngrps
c        write(*,*)'igrp,nprmgrp=',igrp,nprmgrp(igrp)
         if(nprmgrpIN(igrp).eq.0) then
            write(ermess1,*)'nprmgrp(igrp), igrp='
     1       ,nprmgrpIN(igrp),igrp
            ermess2='at least 1 parameter in any group'
            go to 999
         elseif(nprmgrpIN(igrp).gt.nprmsgrp) then
            write(ermess1,*)'nprmgrp(igrp), igrp, nprmsgrp='
     1       ,nprmgrpIN(igrp),igrp, nprmsgrp
            ermess2='too many parameters'
            go to 999
         endif
         nprmgrp(igrp)=nprmgrpIN(igrp)
 
         do iprmgrp=1,nprmgrp(igrp)
            if(iiprmIN(iprmgrp,igrp).lt.1.or.
     1       iiprmIN(iprmgrp,igrp).gt.nprm) then
               write(ermess1,*)'iiprm,iprmgrp,igrp,nprm='
     1          ,iiprmIN(iprmgrp,igrp),iprmgrp,igrp,nprm
               ermess2='iiprm needs to be between 1 and nprm'
               go to 999
            endif
            iiprm(iprmgrp,igrp)=iiprmIN(iprmgrp,igrp)

            iprm=iiprm(iprmgrp,igrp)
            if(iwrapIN(iprm).lt.1.or.iwrapIN(iprm).gt.2) then
               write(ermess1,*)'iwrap(iprm), iprm=',iwrapIN(iprm),iprm
               ermess2='iwrap = 1 (no wrap), or 2 (wrap)'
               go to 999
            endif
         enddo

      enddo igrpLoop

      iprmLoop: do iprm=1,nprm
         do igrp=1,ngrps
            do iprmgrp=1,nprmgrp(igrp)
               if(iprm.eq.iiprm(iprmgrp,igrp)) cycle iprmLoop
            enddo
         enddo
         write(ermess1,*)'iprm=',iprm
         ermess2='iprm does not appear in iiprm'
         go to 999
      enddo iprmLoop

      if(icycle.eq.3) then
         allocate(pfocus(ngrps),pfocusn(ngrps),pfocusc(ngrps))
         if(minval(pfocusIN(1:ngrps)).le.0) then      
            write(ermess1,*)'minval(pfocus).le.0, pfocus='
     1       ,pfocusIN(1:ngrps)
            ermess2='for icycle=3, pfocus needs to be >0'
            go to 999
         endif
         pfocus=pfocusIN(1:ngrps)
         pfocusn=pfocus/sum(pfocus) ! nomalized probability
         do igrp=1,ngrps
            pfocusc(igrp)=sum(pfocusn(:igrp))  ! cumulative probability
         enddo
         pfocusc(ngrps)=1 ! in case it isn't because of roundoff errors
      endif

      if(ialgIN.lt.1.or.ialgIN.gt.3) then
         write(ermess1,*)'ialg=',ialgIN
         ermess2='ialg = 1,2,3 for MRT, B, MC'
         go to 999
      endif
      if(ialgIN.eq.3.and.nhyp.eq.2) then
         write(*,*)'GMCMC: WARNING ialg=',ialgIN
     1    ,' should not be used with nhyp=2'
         write(1,*)'GMCMC: WARNING ialg=',ialgIN
     1    ,' should not be used with nhyp=2'
      endif
      ialg=ialgIN

!  initialize

      do iprm=1,nprm
         if(thetIN(iprm).lt.0.or.thetIN(iprm).gt.1) then
            write(ermess1,*)'thet(iprm), iprm=',thetIN(iprm),iprm
           ermess2='starting theta value required to be between 0 and 1'
            go to 999
         endif

         if(dthetIN(iprm).le.0.or.dthetIN(iprm).gt.1) then
            write(ermess1,*)'iprm,dthet(iprm)='
     1       ,iprm,dthetIN(iprm)
            ermess2='dthet required to be > 0 and <= 1'
            go to 999
         endif
         dthet(iprm)=dthetIN(iprm)

         iwrap(iprm)=iwrapIN(iprm)
         thetu(iprm)=thetIN(iprm)
      enddo

      do igrp=1,ngrps
         proddthet=1
         do iprmgrp=1,nprmgrp(igrp)
            iprm=iiprm(iprmgrp,igrp)
            proddthet=proddthet*min(2*dthet(iprm),1d0)
         enddo
      enddo
      nsumlike=0
      sumlike=0

      write(*,*)' Parameterization:'
      write(1,*)' Parameterization:'
      if(icycle.eq.3) then
         write(*,*)'   iprm    igrp iprmgrp   iwrap'
     1    ,'  pfocus   thet0    dthet'
         write(1,*)'   iprm    igrp iprmgrp   iwrap'
     1    ,'  pfocus   thet0    dthet'
      else
         write(*,*)'   iprm    igrp iprmgrp   iwrap'
     1    ,'  thet0    dthet'
         write(1,*)'   iprm    igrp iprmgrp   iwrap'
     1    ,'  thet0    dthet'
      endif

      do igrp=1,ngrps

         do iprmgrp=1,nprmgrp(igrp)
            iprm=iiprm(iprmgrp,igrp)
            if(icycle.eq.3) then
               write(*,'(4i8,2f8.5,1p2e9.2)')
     1          iprm,igrp,iprmgrp,iwrap(iprm)
     1          ,pfocusn(igrp),thetu(iprm),dthet(iprm)
               write(1,'(4i8,2f8.5,1p2e9.2)')
     1          iprm,igrp,iprmgrp,iwrap(iprm)
     1          ,pfocusn(igrp),thetu(iprm),dthet(iprm)
            else
               write(*,'(4i8,2f8.5,1p2e9.2)')
     1          iprm,igrp,iprmgrp,iwrap(iprm)
     1          ,thetu(iprm),dthet(iprm)
               write(1,'(4i8,2f8.5,1p2e9.2)')
     1          iprm,igrp,iprmgrp,iwrap(iprm)
     1          ,thetu(iprm),dthet(iprm)
            endif
         enddo
      enddo
 
      write(*,*)
      write(1,*)

      ifocus=1
      nvgood=1
      ihypv=1
 
      do iv1=1,nvgood
         bad(iv1)=.false.
         iv1good(iv1)=iv1

         do igrp=1,ngrps
            do iprmgrp=1,nprmgrp(igrp)
               iprm=iiprm(iprmgrp,igrp)
               thetv(iprm,iv1)=thetu(iprm)
c              write(*,*)'GMCMC: iv1,iprm,thetv=',iv1,iprm,thetv(iprm,iv1)
            enddo
         enddo
      enddo
  
      slow=.true.
      write(*,*)'GMCMC: before first call to CALCVX'
      write(1,*)'GMCMC: before first call to CALCVX'
      call calcvx

      iv=1
      write(*,*)'GMCMC: before first call to MOVE'
      write(1,*)'GMCMC: before first call to MOVE'
      call move !  move moves v variables for iv to u variables
c      write(*,*)'after move'

!  restore slow
      slow=slow0

      write(*,*)'GMCMC: first value, Chisq/nData=',c2f*esumu
      write(1,*)'GMCMC: first value, Chisq/nData=',c2f*esumu
      if(esumu.ne.esumu) then
         write(*,*)'GMCMC: Error, Chisq/nData bad, esumu=',esumu
         write(1,*)'GMCMC: Error, Chisq/nData bad, esumu=',esumu
         stop
      endif

      write(*,*)'GMCMC: initial parameter check complete!'
      write(1,*)'GMCMC: initial parameter check complete!'
c     write(*,*)'slow=',slow0

      call chktape(ireturn)
      if(ireturn.eq.1) return
      call runchain
      call finalze
      return

999   continue
      write(*,*)'GMCMC: ERROR ',trim(ermess1),' ',trim(ermess2)
      write(1,*)'GMCMC: ERROR ',trim(ermess1),' ',trim(ermess2)
      stop

      end
!-----------------------------------------------------------------------
      subroutine chktape(ireturn)
! checks for the existence of a tape file
!  if one exists and the data or requested number of chain
!  iterations haven't changed, just read from tape file
!  without recalculating chain
      use GMCMC_MOD, only : burnfrac,dataid,idorun,iseed,xrand,dranr250
     1 ,info,infoid,irun,irundbug,ndebug,nrun,ntape,tapefile
      implicit none

      integer ichk,ireturn,if1,ib1,ibuf1
     1 ,nrunfac,ntapetape,nthintape,nprmtape
      character*80 dataidtape
      character*120 cmd
      character*50 mess
      logical chkfile
      integer*8 nruncalc,nruntape
      common/r250com/if1,ib1,ibuf1(250)

c     write(*,*)
c    1 'CHKTAPE: ntape,nrun,nrun0,nrpts,irpt,nthin,nburn='
c    1 ,ntape,nrun,nrun0,nrpts,irpt,nthin,nburn
      ireturn=0

! set up debug printout
      if(ndebug.gt.0) then
         irundbug=nrun/ndebug
      endif

c      write(*,*)'dataid=',dataid
      if(trim(dataid).ne.'') then
! get identification of data used
         cmd='chksum$ '//trim(dataid)// ' > dataid.txt'
c         write(*,*)'cmd = ',cmd
         call shellcmd(cmd)

         if(chkfile('dataid.txt',infoid)) then
            open(167,file='dataid.txt')
            read(167,'(a80)')dataid
            write(*,*)'CHKTAPE: dataid=',dataid
            write(1,*)'CHKTAPE: dataid=',dataid
            close(167)
         endif
      endif

      idorun=0

      write(*,*)'CHKTAPE: Check for existing ',trim(tapefile)
      write(1,*)'CHKTAPE: Check for existing ',trim(tapefile)

      if(chkfile(trim(tapefile),info)) then
         open(170,file=trim(tapefile),action='read')
         write(*,*)'CHKTAPE: reading file: ',trim(tapefile)
     1    ,' info=',info
         write(1,*)'CHKTAPE: reading file: ',trim(tapefile)
     1    ,' info=',info
         read(170,*)nthintape,ntapetape,nruntape,nprmtape
         read(170,'(a80)')dataidtape
         write(*,*)' dataidtape=',trim(dataidtape)
         write(1,*)' dataidtape=',trim(dataidtape)
         close(170)

         nrunfac=nrun/nruntape
         nruncalc=nrunfac
         nruncalc=nruncalc*nruntape
c        write(*,*)'nrun,nruntape,nrunfac=',nrun,nruntape,nrunfac
c        write(*,*)'nrun,nruncalc=',nrun,nruncalc

         if(dataid.eq.dataidtape.and.nrun.eq.nruntape) then
            ichk=1
            mess=' with same nrun, dataid'
         else
            ichk=0
            mess=', use it'
         endif

         if(ichk.eq.1) then
!  use existing tape, skip runchain
            write(*,*)'CHKTAPE: existing ',trim(tapefile),trim(mess)
            write(1,*)'CHKTAPE: existing ',trim(tapefile),trim(mess)
!  initialize random number generator r250 to correspond to end of chain
            open(169,file=tapefile)
            read(169,*) ! skip two lines
            read(169,*)
100         read(169,*,err=101)irun ! error will be caused by presence of text " rand250 initialization:"
            go to 100
101         continue
!  the is the rand250 info needed for continuation of sequence of random numbers
            read(169,*)if1,ib1,ibuf1
            close(169)
            call finalze
            ireturn=1
            return

         elseif(nrun.eq.nruncalc.and.
     1    dataid.eq.dataidtape) then
!  extend existing tape
            write(*,'('' CHKTAPE: existing '',a9,'', nrun ='',i2
     1       ,'' x nrun tape same dataid'')')trim(tapefile),nrunfac
            write(1,'('' CHKTAPE: existing '',a9,'', nrun ='',i2
     1       ,'' x nrun tape same dataid'')')trim(tapefile),nrunfac
            idorun=2
            return
         else
!  rerun, making new tape
            write(*,*)'CHKTAPE: ',trim(tapefile)
     1       ,' exists but nrun or dataid not right for continuation'
            write(*,*)' on tape nrun = ',nruntape
     1       ,' data id = ',dataidtape
            write(1,*)'CHKTAPE: ',trim(tapefile)
     1       ,' exists but nrun or dataid not right for continuation'
            write(1,*)' on tape nrun = ',nruntape
     1       ,' data id = ',dataidtape
            write(*,*)' deleting ',trim(tapefile),', new run'
            write(1,*)' deleting ',trim(tapefile),', new run'
            cmd='del '//trim(tapefile)
            call shellcmd(cmd)
            idorun=1
            return
         endif

      endif

      idorun=1

      return
      end
!-----------------------------------------------------------------------
      subroutine runchain
! start or continue mcmc run
!  this routine handles the "tape file"
!  deciding whether to start new run or continue old one
      use GMCMC_MOD, only : burnfrac,dataid,idorun
     1 ,ihypu,ihypv,irpt,irun,irundbug,iseed,iv,iv1good
     1 ,nburn,ndebug,nprm,nrpts,nrun,ntape,nthin
     1 ,nv,nvgood,slow,slow0,tapefile,xsec0
      implicit none

      integer time
      character*80 dataidtape
      real*8 xrand,dranr250
      character*120 cmd
      integer i,ibuf1,ib1,iend,if1,iwrt0
     1 ,nrunfac,nrun0,ncopy,nnew,nvsav,nwrt0
      integer*8 irunrpt

      common/r250com/if1,ib1,ibuf1(250)

      if(idorun.eq.1) then
         go to 1
      elseif(idorun.eq.2) then
         go to 2
      else
         write(*,*)'RUNCHAIN: should not come here, idorun=',idorun
         write(1,*)'RUNCHAIN: should not come here, idorun=',idorun
         stop
      endif

1     continue ! come here to start new run

      nthin=nrun/ntape
      nrun=nthin
      nrun=nrun*ntape
      nburn=nint(burnfrac*nrun)
c     write(*,*)'nthin,nrun,nburn=',nthin,nrun,nburn

      write(*,*)'RUNCHAIN: new MCMC run of ',nrun,' iterations'
      write(*,'(''  write on tape'',i5
     1 ,'' records total (one every ''
     1 ,i10,'' iterations)'')')ntape,nthin
      write(1,*)'RUNCHAIN: new MCMC run of ',nrun,' iterations'
      write(1,'(''  write on tape'',i5
     1 ,'' records total (one every ''
     1 ,i10,'' iterations)'')')ntape,nthin

      xsec0=time()
c     write(*,*)'nrpts=',nrpts
      irunrpt=nrun/nrpts
c     write(*,*)'irunrpt,nrpts=',irunrpt,nrpts
      irpt=0
      irun=0

!  open new tape file
      open(168,file='new'//trim(tapefile))
      write(168,*)nthin,ntape,nrun,nprm
      write(168,'(a80)')dataid

!  initialize random number generator r250
      xrand=dranr250(iseed)
      write(*,*)'RUNCHAIN: first random number=',xrand
     1 ,' iseed=',iseed
      write(1,*)'RUNCHAIN: first random number=',xrand
     1 ,' iseed=',iseed

      slow=.true.

!  this do loop is the actual mcmc run

      do i=1,nrun
         call run1
!  first run with slow = true
         if(i.eq.1) slow=slow0

         if(mod(irun,nthin).eq.0) then
!   write record on tape
            call wrt1u(168)
            call increm0
         endif

         if(mod(irun,irunrpt).eq.0) then
!  write one of nrpts reports
            call wrtrpt
         endif
 
      enddo

!   write rand250 initialization info to allow continuation
      write(168,*)'rand250 initialization:'
      write(168,*)if1,ib1,ibuf1
      close(168)

      write(*,*)'RUNCHAIN: copy preliminary new'//trim(tapefile)
     1 //' -> '//trim(tapefile)
      write(1,*)'RUNCHAIN: copy preliminary new'//trim(tapefile)
     1 //' -> '//trim(tapefile)
      cmd='copy new'//trim(tapefile)//' '//trim(tapefile)//' /y'
      call shellcmd(cmd)
      cmd='del new'//trim(tapefile)
      call shellcmd(cmd)

      return

2     continue ! come here if previous run with good tapefile
      open(169,file=tapefile)
      read(169,*)iwrt0,nwrt0,nrun0
      read(169,'(a80)')dataidtape

      write(*,*)'RUNCHAIN: reading file: ',trim(tapefile)
     1 ,' dataid=',trim(dataidtape)
      write(1,*)'RUNCHAIN: reading file: ',trim(tapefile)
     1 ,' dataid=',trim(dataidtape)

!  new run will be some integer nrunfac times as long
      nrunfac=nrun/nrun0
c     write(*,*)'nrunfac,nrunfac*nrun0=',nrunfac,nrunfac*nrun0

      if(nrunfac.gt.1.and.nrun.eq.nrunfac*nrun0.and.
     1 dataid.eq.dataidtape) then
         write(*,*)'RUNCHAIN: ',trim(tapefile),' OK for continuation'
         write(1,*)'RUNCHAIN: ',trim(tapefile),' OK for continuation'
      else
         write(*,*)'RUNCHAIN: ',trim(tapefile)
     1    ,': nrun or dataid not right to allow continuation,'
     1    ,' start new run'
         write(*,*)'  nrun needs to be ',nrunfac*nrun0
         write(*,*)'  data id from tape = ',dataidtape
         write(1,*)'RUNCHAIN: ',trim(tapefile)
     1    ,': nrun or dataid not right to allow continuation,'
     1    ,' start new run'
         write(1,*)'  nrun needs to be ',nrunfac*nrun0
         write(1,*)'  data id from tape = ',dataidtape
         close(169)
         write(*,*)'  deleting ',trim(tapefile)
         write(1,*)'  deleting ',trim(tapefile)
         cmd='del '//trim(tapefile)
         call shellcmd(cmd)
         go to 1
      endif

      nthin=nrun/ntape
      nrun=nthin
      nrun=nrun*ntape
      nburn=nint(burnfrac*nrun)
c     write(*,*)'nthin,nrun,nburn=',nthin,nrun,nburn

      open(168,file='new'//trim(tapefile))
      write(168,*)nthin,ntape,nrun,nprm
      write(168,'(a80)')dataidtape

      ncopy=ntape/nrunfac
      nnew=ntape-ncopy

      write(*,*)'RUNCHAIN: extending mcmc run to ',nrun,' iterations'
      write(*,'(''  write '',i5,'' records total (one every ''
     1 ,i10,'' iterations)'')')ntape,nthin
      write(*,'(''  copy   '',i4,'' records from old tape'')')ncopy
      write(*,'(''  write  '',i4,'' new records'')')nnew
      write(1,*)'RUNCHAIN: extending mcmc run to ',nrun,' iterations'
      write(1,'(''  write '',i5,'' records total (one every ''
     1 ,i10,'' iterations)'')')ntape,nthin
      write(1,'(''  copy   '',i4,'' records from old tape'')')ncopy
      write(1,'(''  write  '',i4,'' new records'')')nnew

      xsec0=time()

      irunrpt=nrun/nrpts
c     write(*,*)'nrun,nrpts,nrpts,irunrpt=',nrun,nrpts,nrpts,irunrpt
      irpt=0

      nvsav=nv
      nv=1
      iv=1
      nvgood=1
      iv1good(1)=1

!  read record from tape
100   read(169,*,err=101)irun ! error will be caused by presence of text " rand250 initialization:"
c     write(*,*)'reading record from old tape'
c     write(*,*)'irun=',irun
!  get full info from tape every nthin records
      if(mod(irun,nthin).eq.0) then
         backspace(169)
         if(irun.eq.nrun0) then
            call rd1tov(169,iend)
            slow=.true.
            call calcvx
            call move
         else
            call rd1tou(169,iend)
         endif
         call wrt1u(168)
         call increm0
      endif
      go to 100

101   continue
!  initialize rand250
!  the is the rand250 info needed for continuation of sequence
!   of random numbers
      read(169,*)if1,ib1,ibuf1
      close(169)

      if(irun.ne.nrun0) then
         write(*,*)'RUNCHAIN: irun.ne.nrun0, irun,nrun0=',irun,nrun0
         write(1,*)'RUNCHAIN: irun.ne.nrun0, irun,nrun0=',irun,nrun0
         stop
      endif

      nv=nvsav

!  continue mcmc run

      do i=nrun0+1,nrun  !  this loop continues old run
         call run1
!   first run slow, because v variables not set for all iv1
!   can't restore for igrp=ifocus
!   after first run restore slow
         if(irun.eq.nrun0+1) slow=slow0

         if(mod(irun,nthin).eq.0) then
!   write record on tape
            call wrt1u(168)
            call increm0
         endif

         if(mod(irun,irunrpt).eq.0) call wrtrpt
      enddo

!   write rand250 initialization info to allow run continuation
      write(168,*)'rand250 initialization:'
      write(168,*)if1,ib1,ibuf1
      close(168)
      close(169)

      write(*,*)'RUNCHAIN: copy preliminary new'//trim(tapefile)
     1 //' -> ',trim(tapefile)
      write(1,*)'RUNCHAIN: copy preliminary new'//trim(tapefile)
     1 //' -> ',trim(tapefile)
      cmd='copy new'//trim(tapefile)//' '//trim(tapefile)
      call shellcmd(cmd)
      cmd='del new'//trim(tapefile)
      call shellcmd(cmd)
      write(*,*)'file ',trim(tapefile),' written'
      write(1,*)'file ',trim(tapefile),' written'
      end
!-----------------------------------------------------------------------
      subroutine increm0
!  increment talleys for wrtrpt
      use GMCMC_MOD, only : accept,c2,c2f,esumu,ihypu,incr
     1 ,incracpt,incrhyp,irun,nburn,nhyp 
      implicit none

      integer ih

c     write(*,*)'INCREM: irun,nburn=',irun,nburn

      incr(1)=incr(1)+1
      if(accept) incracpt(1)=incracpt(1)+1

      ih=ihypu+nhyp-1
c      write(*,*)'INCREM0: ih=',ih
      incrhyp(ih,1)=incrhyp(ih,1)+1
      c2(ih,1)=c2(ih,1)+c2f*esumu

      if(irun.le.nburn) return

      incr(2)=incr(2)+1

      if(accept) incracpt(2)=incracpt(2)+1

      incrhyp(ih,2)=incrhyp(ih,2)+1
      c2(ih,2)=c2(ih,2)+c2f*esumu

c      write(*,*)'INCREM: incr,ihypu,c2/incrhyp=',incr,ihypu,c2/incrhyp
      end
!----------------------------------------------------------------------- 
      subroutine zero0
!  zero built-in tallies
      use GMCMC_MOD, only : c2,incr,incracpt,incrhyp
     1 ,irun,irunacpt,nhyp,irpt
      implicit none

      irun=0
      irunacpt=0
      incr=0
      incracpt=0
      c2=0
      incrhyp=0
      end
!----------------------------------------------------------------------- 
      subroutine wrtrpt
      use GMCMC_MOD, only : c2,incr,incracpt,incrhyp
     1 ,irpt,irun,nhyp,nburn,nrpts,nrun,nv,xsec,xsec0
      implicit none
! write report. This is done nrpts times throughout the entire run
      integer time
      integer ih
      integer,save ::incracpt0,incrhyp0(2),incr10
      real*8 ,save :: c2sum0(2),ttaken
      real*4 tend

      if(irpt.eq.0) then
         write(*,*)'WRTRPT: Number of chain iterations = ',nrun
         write(*,*)' Number of candidates per chain iteration = ',nv
         write(*,*)' Number reports = ',nrpts
         write(1,*)'WRTRPT: Number of chain iterations = ',nrun
         write(1,*)' Number of candidates per chain iteration = ',nv
         write(1,*)' Number reports = ',nrpts
         if(nhyp.eq.2) then
            write(*,'(''      rpt#       irun   Chisq/ND(hyp1,0)  ''
     1       ,''time(min)    #tape''
     1       ,''    acptd NormalHyp    RefHyp'')')
            write(1,'(''      rpt#       irun   Chisq/ND(hyp1,0)  ''
     1       ,''time(min)    #tape''
     1       ,''    acptd NormalHyp    RefHyp'')')
         else
      write(*,'(''      rpt#       irun  Chisq/ND time(min)    #tape''
     1 ,''    acptd'')')
      write(1,'(''      rpt#       irun  Chisq/ND time(min)    #tape''
     1 ,''    acptd'')')
         endif
         tend=0
         incr10=0       
         incrhyp0=0
         incracpt0=0
         c2sum0=0
      endif

      xsec=time() ! system call giving clock time in seconds
      irpt=irpt+1
      tend=(nrun-irun)*(xsec-xsec0)/(60*irun)

c      write(*,*)'irpt,xsec,incr(1),incr10=',irpt,xsec,incr(1),incr10

      if(nhyp.eq.2) then
         write(*,'(i10,i11,1p2g10.3,0pf10.2,2i9,2i10)')
     1    irpt,irun
     1    ,((c2(ih,1)-c2sum0(ih))/max((incrhyp(ih,1)-incrhyp0(ih)),1)
     1    ,ih=nhyp,1,-1)
     1    ,tend
     1    ,incr(1)-incr10,incracpt(1)-incracpt0
     1    ,(incrhyp(ih,1)-incrhyp0(ih),ih=nhyp,1,-1)
         write(1,'(i10,i11,1p2g10.3,0pf10.2,2i9,2i10)')
     1    irpt,irun
     1    ,((c2(ih,1)-c2sum0(ih))/max((incrhyp(ih,1)-incrhyp0(ih)),1)
     1    ,ih=nhyp,1,-1)
     1    ,tend
     1    ,incr(1)-incr10,incracpt(1)-incracpt0
     1    ,(incrhyp(ih,1)-incrhyp0(ih),ih=nhyp,1,-1)
      else
         write(*,'(i10,i11,1p1g10.3,0pf10.2,2i9,2i10)')
     1    irpt,irun
     1    ,((c2(ih,1)-c2sum0(ih))/max((incrhyp(ih,1)-incrhyp0(ih)),1)
     1    ,ih=nhyp,1,-1)
     1    ,tend
     1    ,incr(1)-incr10,incracpt(1)-incracpt0
         write(1,'(i10,i11,1p1g10.3,0pf10.2,2i9,2i10)')
     1    irpt,irun
     1    ,((c2(ih,1)-c2sum0(ih))/max((incrhyp(ih,1)-incrhyp0(ih)),1)
     1    ,ih=nhyp,1,-1)
     1    ,tend
     1    ,incr(1)-incr10,incracpt(1)-incracpt0
      endif

      incracpt0=incracpt(1)

      incrhyp0=incrhyp(:,1)
      c2sum0=c2(:,1)

      incr10=incr(1)

      if(irun.ne.nrun) return

      write(*,*)' Totals (excluding burnfrac):'
      write(1,*)' Totals (excluding burnfrac):'

      if(nhyp.eq.2) then
         write(*,'(i10,i11,1p2g10.3,0pf10.2,2i9,2i10)')
     1    irpt,irun-nburn
     1    ,(c2(ih,2)/max(incrhyp(ih,2),1),ih=nhyp,1,-1)
     1    ,tend
     1    ,incr(2),incracpt(2)
     1    ,(incrhyp(ih,2),ih=nhyp,1,-1)
         write(1,'(i10,i11,1p2g10.3,0pf10.2,2i9,2i10)')
     1    irpt,irun-nburn
     1    ,(c2(ih,2)/max(incrhyp(ih,2),1),ih=nhyp,1,-1)
     1    ,tend
     1    ,incr(2),incracpt(2)
     1    ,(incrhyp(ih,2),ih=nhyp,1,-1)
       else
         write(*,'(i10,i11,1pg10.3,0pf10.2,2i9,2i10)')
     1    irpt,irun-nburn
     1    ,(c2(ih,2)/max(incrhyp(ih,2),1),ih=nhyp,1,-1)
     1    ,tend
     1    ,incr(2),incracpt(2)
         write(1,'(i10,i11,1pg10.3,0pf10.2,2i9,2i10)')
     1    irpt,irun-nburn
     1    ,(c2(ih,2)/max(incrhyp(ih,2),1),ih=nhyp,1,-1)
     1    ,tend
     1    ,incr(2),incracpt(2)
      endif
      xsec=time()
      ttaken=(xsec-xsec0)/60
      write(*,*)' Elapsed time to run chain(min) = ',ttaken
      write(1,*)' Elapsed time to run chain(min) = ',ttaken
      end
!-----------------------------------------------------------------------
      subroutine finalze
!  at end of run the "tape file" is read and tallies incremented
      use GMCMC_MOD, only : burnfrac,dataid,idorun,iv,iv1good
     1 ,irecmax,irun,nburn,irunacpt,ngrps,nhyp,nrun,nsumlike,ntape
     1 ,nthin,nv,nvgood,proddthet,prodsum,slow,sumlike
     1 ,tapefile,thetu,thetv,xsec0,One,esumref
     1 ,incrhyp,esumrefp,esumu,ihypu,ihypv,c2f
      implicit none
      integer time,iend
      real*8 xlike

      if(idorun.ge.1)then
         write(*,*)'FINALZE: Number of moves=',irunacpt
     1    ,' Average acceptance fraction (excluding burnfrac)='
     1    ,One*irunacpt/(nrun-nburn)
         write(1,*)'FINALZE: Number of moves=',irunacpt
     1    ,' Average acceptance fraction (excluding burnfrac)='
     1    ,One*irunacpt/(nrun-nburn)
         if(irunacpt.lt.ntape) then
            write(*,*)'WARNING: NMOVES < NTAPE, increase NRUN'
            write(1,*)'WARNING: NMOVES < NTAPE, increase NRUN'
         endif
         if(ngrps.eq.1) then
            prodsum=sumlike*proddthet/nsumlike
            write(*,*)'FINALZE: Direct estimate of integral of'
     1       ,' likelihood over prior=',prodsum
     1       ,' for',nsumlike,' trials'
            write(*,*)' Equivalent Chisq/nData=',-c2f*log(prodsum)
            write(1,*)'FINALZE: Direct estimate of integral of'
     1       ,' likelihood over prior=',prodsum
     1       ,' for',nsumlike,' trials'
            write(1,*)' Equivalent Chisq/nData=',-c2f*log(prodsum)
         endif
      endif

      nv=1
      iv=1
      nvgood=1
      iv1good(1)=1
      nburn=nint(burnfrac*nrun)
c     write(*,*)'FINALZE: nrun,ntape,burnfrac=',nrun,ntape,burnfrac
      xsec0=time()

!  read and process events from tape
      write(*,*)'FINALZE: reading from ',trim(tapefile)
      write(1,*)'FINALZE: reading from ',trim(tapefile)
      open(168,file=trim(tapefile),action='read')

      call zero0
      call zero
      slow=.true.
!  first record contains nthin,ntape,nrun
      read(168,*)nthin,ntape,nrun
      if(ntape.eq.0) then
         write(*,*)'no records on tape'
         write(1,*)'no records on tape'
         stop
      endif
!  second record contains dataid
      read(168,*)

!  read thru to find irecmax
      irecmax(1)=0
      irecmax(2)=0
100   read(168,*,end=101,err=101)irun
c     write(*,*)'irun,nburn=',irun,nburn
      irecmax(1)=irecmax(1)+1
      if(irun.gt.nburn) irecmax(2)=irecmax(2)+1
      go to 100
101   continue

c     write(*,*)'irecmax=',irecmax
      rewind(168)
      read(168,*)
      read(168,*)

1     call rd1tov(168,iend)
c     write(*,*)'irun=',irun
      if(iend.eq.1) go to 2

c     write(*,*)'before calcvx'
      if(nhyp.eq.2) then
         if(ihypv.eq.0) then
            esumu=esumrefp
            ihypu=ihypv
         else
            call calcvx
            call move ! move information from v to u variables
         endif
      else
         call calcvx
         call move 
      endif

      call increm0

!  increm operates on u variables
      if(irun.gt.nburn) call increm

      go to 1
2     continue

      if(nhyp.eq.2) then
         xlike=(exp(-esumref)*incrhyp(2,2))/incrhyp(1,2)
         write(*,*)'FINALZE: Estimate of integral of likelihood'
     1    ,' over prior from exp(-esumref)*nNormalHyp/nRefHyp'
         write(*,*)' =',xlike,' Equivalent Chisq/nData=',-c2f*log(xlike)
         write(1,*)'FINALZE: Estimate of integral of likelihood'
     1    ,' over prior from exp(-esumref)*nNormalHyp/nRefHyp'
         write(1,*)' =',xlike,' Equivalent Chisq/nData=',-c2f*log(xlike)
      endif
 
      write(*,*)'--------------------end GMCMC----------------------'
      write(1,*)'--------------------end GMCMC----------------------'

      close(168)
      close(1)
      end
!-----------------------------------------------------------------------
      subroutine run1
! generate candidate points in parameter space
!  test for acceptance, if accepted then move to new point
      use GMCMC_MOD, only : accept,icalcvx,irun,irunacpt,irundbug
     1 ,ndebug,nburn
      implicit none

      irun=irun+1

      if(ndebug.gt.0.and.mod(irun,irundbug).eq.0)
     1 write(*,*)'RUN1:irun=',irun

!  generate candidates
      call gencand

c     write(*,*)'RUN1: before calcacpt'
      call calcacpt

      if(accept) then
         if(ndebug.gt.0.and.mod(irun,irundbug).eq.0)
     1     write(*,*)'RUN1: moving'

         call move0
         if(icalcvx.eq.1) then
            call move1
            call movex
         endif
         if(irun.gt.nburn) irunacpt=irunacpt+1
c        write(*,*)'RUN1:accept=.true.,nburn=',nburn
      endif

      return
      end
!-----------------------------------------------------------------------
      subroutine move0
! move basic variable from v to u
      use GMCMC_MOD, only : esumu,esumv,ihypu,ihypv,iv
      implicit none

      ihypu=ihypv
      esumu=esumv(iv)
c      write(*,*)'MOVE0: iv,ihypu,esumu=',iv,ihypu,esumu

      return
      end
!-----------------------------------------------------------------------
      subroutine move1
! move theta variables from v to u
      use GMCMC_MOD, only : ifocus,iiprm,iv,nprm,nprmgrp
     1 ,slow,thetu,thetv
      implicit none
      integer igrp,iprm,iprmgrp

      if(slow) then
         thetu=thetv(:,iv)
      else

         igrp=ifocus

         do iprmgrp=1,nprmgrp(igrp)
            iprm=iiprm(iprmgrp,igrp)
c           write(*,*)'igrp,iprmgrp,iprm=',igrp,iprmgrp,iprm
            thetu(iprm)=thetv(iprm,iv)
c           write(*,*)'MOVE: iprm,thetu=',iprm,thetu(iprm)
         enddo

      endif
      end
!-----------------------------------------------------------------------
      subroutine move
      call move0
      call move1
      call movex
      return
      end
!-----------------------------------------------------------------------
      subroutine gencand
! generate candidates for next point in parameter space denoted by v
! current point is denoted by u
! theta variables are generated
! If there are groups of parameters, candidates generated
! only for group ifocus
      use GMCMC_MOD, only : bad,dthet,esumrefp,esumu,esumv
     1 ,fracrw,ialg,icalcvx,icycle,ifocus
     1 ,ihypu,ihypv,iiprm,irun,irundbug,iv1good,iwrap
     1 ,ndebug,ngrps,nhyp,nprm,nprmgrp,nv,nvgood
     1 ,pfocusc,proddthet,pswitch,ranwalk,slow,thetu,thetv

      implicit none
      logical nogrps
      data nogrps /.true./
      integer icalcd,igrp,igrp1,igrp2,ih,ihyp,iprm,iprmgrp,ivgood
     1 ,iv1
      real*8 xrand,dranr250

      if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) then
         if(nhyp.eq.2) then
            write(*,*)'GENCAND: nv=',nv,' ihypu=',ihypu
         else
            write(*,*)'GENCAND: nv=',nv
         endif
      endif

!  restore starting v values to u values for igrp = ifocus
!   which was moved last iteration
!   if no v values changed last iteration this is not neccessary
!   and icalcvx=0

c      write(*,*)'GENCAND: slow,nogrps,nv,nprm='
c     1 ,slow,nogrps,nv,nprm

      if(slow.or.nogrps) then
!  slow but sure, restore everthing
!  nogrps=.true. everything changed last time

         do ivgood=1,nv
            thetv(:,ivgood)=thetu  !  set to u value
         enddo

      elseif(icalcvx.eq.1) then
         do ivgood=1,nv
!  only parms for ifocus changed
            do iprmgrp=1,nprmgrp(ifocus)
               iprm=iiprm(iprmgrp,ifocus)
!  set to u value
               thetv(iprm,ivgood)=thetu(iprm)
c     write(*,*)'GENCAND: ivgood,iprm,thetv=',ivgood,iprm,thetu(iprm)
            enddo

!  restore dependent vx variables (extra v variables)
            call mvbackvx(ivgood)
         enddo
      endif

c      write(*,*)'thetv=',thetv(:,1)
c      write(*,*)'thetu=',thetu

      if(irun.lt.1.or.ndebug.eq.0) go to 100

!  check restore
c     write(*,*)'check restore, irun=',irun
      do 90 ivgood=1,nv
         do iprm=1,nprm
!  compare to u value
            if(thetv(iprm,ivgood).ne.thetu(iprm)) then
               write(*,*)'GENCAND: restore problem irun=',irun
               write(*,*)'thetv=',thetv(:,ivgood)
               write(*,*)'thetu=',thetu
               go to 90
            endif
         enddo
90    continue

100   continue

!  end check

!  Is random walk to be used?
!   use random walk for fraction fracrw of time
!   otherwise move variables over entire range
!   ranwalk also used in routine CALCPRBD

      if(fracrw.eq.0.) then
         ranwalk=.false.
      elseif(fracrw.eq.1.) then
         ranwalk=.true.
      else
         xrand=dranr250(0)
c        write(*,*)'xrand=',xrand
         if(xrand.lt.fracrw) then
            ranwalk=.true.
         else
            ranwalk=.false.
         endif
      endif

!  change ifocus if grouping used

      if(ngrps.gt.1) then
!  parameter grouping used, calculate ifocus
         if(icycle.eq.1) then
! cycle ifocus: 1, 2, 3 etc
            ifocus=mod(ifocus,ngrps)+1
         elseif(icycle.eq.2) then
! choose ifocus randomly with uniform probability
            xrand=dranr250(0)
            ifocus=ngrps*xrand+1
         else
! choose ifocus randomly with some specified probability distribution
! requires pfocusc, the cumulative probability of igrp
            xrand=dranr250(0)
            ifocus=icalcd(ngrps,pfocusc(1),xrand)
         endif
         igrp1=ifocus
         igrp2=ifocus
      else
         igrp1=1
         igrp2=ngrps
      endif

      nogrps=.false.
!  nogrps used in restore operation at beginning of this subroutine
!   normally=.false., except for nhyp=2 when can switch
!   from grouping with some number ngrps to no grouping

      icalcvx=1
!  icalcvx = 1 means that new energy calculation will be done
!   icalcvx sometimes set to 0  
!  icalcvx=1 indicates that thetv's will be moved and need to be restored
!   at beginning of next iteration

      ih=ihypu+nhyp-1  ! ih = 1, Ref hyp, ihypu=0, ih = 2, normal hyp

!  nhyp = 1, no alternate hypothesis
!  nhyp = 2 alt hyp with no parameters

      if(nhyp.eq.1) goto 1

!  come here for nhyp=2
!  Ref hypothesis, ihypu = 0, has no parameters, esum=esumref
!  Normal hypothesis, ihypu = 1
      if(ihypu.eq.1) then ! normal hypothesis
         if(dranr250(0).lt.pswitch) then
! probability of switching from normal hyp to ref hyp is pswitch
            ihypv=0  !  reference hypothesis
            nvgood=1           
            iv1good(1)=1
            esumv(1)=esumrefp ! corresponds to -log(P0/pswitch)
c            write(*,*)'GENCAND: ihypv=0,esumv,ihypu,esumu='
c     1       ,esumv(1),ihypu,esumu
            icalcvx=0  ! no thetv's moved
            go to 2
         else ! ihyp1 -> ihyp1 normal chain movement
            ihypv=1  !  normal hypothesis
            go to 1
         endif 
      else ! ihypu=0 ref hyp immediately switches ihyp0 -> ihyp1
         ihypv=1  !  normal hypothesis
         ranwalk=.false. ! transitions to entire space 
! ranwalk also used in routine CALCPRBD
         nogrps=.true. ! no grouping used, important for next restore
         igrp1=1
         igrp2=ngrps
         go to 1
      endif

1     continue !  calculate theta's for candidates

      nvgood=0

      do 1000 iv1=1,nv
         bad(iv1)=.false. !  start with clean badness slate

         do igrp=igrp1,igrp2
            do iprmgrp=1,nprmgrp(igrp)
               iprm=iiprm(iprmgrp,igrp)

!  move theta possibly with random walk
!   or possibly chosing from entire space, depending on RANWALK

               if(ranwalk) then
                  if(dthet(iprm).ge.1.) then
                     thetv(iprm,iv1)=dranr250(0)
                  else
                     thetv(iprm,iv1)=thetu(iprm)
     1                +(2*dranr250(0)-1)*dthet(iprm)
                     if(iwrap(iprm).eq.2) then
c                       write(*,*)'before wrap=',thetv(iprm,iv1)
                        thetv(iprm,iv1)=thetv(iprm,iv1)+10
     1                   -int(thetv(iprm,iv1)+10)
c                       write(*,*)'after wrap=',thetv(iprm,iv1)
                     elseif(thetv(iprm,iv1).lt.0.or.
     1                thetv(iprm,iv1).gt.1) then
                        bad(iv1)=.true.
c                       write(*,*)'iprm,thetv=',iprm,thetv(iprm,iv1)
                        if(nv.eq.1) go to 999
                        go to 1000
                     endif
                  endif
               else
                  thetv(iprm,iv1)=dranr250(0)
               endif
            enddo
         enddo

c     write(*,*)'iv1,iprm,thetv,thetu,dthet,deviation=',iv1,iprm
c    1 ,thetv(iprm,iv1),thetu(iprm),dthet(iprm)
c    1 ,(thetv(iprm,iv1)-thetu(iprm))/dthet(iprm)
 
         nvgood=nvgood+1
         iv1good(nvgood)=iv1

1000  continue

2     continue !  skip normal candidate calculation above

999   continue
      if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) then
       write(*,*)'iv1good     iv1    ihyp    igrp    iprm       thetv'
     1    ,'       dthet  deviation    ranwalk       bad'
         ihyp=ihypv
         do ivgood=1,nvgood
            iv1=iv1good(ivgood)
c           write(*,*)'iv1,ngrps,igrp,nprmgrp='
c    1       ,iv1,ngrps,ifocus,nprmgrp(igrp)
            igrp=ifocus

            do igrp=igrp1,igrp2
               do iprmgrp=1,nprmgrp(igrp)
                  iprm=iiprm(iprmgrp,igrp)
                  write(*,'(5i8,2f12.8,g12.4,2l10)')
     1             iv1,iv1,ihyp,igrp,iprm
     1             ,thetv(iprm,iv1),dthet(iprm)
     1             ,(thetv(iprm,iv1)-thetu(iprm))/dthet(iprm)
     1             ,ranwalk,bad(iv1)
               enddo
            enddo
         enddo
c        write(*,*)'GENCAND: before call calcvx'
      endif

      if(icalcvx.eq.1) call calcvx
      end
!-----------------------------------------------------------------------
      subroutine calcacpt
!  call calcprbd to calculate the unnormalized probability prbnd of
!   iv1good=1,nvgood new configurations
!   including the current configuration as iv1good=0
!   These are stored in prbcd(i), i=1,nvgood+1
!  Generate ivgoodnew from prbcd, which determines whether
!   to stay (ivgoodnew=0)
!   or move to a new configuration
      use GMCMC_MOD, only : accept,ihypu,ihypv,irun,irundbug
     1 ,iv,iv1good,ndebug,nprm,nvgood,prbcd,thetv
      implicit none
      integer i,igoodnew,ihyp,ip,iprm,nprb
      real*8 xrand,dranr250

!  ivgood=1,nvgood are new configurations with bad.ne..true.

      if(nvgood.eq.0) then
c        write(*,*)'CALCACPT: nvgood=',nvgood
         igoodnew=0
      else
         call calcprbd
         nprb=nvgood+1
!        igoodnew=igenpd(nprb,prbcd)-1
         xrand=dranr250(0)
c        write(*,*)'IGENPD:xrand=',xrand
         do ip=1,nprb
            if(xrand.le.prbcd(ip)) go to 10
         enddo
         ip=nprb
         write(*,*)'CALCACPT ERROR: should not come here'
     1    ,', irun=',irun,', xrand,prbcd=',xrand,(prbcd(i),i=1,nprb)
10       continue
         igoodnew=ip-1
      endif

      if(igoodnew.gt.0) then
         accept=.true.
         iv=iv1good(igoodnew)
      else
         accept=.false.
         iv=0
      endif

      if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) then
         if(iv.eq.0) then
            ihyp=ihypu
         else
            ihyp=ihypv
         endif

         write(*,*)'CALCACPT: accept=',accept,' iv=',iv
     1    ,' ihyp=',ihyp
         if(accept) then
            write(*,*)'   iprm       thetv'
            do 12 iprm=1,nprm
            write(*,'(i8,f12.8)')iprm,thetv(iprm,iv)
12          continue
         endif
      endif
      end
!-----------------------------------------------------------------------
      subroutine calcprbd
! calculates relative probability of new points in parameter space
      use GMCMC_MOD, only : bad,eminv,esumu,esumv
     1 ,ialg,ifocus,ihypu,ihypv,iprbmax,irun,irundbug,iv1good,nburn
     1 ,ndebug,ngrps,nhyp,nsumlike ,nvgood,prbcd,prbnd,pswitch
     1 ,sumlike,ranwalk
      implicit none
      integer igood,iiv,iprb,ivgood,ivgoodmin,ivmax,ivmin,iv1
      real*8 eminv0,esum,sum,sumc
      logical dosum

      ivmin=iv1good(1)
      ivgoodmin=1

!   find minimum energy of new candidates
      do ivgood=1,nvgood
         iv1=iv1good(ivgood)

         if(esumv(iv1).lt.esumv(ivmin)) then
c           write(*,*)'iv1,esumv(iv1),esumv(ivmin)='
c    1       ,iv1,esumv(iv1),esumv(ivmin)
            ivgoodmin=ivgood
            ivmin=iv1
         endif
      enddo

      eminv=esumv(ivmin)

      if(esumu.lt.eminv) then
         iprbmax=1
         eminv0=esumu
      else
         iprbmax=ivgoodmin+1
         eminv0=eminv
      endif

      sum=0

      do iprb=2,nvgood+1
         ivgood=iprb-1
         iv1=iv1good(ivgood)

         if(ialg.eq.1) then
            if(esumv(iv1).lt.esumu) then
               prbnd(iprb)=1.d0/nvgood
            else
               prbnd(iprb)=dexp(esumu-esumv(iv1))/nvgood
            endif
         elseif(ialg.eq.2) then
            prbnd(iprb)=1/((1+dexp((esumv(iv1)-esumu)))*nvgood)
         else
            prbnd(iprb)=dexp(eminv0-esumv(iv1))
         endif
         sum=sum+prbnd(iprb)
      enddo

      if(ialg.eq.3) then
!  multiple candidate algorithm
         prbnd(1)=dexp(eminv0-esumu)
         sum=sum+prbnd(1)
         do iprb=1,nvgood+1
            prbnd(iprb)=prbnd(iprb)/sum
         enddo
      else
!  calculate probability for iprb=1, chain stays put
         prbnd(1)=1-sum
      endif

      sumc=0
      dosum=ranwalk.and.ngrps.eq.1.and.ihypv.eq.1
     1    .and.irun.gt.nburn

      do iprb=1,nvgood+1
         sumc=sumc+prbnd(iprb)
!  calculate the cummulative probability over the candidates
         prbcd(iprb)=sumc

!  This code is used to calculate sumlike, which can be used
!   to directly estimate the normalization integral
!   only works for ranwalk=.true. and ngrps=1
         if(dosum.and.iprb.gt.1) then

            ivgood=iprb-1
            iv1=iv1good(ivgood)
            esum=esumv(iv1)

            sumlike=sumlike+dexp(-esum)
            nsumlike=nsumlike+1
         endif
      enddo

      if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) then
         if(iprbmax.eq.1) then
            ivmax=0
         else
            ivmax=iv1good(iprbmax-1)
         endif
         write(*,*)'CALCPRBD: esumu,eminv,ivMaxProb='
     1    ,esumu,eminv,ivmax
         write(*,*)'     iv   igood        esum         prb'

         do iprb=1,nvgood+1
            if(iprb.eq.1) then
               esum=esumu
               iiv=0
               igood=0
            else
               igood=iprb-1
               iiv=iv1good(igood)
               esum=esumv(iiv)
            endif
            write(*,'(2i8,g12.3,f12.8)')iiv,igood,esum,prbnd(iprb)
         enddo

      endif
      end
!-----------------------------------------------------------------------
      function igenpd(n,pc)
! generate a random integer from 1 to n
!  with probability given by the cumulative probability pc
      implicit none
      integer i,igenpd,n
      real*8 xrand,dranr250
      real*8 pc(n)

      xrand=dranr250(0)
c     write(*,*)'IGENPD:xrand=',xrand
      do igenpd=1,n
         if(xrand.le.pc(igenpd)) return
      enddo

      igenpd=n
      write(*,*)'IGENPD ERROR: should not come here'
     1 ,', xrand,pc=',xrand,(pc(i),i=1,n)
      end
!-----------------------------------------------------------------------
      function icalcd(n,pc,theta)
! given cumulative probability distribution pc(i) over i=1,n
!  calculate the integer icalcd that has cumulative probability theta
      implicit none
      integer icalcd,n
      real*8 theta,pc(n)

      do icalcd=1,n
         if(theta.le.pc(icalcd)) return
      enddo

      icalcd=n
      write(*,'('' ICALCD ERROR: should not come here''
     1 ,'' theta,pc(n)='',2z17.16)') theta,pc(n)
      write(1,'('' ICALCD ERROR: should not come here''
     1 ,'' theta,pc(n)='',2z17.16)') theta,pc(n)
      end
!-----------------------------------------------------------------------
      subroutine wrt1u(iunit)
!  write one u record on "tape"
      use GMCMC_MOD, only:cesum,cthet,esumu,ifocus
     1 ,ihypu,irun,nprm,thetu
      implicit none
      integer iunit

      write(cesum,'(z16)') esumu
      write(cthet,'(z16)') thetu

      write(iunit,*)irun
     1 ,ihypu,cesum,ifocus,cthet

      end
!-----------------------------------------------------------------------
      subroutine rd1tou(iunit,iend)
!  read one record from "tape" to u arrays
      use GMCMC_MOD, only:cesum,cthet,esumu,ifocus,ihypu,irun
     1 ,ngrps,nprm,thetu
      implicit none
      integer iend,iunit

      iend=0

      read(iunit,*,end=99)irun
     1 ,ihypu,cesum,ifocus,cthet
      read(cesum,'(z16)') esumu
      read(cthet,'(z16)') thetu

      if(ngrps.eq.1) ifocus=1

      return

99    continue
      iend=1
      end
!-----------------------------------------------------------------------
      subroutine rd1tov(iunit,iend)
!  read one record from "tape" to v arrays
      use GMCMC_MOD, only:cesum,cthet,esumv,ifocus,ihypv,irun
     1 ,iv,ngrps,nprm,thetv
      implicit none
      integer iend,iunit

      iend=0
      read(iunit,*,end=99,err=99)irun
     1 ,ihypv,cesum,ifocus,cthet
      read(cesum,'(z16)') esumv(iv)
      read(cthet,'(z16)') thetv(:,iv)

      if(ngrps.eq.1) ifocus=1

c     write(*,*)'RD1TOV: iv,ihypv,thetv(1,iv)='
c    1 ,iv,ihypv,thetv(1,iv)

      return

99    continue
      iend=1
      end
!-----------------------------------------------------------------------
      logical function chkfile(filename,info)
! check if a file with path filename of length >0 exists.
      implicit none
      integer i
      character*(*) filename,info
      character*120 cmd

!  use dir command to get file info, put in file FILEINFO.TXT
      cmd='dir '//trim(filename)//' > fileinfo.txt'
c      write(*,*)'cmd=',cmd

      call shellcmd(cmd)

!  if file exists, file info is in 6th record
c      write(*,*)'before open'
      open(221,file='fileinfo.txt')
c      write(*,*)'before do i=1,6'
      do i=1,6
         read(221,'(a39)',end=99,err=99) info
c         write(*,*)'i,info=',i,info
      enddo
      close(221)

      chkfile=.true.
!  check for 0 file size
      if(info(37:38).eq.' 0') chkfile=.false.
c      write(*,*)'chkfile,info=',chkfile,info
      return

99    continue
      chkfile=.false.
c      write(*,*)'ISTHERE: chkfile,info=',chkfile,info
      close(221)
      end
!-----------------------------------------------------------------------
      subroutine shellcmd(cmd)
      implicit none
      character*(*) cmd
      integer system,iexit

      cmd='cmd /c '//trim(cmd)
      iexit=system(trim(cmd))
c     write(*,*)'iexit=',iexit
      if(iexit.ne.0.and.iexit.ne.1) then
         write(*,*)'SHELLCMD: abnormal ending, iexit='
     1   ,iexit,', cmd=',cmd
         write(1,*)'SHELLCMD: abnormal ending, iexit='
     1   ,iexit,', cmd=',cmd
      endif
      return
      end
!------------------------------Begin r250-------------------------------
      real*8 function dranr250(ir)
      if(ir.gt.0) call r250init(ir)
      dranr250=irand250()/2147483648.0 ! largest 31 bit integer + 1 = 2**31
      return
      end
!-----------------------------------------------------------------------
      subroutine rtest
!     Test the random number generators

      Print *,'The EXPECTED values are in parenthesis'

!     First test the LCM generator
      ix = lcmrand( 1 )
      Do 10 I = 1, 9999
         ix = lcmrand( -1 )
   10 EndDo
      Print *,'Final LCM value: (1043618065)', ix

!     Now initialize and test R250
      Call R250Init( 1 )
      Do 20 I = 1, 10000
        ix = irand250()
   20 EndDo
      Print *,'Final R250 value: (267771767)', ix
 
      END
!-----------------------------------------------------------------------
! R250.F77     The R250 Pseudo-random number generator

! algorithm from:
! Kirkpatrick, S., and E. Stoll, 1981; A Very Fast Shift-Register
! Sequence Random Number Generator, Journal of Computational Physics,
! V. 40. p. 517

! see also:
! Maier, W.L., 1991; A Fast Pseudo Random Number Generator,
!                    Dr. Dobb's Journal, May, pp. 152 - 157

! Uses the Linear Congruential Method,
! the "minimal standard generator"
! Park & Miller, 1988, Comm of the ACM, 31(10), pp. 1192-1201
! for initialization

! For a review of BOTH of these generators, see:
! Carter, E.F, 1994; Generation and Application of Random Numbers,
! Forth Dimensions, Vol. XVI, Numbers 1,2 May/June, July/August

! $Author:   skip  $
! $Workfile:   r250.f  $
! $Revision:   1.1  $
! $Date:   07 Nov 1996 01:23:06  $
! ===================================================================
      Function lcmrand(ix)
!     The minimal standard PRNG for 31 bit unsigned integers
!     designed with automati! overflow protection
!     uses ix as the seed value if it is greater than zero
!     otherwise it is ignored
      Integer*4 ix
      Integer*4 a, b, m, q, r
      Integer*4 hi, lo, test
      Integer*4 x
      SAVE x
      Parameter (a = 16807, b = 0, m = 2147483647)
      Parameter (q = 127773, r = 2836)

      If ( ix .gt. 0 ) x = ix
 
      hi = x / q
      lo = mod( x, q )
      test = a * lo - r * hi
      if ( test .gt. 0 ) then
          x = test
      else
          x = test + m
      endif
 
      lcmrand = x
      return
      End
 
! ===================================================================
C
!  R250, call R250Init with the desired initial seed BEFORE
!  the first invocation of IRAND250()
C
! ===================================================================
 
      Subroutine R250Init(iseed)
      Integer k, mask, msb
      Integer indexf, indexb, buffer(250)
      Common/R250COM/indexf,indexb,buffer
      Integer ms_bit, all_bits, half_range, step
      Parameter ( ms_bit = Z'40000000')
      Parameter ( half_range = Z'20000000' )
      Parameter ( all_bits = Z'7FFFFFFF' )
      Parameter ( step = 7 )

      indexf = 1
      indexb = 104
      k = iseed
      Do 10 i = 1, 250
	  buffer(i) = lcmrand( k )
	  k = -1
  10  EndDo
      Do 20 i = 1, 250
	 if ( lcmrand( -1 ) .gt. half_range ) then
	     buffer(i) = ior( buffer(i), ms_bit )
	 endif
 20   EndDo
 
      msb = ms_bit
      mask = all_bits
 
      Do 30 i = 0,30
	  k = step * i + 4
	  buffer(k) = iand( buffer(k), mask )
	  buffer(k) = ior( buffer(k), msb )
	  msb = msb / 2
          mask = mask / 2
  30  EndDo
 
      Return
      END
!----------------------------------------------------------------------- 
      Function irand250()
!     R250 PRNG, run after R250_Init
      Integer newrand
      Integer indexf, indexb, buffer(250)
      Common/R250COM/indexf, indexb,buffer
 
      newrand = ieor( buffer(indexf), buffer(indexb) )
      buffer(indexf) = newrand
 
      indexf = indexf + 1
      if ( indexf .gt. 250 ) indexf = 1
 
      indexb = indexb + 1
      if ( indexb .gt. 250 ) indexb = 1
	
 
      irand250= newrand
      return
      End
!----------------------------end r250-----------------------------------
!###############################END GMCMC###############################
!-----------------------------------------------------------------------
!########################test program follows###########################
!-----------------------------------------------------------------------
      module testMCMC_mod
      integer, save :: ndim,nsum
      real*8, save :: x0,sig0
      real*8, save, dimension(:), allocatable :: xu
      real*8, save, dimension(:,:), allocatable :: xv,xsum
      end module testMCMC_mod
!-----------------------------------------------------------------------
      program testMCMC
! test of MCMC basic programs
! Markov Chain Monte Carlo, writing tape file of events
! In the core MCMC prog, all parameters are expressed
!   in terms of theta variables, which run from 0 to 1
! All priors on theta variables are independent and uniform
!   The actual variables ("x" variables) are expressed as
!   functions of the theta variables, where the function is such that
!   actual variable has the desired prior probability distribution when
!   theta variable is unifomly distributed from 0 to 1
! set up call to gmcmc
      use GMCMC_MOD, only: xsec,xsec0
      use testMCMC_mod
      implicit none
      integer ialg,icycle,igrp,iseed,j,iprmgrp,iprm,iwrap
     1 ,mdata,nburn,ndebug,ngrps,nprm
     1 ,nprmsgrp,nrpts,ntape,nv,time
      real*8 dthetrw,fracrw,thet0,xliketrue,avg

      integer, dimension(:), allocatable :: iwrap1,nprmgrp
      integer, dimension(:,:), allocatable :: iiprm
      real*8, dimension(:), allocatable:: pfocus,thet
      real*8, dimension(:,:), allocatable:: dthet
      real*8 burnfrac
      character*70 datafile
      data datafile/' '/
      character*10 tapefile
      logical slow,exist
      data slow/.false./
      integer*8 nrun

      namelist/l/iseed,thet0,nrun,ntape,nrpts,ndebug,burnfrac
     1 ,fracrw,dthetrw,iwrap,icycle,ngrps,nv,ialg,ndim,x0,sig0

      tapefile='tape.txt'  !  name of tape file 
      ntape=100 ! number of iterations written on tape file 
      ialg=1  ! MCMC algorithm: 1 for MRT, 2 for B, 3 for MC
      nv=1  ! number of multiple candidates
      nrun=1000000 ! number of chain iterations in MCMC run
      iwrap=1  !  wrap variable if random walk goes outside of (0,1) range 
      icycle=1  !  cycling of parameter groups that are varied: 1 for cyclic, 2 random,
!                   3 specified probability supplied in pfocus
      ngrps=6  ! number of parameter groups
      thet0=0  ! starting theta value (for all dimensions)
      burnfrac=.05d0  !  initial discarded fraction of chain iterations
      dthetrw=.001d0 ! random walk Delta theta
      fracrw=0.5d0  ! fraction of time random walk is used rather than generating candidates over entire space
      nrpts=10  !  number of reports during course of run
      iseed=1  ! random number generator seed
      ndebug=0  ! >0 for debug printout

!  Specification of desired steady state distribution:
      ndim=6  ! dimensionality
      x0=0.5  !  location of mode(s)
      sig0=1.d-6 ! width of mode(s)

      inquire(file='gmcmc.in',exist=exist)
      if(.not.exist)goto 998
      open(10,file='gmcmc.in')
      read(10,l,err=999)
      close(10)
998   continue
      write(*,l)

      nprm=ndim
      nprmsgrp=ndim

      allocate (iiprm(nprmsgrp,ngrps),nprmgrp(ngrps),iwrap1(nprm))
      allocate (pfocus(nprm),thet(nprm),dthet(nprm,2))
      allocate (xv(ndim,nv),xu(ndim),xsum(3,ndim))

      if(ngrps.eq.1) then
         nprmgrp(1)=ndim
         igrp=1
         pfocus(igrp)=1
         do iprmgrp=1,nprmgrp(igrp)
            iiprm(iprmgrp,igrp)=iprmgrp
         enddo
      elseif(ngrps.eq.ndim) then
         do igrp=1,ngrps
            pfocus(igrp)=1
            nprmgrp(igrp)=1
            iiprm(1,igrp)=igrp
         enddo
      else
         write(*,*)'bad ngrps=',ngrps
         stop
      endif

      nburn=nint(burnfrac*nrun)

      do igrp=1,ngrps
         do iprmgrp=1,nprmgrp(igrp)
            iprm=iiprm(iprmgrp,igrp)
            iwrap1(iprm)=iwrap
            dthet(iprm,1)=dthetrw
            dthet(iprm,2)=dthet(iprm,1)
            thet(iprm)=thet0
         enddo
      enddo

      xliketrue=1
      do j=1,ndim
         xliketrue=xliketrue*sqrt(2*3.14159256)*sig0
      enddo
      write(*,*)'integral of likelihood function over prior='
     1 ,xliketrue
 
      mdata=ndim

      call GMCMC(nrun,ntape,nrpts,ndebug
     1 ,mdata,iseed,icycle,nv,ialg
     1 ,slow,burnfrac,fracrw
     1 ,nprm,ngrps,nprmsgrp,nprmgrp,iiprm
     1 ,pfocus,iwrap1,thet,dthet
     1 ,datafile)

      write(*,*)'After GMCMC: nsum=',nsum
      write(*,*)'j xavg xsd'
      do j=1,ndim
         avg=xsum(1,j)/nsum
         write(*,*)j,avg,sqrt(xsum(2,j)/nsum-avg**2)
      enddo

      write(*,*)

      xsec=time()
      write(*,*)'Elapsed time for post processing (sec) ='
     c ,(xsec-xsec0)
      stop

999   write(*,*)'ERROR involving namelist read of GMCMC.IN'
      go to 998
      end
!-----------------------------------------------------------------------
      subroutine zero
      use testMCMC_mod, only:xsum,nsum
      implicit none
!  zero tallies
      xsum=0
      nsum=0
      end
!----------------------------------------------------------------------- 
      subroutine increm
!  increment talleys
      use testMCMC_mod, only:xsum,xu,ndim,nsum
      implicit none
      integer j
      do j=1,ndim
         xsum(1,j)=xsum(1,j)+xu(j)
         xsum(2,j)=xsum(2,j)+xu(j)**2
      enddo
      nsum=nsum+1
      end
!-----------------------------------------------------------------------
      subroutine calcvx
! calculate dependent v variables from theta variables
      use GMCMC_MOD, only:esumv,iv1good,nvgood,thetv
      use testMCMC_mod, only:ndim,xv,x0,sig0
      implicit none
      integer iv1,ivgood,iprm,j
      real*8 chi

!  here is the calculation of the likelihood function
!   normally done in two possible ways, slow but sure and as fast as possible
!    but hopefully also sure!
!  This simulates an energy calculation in terms of the actual parameters
!    rather than the theta parameters

      do ivgood=1,nvgood
         iv1=iv1good(ivgood)
         esumv(iv1)=0
         do j=1,ndim
            iprm=j
            xv(j,iv1)=thetv(iprm,iv1) ! normally there is some general functional relationship between the problem
!  variables xv and the MCMC theta variables (theta from 0 to 1 with uniform prior probability)
            chi=(xv(j,iv1)-x0)/sig0 ! normally some general functional relationship between problem variables and chi 
            esumv(iv1)=esumv(iv1)+chi**2/2 
         enddo
      enddo
      end
!-----------------------------------------------------------------------
      subroutine mvbackvx(iv1)
!  restore extra dependent v variables
      use testMCMC_mod, only:ndim,xu,xv
      implicit none
      integer iv1,j

      do j=1,ndim
         xv(j,iv1)=xu(j)
      enddo
      end
!-----------------------------------------------------------------------
      subroutine movex
! If a candidate Xv is accepted, chain is moved and Xv->Xv
!  This routine moves v variables to u variables
!  if ngrps > 0 for the group ifocus only, the minimum required
!    after gencand, calcacpt, because other variables have not changed
      use GMCMC_MOD, only: iv
      use testMCMC_mod, only:ndim,xu,xv
      implicit none
      integer j

!  move other dependent variables

      do j=1,ndim
         xu(j)=xv(j,iv)
      enddo
      end
!-----------------------------------------------------------------------
