text:
cc Copyright (C) 2009-2012: Leslie Greengard and Zydrunas Gimbutas
cc Contact: greengard@cims.nyu.edu
cc
cc This software is being released under a modified FreeBSD license
cc (see COPYING in home directory).
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c    $Date: 2011-06-21 12:54:41 -0400 (Tue, 21 Jun 2011) $
c    $Revision: 2105 $
c
c
c    Subroutine library for shifting via projection
c
C***********************************************************************
      subroutine h3dprojloc(nterms,ldl,nquad,xnodes,wts,phival,local,
     1           ynm,w,lw,ier)
C***********************************************************************
C     Purpose:
C
C           This subroutine projects a function given on the unit sphere
C           onto spherical harmonics.
C           It is intended for use with a library of special
C           functions in which Ynm does NOT include the
C           1/sqrt(4*pi) normalization. Thus, the Ynm are orthogonal
C           but not orthonormal -> the inner product value is 4*pi.
c
C           phival(theta,phi) tabulated at nquad*nquad points:
C           Legendre nodes x_j = cos(theta_j) in polar angle and
C           equispaced in the azimuthal direction.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl =    defines leading dimension of local
C           nquad  = number of quadrature nodes in each direction.
C           xnodes = Gauss-Legendre nodes x_j = cos theta_j
C           wts    = Gauss quadrature weights
C           phival = tabulated function
C                    phival(i,j) = phi(sin theta_j cos phi_i,
C                                      sin theta_j sin phi_i,
C                                      cos theta_j).
C           ynm     = workspace of dimension (0:nterms,0:nterms)
c                     used to get Ynm values.
C           w       = workspace of length lw
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local = coefficients of s.h. expansion
C           ier = error flag
C                 0 normal return
C                 8 insufficient memory in w.
C
C***********************************************************************
      implicit real *8 (a-h,o-z)
      integer nterms,nquad,ier,ldl
      integer l,m,jj,kk
      real *8 w(lw)
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk,phival(nquad,nquad)
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 ephi,imag,emul,sum,zmul
      data imag/(0.0d0,1.0d0)/
C
      ier = 0
      imarray = 1
      lmarray = 2*nquad*(2*nterms+1)+7
C
      lused = imarray+lmarray
C
      if (lused.gt.lw) then
         ier=6
	 return
      endif
c
      call h3dprojloc0(nterms,ldl,nquad,xnodes,wts,phival,local,
     1           w(imarray),ynm)
      return
      end
c
C***********************************************************************
      subroutine h3dprojloc0(nterms,ldl,nquad,xnodes,wts,phival,local,
     1           marray,ynm)
C***********************************************************************
C
C     Usage:
C
C           compute spherical harmonic expansion on unit sphere
C           of function tabulated at nquad*nquad grid points.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl    = dimension parameter for local expansion
C           nquad  = number of quadrature nodes in each direction.
C           xnodes = Gauss-Legendre nodes x_j = cos theta_j
C           wts    = Gauss quadrature weights
C           phival = tabulated function
C                    phival(i,j) = phi(sin theta_j cos phi_i,
C                                      sin theta_j sin phi_i,
C                                      cos theta_j).
C
C           marray  = workspace of dimension (nquad,-nterms:nterms)
C           w       = workspace of length lw
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local = coefficients of s.h. expansion
C
C    NOTE:
C
C    yrecursion.f produces Ynm with a nonstandard scaling:
C    (without the 1/sqrt(4*pi)). Thus the orthogonality relation
C    is
C             \int_S  Y_nm Y_n'm'*  dA = delta(n) delta(m) * 4*pi.
C
C   In the first loop below, you see
C
Cccc	    marray(jj,m) = sum*2*pi/nquad
C	    marray(jj,m) = sum/(2*nquad)
C
C   The latter has incorporated the 1/(4*pi) normalization factor
C   into the azimuthal quadrature weight (2*pi/nquad).
C
C***********************************************************************
      implicit real *8 (a-h,o-z)
      integer nterms,nquad,ier
      integer l,m,jj,kk
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk,phival(nquad,nquad)
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 marray(nquad,-nterms:nterms)
      complex *16 ephi,imag,emul,sum,zmul,emul1
      data imag/(0.0d0,1.0d0)/
C
      pi = 4.0d0*datan(1.0d0)
c
c     initialize local exp to zero
c
      do l = 0,ldl
         do m = -l,l
	    local(l,m) = 0.0d0
         enddo
      enddo
c
c     create marray (intermediate array)
c
      emul1 = cdexp(imag*2*pi/nquad)
      emul = cdexp(-imag*2*nterms*pi/nquad)
      do m=-nterms,nterms
cc	    emul = cdexp(imag*m*2*pi/nquad)
	 do jj=1,nquad
	    sum = 0
	    ephi = emul
	    do kk = 1,nquad
               sum = sum + phival(jj,kk)*dconjg(ephi)
	       ephi = ephi*emul
            enddo
ccc	    marray(jj,m) = sum*2*pi/nquad
	    marray(jj,m) = sum/(2*nquad)
         enddo
	 emul = emul*emul1
      enddo
c
c     get local exp
c
      do jj=1,nquad
	 cthetaj = xnodes(jj)
	 call ylgndr(nterms,cthetaj,ynm)
