c
c     ==========================================
      subroutine mapc2p(x1,y1,z1,xp,yp,zp)
c     ==========================================
c
c     Map computational point (x1,y1,z1) to the point (xp,yp,zp)
c     in the spherical shell with inner radius Rinner and depth dR.
c     (x1,y1) gets mapped to surface of sphere and then z1*dR is distance
c     outwards from inner sphere,  0 <= z1 <= 1.
c     
c     (x1,y1) points in [-1,1]x[-1,1]  are mapped to the upper hemisphere
c     (x1,y1) points in [-3,-1]x[-1,1] are mapped to the lower hemisphere
c
      implicit real*8(a-h,o-z)
      common /comsphere/ Rinner,dR   !# must be set elsewhere or here

      xc = x1
      yc = y1

c     # ghost cell values outside of [-3,1]x[-1,1] get mapped to other
c     # hemisphere:
      if (xc.gt.1.d0) xc = xc-4.d0
      if (xc.lt.-3.d0) xc = xc+4.d0
      if (yc.gt.1.d0) then 
         yc = 2.d0 - yc
         xc = -2.d0 - xc
         endif
      if (yc.lt.-1.d0) then
         yc = -2.d0 - yc
         xc = -2.d0 - xc
         endif


      if (xc.lt.-1.d0) then
c         # points in [-3,-1] map to lower hemisphere - reflect about x=-1
c         # to compute x,y mapping and set sgnz appropriately:
          xc = -2.d0 - xc
          sgnz = -1.d0
         else
          sgnz = 1.d0
         endif
      sgnxc = dsign(1.d0,xc)
      sgnyc = dsign(1.d0,yc)

      xc1 = dabs(xc)
      yc1 = dabs(yc)
      d = dmax1(dmax1(xc1,yc1), 1.d-10)

      d1 = Rinner*d*(2.d0 - d) / dsqrt(2.d0)

      center = d1 - dsqrt(dmax1(Rinner**2 - d1**2, 0.d0))
      xp = d1/d * xc1
      yp = d1/d * yc1

      if (yc1 .gt. xc1) then
          yp = center + dsqrt(dmax1(Rinner**2 - xp**2, 0.d0))
        else
          xp = center + dsqrt(dmax1(Rinner**2 - yp**2, 0.d0))
        endif

      zp = dsqrt(dmax1(Rinner**2 - (xp**2 + yp**2), 0.d0))
      r2 = (Rinner + z1*dR) / Rinner

      xp = r2*xp*sgnxc
      yp = r2*yp*sgnyc
      zp = r2*zp*sgnz
c      
      return
      end

