      program su3datcollect ! Berg Jan 11 2009. Gather (g) SU3 data
C  rom sublattices into two data file for nltat=1 and 2 lattices.
      include '../../Libs/Fortran/implicit.08'
      include 'mpif.h'
      character cd*1,cf*1,cdmpi*1,clx*3,clt*3,cmy*4,clat*1
      include '../../Libs/Fortran/constants.08'
      include 'latmpi.par'
      include 'mc.par'
      include '../../Libs/MPISU3/common_cbsu3.f'
      dimension tsag(nmeas),tswsg(nmeas),tswtg(nmeas),tsprg(nmeas),
     &  tspig(nmeas),tsk1rg(nmeas),tsk1ig(nmeas),aacpt(msmpi),nla(nd),
     &  nla_mpi(ndmpi)
      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 "su3datcollect: error 1-processor rusn only!"
      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)') nl1
      write(clt,'(I3.3)') nl4
      write(iuo,'(/,"su3datcollect: msmpi =",I9,/)') msmpi
C
C Read headers:
C
      do ilat=1,nlat ! nlat=2 for double-layered torus. Otherwise 1.
        write(iuo,'(2X,/,"Read and sum data for ilat =",I6,":",/)') ilat
        if(nlat.eq.1) then
          open(iud1,file=cfln//cd//"d"//cf//"f"//cdmpi//"d"//clx//"x"//
     &                   clt//"t.D",form='unformatted',status='unknown')
        else
          write(clat,'(I1.1)') ilat
          open(iud1,file=cfln//cd//"d"//cf//"f"//cdmpi//"d"//clx//"x"//
     &         clt//"t"//clat//".D",form='unformatted',status='unknown')
        endif
        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
          write(iuo,'(2X,"data filename: ",2A,".D")') cfln,cfile
          iud=iud2+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
          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
          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
        write(iud1) 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
C     
C Read and sum data:
C
        write(iuo,'(/,2x,"irpt,tsag(1),twswg,tswtg,abs(tspg(1)),")')
        write(iuo,'(2x,"tsprg(1),tspig(1),tsk1rg(1),tsk1ig(1):",/)')
        do irpt=1,nrpt
          call razero(tsag,ione,nmeas)
          call razero(tswsg,ione,nmeas)
          call razero(tswtg,ione,nmeas)
          call razero(tsprg,ione,nmeas)
          call razero(tspig,ione,nmeas)
          call razero(tsk1rg,ione,nmeas)
          call razero(tsk1ig,ione,nmeas)
          acptg=zero
          aming=zero
          amaxg=zero
          do ismpi=1,msmpi ! Read data:
            iud=iud2+ismpi-1+(ilat-1)*msmpi
            read(iud) tsa,tsws,tswt,tspr,tspi,tsk1r,tsk1i,
     &                acpt,a_min,a_max,i_in
            if(irpt.eq.1) then
              aacpt(ismpi)=acpt
              if(ismpi.eq.msmpi) 
     &          write(iuo,'("  acpt rates:",4G12.3)') aacpt
            endif
            if(i_in.ne.irpt) then
              write(iuo,'(" su2datcollect: i_in.ne.irpt.")')
              ltest=.true.
            endif
      if(ltest) then ! -----------------------------------------
        write(iuo,'(" su2datcollect: ltest STOP.")')
        call mpi_finalize(ierr)
      stop
      endif ! --------------------------------------------------
            do imeas=1,nmeas
              tsag(imeas)=tsag(imeas)+tsa(imeas)/msmpi
              tswsg(imeas)=tswsg(imeas)+tsws(imeas)/msmpi
              tswtg(imeas)=tswtg(imeas)+tswt(imeas)/msmpi
              tsprg(imeas)=tsprg(imeas)+tspr(imeas)/msmpi
              tspig(imeas)=tspig(imeas)+tspi(imeas)/msmpi
              tsk1rg(imeas)=tsk1rg(imeas)+tsk1r(imeas)/msmpi
              tsk1ig(imeas)=tsk1ig(imeas)+tsk1i(imeas)/msmpi
              acptg=acptg+acpt/msmpi
              aming=aming+a_min/msmpi
              amaxg=amaxg+a_max/msmpi
            enddo
c         write(iuo,'(I6,4G15.6)')iud,tsa(1),tswt(1),tspr(1),tsk1r(1)
          end do ! Write summed data:
          atspg=sqrt(tsprg(1)**2+tspig(1)**2) 
          write(iuo,'(I6,4G15.6)') irpt,tsag(1),tswsg(1),tswtg(1),atspg
          write(iuo,'(6X,4G15.6)') tsprg(1),tspig(1),tsk1rg(1),tsk1ig(1)
          write(iud1)tsag,tswsg,tswtg,tsprg,tspig,tsk1rg,tsk1ig,
     &               acptg,aming,amaxg,irpt
        end do
        close(iud1)
      enddo
C
      call mpi_finalize(ierr)
      stop "su3datcollect: All done."
      end

      include '../../Libs/Fortran/razero.f'