ccc	 call prinf(' nterms is *',nterms,1)
ccc	 call prinf(' ldl is *',ldl,1)
ccc	 call prinf(' l is *',l,1)
         do m=-nterms,nterms
	    zmul = marray(jj,m)*wts(jj)
            do l=abs(m),nterms
               local(l,m) = local(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
      return
      end
c
C***********************************************************************
      subroutine h3drescale(nterms,local,radius,zk0,scale,w,lw,ier)
C***********************************************************************
C
C     This subroutine rescales a spherical harmonic expansion
C     on the surface by j_n(zk0 * radius), converting
C     a surface function to the corresponding Bessel expansion
C     in the interior.
C     It is intended for use with a specific library of special
C     functions in which local expansions are also scaled by
C     i* zk0, so this factor is scaled out here for consistency.
c
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           local = coefficients of s.h. expansion
C           radius = sphere radius
C           zk0 = Helmholtz parameter
C           scale  = scale parameter for J expansion.
C           w       = workspace of length lw
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local = rescaled by 1/(j_n(zk0* radius)*i*zk0)
C           ier = error flag
C                 0 normal return
C                 4 insufficient memory in w.
C                 8 lwfjs insufficient in calling jfuns3d.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,ier,l,m,jj,kk,lwfjs
      real *8 w(lw)
      complex *16 local(0:nterms,-nterms:nterms)
      complex *16 ephi,imag,emul,sum,zmul
      complex *16 zk0,z
      data imag/(0.0d0,1.0d0)/
C
      ier = 0
      lwfjs = nterms+1000
c
      ifjs = 1
      lfjs = 2*(lwfjs+1)+7
C
      ifjder = ifjs+lfjs
      lfjder = 2*(nterms+1)+7
C
      iiscale = ifjder+lfjder
      liscale = (lwfjs+1) +7
C
      lused = iiscale+liscale
      if (lused.gt.lw) then
         ier=4
	 return
      endif
c
      z = zk0*radius
      ifder = 0
      call jfuns3d(ier1,nterms,z,scale,w(ifjs),
     1             ifder,w(ifjder),lwfjs,w(iiscale),ntop)
      if (ier1.eq.8) then
         ier = 8
	 return
      endif
ccc      call prin2(' fjs is *',w(ifjs),2*nterms)
      do l=0,nterms
         do m=-l,l
	    zmul = dcmplx(w(ifjs+2*l),w(ifjs+2*l+1))
	    local(l,m) = local(l,m)/zmul
         enddo
      enddo
      return
      end
c
C***********************************************************************
      subroutine h3dmpevalsphereslow(mpole,zk,scale,phival,zshift,
     1           radius,nterms,nquad,nquadm,xnodes)
C***********************************************************************
C
C     This subroutine evaluates a multipole expansion on an
C     nquad x nquadm spectral grid on a target sphere
C     along the z axis a distance zshift.
C
C     TESTING CODE  ---- NOT USED
C
C---------------------------------------------------------------------
C     INPUT:
C
C     mpole    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : mpole scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     nquad    : number of quadrature nodes in theta direction
C     nquadm   : number of quadrature nodes in phi direction
C
C                total number of nodes on target sphere is nquad*nquadm
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phival   : value of potential on tensor product
C                              mesh on target sphere.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      complex *16 mpole(0:nterms,-nterms:nterms)
      complex *16 phival(nquad,nquadm)
      complex *16 pot,fld(3), zk
      real *8 w(100000)
C
      integer l,m,jnew,knew
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      lw = 100000
      iffld = 1
cc      pot = dcmplx(1.0d0,1.0d0)
      do jj=1,nquad
      do kk=1,nquadm
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
	 phi = 2*pi*kk/nquadm
	 cosphi = dcos(phi)
	 sinphi = dsin(phi)
	 targ(1) = radius*stheta*cosphi
	 targ(2) = radius*stheta*sinphi
	 targ(3) = zshift + radius*ctheta
	 rn1 = stheta*cosphi
	 rn2 = stheta*sinphi
	 rn3 = ctheta
         call h3dmpeval(zk,scale,center,mpole,nterms,targ,
     1        pot,iffld,fld,w,lw,lused,ier)
         phival(jj,kk) = pot
ccc         phival(jj,kk) = -(fld(1)*rn1+fld(2)*rn2+fld(3)*rn3)
      enddo
      enddo
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dmpevalsphere(mpole,zk,scale,phival,zshift,radius,
     1           nterms,lmp,ynm,phitemp,nquad,xnodes)
C***********************************************************************
C
C     This subroutine evaluates a multipole expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C     TESTING CODE  ---- NOT USED
C
C---------------------------------------------------------------------
C     INPUT:
C
C     mpole    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : mpole scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     lmp      : dimension parameter for mpole
C     ynm      : storage for Ynm
C     phitemp  : intermediate storage for p^3 method
C     nquad    : number of quadrature nodes
C                              on target sphere is nquad*nquad.
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phival   : value of potential on tensor product
C                              mesh on target sphere.
C
C
C     NOTE: fhs and fhder hardwired at 2000 here.
C
C***********************************************************************
      implicit real *8 (a-h,o-z)
      integer nterms
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 phival(nquad,nquad)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z
      complex *16 ephi,fhs(0:2000),fhder(0:2000)
      complex *16 ephi1,ephik
C
      integer l,m,jnew,knew
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 0
      ifder = 0
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
      enddo
      enddo
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 z = zk*rj
	 call ylgndr(nterms,cthetaj,ynm)
	 call h3dall(nterms,z,scale,fhs,ifder,fhder)
         do m=-nterms,nterms
	    mabs = abs(m)
	    do n=mabs,nterms
	       phitemp(jj,m) = phitemp(jj,m) +
     1                mpole(n,m)*fhs(n)*ynm(n,mabs)
	    enddo
	 enddo
      enddo
      ephi1 = cdexp(2*pi*imag/nquad)
      ephik = 1.0d0
      do jj = 1,nquad
      do kk = 1,nquad
         phival(jj,kk) = 0.0d0
         ephik = ephik*ephi1
         ephi = ephik**(-nterms)
	 do m = -nterms,nterms
cc	    phi = 2*pi*kk*m/nquad
cc	    ephi = cdexp(imag*phi)
            phival(jj,kk) = phival(jj,kk) + phitemp(jj,m)*ephi
	    ephi = ephi*ephik
         enddo
      enddo
      enddo
      return
      end
C
C
C
C
C
C***********************************************************************
      subroutine h3dlocevalsphere(local,zk,scale,phival,zshift,radius,
     1           nterms,nterms2,lmp,ynm,phitemp,nquad,xnodes)
C***********************************************************************
C
C     This subroutine evaluates a local expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C     TESTING CODE - OBSOLETE.
C
C---------------------------------------------------------------------
C     INPUT:
C
C     local    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     nterms2  : number of terms in exp on target sphere
C     lmp      : dimension param for local expansion
C     ynm      : storage for Ynm out to nterms.
C     phitemp  : storage for phitemp in p^3 algorithm.
C     nquad    : number of quadrature nodes
C                              on target sphere is nquad*nquad.
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phival   : value of potential on tensor product
C                              mesh on target sphere.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer iscale(0:2000)
      integer nterms
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 local(0:lmp,-lmp:lmp)
      complex *16 phival(nquad,nquad)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z
      complex *16 ephi1,ephik,ephi,fjs(0:2000),fjder(0:2000)
C
      integer l,m,jnew,knew
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 0
      ifder = 0
      lwfjs = 2000
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
      enddo
      enddo
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 z = zk*rj
	 call ylgndr(nterms,cthetaj,ynm)
	 call jfuns3d(jer,nterms,z,scale,fjs,ifder,fjder,
     1        lwfjs,iscale,ntop)
         do m=-nterms2,nterms2
	    mabs = abs(m)
	    do n=mabs,nterms
	       phitemp(jj,m) = phitemp(jj,m) +
     1                local(n,m)*fjs(n)*ynm(n,mabs)
	    enddo
	 enddo
      enddo
c
      ephi1 = cdexp(2*pi*imag/nquad)
      ephik = 1.0d0
      do jj = 1,nquad
      do kk = 1,nquad
         phival(jj,kk) = 0.0d0
         ephik = ephik*ephi1
         ephi = ephik**(-nterms2)
	 do m = -nterms2,nterms2
cc	    phi = 2*pi*kk*m/nquad
cc	    ephi = cdexp(imag*phi)
            phival(jj,kk) = phival(jj,kk) + phitemp(jj,m)*ephi
	    ephi = ephi*ephik
         enddo
      enddo
      enddo
      return
      end
c
C
C***********************************************************************
      subroutine h3drescalemp(nterms,lmp,mpole,
     1           radius,zk0,scale,fhs,fhder)
C***********************************************************************
C
C     This subroutine rescales a spherical harmonic expansion
C     on the surface by h_n(zk0 * radius), converting
C     a surface function to the corresponding Hankel expansion
C     in the exterior.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           mpole = coefficients of s.h. expansion
C           radius = sphere radius
C           zk0 = Helmholtz parameter
C           scale = scale parameter for expansions
C           w       = workspace of length lw
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           mpole = rescaled by 1/h_n(zk0* radius)
C           ier = error flag
C                 0 normal return
C                 4 insufficient memory in w.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,ier
      integer l,m,jj,kk
      integer lwfhs
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 ephi,imag,emul,sum,zmul
      complex *16 fhs(0:nterms),fhder(0:nterms)
      complex *16 zk0,z
      data imag/(0.0d0,1.0d0)/
C
      z = zk0*radius
      ifder = 0
      call h3dall(nterms,z,scale,fhs,
     1             ifder,fhder)
ccc      call prin2(' fhs is *',w(ifhs),2*nterms)
      do l=0,nterms
         do m=-l,l
	    zmul = fhs(l)
	    mpole(l,m) = mpole(l,m)/zmul
         enddo
      enddo
      return
      end
c
C
C***********************************************************************
      subroutine h3drescalempstab(nterms,lmp,mpole,mpn,
     1           radius,zk0,scale,fhs,fhder)
C***********************************************************************
C
C     This subroutine takes as input the potential and its normal
C     derivative on a sphere of radius RADIUS and returns the
C     h-expansion coefficients consist with the data (in a least
C     squares sense). That is
C
C           phi  = sum mpole(n,m) h_n Y_n^m  ->
C           phin = sum mpole(n,m) zk0 *h_n' Y_n^m and
C
C           mpole(n,m) * h_n        = phi_n^m
C           mpole(n,m) * h_n' *zk0  = phin_n^m
C
C     The 1x1 normal equations are:
C
C           mpole(n,m)* ( h_n**2 + (zk0* h_n')**2) =
C
C                          h_n * phi_n^m + zk0 * h_n' * phin_n^m.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           mpole = coefficients of s.h. expansion for phi
C           mpn   = coefficients of s.h. expansion for dphi/dn
C           radius = sphere radius
C           zk0 = Helmholtz parameter
C           scale = scale parameter for expansions
C           w       = workspace of length lw
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           mpole = computed as above.
C           ier = error flag
C                 0 normal return
C                 4 insufficient memory in w.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,ier
      integer l,m,jj,kk
      integer lwfhs
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 mpn(0:lmp,-lmp:lmp)
      complex *16 ephi,imag,emul,sum,zmul
      complex *16 fhs(0:nterms),fhder(0:nterms)
      complex *16 zk0,z,zh,zhn
      data imag/(0.0d0,1.0d0)/
C
      z = zk0*radius
      ifder = 1
      call h3dall(nterms,z,scale,fhs,
     1             ifder,fhder)
      do l=0,nterms
         do m=-l,l
	    zh = fhs(l)
	    zhn = fhder(l)
	    zmul = zh*zh + zhn*zhn
	    mpole(l,m) = (zh*mpole(l,m) + zhn*mpn(l,m))/zmul
         enddo
      enddo
      return
      end
c
C
C***********************************************************************
      subroutine h3drescalestab(nterms,lmp,local,localn,
     1           radius,zk0,scale,fjs,fjder,iscale,lwfjs,ier)
C***********************************************************************
C
C     This subroutine takes as input the potential and its normal
C     derivative on a sphere of radius RADIUS and returns the
C     j-expansion coefficients consist with the data (in a least
C     squares sense). That is
C
C           phi  = sum local(n,m) j_n Y_n^m  ->
C           phin = sum local(n,m) zk0 *j_n' Y_n^m and
C
C           local(n,m) * j_n        = phi_n^m
C           local(n,m) * j_n' *zk0  = phin_n^m
C
C     The 1x1 normal equations are:
C
C           local(n,m)* ( j_n**2 + (zk0* j_n')**2) =
C
C                          j_n * phi_n^m + zk0 * j_n' * phin_n^m.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           local = coefficients of s.h. expansion of phi
C           localn = coefficients of s.h. expansion of dphi/dn
C           radius = sphere radius
C           zk0 = Helmholtz parameter
C           scale = scale parameter for expansions
C           w       = workspace of length lw
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local = computed as above.
C           ier = error flag
C                 0 normal return
C                 4 insufficient memory in w.
C                 8 lwfjs insufficient in calling jfuns3d.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,ier
      integer l,m,jj,kk
      integer lwfhs
      complex *16 fjs(0:lwfjs)
      complex *16 fjder(0:nterms)
      complex *16 local(0:lmp,-lmp:lmp)
      complex *16 localn(0:lmp,-lmp:lmp)
      complex *16 ephi,imag,emul,sum,zmul
      complex *16 zk0,z,zh,zhn
      data imag/(0.0d0,1.0d0)/
C
C
      z = zk0*radius
      ifder = 1
      call jfuns3d(ier1,nterms,z,scale,fjs,
     1             ifder,fjder,lwfjs,iscale,ntop)
      if (ier1.eq.8) then
         ier = 8
	 return
      endif
c
      do l=0,nterms
         do m=-l,l
	    zh = fjs(l)
	    zhn = fjder(l)*zk0
	    zmul = zh*zh + zhn*zhn
	    local(l,m) = (zh*local(l,m) + zhn*localn(l,m))/zmul
         enddo
      enddo
      return
      end
c
C
C
C
C
C***********************************************************************
      subroutine h3dlocevalspherestab(local,zk,scale,zshift,radius,
     1           nterms,nterms2,lmp,ynm,ynmd,phitemp,phitempn,
     2           nquad,xnodes,iscale,fjs,fjder,lwfjs,ier)
C***********************************************************************
C
C     This subroutine evaluates a local expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     local    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     nterms2  : number of terms in exp on target sphere
C     lmp      : dimension param for local expansion
C     ynm      : storage for Ynm out to nterms.
C     nquad    : number of quadrature nodes
C                              on target sphere is nquad*nquad.
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phitemp(i,j)  : jth mode of phi at ith quad node.
C     phitempn(i,j) : jth mode of phi at ith quad node.
C
C***********************************************************************
      implicit real *8 (a-h,o-z)
      integer iscale(0:lwfjs)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      real *8 ynmd(0:nterms,0:nterms)
      complex *16 local(0:lmp,-lmp:lmp)
ccc      complex *16 phitemp(nquad,-nterms:nterms)
ccc      complex *16 phitempn(nquad,-nterms:nterms)
      complex *16 phitemp(nquad,-nterms2:nterms2)
      complex *16 phitempn(nquad,-nterms2:nterms2)
      complex *16 imag,pot,fld(3), zk,z,uval,unval,ur,utheta
      complex *16 ephi1,ephik,ephi,fjs(0:lwfjs),fjder(0:lwfjs),ztmp1
      complex *16 ut1,ut2,ut3
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      ier = 0
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 1
      ifder = 1
ccc      call prinf(' nterms is *',nterms,1)
ccc      call prinf(' nterms2 is *',nterms2,1)
      do jj=1,nquad
ccc      do m=-nterms,nterms
      do m=-nterms2,nterms2
         phitemp(jj,m) = 0.0d0
         phitempn(jj,m) = 0.0d0
      enddo
      enddo
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 sthetaj = dsqrt(1.0d0-cthetaj**2)
	 rn = sthetaj*stheta + cthetaj*ctheta
	 thetan = (cthetaj*stheta - sthetaj*ctheta)/rj
	 z = zk*rj
	 call ylgndr2s(nterms,cthetaj,ynm,ynmd)
	 call jfuns3d(jer,nterms,z,scale,fjs,ifder,fjder,
     1        lwfjs,iscale,ntop)
         if (jer.eq.8) then
            ier = 8
	    return
         endif
	 do n = 0,nterms
	    fjder(n) = fjder(n)*zk
         enddo
	 do n = 1,nterms
	    do m = 1,n
	       ynm(n,m) = ynm(n,m)*sthetaj
            enddo
         enddo
	 phitemp(jj,0) = local(0,0)*fjs(0)
	 phitempn(jj,0) = local(0,0)*fjder(0)*rn
ccc         do n=1,nterms2
c
c  4/4/2010 - restoring n=1,nterms
c
         do n=1,nterms
	    phitemp(jj,0) = phitemp(jj,0) +
     1                local(n,0)*fjs(n)*ynm(n,0)
	    ut1 = fjder(n)*rn
	    ut2 = fjs(n)*thetan
	    ut3 = ut1*ynm(n,0)-ut2*ynmd(n,0)*sthetaj
	    phitempn(jj,0) = phitempn(jj,0)+ut3*local(n,0)
ccc            ur = fjder(n)*ynm(n,0)*local(n,0)
ccc            utheta = -local(n,0)*fjs(n)*ynmd(n,0)*sthetaj
ccc	    phitempn(jj,0) = phitempn(jj,0)+ur*rn+utheta*thetan
ccc	    do m=1,n
c
c  4/4/2010 - only modes up to nterms2 are needed on smaller sphere.
c
	    do m=1,min(n,nterms2)
c
	       ztmp1 = fjs(n)*ynm(n,m)
	       phitemp(jj,m) = phitemp(jj,m) +
     1                local(n,m)*ztmp1
	       phitemp(jj,-m) = phitemp(jj,-m) +
     1                local(n,-m)*ztmp1
	       ut3 = ut1*ynm(n,m)-ut2*ynmd(n,m)
	       phitempn(jj,m) = phitempn(jj,m)+ut3*local(n,m)
	       phitempn(jj,-m) = phitempn(jj,-m)+ut3*local(n,-m)
ccc               ur = fjder(n)*ynm(n,m)*local(n,m)
ccc               utheta = -local(n,m)*fjs(n)*ynmd(n,m)
ccc	       phitempn(jj,m) = phitempn(jj,m)+ur*rn+utheta*thetan
ccc               ur = fjder(n)*ynm(n,m)*local(n,-m)
ccc               utheta = -local(n,-m)*fjs(n)*ynmd(n,m)
ccc	       phitempn(jj,-m) = phitempn(jj,-m)+ur*rn+utheta*thetan
	    enddo
	 enddo
      enddo
c
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dlocevalspherestab_fast(local,zk,scale,zshift,radius,
     1           nterms,nterms2,lmp,ynm,ynmd,phitemp,phitempn,
     2           nquad,xnodes,iscale,fjs,fjder,rat1,rat2,lwfjs,ier)
C***********************************************************************
C
C     This subroutine evaluates a local expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     local    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     nterms2  : number of terms in exp on target sphere
C     lmp      : dimension param for local expansion
C     ynm      : storage for Ynm out to nterms.
C     nquad    : number of quadrature nodes
C                              on target sphere is nquad*nquad.
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phitemp(i,j)  : jth mode of phi at ith quad node.
C     phitempn(i,j) : jth mode of phi at ith quad node.
C
C***********************************************************************
      implicit real *8 (a-h,o-z)
      integer iscale(0:lwfjs)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      real *8 ynmd(0:nterms,0:nterms)
      complex *16 local(0:lmp,-lmp:lmp)
ccc      complex *16 phitemp(nquad,-nterms:nterms)
ccc      complex *16 phitempn(nquad,-nterms:nterms)
      complex *16 phitemp(nquad,-nterms2:nterms2)
      complex *16 phitempn(nquad,-nterms2:nterms2)
      complex *16 imag,pot,fld(3), zk,z,uval,unval,ur,utheta
      complex *16 ephi1,ephik,ephi,fjs(0:lwfjs),fjder(0:lwfjs),ztmp1
      complex *16 ut1,ut2,ut3
      real *8 rat1(0:nterms,0:nterms),rat2(0:nterms,0:nterms)
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      ier = 0
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 1
      ifder = 1
ccc      call prinf(' nterms is *',nterms,1)
ccc      call prinf(' nterms2 is *',nterms2,1)
      do jj=1,nquad
ccc      do m=-nterms,nterms
      do m=-nterms2,nterms2
         phitemp(jj,m) = 0.0d0
         phitempn(jj,m) = 0.0d0
      enddo
      enddo
      call ylgndrini(nterms,rat1,rat2)
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 sthetaj = dsqrt(1.0d0-cthetaj**2)
	 rn = sthetaj*stheta + cthetaj*ctheta
	 thetan = (cthetaj*stheta - sthetaj*ctheta)/rj
	 z = zk*rj
	 call ylgndr2sf(nterms,cthetaj,ynm,ynmd,rat1,rat2)
	 call jfuns3d(jer,nterms,z,scale,fjs,ifder,fjder,
     1        lwfjs,iscale,ntop)
         if (jer.eq.8) then
            ier = 8
	    return
         endif
	 do n = 0,nterms
	    fjder(n) = fjder(n)*zk
         enddo
	 do n = 1,nterms
	    do m = 1,n
	       ynm(n,m) = ynm(n,m)*sthetaj
            enddo
         enddo
	 phitemp(jj,0) = local(0,0)*fjs(0)
	 phitempn(jj,0) = local(0,0)*fjder(0)*rn
ccc         do n=1,nterms2
c
c  4/4/2010 - restoring n=1,nterms
c
         do n=1,nterms
	    phitemp(jj,0) = phitemp(jj,0) +
     1                local(n,0)*fjs(n)*ynm(n,0)
	    ut1 = fjder(n)*rn
	    ut2 = fjs(n)*thetan
	    ut3 = ut1*ynm(n,0)-ut2*ynmd(n,0)*sthetaj
	    phitempn(jj,0) = phitempn(jj,0)+ut3*local(n,0)
ccc            ur = fjder(n)*ynm(n,0)*local(n,0)
ccc            utheta = -local(n,0)*fjs(n)*ynmd(n,0)*sthetaj
ccc	    phitempn(jj,0) = phitempn(jj,0)+ur*rn+utheta*thetan
c
ccc	    do m=1,n
c
c  4/4/2010 - only modes up to nterms2 are needed on smaller sphere.
c
	    do m=1,min(n,nterms2)
c
	       ztmp1 = fjs(n)*ynm(n,m)
	       phitemp(jj,m) = phitemp(jj,m) +
     1                local(n,m)*ztmp1
	       phitemp(jj,-m) = phitemp(jj,-m) +
     1                local(n,-m)*ztmp1
	       ut3 = ut1*ynm(n,m)-ut2*ynmd(n,m)
	       phitempn(jj,m) = phitempn(jj,m)+ut3*local(n,m)
	       phitempn(jj,-m) = phitempn(jj,-m)+ut3*local(n,-m)
ccc               ur = fjder(n)*ynm(n,m)*local(n,m)
ccc               utheta = -local(n,m)*fjs(n)*ynmd(n,m)
ccc	       phitempn(jj,m) = phitempn(jj,m)+ur*rn+utheta*thetan
ccc               ur = fjder(n)*ynm(n,m)*local(n,-m)
ccc               utheta = -local(n,-m)*fjs(n)*ynmd(n,m)
ccc	       phitempn(jj,-m) = phitempn(jj,-m)+ur*rn+utheta*thetan
	    enddo
	 enddo
      enddo
c
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dlocevalspherestabold(local,zk,scale,phival,phivaln,
     1           zshift,radius,nterms,nterms2,lmp,ynm,ynmd,
ccc     2           phitemp,phitempn,nquad,xnodes)
     2         phitemp,phitempn,nquad,xnodes,iscale,fjs,fjder,lwfjs,ier)
C***********************************************************************
C
C     This subroutine evaluates a local expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     local    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     nterms2  : number of terms in exp on target sphere
C     lmp      : dimension param for local expansion
C     ynm      : storage for Ynm out to nterms.
C     phitemp  : intermediate storage in O(p^3) algorithm.
C     phitempn : intermediate storage in O(p^3) algorithm.
C     nquad    : number of quadrature nodes
C                              on target sphere is nquad*nquad.
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phival   : value of potential on tensor product
C                              mesh on target sphere.
C     phivanl  : value of normal deriv of potential on tensor product
C                              mesh on target sphere.
C
C***********************************************************************
      implicit real *8 (a-h,o-z)
      integer iscale(0:lwfjs)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      real *8 ynmd(0:nterms,0:nterms)
      complex *16 local(0:lmp,-lmp:lmp)
      complex *16 phival(nquad,nquad)
      complex *16 phivaln(nquad,nquad)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 phitempn(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z,uval,unval,ur,utheta
      complex *16 ephi1,ephik,ephi,fjs(0:lwfjs),fjder(0:lwfjs),ztmp1
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      ier = 0
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 1
      ifder = 1
ccc      call prinf(' nterms is *',nterms,1)
ccc      call prinf(' nterms2 is *',nterms2,1)
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
         phitempn(jj,m) = 0.0d0
      enddo
      enddo
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 sthetaj = dsqrt(1.0d0-cthetaj**2)
	 rn = sthetaj*stheta + cthetaj*ctheta
	 thetan = (cthetaj*stheta - sthetaj*ctheta)/rj
	 z = zk*rj
	 call ylgndr2s(nterms,cthetaj,ynm,ynmd)
	 call jfuns3d(jer,nterms,z,scale,fjs,ifder,fjder,
     1        lwfjs,iscale,ntop)
         if (jer.eq.8) then
            ier = 8
	    return
         endif
	 do n = 0,nterms
	    fjder(n) = fjder(n)*zk
         enddo
	 do n = 1,nterms
	    do m = 1,n
	       ynm(n,m) = ynm(n,m)*sthetaj
            enddo
         enddo
	 phitemp(jj,0) = local(0,0)*fjs(0)
	 phitempn(jj,0) = local(0,0)*fjder(0)*rn
ccc         do n=1,nterms2
c
c   not tested, but same fix as for h3dlocevalspherestab_fast
c
         do n=1,nterms
	    phitemp(jj,0) = phitemp(jj,0) +
     1                local(n,0)*fjs(n)*ynm(n,0)
            ur = fjder(n)*ynm(n,0)*local(n,0)
            utheta = -local(n,0)*fjs(n)*ynmd(n,0)*sthetaj
	    phitempn(jj,0) = phitempn(jj,0)+ur*rn+utheta*thetan
ccc	    do m=1,n
	    do m=1,min(n,nterms2)
	       ztmp1 = fjs(n)*ynm(n,m)
	       phitemp(jj,m) = phitemp(jj,m) +
     1                local(n,m)*ztmp1
	       phitemp(jj,-m) = phitemp(jj,-m) +
     1                local(n,-m)*ztmp1
               ur = fjder(n)*ynm(n,m)*local(n,m)
               utheta = -local(n,m)*fjs(n)*ynmd(n,m)
	       phitempn(jj,m) = phitempn(jj,m)+ur*rn+utheta*thetan
               ur = fjder(n)*ynm(n,m)*local(n,-m)
               utheta = -local(n,-m)*fjs(n)*ynmd(n,m)
	       phitempn(jj,-m) = phitempn(jj,-m)+ur*rn+utheta*thetan
	    enddo
	 enddo
      enddo
c
      do jj = 1,nquad
         ephi1 = cdexp(2*pi*imag/nquad)
         ephik = 1.0d0
         do kk = 1,nquad
         phival(jj,kk) = phitemp(jj,0)
         phivaln(jj,kk) = phitempn(jj,0)
         ephik = ephik*ephi1
         ephi = ephik
	 do m = 1,nterms2
ccc	 do m = 1,nterms
cc	    phi = 2*pi*kk*m/nquad
cc	    ephi = cdexp(imag*phi)
            phival(jj,kk)=phival(jj,kk)+phitemp(jj,m)*ephi
            phival(jj,kk)=phival(jj,kk)+phitemp(jj,-m)*dconjg(ephi)
            phivaln(jj,kk)=phivaln(jj,kk)+phitempn(jj,m)*ephi
            phivaln(jj,kk)=phivaln(jj,kk)+phitempn(jj,-m)*dconjg(ephi)
	    ephi = ephi*ephik
         enddo
      enddo
      enddo
      return
      end
C
C
C
C
C
C
C***********************************************************************
      subroutine h3dprojlocnmsep(nterms,ldl,nquadn,ntold,xnodes,wts,
     1           phitemp,local,ynm)
C***********************************************************************
C
C     compute spherical harmonic expansion on unit sphere
C     of function tabulated at nquadn*nquadm grid points.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl = dimension param for local expansion
C           nquadn = number of quadrature nodes in polar angle.
C           ntold =  number of azimuthal (m) modes in phitemp
C           xnodes = quad nodes in theta (nquadn of them)
C           wts = quad weights in theta (nquadn of them)
C           phitemp = tabulated function
C                    phitemp(i,j) = jth mode in phi at ith quad node
C                                   in theta
C           ynm     = workspace for assoc Legendre functions
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local = coefficients of s.h. expansion
C
C    NOTE:
C
C    yrecursion.f produces Ynm with a nonstandard scaling:
C    (without the 1/sqrt(4*pi)). Thus the orthogonality relation
C    is
C             \int_S  Y_nm Y_n'm'*  dA = delta(n) delta(m) * 4*pi.
C
C   This accounts for factor (1/2) = (2*pi) * 1/(4*pi) in zmul below.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,nquadn,nquadm,ier
      integer l,m,jj,kk
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk,phitemp(nquadn,-ntold:ntold)
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 ephi,imag,emul,sum,zmul,emul1
      data imag/(0.0d0,1.0d0)/
C
      pi = 4.0d0*datan(1.0d0)
c
c     initialize local exp to zero
c
      do l = 0,ldl
         do m = -l,l
	    local(l,m) = 0.0d0
         enddo
      enddo
c
c     get local exp
c
      do jj=1,nquadn
	 cthetaj = xnodes(jj)
	 call ylgndr(nterms,cthetaj,ynm)
         do m=-ntold,ntold
	    zmul = phitemp(jj,m)*wts(jj)/2
            do l=abs(m),nterms
               local(l,m) = local(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
      return
      end
c
c
c
C***********************************************************************
      subroutine h3dprojlocnmsep_fast
     $     (nterms,ldl,nquadn,ntold,xnodes,wts,
     1           phitemp,local,ynm,rat1,rat2)
C***********************************************************************
C
C     compute spherical harmonic expansion on unit sphere
C     of function tabulated at nquadn*nquadm grid points.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl = dimension param for local expansion
C           nquadn = number of quadrature nodes in polar angle.
C           ntold =  number of azimuthal (m) modes in phitemp
C           xnodes = quad nodes in theta (nquadn of them)
C           wts = quad weights in theta (nquadn of them)
C           phitemp = tabulated function
C                    phitemp(i,j) = jth mode in phi at ith quad node
C                                   in theta
C           ynm     = workspace for assoc Legendre functions
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local = coefficients of s.h. expansion
C
C    NOTE:
C
C    yrecursion.f produces Ynm with a nonstandard scaling:
C    (without the 1/sqrt(4*pi)). Thus the orthogonality relation
C    is
C             \int_S  Y_nm Y_n'm'*  dA = delta(n) delta(m) * 4*pi.
C
C   This accounts for factor (1/2) = (2*pi) * 1/(4*pi) in zmul below.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,nquadn,nquadm,ier
      integer l,m,jj,kk
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk,phitemp(nquadn,-ntold:ntold)
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 ephi,imag,emul,sum,zmul,emul1
      real *8 rat1(0:nterms,0:nterms),rat2(0:nterms,0:nterms)
      data imag/(0.0d0,1.0d0)/
C
      pi = 4.0d0*datan(1.0d0)
c
c     initialize local exp to zero
c
      do l = 0,ldl
         do m = -l,l
	    local(l,m) = 0.0d0
         enddo
      enddo
c
c     get local exp
c
      call ylgndrini(nterms,rat1,rat2)
      do jj=1,nquadn
	 cthetaj = xnodes(jj)
	 call ylgndrf(nterms,cthetaj,ynm,rat1,rat2)
         do m=-ntold,ntold
	    zmul = phitemp(jj,m)*wts(jj)/2
            do l=abs(m),nterms
               local(l,m) = local(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
      return
      end
c
C***********************************************************************
      subroutine h3dprojlocnmold(nterms,ldl,nquadn,nquadm,xnodes,wts,
     1           phival,local,ynm,w,lw,ier)
C***********************************************************************
C
C     This subroutine projects a function given on the unit sphere
C     onto spherical harmonics.
C     It is intended for use with a library of special
C     functions in which Ynm does NOT include the
C     1/sqrt(4*pi) factor. Thus, the Ynm are orthogonal
C     but not orthonormal -> the inner product value is 4*pi.
C     (There is no universal definition - physics, spectral analysis,
C     and magnetics communities all use different normalizations.)
c
C     phival(theta,phi) is the function tabulated at nquadn*nquadm pts:
C     Legendre nodes x_j = cos(theta_j) in polar angle and
C     equispaced in the azimuthal direction.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl =    defines leading dimension of local
C           nquadn = number of quadrature nodes in polar angle.
C           nquadm = number of quadrature nodes in azimuthal direction.
C           xnodes = Gauss-Legendre nodes x_j = cos theta_j
C           wts    = Gauss quadrature weights
C           phival = tabulated function
C                    phival(i,j) = phi(sin theta_j cos phi_i,
C                                      sin theta_j sin phi_i,
C                                      cos theta_j).
C           ynm    = work array to hold spherical harmonics
C           w       = workspace of length lw
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local = coefficients of s.h. expansion
C           ier = error flag
C                 0 normal return
C                 6 insufficient memory in w.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,nquadn,nquadm,ier,ldl
      integer l,m,jj,kk
      real *8 w(lw)
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk,phival(nquadn,nquadm)
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 ephi,imag,emul,sum,zmul
      data imag/(0.0d0,1.0d0)/
C
      ier = 0
      imarray = 1
      lmarray = 2*nquadn*(2*nterms+1)+7
C
      lused = imarray+lmarray
C
      if (lused.gt.lw) then
         ier=6
	 return
      endif
c
      call h3dprojloc0nm(nterms,ldl,nquadn,nquadm,xnodes,wts,phival,
     1           local,w(imarray),ynm)
      return
      end
c
C***********************************************************************
      subroutine h3dprojloc0nm(nterms,ldl,nquadn,nquadm,xnodes,wts,
     1           phival,local,marray,ynm)
C***********************************************************************
C
C     compute spherical harmonic expansion on unit sphere
C     of function tabulated at nquadn*nquadm grid points.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl = dimension param for local expansion
C           nquadn = number of quadrature nodes in polar angle.
C           nquadm = number of quadrature nodes in azimuthal direction.
C           xnodes = quad nodes in theta (nquadn of them)
C           wts = quad weights in theta (nquadn of them)
C           phival = tabulated function
C                    phival(i,j) = phi(sin theta_j cos phi_i,
C                                      sin theta_j sin phi_i,
C                                      cos theta_j).
C
C           marray  = workspace for intermediate values
C           ynm     = workspace for assoc Legendre functions
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local = coefficients of s.h. expansion
C
C    NOTE:
C
C    yrecursion.f produces Ynm with a nonstandard scaling:
C    (without the 1/sqrt(4*pi)). Thus the orthogonality relation
C    is
C             \int_S  Y_nm Y_n'm'*  dA = delta(n) delta(m) * 4*pi.
C
C   In the first loop below, you see
C
Cccc	    marray(jj,m) = sum*2*pi/nquad
C	    marray(jj,m) = sum/(2*nquad)
C
C   The latter has incorporated the 1/(4*pi) normalization factor
C   into the azimuthal quadrature weight (2*pi/nquad).
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,nquadn,nquadm,ier
      integer l,m,jj,kk
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk,phival(nquadn,nquadm)
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 marray(nquadn,-nquadm/2:nquadm/2)
      complex *16 ephi,imag,emul,sum,zmul,emul1
      data imag/(0.0d0,1.0d0)/
C
      pi = 4.0d0*datan(1.0d0)
c
c     initialize local exp to zero
c
      do l = 0,ldl
         do m = -l,l
	    local(l,m) = 0.0d0
         enddo
      enddo
c
c     create marray (intermediate array)
c
      emul1 = cdexp(imag*2*pi/nquadm)
ccc      emul = cdexp(-imag*2*nterms*pi/nquadm)
      emul = cdexp(-imag*2*(nquadm/2)*pi/nquadm)
      do m=-nquadm/2,nquadm/2
cc	    emul = cdexp(imag*m*2*pi/nquadm)
	 do jj=1,nquadn
	    sum = 0
	    ephi = emul
	    do kk = 1,nquadm
               sum = sum + phival(jj,kk)*dconjg(ephi)
	       ephi = ephi*emul
            enddo
ccc	    marray(jj,m) = sum*2*pi/nquad
	    marray(jj,m) = sum/(2*nquadm)
         enddo
	 emul = emul*emul1
      enddo
c
c     get local exp
c
      do jj=1,nquadn
	 cthetaj = xnodes(jj)
	 call ylgndr(nterms,cthetaj,ynm)
         do m=-nquadm/2,nquadm/2
	    zmul = marray(jj,m)*wts(jj)
            do l=abs(m),nterms
               local(l,m) = local(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
      return
      end
c
C
C
C
C***********************************************************************
      subroutine h3dmpevalspherenm(mpole,zk,scale,zshift,radius,
     1           nterms,lmp,ynm,phitemp,nquad,xnodes,fhs,fhder)
C***********************************************************************
C
C     This subroutine evaluates a multipole expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     mpole    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : mpole scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     lmp      : dimension param for mpole array
C     ynm      : storage for assoc Legendre functions
C     phitemp  : storage for temporary array in O(p^3) scheme
C     nquad    : number of quadrature nodes in theta
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phival   : value of potential on tensor product
C                              mesh on target sphere.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z
      complex *16 ephi1,ephik,ephi,fhs(0:nterms),fhder(0:nterms)
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 0
      ifder = 0
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
      enddo
      enddo
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 z = zk*rj
	 call ylgndr(nterms,cthetaj,ynm)
	 call h3dall(nterms,z,scale,fhs,ifder,fhder)
         do m=-nterms,nterms
	    mabs = abs(m)
	    do n=mabs,nterms
	       phitemp(jj,m) = phitemp(jj,m) +
     1                mpole(n,m)*fhs(n)*ynm(n,mabs)
	    enddo
	 enddo
      enddo
c
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dmpevalspherenm_fast(mpole,zk,scale,zshift,radius,
     1           nterms,lmp,ynm,phitemp,nquad,xnodes,fhs,fhder,
     2           rat1,rat2)
C***********************************************************************
C
C     This subroutine evaluates a multipole expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     mpole    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : mpole scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     lmp      : dimension param for mpole array
C     ynm      : storage for assoc Legendre functions
C     phitemp  : storage for temporary array in O(p^3) scheme
C     nquad    : number of quadrature nodes in theta
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phival   : value of potential on tensor product
C                              mesh on target sphere.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z
      complex *16 ephi1,ephik,ephi,fhs(0:nterms),fhder(0:nterms)
      real *8 rat1(0:nterms,0:nterms),rat2(0:nterms,0:nterms)
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 0
      ifder = 0
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
      enddo
      enddo
      call ylgndrini(nterms,rat1,rat2)
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 z = zk*rj
	 call ylgndrf(nterms,cthetaj,ynm,rat1,rat2)
	 call h3dall(nterms,z,scale,fhs,ifder,fhder)
         do m=-nterms,nterms
	    mabs = abs(m)
	    do n=mabs,nterms
	       phitemp(jj,m) = phitemp(jj,m) +
     1                mpole(n,m)*fhs(n)*ynm(n,mabs)
	    enddo
	 enddo
      enddo
c
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dmpevalspherenmold(mpole,zk,scale,phival,zshift,
     1     radius,nterms,lmp,ynm,phitemp,nquad,nquadm,xnodes,fhs,fhder)
C***********************************************************************
C
C     This subroutine evaluates a multipole expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     mpole    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : mpole scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     lmp      : dimension param for mpole array
C     ynm      : storage for assoc Legendre functions
C     phitemp  : storage for temporary array in O(p^3) scheme
C     nquad    : number of quadrature nodes in theta
C     nquadm   : number of quadrature nodes in phi
C                   total on target sphere is nquad*nquadm.
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phival   : value of potential on tensor product
C                              mesh on target sphere.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 phival(nquad,nquadm)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z
      complex *16 ephi1,ephik,ephi,fhs(0:nterms),fhder(0:nterms)
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 0
      ifder = 0
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
      enddo
      enddo
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 z = zk*rj
	 call ylgndr(nterms,cthetaj,ynm)
	 call h3dall(nterms,z,scale,fhs,ifder,fhder)
         do m=-nterms,nterms
	    mabs = abs(m)
	    do n=mabs,nterms
	       phitemp(jj,m) = phitemp(jj,m) +
     1                mpole(n,m)*fhs(n)*ynm(n,mabs)
	    enddo
	 enddo
      enddo
c
      do jj = 1,nquad
      ephi1 = cdexp(2*pi*imag/nquadm)
      ephik = 1.0d0
      do kk = 1,nquadm
         phival(jj,kk) = 0.0d0
         ephik = ephik*ephi1
         ephi = ephik**(-nterms)
	 do m = -nterms,nterms
ccc	    phi = 2*pi*kk*m/nquadm
ccc	    ephi = cdexp(imag*phi)
            phival(jj,kk) = phival(jj,kk) + phitemp(jj,m)*ephi
	    ephi = ephi*ephik
         enddo
      enddo
      enddo
c
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dmpevalspherenmstab(mpole,zk,scale,zshift,
     1           radius,nterms,lmp,ynm,ynmd,phitemp,phitempn,
     2           nquad,xnodes,fhs,fhder)
C***********************************************************************
C
C     This subroutine evaluates a multipole expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     mpole    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : mpole scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     lmp      : dimension param for mpole
C     ynm      : storage for assoc Legendre functions
C     ynmd     : storage for derivs of assoc Legendre functions
C     nquad    : number of quadrature nodes in theta
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phitemp  : nquad by (-nterms,nterms)
C     phitempn : nquad by (-nterms,nterms)
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      real *8 ynmd(0:nterms,0:nterms)
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 phitempn(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z,uval,unval,ur,utheta,ut1,ut2
      complex *16 ut3,ephi1,ephik,ephi,fhs(0:nterms),fhder(0:nterms)
      complex *16 ztmp1
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 1
      ifder = 1
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
         phitempn(jj,m) = 0.0d0
      enddo
      enddo
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 sthetaj = dsqrt(1.0d0 - cthetaj**2)
	 rn = sthetaj*stheta + cthetaj*ctheta
	 thetan = (cthetaj*stheta - ctheta*sthetaj)/rj
	 z = zk*rj
	 call ylgndr2s(nterms,cthetaj,ynm,ynmd)
	 call h3dall(nterms,z,scale,fhs,ifder,fhder)
         do i = 0,nterms
	    fhder(i) = fhder(i)*zk
         enddo
         do n=1,nterms
	    do m = 1,n
  	       ynm(n,m) = ynm(n,m)*sthetaj
              enddo
         enddo
	 phitemp(jj,0) = mpole(0,0)*fhs(0)
	 phitempn(jj,0) = mpole(0,0)*fhder(0)*rn
	 do n=1,nterms
	    phitemp(jj,0) = phitemp(jj,0) +
     1                mpole(n,0)*fhs(n)*ynm(n,0)
	    ut1 = fhder(n)*rn
	    ut2 = fhs(n)*thetan
	    ut3 = ut1*ynm(n,0)-ut2*ynmd(n,0)*sthetaj
	    phitempn(jj,0) = phitempn(jj,0)+ut3*mpole(n,0)
ccc	    ur = fhder(n)*ynm(n,0)*mpole(n,0)
ccc	    utheta = -mpole(n,0)*fhs(n)*ynmd(n,0)*sthetaj
ccc	    phitempn(jj,0) = phitempn(jj,0)+ur*rn+utheta*thetan
            do m=1,n
	       ztmp1 = fhs(n)*ynm(n,m)
	       phitemp(jj,m) = phitemp(jj,m) +
     1                mpole(n,m)*ztmp1
	       phitemp(jj,-m) = phitemp(jj,-m) +
     1                mpole(n,-m)*ztmp1
	       ut3 = ut1*ynm(n,m)-ut2*ynmd(n,m)
	       phitempn(jj,m) = phitempn(jj,m)+ut3*mpole(n,m)
	       phitempn(jj,-m) = phitempn(jj,-m)+ut3*mpole(n,-m)
cc	       ur = fhder(n)*ynm(n,m)*mpole(n,m)
cc	       utheta = -fhs(n)*ynmd(n,m)*mpole(n,m)
cc	       phitempn(jj,m) = phitempn(jj,m)+ur*rn+utheta*thetan
cc	       ur = fhder(n)*ynm(n,m)*mpole(n,-m)
cc	       utheta = -fhs(n)*ynmd(n,m)*mpole(n,-m)
cc	       phitempn(jj,-m) = phitempn(jj,-m)+ur*rn+utheta*thetan
	    enddo
	 enddo
      enddo
c
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dmpevalspherenmstab_fast(mpole,zk,scale,zshift,
     1           radius,nterms,lmp,ynm,ynmd,phitemp,phitempn,
     2           nquad,xnodes,fhs,fhder,rat1,rat2)
C***********************************************************************
C
C     This subroutine evaluates a multipole expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     mpole    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : mpole scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     lmp      : dimension param for mpole
C     ynm      : storage for assoc Legendre functions
C     ynmd     : storage for derivs of assoc Legendre functions
C     nquad    : number of quadrature nodes in theta
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phitemp  : nquad by (-nterms,nterms)
C     phitempn : nquad by (-nterms,nterms)
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      real *8 ynmd(0:nterms,0:nterms)
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 phitempn(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z,uval,unval,ur,utheta,ut1,ut2
      complex *16 ut3,ephi1,ephik,ephi,fhs(0:nterms),fhder(0:nterms)
      complex *16 ztmp1
      real *8 rat1(0:nterms,0:nterms),rat2(0:nterms,0:nterms)
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 1
      ifder = 1
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
         phitempn(jj,m) = 0.0d0
      enddo
      enddo
      call ylgndrini(nterms,rat1,rat2)
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 sthetaj = dsqrt(1.0d0 - cthetaj**2)
	 rn = sthetaj*stheta + cthetaj*ctheta
	 thetan = (cthetaj*stheta - ctheta*sthetaj)/rj
	 z = zk*rj
	 call ylgndr2sf(nterms,cthetaj,ynm,ynmd,rat1,rat2)
	 call h3dall(nterms,z,scale,fhs,ifder,fhder)
         do i = 0,nterms
	    fhder(i) = fhder(i)*zk
         enddo
         do n=1,nterms
	    do m = 1,n
  	       ynm(n,m) = ynm(n,m)*sthetaj
              enddo
         enddo
	 phitemp(jj,0) = mpole(0,0)*fhs(0)
	 phitempn(jj,0) = mpole(0,0)*fhder(0)*rn
	 do n=1,nterms
	    phitemp(jj,0) = phitemp(jj,0) +
     1                mpole(n,0)*fhs(n)*ynm(n,0)
	    ut1 = fhder(n)*rn
	    ut2 = fhs(n)*thetan
	    ut3 = ut1*ynm(n,0)-ut2*ynmd(n,0)*sthetaj
	    phitempn(jj,0) = phitempn(jj,0)+ut3*mpole(n,0)
ccc	    ur = fhder(n)*ynm(n,0)*mpole(n,0)
ccc	    utheta = -mpole(n,0)*fhs(n)*ynmd(n,0)*sthetaj
ccc	    phitempn(jj,0) = phitempn(jj,0)+ur*rn+utheta*thetan
            do m=1,n
	       ztmp1 = fhs(n)*ynm(n,m)
	       phitemp(jj,m) = phitemp(jj,m) +
     1                mpole(n,m)*ztmp1
	       phitemp(jj,-m) = phitemp(jj,-m) +
     1                mpole(n,-m)*ztmp1
	       ut3 = ut1*ynm(n,m)-ut2*ynmd(n,m)
	       phitempn(jj,m) = phitempn(jj,m)+ut3*mpole(n,m)
	       phitempn(jj,-m) = phitempn(jj,-m)+ut3*mpole(n,-m)
cc	       ur = fhder(n)*ynm(n,m)*mpole(n,m)
cc	       utheta = -fhs(n)*ynmd(n,m)*mpole(n,m)
cc	       phitempn(jj,m) = phitempn(jj,m)+ur*rn+utheta*thetan
cc	       ur = fhder(n)*ynm(n,m)*mpole(n,-m)
cc	       utheta = -fhs(n)*ynmd(n,m)*mpole(n,-m)
cc	       phitempn(jj,-m) = phitempn(jj,-m)+ur*rn+utheta*thetan
	    enddo
	 enddo
      enddo
c
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dmpevalspherenmstabold(mpole,zk,scale,phival,phivaln,
     1           zshift,radius,nterms,lmp,ynm,ynmd,phitemp,phitempn,
cc     2           nquad,nquadm,xnodes)
     2           nquad,nquadm,xnodes,fhs,fhder)
C***********************************************************************
C
C     This subroutine evaluates a multipole expansion on a target
C     sphere at a distance (0,0,zshift) from the origin of radius
C     "radius".
C
C---------------------------------------------------------------------
C     INPUT:
C
C     mpole    : coefficients of original multipole exp.
C     zk       : Helmholtz parameter
C     scale    : mpole scaling parameter
C     zshift   : shift distance along z-axis.
C     radius   : radius of sphere about (0,0,zshift)
C                              where phival is computed.
C     nterms   : number of terms in the orig. expansion
C     lmp      : dimension param for mpole
C     ynm      : storage for assoc Legendre functions
C     ynmd     : storage for derivs of assoc Legendre functions
C     phitemp  : storage for temporary array in O(p^3) scheme
C     phitempn : storage for temporary array in O(p^3) scheme
C     nquad    : number of quadrature nodes in theta
C     nquadm   : number of quadrature nodes in phi
C                       total on target sphere is nquad*nquadm.
C     xnodes   : Legendre nodes x_j = cos theta_j.
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C     phival   : value of potential on tensor product
C                              mesh on target sphere.
C     phivaln  : value of normal deriv of potential on tensor product
C                              mesh on target sphere.
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms
      integer l,m,jnew,knew
      real *8 zshift, targ(3), center(3)
      real *8 xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      real *8 ynmd(0:nterms,0:nterms)
      complex *16 mpole(0:lmp,-lmp:lmp)
      complex *16 phival(nquad,nquadm)
      complex *16 phivaln(nquad,nquadm)
      complex *16 phitemp(nquad,-nterms:nterms)
      complex *16 phitempn(nquad,-nterms:nterms)
      complex *16 imag,pot,fld(3), zk,z,uval,unval,ur,utheta,ut1,ut2
      complex *16 ut3,ephi1,ephik,ephi,fhs(0:nterms),fhder(0:nterms)
      complex *16 ztmp1
      data imag/(0.0d0,1.0d0)/
C
C----- shift along z-axis.
C      note that everything is scaled.
C
      pi = 4.0d0*datan(1.0d0)
      center(1) = 0.0d0
      center(2) = 0.0d0
      center(3) = 0.0d0
      iffld = 1
      ifder = 1
      do jj=1,nquad
      do m=-nterms,nterms
         phitemp(jj,m) = 0.0d0
         phitempn(jj,m) = 0.0d0
      enddo
      enddo
      do jj=1,nquad
	 ctheta = xnodes(jj)
	 stheta = dsqrt(1.0d0 - ctheta**2)
         rj = (zshift+ radius*ctheta)**2 + (radius*stheta)**2
         rj = dsqrt(rj)
	 cthetaj = (zshift+radius*ctheta)/rj
	 sthetaj = dsqrt(1.0d0 - cthetaj**2)
	 rn = sthetaj*stheta + cthetaj*ctheta
	 thetan = (cthetaj*stheta - ctheta*sthetaj)/rj
	 z = zk*rj
	 call ylgndr2s(nterms,cthetaj,ynm,ynmd)
	 call h3dall(nterms,z,scale,fhs,ifder,fhder)
         do i = 0,nterms
	    fhder(i) = fhder(i)*zk
         enddo
         do n=1,nterms
	    do m = 1,n
  	       ynm(n,m) = ynm(n,m)*sthetaj
              enddo
         enddo
	 phitemp(jj,0) = mpole(0,0)*fhs(0)
	 phitempn(jj,0) = mpole(0,0)*fhder(0)*rn
	 do n=1,nterms
	    phitemp(jj,0) = phitemp(jj,0) +
     1                mpole(n,0)*fhs(n)*ynm(n,0)
	    ut1 = fhder(n)*rn
	    ut2 = fhs(n)*thetan
	    ut3 = ut1*ynm(n,0)-ut2*ynmd(n,0)*sthetaj
	    phitempn(jj,0) = phitempn(jj,0)+ut3*mpole(n,0)
ccc	    ur = fhder(n)*ynm(n,0)*mpole(n,0)
ccc	    utheta = -mpole(n,0)*fhs(n)*ynmd(n,0)*sthetaj
ccc	    phitempn(jj,0) = phitempn(jj,0)+ur*rn+utheta*thetan
            do m=1,n
	       ztmp1 = fhs(n)*ynm(n,m)
	       phitemp(jj,m) = phitemp(jj,m) +
     1                mpole(n,m)*ztmp1
	       phitemp(jj,-m) = phitemp(jj,-m) +
     1                mpole(n,-m)*ztmp1
	       ut3 = ut1*ynm(n,m)-ut2*ynmd(n,m)
	       phitempn(jj,m) = phitempn(jj,m)+ut3*mpole(n,m)
	       phitempn(jj,-m) = phitempn(jj,-m)+ut3*mpole(n,-m)
cc	       ur = fhder(n)*ynm(n,m)*mpole(n,m)
cc	       utheta = -fhs(n)*ynmd(n,m)*mpole(n,m)
cc	       phitempn(jj,m) = phitempn(jj,m)+ur*rn+utheta*thetan
cc	       ur = fhder(n)*ynm(n,m)*mpole(n,-m)
cc	       utheta = -fhs(n)*ynmd(n,m)*mpole(n,-m)
cc	       phitempn(jj,-m) = phitempn(jj,-m)+ur*rn+utheta*thetan
	    enddo
	 enddo
      enddo
c
      do jj = 1,nquad
         ephi1 = cdexp(2*pi*imag/nquadm)
         ephik = 1.0d0
         do kk = 1,nquadm
            phival(jj,kk) = phitemp(jj,0)
            phivaln(jj,kk) = phitempn(jj,0)
            ephik = ephik*ephi1
            ephi = ephik
	    do m = 1,nterms
cc	       phi = 2*pi*kk*m/nquadm
cc	       ephi = exp(imag*phi)
               phival(jj,kk) = phival(jj,kk) + phitemp(jj,m)*ephi +
     1      	                    phitemp(jj,-m)*dconjg(ephi)
               phivaln(jj,kk) = phivaln(jj,kk) + phitempn(jj,m)*ephi +
     1                              phitempn(jj,-m)*dconjg(ephi)
	       ephi = ephi*ephik
            enddo
         enddo
      enddo
c
      return
      end
C
C
C
C
C***********************************************************************
      subroutine h3dprojlocsepstab(nterms,ldl,nquadn,ntold,xnodes,wts,
     1           phitemp,phitempn,local,local2,ynm)
C***********************************************************************
C
C     compute spherical harmonic expansion on unit sphere
C     of function tabulated at nquadn*nquadm grid points.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl = order of spherical harmonic expansion
C           nquadn = number of quadrature nodes in polar angle (theta).
C           ntold = number of modes in azimuthal direction.
C           xnodes = quad nodes in theta (polar angle) - nquadn of them
C           wts  = quad weights in theta (polar angle)
C           phitemp, phitempn = tabulated function and normal deriv
C                    phitemp(i,j) = jth mode of phi at ith quad node.
C                    phivaln(i,j) = jth mode of phi at ith quad node.
C
C           marray  = workspace
C           ynm     = workspace for assoc Legendre functions
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local  = coefficients of s.h. expansion for phi
C           local2 = coefficients of s.h. expansion for dphi/dn
C
C    NOTE:
C
C    yrecursion.f produces Ynm with a nonstandard scaling:
C    (without the 1/sqrt(4*pi)). Thus the orthogonality relation
C    is
C             \int_S  Y_nm Y_n'm'*  dA = delta(n) delta(m) * 4*pi.
C
C   In the first loop below, you see
C
Cccc	    marray(jj,m) = sum*2*pi/nquad
C	    marray(jj,m) = sum/(2*nquad)
C
C   The latter has incorporated the 1/(4*pi) normalization factor
C   into the azimuthal quadrature weight (2*pi/nquad).
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,nquadn,nquadm,ier
      integer l,m,jj,kk
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 local2(0:ldl,-ldl:ldl)
      complex *16 phitemp(nquadn,-ntold:ntold)
      complex *16 phitempn(nquadn,-ntold:ntold)
      complex *16 ephi,imag,emul,sum,zmul,emul1
      data imag/(0.0d0,1.0d0)/
C
      pi = 4.0d0*datan(1.0d0)
c
c     initialize local exp to zero
c
      do l = 0,ldl
         do m = -l,l
	    local(l,m) = 0.0d0
	    local2(l,m) = 0.0d0
         enddo
      enddo
c
c     get local exp
c
      do jj=1,nquadn
	 cthetaj = xnodes(jj)
	 call ylgndr(nterms,cthetaj,ynm)
         do m=-ntold,ntold
	    zmul = phitemp(jj,m)*wts(jj)/2.0d0
            do l=abs(m),nterms
               local(l,m) = local(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
c
c     get local exp
c
      do jj=1,nquadn
	 cthetaj = xnodes(jj)
	 call ylgndr(nterms,cthetaj,ynm)
         do m=-ntold,ntold
	    zmul = phitempn(jj,m)*wts(jj)/2.0d0
            do l=abs(m),nterms
               local2(l,m) = local2(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
      return
      end
c
c
c
C
C
C***********************************************************************
      subroutine h3dprojlocsepstab_fast
     $          (nterms,ldl,nquadn,ntold,xnodes,wts,
     1           phitemp,phitempn,local,local2,ynm,rat1,rat2)
C***********************************************************************
C
C     compute spherical harmonic expansion on unit sphere
C     of function tabulated at nquadn*nquadm grid points.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl = order of spherical harmonic expansion
C           nquadn = number of quadrature nodes in polar angle (theta).
C           ntold = number of modes in azimuthal direction.
C           xnodes = quad nodes in theta (polar angle) - nquadn of them
C           wts  = quad weights in theta (polar angle)
C           phitemp, phitempn = tabulated function and normal deriv
C                    phitemp(i,j) = jth mode of phi at ith quad node.
C                    phivaln(i,j) = jth mode of phi at ith quad node.
C
C           marray  = workspace
C           ynm     = workspace for assoc Legendre functions
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local  = coefficients of s.h. expansion for phi
C           local2 = coefficients of s.h. expansion for dphi/dn
C
C    NOTE:
C
C    yrecursion.f produces Ynm with a nonstandard scaling:
C    (without the 1/sqrt(4*pi)). Thus the orthogonality relation
C    is
C             \int_S  Y_nm Y_n'm'*  dA = delta(n) delta(m) * 4*pi.
C
C   In the first loop below, you see
C
Cccc	    marray(jj,m) = sum*2*pi/nquad
C	    marray(jj,m) = sum/(2*nquad)
C
C   The latter has incorporated the 1/(4*pi) normalization factor
C   into the azimuthal quadrature weight (2*pi/nquad).
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,nquadn,nquadm,ier
      integer l,m,jj,kk
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 local2(0:ldl,-ldl:ldl)
      complex *16 phitemp(nquadn,-ntold:ntold)
      complex *16 phitempn(nquadn,-ntold:ntold)
      complex *16 ephi,imag,emul,sum,zmul,emul1
      real *8 rat1(0:nterms,0:nterms),rat2(0:nterms,0:nterms)
      data imag/(0.0d0,1.0d0)/
C
      pi = 4.0d0*datan(1.0d0)
c
c     initialize local exp to zero
c
      do l = 0,ldl
         do m = -l,l
	    local(l,m) = 0.0d0
	    local2(l,m) = 0.0d0
         enddo
      enddo
c
c     get local exp
c
c
      call ylgndrini(nterms,rat1,rat2)
      do jj=1,nquadn
	 cthetaj = xnodes(jj)
	 call ylgndrf(nterms,cthetaj,ynm,rat1,rat2)
         do m=-ntold,ntold
	    zmul = phitemp(jj,m)*wts(jj)/2.0d0
            do l=abs(m),nterms
               local(l,m) = local(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
	    zmul = phitempn(jj,m)*wts(jj)/2.0d0
            do l=abs(m),nterms
               local2(l,m) = local2(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
c
      return
      end
c
c
c
C
C
C***********************************************************************
      subroutine h3dprojloc0nmstab(nterms,ldl,nquadn,nquadm,xnodes,wts,
     1           phival,phivaln,local,local2,marray,ynm)
C***********************************************************************
C
C     compute spherical harmonic expansion on unit sphere
C     of function tabulated at nquadn*nquadm grid points.
C
C---------------------------------------------------------------------
C     INPUT:
C
C           nterms = order of spherical harmonic expansion
C           ldl = order of spherical harmonic expansion
C           nquadn = number of quadrature nodes in polar angle (theta).
C           nquadm = number of quadrature nodes in azimuthal direction.
C           xnodes = quad nodes in theta (polar angle) - nquadn of them
C           wts  = quad weights in theta (polar angle)
C           phival = tabulated function
C                    phival(i,j) = phi(sin theta_j cos phi_i,
C                                      sin theta_j sin phi_i,
C                                      cos theta_j).
C           phivaln = tabulated normal deriv of function
C                    phivaln(i,j) = phi(sin theta_j cos phi_i,
C                                      sin theta_j sin phi_i,
C                                      cos theta_j).
C
C           marray  = workspace
C           ynm     = workspace for assoc Legendre functions
C
C---------------------------------------------------------------------
C     OUTPUT:
C
C           local  = coefficients of s.h. expansion for phi
C           local2 = coefficients of s.h. expansion for dphi/dn
C
C    NOTE:
C
C    yrecursion.f produces Ynm with a nonstandard scaling:
C    (without the 1/sqrt(4*pi)). Thus the orthogonality relation
C    is
C             \int_S  Y_nm Y_n'm'*  dA = delta(n) delta(m) * 4*pi.
C
C   In the first loop below, you see
C
Cccc	    marray(jj,m) = sum*2*pi/nquad
C	    marray(jj,m) = sum/(2*nquad)
C
C   The latter has incorporated the 1/(4*pi) normalization factor
C   into the azimuthal quadrature weight (2*pi/nquad).
C
C---------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      integer nterms,nquadn,nquadm,ier
      integer l,m,jj,kk
      real *8 wts(1),xnodes(1)
      real *8 ynm(0:nterms,0:nterms)
      complex *16 zk,phival(nquadn,nquadm)
      complex *16 phivaln(nquadn,nquadm)
      complex *16 local(0:ldl,-ldl:ldl)
      complex *16 local2(0:ldl,-ldl:ldl)
      complex *16 marray(nquadn,-nquadm/2:nquadm/2)
      complex *16 ephi,imag,emul,sum,zmul,emul1
      data imag/(0.0d0,1.0d0)/
C
      pi = 4.0d0*datan(1.0d0)
c
c     initialize local exp to zero
c
      do l = 0,ldl
         do m = -l,l
	    local(l,m) = 0.0d0
	    local2(l,m) = 0.0d0
         enddo
      enddo
c
c     create marray for potential (intermediate array)
c
      emul1 = cdexp(imag*2*pi/nquadm)
      emul = cdexp(-imag*2*(nquadm/2)*pi/nquadm)
      do m=-nquadm/2,nquadm/2
	 do jj=1,nquadn
	    sum = 0
	    ephi = emul
	    do kk = 1,nquadm
               sum = sum + phival(jj,kk)*dconjg(ephi)
	       ephi = ephi*emul
            enddo
	    marray(jj,m) = sum/(2*nquadm)
         enddo
	 emul = emul*emul1
      enddo
c
c     get local exp
c
      do jj=1,nquadn
	 cthetaj = xnodes(jj)
	 call ylgndr(nterms,cthetaj,ynm)
         do m=-nquadm/2,nquadm/2
	    zmul = marray(jj,m)*wts(jj)
            do l=abs(m),nterms
               local(l,m) = local(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
c
c     create marray for normal derivative (intermediate array)
c
      emul1 = cdexp(imag*2*pi/nquadm)
      emul = cdexp(-imag*2*(nquadm/2)*pi/nquadm)
      do m=-nquadm/2,nquadm/2
	 do jj=1,nquadn
	    sum = 0
	    ephi = emul
	    do kk = 1,nquadm
               sum = sum + phivaln(jj,kk)*dconjg(ephi)
	       ephi = ephi*emul
            enddo
	    marray(jj,m) = sum/(2*nquadm)
         enddo
	 emul = emul*emul1
      enddo
c
c     get local exp
c
      do jj=1,nquadn
	 cthetaj = xnodes(jj)
	 call ylgndr(nterms,cthetaj,ynm)
         do m=-nquadm/2,nquadm/2
	    zmul = marray(jj,m)*wts(jj)
            do l=abs(m),nterms
               local2(l,m) = local2(l,m) +
     1   	       zmul*ynm(l,abs(m))
            enddo
         enddo
      enddo
      return
      end
c
c
c
c

ext: f


program strided_get
  use iso_c_binding, only : DPC=>c_double
  implicit none

  integer :: i,me,np
  integer,allocatable :: a(:,:,:,:)[:],b(:,:,:,:)
  complex(kind=DPC),allocatable :: ac(:,:,:,:)[:],bc(:,:,:,:)

  me = this_image()
  np = num_images()

  allocate(ac(0:11,-10:-5,-1:0,-1:5)[*],bc(6,6,2,7))

  ac = me
  bc = me

  sync all

  if(me == 2) then
    bc(1:2,:,:,:) = ac(0:1,:,:,:)[me-1]
    if(any(bc(1:2,:,:,:) /= 1)) error stop "strided get test failed"
  endif

  sync all

  if (me == 2) write(*,*) 'Test passed.'
end program


!     -------------------------------------------------------------
!      T W O 1 1
!     -------------------------------------------------------------
!                                                                  *
!                                                                  *
!     THIS PACKAGE DETERMINES THE VALUES OF MATRIX ELEMENTS        *
!     OF ONE PARTICLE OPERATOR IN CASE :           N'1 = N1        *
!                                                                  *
!                                                  1111            *
!                                                                  *
!                                                                  *
!     Written by G. Gaigalas,                                      *
!     Vanderbilt University,  Nashville             October 1996   *
!
      SUBROUTINE TWO11(KL1, KS1, KL2, KS2, K, IA, C1)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE MEDEFN_C
      USE DIAGNL_C
      USE CONSTS_C
      USE TRK_C
      USE PERMAT_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  12:23:13  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE ixjtik_I
      USE w1_I
      USE sixj_I
      USE wwpls1_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER  :: KL1
      INTEGER  :: KS1
      INTEGER  :: KL2
      INTEGER  :: KS2
      INTEGER  :: K
      INTEGER  :: IA
      REAL(DOUBLE) , INTENT(OUT) :: C1
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      REAL(DOUBLE) :: W, SI1, SI2, WW
!-----------------------------------------------
      C1 = ZERO
      W = ZERO
      IF (IXJTIK(2*KL1,2*KL2,2*K,2*ID1(3),2*ID1(3),2*ID1(3)) /= 0) THEN
         IF (IXJTIK(2*KS1,2*KS2,2*K,1,1,1) /= 0) THEN
            CALL W1 (IK1, BK1, ID1, BD1, K, K, HALF, (-HALF), W)
            IF (DABS(W) > EPS) THEN
               CALL SIXJ (2*KL1, 2*KL2, 2*K, 2*ID1(3), 2*ID1(3), 2*ID1(3), 0, &
                  SI1)
               CALL SIXJ (1, 1, 2*K, 2*KS1, 2*KS2, 1, 0, SI2)
               W = W*SI1*SI2
            ENDIF
         ENDIF
      ENDIF
      CALL WWPLS1 (KL1, KS1, KL2, KS2, K, K, HALF, (-HALF), HALF, (-HALF), WW)
      WW = WW/DSQRT(DBLE((2*KL1 + 1)*(2*KL2 + 1)*(2*KS1 + 1)*(2*KS2 + 1)))
      C1 = HALF*RS(1,1)*RL(1,1)*(WW - W)
      RETURN
      END SUBROUTINE TWO11


text:
! RUN: %python %S/test_folding.py %s %flang_fc1
! Ensure that lower bounds are accounted for in intrinsic folding;
! this is a regression test for a bug in which they were not
module m
  real, parameter :: a(-1:-1) = 1.
  real, parameter :: b(-1:-1) = log(a)
  integer, parameter :: c(-1:1) = [33, 22, 11]
  integer, parameter :: d(1:3) = [33, 22, 11]
  integer, parameter :: e(-2:0) = ([33, 22, 11])
  logical, parameter :: test_1 = lbound((a),1)==-1 .and. lbound(b,1)==-1 .and. &
                               lbound(log(a),1)==1 .and. all(b==0)
  logical, parameter :: test_2 = all(c .eq. d)
  logical, parameter :: test_3 = all(c .eq. e)
  logical, parameter :: test_4 = all(d .eq. e)
end

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
Claire OBrien (born 25 May 1903, Sacramento, CA; died 23 Aug. 1996, Sacramento, CA) was a former teacher and viceprincipal at Davis Senior High School, as well as its principal from 1950 to 1959.  OBrien received both graduate and undergraduate degrees from UC Berkeley.  She became one of the first teachers at Davis High School when it opened in 1927, teaching classes in typing and shorthand.  After leaving the principals position, she became assistant chief in the Bureau of Business Education in Sacramento until she retired in 1980.  She lived in the Land Park neighborhood of Sacramento.

Source:

1996 Aug. 27, Obituary: Claire OBrien, Davis Enterprise

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
        MODULE MPI_CONSTANTS
        IMPLICIT NONE
        INCLUDE 'mpifnoext.h'
        END MODULE MPI_CONSTANTS

max_stars_repo_name: raffenet/mpich-CVS

max_issues_repo_name: raffenet/mpich-CVS

max_forks_repo_name: raffenet/mpich-CVS