*-----------------------------------------------------------------------
*     uvVECT
*-----------------------------------------------------------------------
      SUBROUTINE uvVECT(U,MU,V,MV,NX,NY)

      parameter ( iundef = 9999 )

      REAL      U(MU,*),V(MV,*)

      real      uxp(1000), uyp(1000)

      CHARACTER CMSG*80
      LOGICAL   LMISS,LNRMAL,LSETX,LSETY,LMSG,LEQRAT,
     +          LMISSP,LSMALL,LUNIT,LUMSG
      logical   ldeg

      logical   linvis

      equivalence ( vrcnst, ivcnst )

      SAVE

*-----------------------------------------------------------------------
*     / GET INTERNAL PARAMETERS /
*-----------------------------------------------------------------------

      CALL GLPGET('LMISS   ',LMISS )
      CALL GLPGET('RMISS   ',RMISS )

      call glpget('REALMAX ',REMAX )
      call glpget('REALMIN ',REMIN )

      CALL uvPGET('INDEX   ',INDEX )
      CALL uvPGET('LNRMAL  ',LNRMAL)
      CALL uvPGET('LEQRAT  ',LEQRAT)
      CALL uvPGET('LMSG    ',LMSG  )
      CALL uvPGET('LUNIT   ',LUNIT )
      CALL uvPGET('LUMSG   ',LUMSG )
      CALL uvPGET('ICENT   ',ICENT )
      CALL uvPGET('LMISSP  ',LMISSP)
      CALL uvPGET('ITYPE1  ',ITYPE1)
      CALL uvPGET('LSMALL  ',LSMALL)
      CALL uvPGET('RSMALL  ',RSMALL)
      CALL uvPGET('ITYPE2  ',ITYPE2)
      CALL uvPGET('RSIZEM  ',RSIZEM)
      CALL uvPGET('RSIZET  ',RSIZET)
      CALL uvPGET('XTTL    ',XTTL  )
      CALL uvPGET('IXINT   ',IXINT )
      CALL uvPGET('IYINT   ',IYINT )

      call uvpget('linvis  ',linvis )
      call uvpget('vvpmax  ',vvpmax )
      call uvpget('iarlen  ',iarlen )
      call uvpget('vrcnst  ',vrcnst )

*-----------------------------------------------------------------------
*     / get info. on coordinate and mapping /
*-----------------------------------------------------------------------

      call sgpget ( 'LDEG', ldeg )

      if ( ldeg ) then
         cp = 3.14159265 / 180.0
      else
         cp = 1.0
      end if

      call sgqtrn ( itr )
      call sgqwnd ( uxmn, uxmx, uymn, uymx )
      call sgqvpt ( vxmn, vxmx, vymn, vymx )

      if ( itr .ge. 5  .and.  itr .le. 7 ) then
         call sgqsim ( aotg, vxoff, vyoff )
      else if ( itr .ge. 10   .and.  itr .le. 13 ) then
         call sgqsim ( radius, vxoff, vyoff )
         call sgqmpl ( xlon0, ylat0, rot0 )
      else if ( itr .ge. 20   .and.  itr .le. 23 ) then
         call sgqsim ( radius, vxoff, vyoff )
         call sgqmpl ( xlon0, ylat0, rot0 )
      else if ( itr .ge. 30   .and.  itr .le. 33 ) then
         call sgqsim ( radius, vxoff, vyoff )
         call sgqmpl ( xlon0, ylat0, rot0 )
      end if

      if ( itr .ge. 2  .and.  itr .le. 4 ) then
         CMSG='LOG coodinate not supported'
         CALL MSGDMP('E','uvVECT',CMSG)
      end if

      if ( itr .eq. 1   .and.  iarlen .eq. 3 ) then
         CMSG=' iarlen=3 is not good with itr=1 '
         CALL MSGDMP('W','uvVECT',CMSG)
      end if

         
*-----------------------------------------------------------------------
*     / get ( and set if necessary ) info. on grid points /
*-----------------------------------------------------------------------

*     / SET GRID ATTRIBUTE IF IT HAS NOT BEEN SET YET /

      CALL UWQGXZ(LSETX)
      CALL UWQGYZ(LSETY)

      IF (.NOT.LSETX) THEN
         CALL UWSGXB(UXMN,UXMX,NX)
         CALL UWSGXZ(.FALSE.)
      END IF

      IF (.NOT.LSETY) THEN
         CALL UWSGYB(UYMN,UYMX,NY)
         CALL UWSGYZ(.FALSE.)
      END IF

*     / get grid position on WC /

      call uwqgxa ( uxp, nx )
      call uwqgya ( uyp, ny )

