      program su3_tshb ! Alexe Bazavov / Bernd Berg Mar 14 2008.
c SU3 MCMC time series program based on Cabibbo-Marinari heatbath.
      include '../../Libs/Fortran/implicit.sta'
      include '../../Libs/Fortran/constants.par'
      include 'lat.par'
      include 'mc.par'
      include '../../Libs/Fortran_par/common_lat.f'
      include '../../Libs/SU3_par/common_su3.f'
      include 'lat.dat'
      ltest=.true.
      ltest=.false.
c
      write(iuo,'(/,"su3_tshb: beta =",F20.10,/)') beta
      call su3_init ! Initialization.
c
      do iequi=1,nequi ! Equilibration:
        call su3_mchb  ! CM heatbath.
      end do
c 
c Write header information into time series file:
      open(iud1,file=cfln//cfile//".D",form="unformatted",
     &                                 status="unknown")
      write(iud1) beta,nd,nla,mlink,nequi,nrpt,nmeas,nsw
      close(iud1)
      write(iuo,'(2X,"Equilibration done act/np =",G16.7)') act/np
c
      write(iuo,'(/,3X,"irpt",4x,"act/np",7x,"(act-actdirect)/np",
     &              4X,"acpt",/)')
      do irpt=1,nrpt ! nrpt repititions of nmeas measurement sweeps.
        iact=nint(act)
        call write_progress(iud1,"irpt,i(act),acpt:",irpt,iact,acpt)
        acpt=zero
        a_min=act
        a_max=act
        do imeas=1,nmeas
          do isw=1,nsw
            call su3_mchb ! CM heatbath.
          enddo
          tsa(imeas)=act/np      ! Action per plaquette.
          call su3_wloops(imeas) ! Record space and time plqs.
          call su3_ploop(imeas)  ! Polyakov loop measurements
        end do
c Output data:
        open(unit=iud1,file=cfln//cfile//".D",form="unformatted",
     &       status="old",access="append")
        write(iud1) tsa,tsws,tswt,tspr,tspi,acpt,a_min,a_max,irpt
        close(iud1)
        call su3_act(actdrct) ! Check action by direct calculation.
        actdif=abs(act-actdrct)/np
        write(iuo,'(I6,F16.11,G16.4,F13.3)') irpt,act/np,actdif,acpt
        act=actdrct ! Continue with directly measured action.
      end do

      stop "su3_tshb: All done."
      end


c lattice layout:
      include '../../Libs/Fortran/ixcor.f'
      include '../../Libs/Fortran/lat_init.f'
      include '../../Libs/Fortran/ipointer.f'
      include '../../Libs/Fortran/nsfun.f'
      include '../../Libs/Fortran/isfun.f'
c random number generator:
      include '../../Libs/Fortran/ranmar.f'
      include '../../Libs/Fortran/rmaset.f'
      include '../../Libs/Fortran/write_progress.f'
c needed for su3:
c     include "../../Libs/Fortran/dim_reduced.f"
      include "../../Libs/SU2/su2_a0.f"
      include "../../Libs/SU3/su3reunit.f"
c matrix operations
      include "../../Libs/SU3/su3add_m_m.f"
      include "../../Libs/SU3/su3copy_m_m.f"
      include "../../Libs/SU3/su3init0.f"
      include "../../Libs/SU3/su3init1.f"
      include "../../Libs/SU3/su3mult_m_m_m.f"
      include "../../Libs/SU3/su3mult_m_mh_m.f"
      include "../../Libs/SU3/su3mult_mh_mh_m.f"

c needed for su3, including parameters:
      include "../../Libs/SU3_par/su3_init.f"
      include "../../Libs/SU3_par/su3_act.f"
      include "../../Libs/SU3_par/su3_wloops.f"
      include "../../Libs/SU3_par/su3_ran.f"
      include "../../Libs/SU3_par/su3_order.f"
      include "../../Libs/SU3_par/su3_mchb.f"
      include "../../Libs/SU3_par/su3_getstaple.f"
      include "../../Libs/SU3_par/su3_update_su2.f"
c Polyakov loop measurements:
      include "../../Libs/SU3_par/su3_ploop.f"
