subroutine delsq(ap,rho,ibc)
c
      implicit none
      include 'itg.par'
c
      real ap(lz,mz,nz,nspecz),rho(lz,mz,nz,nspecz),
     * den,det,r32,r22
      integer ibc,mr,mr2,lp,lm,l,m,n,ldbb,i
      include 'itg.cmn'
c
c     evaluate	rho = grad**2 ap
c		rho = ((d/dx)(d/dx) - (m/y0)**2) ap
c
c Boundary conditions at r=0 and r=rmax=r(ld) are:
c
c ibc=-1	rho(0)=rho(rmax)=0.0
c ibc=1		rho(rmax)=0, rho even about r=0, (i.e., d/dr rho = 0)
c ibc=2		periodic B.C.'s, rho(0)=rho(rmax)
c ibc=-2	boundary values are extrapolated from interior points.
c

c do the calculation first on the internal points:
      do 10 i=1,nspecies
      do 10 n=1,nd
      do 11 m=1,md
      mr=mrr(m,n)
      mr2=mr**2
      do 1 l=2,ldb
      lp=l+1
      lm=l-1
      rho(l,m,n,i)=eb(l)*(ap(lp,m,n,i)-ap(l,m,n,i))
     1 -fb(l)*(ap(l,m,n,i)-ap(lm,m,n,i))
     1 -mr2*ap(l,m,n,i)/(y0**2)
1     continue
      rho(ld,m,n,i)=0.0
11    continue
10    continue

c now do the boundaries:
      do 30 i=1,nspecies
      do 30 n=1,nd
      do 30 m=1,md
      mr=mrr(m,n)
      mr2=mr*mr
c
c ibc=1		rho(rmax)=0, rho even about r=0, (i.e., d/dr rho = 0)
c
      if(ibc.eq.1) then
      den=1.+.5*(mr*r(2)/y0)**2
      rho(1,m,n,i)=(ap(2,m,n,i)-den*ap(1,m,n,i))*2./(r(2)**2)

c ibc=2		periodic B.C.'s, rho(0)=rho(rmax)

      else if(ibc.eq.2) then
      rho(1,m,n,i)=eb(1)*(ap(2,m,n,i)+ap(ldb,m,n,i)-2.0*ap(1,m,n,i))
     *          -mr2*ap(1,m,n,i)/(y0*y0)
      rho(ld,m,n,i)=rho(1,m,n,i)

c ibc=-2	boundary values are extrapolated from interior points.

      else if(ibc.eq.-2) then
      r32=r(3)*r(3)
      r22=r(2)*r(2)
      det=r32-r22
      rho(1,m,n,i)=(r32*ap(2,m,n,i)-r22*ap(3,m,n,i))/det
      ldbb=ld-2
      r32=(r(ldbb)-r(ld))*(r(ldbb)-r(ld))
      r22=(r(ldb)-r(ld))*(r(ldb)-r(ld))
      det=r32-r22
      rho(ld,m,n,i)=(r32*ap(ldb,m,n,i)-r22*ap(ldbb,m,n,i))/det

c ibc=-1	rho(0)=rho(rmax)=0.0

      else
      rho(1,m,n,i)=0.0
      endif

30    continue 
      return 
      end