subroutine xy_map(stcprm, x,y, map)
      real stcprm(15)
      real map(3)
      real lnabva
      parameter (REARTH=6371.2)
        xi0 = (x - stcprm(11)) * stcprm(15) / REARTH
        eta0 = (y - stcprm(12)) * stcprm(15) / REARTH
        xi = xi0 * stcprm(13) - eta0 * stcprm(14)
        eta = xi0 * stcprm(14) + eta0 * stcprm(13)
        vsq = stcprm(1) * (xi*xi + eta*eta) - 2. * eta
        if (stcprm(1) * vsq + 1. .le. 3.e-8) then
C/* Case of point at Projection Pole */
          map(1) = 0.
          map(2) = 0.
          map(3) = sign(1.,stcprm(1))
        else
          if ( stcprm(1) .ge. 1.) then
C/* Stereographic Case, std */
            fact = 2. / (2. + vsq)
            map(2) = xi * fact
            map(1) = (1. - eta) * fact
            map(3) = fact - 1.
          else if (stcprm(1) .le. -1.) then
C/* Stereographic Case, centered at antipodes */
            fact = -2. / (2. - vsq)
            map(2) = xi * fact
            map(1) = (- 1. - eta) * fact
            map(3) = fact + 1.
          else
C/* Mercator or Lambert Case */
            ymerc = -.5 * lnabva(stcprm(1), vsq)
            if ( ymerc .gt. 0.) then
              fact = exp(-ymerc)
              cosphi = fact / (1. + fact * fact)
              cosphi = cosphi + cosphi
              map(3) = 1. - fact * cosphi
            else
              fact = exp(ymerc)
              cosphi = fact / (1. + fact * fact)
              cosphi = cosphi + cosphi
              map(3) = fact * cosphi - 1.
            endif
C/*      map.v[2] = 1. - fact * cosphi; */
            if (abs(map(3)) .lt. 1.) then
              theta = atnabv(stcprm(1), xi, eta)
              map(1) = cosphi * cos(theta)
              map(2) = cosphi * sin(theta)
            else
              map(1) = 0.
              map(2) = 0.
            endif
          endif
        endif
      return
      end