      program scalarhb ! Bernd Berg Mar 28 2009.
C Gaussian heatbath for free scalar field. 
C Action and mean field measurements.
      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/Scalar_par/common_scalar.f'
      include 'lat.dat'
      ltest=.true.
      ltest=.false.
c     if(ltest) stop "scalarhb: ltest."
      write(iuo,'(/,"scalarhb test program:")')
      call scalar_init ! Initialization.
c
      write(iuo,'(" Equilibration nreq,nequi =",2I9)') nreq,nequi
      do ireq=1,nreq
        do iequi=1,nequi   ! Sweeps for reaching equilibrium.
          call freehb_ssweep ! Free scalar field.
c         write(iuo,'(I9,1G15.6)') iequi,phisum/ms
        end do ! Check action and field:
        call scalar_actdif(sact,phisum,phi,xmfactor,xla,ipf,ms,nd,da,dp)
      end do
      write(iuo,'(" After equilib phim =",G15.6)') phisum/ms
      write(iuo,'(" sactnn/ms,sact2/ms =",2G15.6,/)') sactnn/ms,sact2/ms
c
c Writing header information into file time series (action)
      open(iud1,file=cfln//cfile//".D",form='unformatted',
     &                                   status='unknown')
      write(iud1) xm2,xla,nd,nla,mlink,nequi,nrpt,nmeas,nsw
      close(iud1)
      write(iuo,'(" irpt",6X,"action",9X,"field")')
      do irpt=1,nrpt ! Repetitions.
        call write_progress(iud1,"irpt,nrpt,sact:",irpt,nrpt,sact)
        do imeas=1,nmeas   ! Measurements Loop.
          do isw=1,nsw
            call freehb_ssweep ! Free scalar field.
          end do
          tsann(imeas)=sactnn/ms
          tsa2(imeas)=sact2/ms
          tsphi(imeas)=phisum/ms
          call scalar_cor(imeas)  ! Measure correlations.
        end do
        tsa0=half+sact/ms
        call scalar_actdif(sact,phisum,phi,xmfactor,xla,ipf,ms,nd,da,dp)
        write(iuo,'(I4,4X,2G15.6,F11.3)') irpt,tsa0,phisum/ms
        open(iud1,file=cfln//cfile//".D",form='unformatted',
     &                           status='old',access='append')
        write(iud1) tsann,tsa2,tsphi,one,phi0sum,cor0,cor1r,cor1i,irpt
        close(iud1)
      end do

      stop
      end

c General layout includes
      include '../../Libs/Fortran/ipointer.f'
      include '../../Libs/Fortran/isfun.f'
      include '../../Libs/Fortran/ixcor.f'
      include '../../Libs/Fortran/lat_init.f'
      include '../../Libs/Fortran/nsfun.f'
      include '../../Libs/Fortran/write_progress.f'
      include '../../Libs/Fortran/razero.f'
      include '../../Libs/Fortran/ranmar.f'
      include '../../Libs/Fortran/rmafun.f'
      include '../../Libs/Fortran/rmagau.f'
      include '../../Libs/Fortran/rmaset.f'
      include '../../Libs/Fortran/rmasave.f'

c Scalar routines:
      include '../../Libs/Scalar/sactphi1.f'     ! Action at one site.
      include '../../Libs/Scalar/freehb_supdt.f' ! Free scalar update.
      include '../../Libs/Scalar/scalar_act.f' ! Action.
      include '../../Libs/Scalar/scalar_actdif.f' ! Action difference.
      include '../../Libs/Scalar/scalar_start.f'  ! Start configuration.
      include '../../Libs/Scalar_par/scalar_cor.f'  ! Correlations.
      include '../../Libs/Scalar_par/scalar_init.f' ! Initialization.
      include '../../Libs/Scalar_par/freehb_ssweep.f' ! Free scalar sweep.
