subroutine bremm(t,rho,mue,z,qbrem)

      implicit double precision (a-h,o-z)

      dimension b(3),c(4)
      double precision mue

      a=5.167d-1
      b(1)=7.577d0
      b(2)=-2.801d0
      b(3)=0.267d0
      c(1)=0.4775d0
      c(2)=1.707d0
      c(3)=-9.090d0
      c(4)=5.250d0
      x=dlog10(rho/mue)
      y=dlog10(t/1.d8)
      g=c(1)*atan(c(2)*(x-y)+c(3))+c(4)
      f=y*g+a*x-((b(1)+b(2)*x+(b(3)*(x**2)))**(0.5d0))
      eps=7.870d-3*(2.d0/mue)*z*(10.d0**(f))
      qbrem=eps*rho

      return
      end

************************************************************************