      program ana2sublatsu3 ! Berg Dec 10 2008.
      include '../../Libs/Fortran/implicit.sta'
      include 'mpif.h'
      character cd*1,cf*1,cdmpi*1,clx*3,clt*3,cmy*4,clat*1
C Analyze SU3 data from simulation on a double-layered torus.
C SU3   average action and integrated autocorrelation time.
      include '../../Libs/Fortran/constants.par'
      include 'latmpi.par'
      include 'mc.par'
      include '../../Libs/MPI_par/common_cblat.f'
      include '../../Libs/MPISU3/common_cbsu3.f'
      dimension act0(nrpt,msmpi,nlat),act0m(msmpi,nlat),
     &  act0me(msmpi,nlat),act(nrpt,nlat),actm(nlat),
     &  actme(nlat),acpt(nrpt)
      include 'latmpi.dat'
      ltest=.true.
      ltest=.false.
C
C Initializations:
C
      call mpi_init(ierr)
      call mpi_comm_size(mpi_comm_world,n1proc,ierr)
      if(n1proc.ne.1) then
        call mpi_finalize(ierr)
      stop "ana2sublatsu3: error 1-processor run required!"
      endif
      write(cd,'(I1.1)') nd ! Prepare filename to read the data:
      write(cf,'(I1.1)') mpifactor
      write(cdmpi,'(I1.1)') ndmpi
      write(clx,'(I3.3)') nla(1)
      write(clt,'(I3.3)') nla(nd)
      write(iuo,'(/,2X,"ana2sublatsu3 nlat,msmpi,lat2test =",
     &   2I9,L7,/)') nlat,msmpi,lat2test
C
C Read headers:
C
      do ilat=1,nlat
        call razero(act(1,ilat),ione,nrpt)
        do ismpi=1,msmpi
          my_id=ismpi-1+(ilat-1)*msmpi
          write(cmy,'(I4.4)') my_id
          write(cfile,'(A,"d",A,"f",A,"d",A,"x",A,"t",A)') 
     &          cd,cf,cdmpi,clx,clt,cmy
          iud=iud1+my_id
          open(iud,file=cfln//cfile//".D",
     &         form="unformatted",status="old")
          read(iud) b0_in,b1_in,nd_in,nlat_in,ns,nlink,nla,ndmpi_in,
     &      msmpi_in,nla_mpi,nreq_in,nequi_in,nrpt_in,nmeas_in,nsw_in
          if(my_id.eq.0) then
            write(iuo,'(2X,"beta0,beta1,nlink,nd:",2F15.10,I10,I3,/)')
     &                      beta0,beta1,nlink,nd_in
            write(iuo,'(2X,"nla:",20X,4I5)') nla
            write(iuo,'(2X,"nla_mpi:",16X,4I5)') nla_mpi
            write(iuo,'(2X,"msmpi,ndmpi:",I17,I5)') msmpi_in,ndmpi_in
            write(iuo,'(2X,"nreq,nequi:",I18,I10)') nreq_in,nequi_in
            write(iuo,'(2X,"nrpt,nmeas,nsw:",I14,2I10)')
     &                      nrpt_in,nmeas_in,nsw_in
          endif
          write(iuo,'(2X,"Data filename: ",2A,".D;  ilat,ismpi,my_id =",
     &          I2,2I6)') cfln,cfile,ilat,ismpi,my_id
          if(nd.ne.nd_in) stop "nd.ne.nd_in."
          if(nlat.ne.nlat_in) stop "nlat.ne.nlat_in."
          if(ndmpi.ne.ndmpi_in) stop "ndmpi.ne.ndmpi_in."
          if(nreq.ne.nreq_in) stop "nreq.ne.nreq_in."
          if(nequi.ne.nequi_in) stop "nequi.ne.nequi_in."
          if(nrpt.ne.nrpt_in) stop "nrpt.ne.nrpt_in."
          if(nmeas.ne.nmeas_in) stop "nmeas.ne.nmeas_in."
          if(nsw.ne.nsw_in) stop "nsw.ne.nsw_in."
        end do
      end do
      write(iuo,'(2X,"All data files open.")')
C     
C Read and analyse data:
C     
      do ilat=1,nlat
        write(iuo,'(/,"Analysis for ilat =",I6,":",/)') ilat
        idat=0 
        itime=0
        write(clat,'(I1.1)') ilat
        do ismpi=1,msmpi                           ! Read data:
          iud=iud1+ismpi-1+(ilat-1)*msmpi
          acpt(ismpi)=zero
          do irpt=1,nrpt
            read(iud) tsa,acpt(irpt) ! ,i_in
            act0(irpt,ismpi,ilat)=stmean(nmeas,tsa)
            act(irpt,ilat)=act(irpt,ilat)+act0(irpt,ismpi,ilat)/msmpi
          enddo
          call steb0(nrpt,act0(1,ismpi,ilat),act0m(ismpi,ilat),
     &               av,act0me(ismpi,ilat))
          write(iuo,'(" ismpi, actm =",I6,F12.8,2X,"+/-",F12.8)') 
     &                  ismpi,act0m(ismpi,ilat),act0me(ismpi,ilat)
          close(iud)
        enddo
        call steb0(nrpt,act(1,ilat),actm(ilat),av,actme(ilat))
          write(iuo,'(/,8X,"actm =",F12.8,2X,"+/-",F12.8)') 
     &                      actm(ilat),actme(ilat)
      enddo
      if(ltest) stop "ana2sublatsu3.f: ltest."      
C
      if(nlat.eq.2) then ! Combine both layers.
        do irpt=1,nrpt
          act(irpt,1)=half*(act(irpt,1)+act(irpt,nlat))
        enddo
        call steb0(nrpt,act(1,1),actm(1),actv,actme(1))
        write(iuo,'(/," Average from both layers:",/,/,8X,
     &    "actm =",F12.8,2X,"+/-",F12.8,/)') actm(1),actme(1)
      endif
C
      call mpi_finalize(ierr)
      stop
      end
      
      include '../../Libs/Fortran/razero.f'
      include '../../Libs/Fortran/steb0.f'
      include '../../Libs/Fortran/stmean.f'