*-----------------------------------------------------------------------
*     / CHECK INPUT DATA ( missing, maximum, minimum ) /
*-----------------------------------------------------------------------

      ndata = 0
      rminu = REMAX
      rminv = REMAX
      rmaxu = REMIN
      rmaxv = REMIN

      DO 10 J=1,NY,IYINT
      DO 10 I=1,NX,IXINT
         if( .NOT. ( LMISS .AND. 
     &              (U(I,J).EQ.RMISS .OR. V(I,J).EQ.RMISS)) ) then
            ndata = ndata + 1
            RMAXU = max ( U(I,J), rmaxu )
            RMAXV = max ( V(I,J), rmaxv )
            RMINU = min ( U(I,J), rminu )
            RMINV = min ( V(I,J), rminv )
         end if
   10 CONTINUE

*     / WRITE MESSAGE and return IF MISSING FIELD /

      IF (  ndata .eq. 0 ) then
         CMSG='MISSING FIELD.'
         CALL MSGDMP('W','uvVECT',CMSG)
         IF (LMSG) THEN
            call uvwmsg( cmsg )
         END IF
         RETURN
*                                                     -------> exit
      END IF

*     / WRITE MESSAGE and return ZERO FIELD /

      IF (          RMINU.EQ.0 .AND. RMAXU.EQ.0
     +        .AND. RMINV.EQ.0 .AND. RMAXV.EQ.0 ) THEN

         CMSG='ZERO FIELD.'
         CALL MSGDMP('W','uvVECT',CMSG)
         IF (LMSG) THEN
            call uvwmsg( cmsg )
         END IF
         RETURN
*                                                     -------> exit
      END IF

*-----------------------------------------------------------------------
*     / CALCULATE NORMALIZATION FACTOR if necessary /
*-----------------------------------------------------------------------

      IF (LNRMAL) THEN

         UE=(VXMX-VXMN)/(NX/IXINT)
         VE=(VYMX-VYMN)/(NY/IYINT)

         if ( itr .le. 4 ) then
            AU=RGNLE(UE/MAX(ABS(RMINU),ABS(RMAXU)))
            AV=RGNLE(VE/MAX(ABS(RMINV),ABS(RMAXV)))
            IF (LEQRAT) THEN
               AU=MIN(AU,AV)
               AV=au
            END IF
         else if ( itr.ge.5  .and.  itr.le.7 ) then
            if ( leqrat ) then
               rumax = max( max( abs(rminu), abs(rmaxu) ),
     &                      max( abs(rminv), abs(rmaxv) )  )
               aa = sqrt ( ue * ve ) / rumax
               au = rgnle(aa)
               av = au
            else
               cmsg = 'LEQRAT must be TRUE if 5 <= ITR <= 7 '
               call msgdmp('W','uvVECT',cmsg)
               return
*                                                     -------> exit
            end if
         else if ( itr.ge.10 .and. itr.le.33 ) then
            if ( leqrat ) then
               rumax = max( max( abs(rminu), abs(rmaxu) ),
     &                      max( abs(rminv), abs(rmaxv) )  )
               aa = 3.1416 * sqrt ( ue * ve ) * radius / rumax
               au = rgnle(aa)
               av = au
            else
               cmsg = 'LEQRAT must be TRUE if 10<=ITR<=33 '
               call msgdmp('W','uvVECT',cmsg)
               return
*                                                     -------> exit
            end if
         end if
      ELSE
         CALL uvPGET('XFACT1  ',XFACT1)
         CALL uvPGET('YFACT1  ',YFACT1)
         AU=XFACT1
         AV=YFACT1
      END IF

      CALL uvPSET('XFACT2  ',AU)
      CALL uvPSET('YFACT2  ',AV)

*-----------------------------------------------------------------------
*     / output info. on scaling  /
*-----------------------------------------------------------------------

*     / determine length of the UNIT VECTOR and draw them /

      IF (LUNIT) THEN

         CALL uvUNIT

*     / WRITE UNIT VALUE /

         IF (LUMSG) THEN
            if ( iarlen .eq. 1  .or.  iarlen .eq. 3 ) then
               CALL uvPGET('UXUNIT  ',UXUNIT)
               CALL uvPGET('UYUNIT  ',UYUNIT)
               CMSG='XUNIT =##########, YUNIT =##########'
               WRITE(CMSG( 8:17),'(1P,E10.3)') UXUNIT
               WRITE(CMSG(27:36),'(1P,E10.3)') UYUNIT
            else if ( iarlen .eq. 2 ) then
               cmsg = ' XUNIT = null, YUNIT = null '
            end if
            call uvwmsg( cmsg )
         END IF

*     / DRAW TITLE FOR UNIT VECTOR /

         CALL uvDUT

      END IF

*     / WRITE SCALING FACTOR /

      IF (LMSG) THEN
         if ( iarlen .eq. 1 .or. iarlen .eq. 3 ) then
            CMSG='XFACT =##########, YFACT =##########'
            WRITE(CMSG( 8:17),'(1P,E10.3)') AU
            WRITE(CMSG(27:36),'(1P,E10.3)') AV
         else if ( iarlen .eq. 2 ) then
            cmsg = 'XFACT = NULL , YFACT = NULL'
         end if
         call uvwmsg ( cmsg )
      END IF

