subroutine proj_3d(stcprm, point, vect, enx,eny,enz)
C/*
C *  At a given point, resolves the components of a vector vect in the
C *  local coordinate system of the map projection.  It is assumed that
C *  vect and point is given in the 3-dimensional geocentric _map_
C *  coordinate system, rather than the North centered _geo_ system.
C *  returns these components as enx, eny, and enz.  Replaces vect with
C *  its projection on the tangent plane at point.
C */
      real stcprm(15)
      real point(3),vect(3)
      dot_pr = 0.
      do k=1,3
        dot_pr = dot_pr + point(k)*vect(k)
      enddo
C/*
C *  dot_prod is the local vertical component of vect.  Note, point is
C *  assumed to be a unit vector.
C */
      enz = dot_pr
      do k=1,3
        vect(k) = vect(k) - dot_pr * point(k)
      enddo
C/*  vect has now been projected to a plane tangent to point.
C */
      fact = 1. + point(3)
      xi = vect(1) - vect(3) * point(1) / fact
      eta = vect(2) - vect(3) * point(2) / fact
C/*
C *  xi, eta represent components of vector on the projection plane,
C *  i.e. the 2-dimensional map system.
C */
      fact = stcprm(1) -1.
      if (fact .lt. 0.) then
C/*
C *  For Lambert Conformal and Mercator projections (gamma < 1.0) ,
C *  a rotation of the vector components is needed.
C */
        if ( abs(point(3)) .lt. 1.) then
          theta = fact * atan2(point(2),point(1))
          cgthta = cos(theta)
          sgthta = sin(theta)
          fact = xi * cgthta - eta * sgthta
          eta = eta * cgthta + xi * sgthta
          xi = fact
        endif
      endif
C/*
C *  Now rotate xi, eta to the final map projection direction.
C */
      enx = eta * stcprm(13) - xi * stcprm(14)
      eny = - xi * stcprm(13) - eta * stcprm(14)
      return
      end