      subroutine ! bernd berg oct 17 2001.
     & rmag_steb0(nrpt,iud,n0,ms,nqm1,mlink,ns,hm,hmm,hme,qm,ha)
c extends read_steb0 to include magnetization.
c reads nrpt times the the magnetization arrays hm() from unit iud and 
c calculates the histogram mean values hmm() and their error bars hme().
c further, each of the nq magnetizations is calculated nrpt times in qm().
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      dimension hm(n0:ms,n0:nqm1),hmm(n0:ms,n0:nqm1),hme(n0:ms,n0:nqm1)
      dimension qm(nrpt,n0:nqm1),ha(n0:mlink)
c
      do is=n0,ns
      do iq=n0,nqm1
        hmm(is,iq)=zero
      end do
      end do
      do irpt=1,nrpt
        read(iud) ha,acpt,irpt_in,nstate0,hm
        do iq=0,nqm1 ! potts_actl works also for the magnetization!
          qm(irpt,iq)=potts_actl(ns,hm(n0,iq))
        end do
        do iq=n0,nqm1
        do is=n0,ns
          hmm(is,iq)=hmm(is,iq)+hm(is,iq)
        end do
        end do
      end do
      do is=n0,ns
      do iq=n0,nqm1
        hmm(is,iq)=hmm(is,iq)/nrpt
      end do
      end do
c
      do irpt=1,nrpt
        backspace iud
      end do
c
      cut=half/nrpt
      do iq=n0,nqm1
      do is=n0,ns
        hme(is,iq)=zero
      end do
      end do
      do irpt=1,nrpt
        read(iud) ha,irpt_in,nstate0,hm
        do iq=n0,nqm1
        do is=n0,ns
        if(hmm(is,iq).gt.cut) then ! prevents underflows.
          hme(is,iq)=hme(is,iq)+(hm(is,iq)-hmm(is,iq))**2
        end if
        end do
        end do
      end do
      do iq=n0,nqm1
      do is=n0,ns
        hme(is,iq)=sqrt(hme(is,iq)/(nrpt*(nrpt-1)))
      end do
      end do
c
      return
      end
