*----------------------------------------------------------------------
      subroutine uvtang
     I  ( x, y, itr, cp, 
     O    veux, veuy, vevx, vevy  )
*----------------------------------------------------------------------
*     / calculate d(NDC)/d(WC) by numerical differenciation /
*       in case of map projection, radius is assumed to be 1.0 
*
      character cmsg*32

*     /  components of unit tangent vector in WC /
      
      call glpget ( 'rmiss', rmiss )

      if ( itr .eq. 1 ) then                         !" cartesian
         uex = 1.0
         uey = 1.0

      else if ( itr .eq. 5 ) then                    !" polar 
         if ( x .gt. 0.0 ) then
            uex = 1.0
            uey = 1.0 / (cp*x)
         else
            veux = rmiss
            veuy = rmiss
            vevx = rmiss
            vevy = rmiss
            return
*                     ----------------------------> invalid value
         end if

      else if (       itr .ge. 10
     &          .and. itr .le. 33  ) then            !" spherical-shell
         cspy = cos(y*cp)
         if( cspy .gt. 0.0 ) then
            uex = 1.0 / cspy / cp
            uey = 1.0 /        cp
         else 
            veux = rmiss
            veuy = rmiss
            vevx = rmiss
            vevy = rmiss
            return
*                     ----------------------------> invalid value
         end if

      else
         CMSG='itr not implemented'
         CALL MSGDMP('E','uetang',CMSG)
      end if

*     /  get d(NDC)/d(WC) by numerical differenciation /
*        uvecta is a small factor.

      call uvpget ( 'uvecta', uvecta )

      UDX = uvecta * uex
      UDY = uvecta * uey

      CALL STFTRF( x    , y    , VX0 , VY0  )

      CALL STFTRF( x+udx, y    , VXn , VYn  )
      CALL STFTRF( x-udx, y    , VXs , VYs  )
      CALL STFTRF( x    , y+udy, VXe , VYe  )
      CALL STFTRF( x    , y-udy, VXw , VYw  )

      d2n = ( vxn - vx0 )**2 + ( vyn - vy0 )**2
      d2s = ( vxs - vx0 )**2 + ( vys - vy0 )**2
      d2e = ( vxe - vx0 )**2 + ( vye - vy0 )**2
      d2w = ( vxw - vx0 )**2 + ( vyw - vy0 )**2

      if ( d2n .le. d2s ) then
         veux = ( vxn - vx0 ) / uvecta
         veuy = ( vyn - vy0 ) / uvecta
      else
         veux = - ( vxs - vx0 ) / uvecta
         veuy = - ( vys - vy0 ) / uvecta
      end if

      if ( d2e .le. d2w ) then
         vevx = ( vxe - vx0 ) / uvecta
         vevy = ( vye - vy0 ) / uvecta
      else
         vevx = - ( vxw - vx0 ) / uvecta
         vevy = - ( vyw - vy0 ) / uvecta
      end if

      return
      end