      subroutine rhoiter(loop,nr,rhoin,rhoout,phi,f0,f1,
     x a,b,rho0,tmp,ifsav,iscr,miter,den)
c Broyden's iteration for density saved on unit ifsav;
c uses scratch file iscr
      implicit real*8 (a-h,o-z)
      real*8 rhoin(nr),rhoout(nr),f0(nr),f1(nr),a(nr),
     x b(nr),rho0(nr),tmp(nr)
      write(6,*) 'entering convrg with miter=',miter
      write(6,*) 'entering convrg with loop',loop
      if (loop.le.0) then
      do 20 ig=1,nr
      rhoin(ig)=rhoin(ig)+phi*(rhoout(ig)-rhoin(ig))
   20 continue
      return
      endif
      if (loop.gt.0) then
      rewind ifsav
      read(ifsav) (rhoin(ig),ig=1,nr)
      read(ifsav) (rhoout(ig),ig=1,nr)
      do 30 ig=1,nr
      f0(ig)=rhoout(ig)-rhoin(ig)
      rho0(ig)=rhoin(ig)
   30 continue
      rewind iscr
      do 100 il=1,loop
      read(ifsav) (rhoin(ig),ig=1,nr)
      read(ifsav) (rhoout(ig),ig=1,nr)
      do 40 ig=1,nr
      f1(ig)=rhoout(ig)-rhoin(ig)
      tmp(ig)=(f1(ig)-f0(ig))
      b(ig)=tmp(ig)
      a(ig)=(rhoin(ig)-rho0(ig))+phi*tmp(ig)
   40 continue
      if (il.gt.1) then
      im=il-1
      rewind iscr
      do 60 ij=1,im
      read(iscr) (tmp(i),i=1,nr)
      xxx=dot_product(tmp,b)/dot_product(tmp,tmp)
      if (il-ij.lt.miter) then
      do 53 ig=1,nr
      b(ig)=b(ig)-tmp(ig)*xxx
   53 continue
      endif
      read(iscr) (tmp(i),i=1,nr)
      do 50 ig=1,nr
      a(ig)=a(ig)-tmp(ig)*xxx
   50 continue
   60 continue
      endif
      den=dot_product(b,b)
      write(6,*) 'den=',il,den
      do 70 ig=1,nr
      f0(ig)=f1(ig)
      rho0(ig)=rhoin(ig)
   70 continue
      write(iscr) (b(i),i=1,nr)
      write(iscr) (a(i),i=1,nr)
  100 continue
      do 110 ig=1,nr   
      rhoin(ig)=rhoin(ig)+phi*f1(ig)
  110 continue
      rewind iscr
      do 120 il=1,loop
      read(iscr) (b(i),i=1,nr)
      read(iscr) (a(i),i=1,nr)
      cc=dot_product(f1,b)/dot_product(b,b)
      do 115 ig=1,nr
      rhoin(ig)=rhoin(ig)-cc*a(ig)
  115 continue
  120 continue
      endif
      return
      end