*-----------------------------------------------------------------------
*     / DRAW VECTORS /
*-----------------------------------------------------------------------
*     / SET line type of arrow /

*     / SET the length of arrow for direction-only-plot /

      if ( ivcnst .eq. iundef ) then
         vrcnst = 0.5  / sqrt( (real(nx)/ixint)*(real(ny)/iyint) )
      else
         call uvpget ( 'VRCNST', vrcnst )
      end if
*
      DO 20 J=1,NY,IYINT
      DO 20 I=1,NX,IXINT
*
         CALL STFTRF( uxp(i), uyp(j), VXSTRT, VYSTRT )

*     / skip missing value and mark the point if required /

         if ( lmiss .and. 
     &      (u(i,j).eq.rmiss .or. v(i,j).eq.rmiss) ) then

            IF (LMISSP) THEN
               if ( vxstrt .ne. rmiss .and.
     &              vystrt .ne. rmiss       ) then
                  CALL SGPMZV
     &                (1,VXstrt,VYstrt,ITYPE1,INDEX,RSIZEM)
               end if
            END IF
            go to 20
         end if

*     / skip if the grid is on singular point of mapping /

         IF ( LMISS )THEN
            CALL STITRF( VXSTRT, VYSTRT, UXSTRT, UYSTRT )
            IF ( ( UXSTRT .EQ. RMISS ) .OR. ( UYSTRT .EQ. RMISS ) )THEN
               go to 20
            ENDIF
         end if

*     / skip far-side point if orthographic projection /

         if ( itr.eq.30 ) then
            xradc = xlon0 * cp
            yradc = 3.14159265/2.0 - ylat0 * cp
            xradp = uxp(i) * cp
            yradp = 3.14159265/2.0 - uyp(j) * cp
            rtemp = 1.0
            call ct3sc( rtemp, yradc, xradc, xc, yc, zc )
            call ct3sc( rtemp, yradp, xradp, xp, yp, zp )
            zz = xc*xp + yc*yp + zc*zp
            if ( zz .lt. 0.0 ) then
               go to 20
            end if
         end if

*     / get unit (x,y)-tangent vector in NDC /

         call uvtang
     I  ( uxp(i), uyp(j), itr, cp, 
     O    veux, veuy, vevx, vevy )

*     / skip singular point /

         if ( lmiss .and. 
     &      ( veux .eq. rmiss .or. veuy .eq. rmiss .or.
     &        vevx .eq. rmiss .or. vevy .eq. rmiss ) ) then
            go to 20
         end if

*     / vector arrow component in NDC /

         if ( itr.ge.10  .and.  itr.le.13 ) then
            vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) / radius
            vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) / radius
         else if ( itr.ge.20  .and.  itr.le.23 ) then
            vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) / radius
            vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) / radius
         else if ( itr.ge.30  .and.  itr.le.33 ) then
            vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) / radius
            vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) / radius
         else if ( itr.ge.5  .and.  itr.le.7 ) then
            vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) / aotg
            vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) / aotg
         else 
            vdx = ( u(i,j) * au * veux + v(i,j) * av * vevx ) 
            vdy = ( u(i,j) * au * veuy + v(i,j) * av * vevy ) 
         end if

         if      ( iarlen .eq. 1 ) then
            tabs = sqrt( (au*u(i,j))**2 + (av*v(i,j))**2 )
            vabs = sqrt( vdx**2 + vdy**2 )
            vdx = tabs * vdx / vabs
            vdy = tabs * vdy / vabs
         else if ( iarlen .eq. 2 ) then
            vabs = sqrt(vdx**2 + vdy**2)
            if ( vabs .ne. 0 ) then
               vdx = vrcnst * vdx / vabs
               vdy = vrcnst * vdy / vabs
            end if
         else if ( iarlen .eq. 3 ) then
            continue
         else
            CMSG='iarlen not implemented'
            CALL MSGDMP('E','uvvect',CMSG)
         end if

*     / skip if vector is too long or too short /

         rlvec2 = vdx**2 + vdy**2

         if (  rlvec2 .gt. vvpmax**2 ) then
            CALL SGPMZV(1,VXstrt,VYstrt,ITYPE2,INDEX,RSIZEM)
            go to 20
         else IF (LSMALL .AND. rlvec2.LE.RSMALL**2) THEN
            CALL SGPMZV(1,VXstrt,VYstrt,ITYPE2,INDEX,RSIZEM)
            go to 20
         END IF

*     / do centering of arrow, and plot it ! /

      VX1 = VXstrt - vdx*(ICENT+1)/2.0
      VX2 = VXstrt - vdx*(ICENT-1)/2.0
      VY1 = VYstrt - vdy*(ICENT+1)/2.0
      VY2 = VYstrt - vdy*(ICENT-1)/2.0

c$$$      CALL UFLNZV(VX1,VY1,VX2,VY2,INDEX)
      CALL SGLAZV(VX1,VY1,VX2,VY2,1,INDEX)

   20 CONTINUE

      return
      END

