      SUBROUTINE GDSWIZ05(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
     &                    LROT,CROT,SROT)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  GDSWIZ05   GDS WIZARD FOR POLAR STEREOGRAPHIC AZIMUTHAL
C   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB GRID DESCRIPTION SECTION
C           (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63)
C           AND RETURNS ONE OF THE FOLLOWING:
C             (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES
C             (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES
C           FOR POLAR STEREOGRAPHIC AZIMUTHAL PROJECTIONS.
C           IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT
C           BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT
C           OUTPUT ELEMENTS ARE SET TO FILL VALUES.
C           THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO.
C           CALCULATIONS FOR A SPHERICAL OR ELLIPTICAL EARTH
C           ARE INVOKED ACCORDING TO THE INFORMATION IN 
C           SECTION 1, OCTET 17, BIT 2.  NOTE: ELLIPTICAL EARTH 
C           CALCULATIONS ARE BASED ON THE WGS84 DATUM.
C
C PROGRAM HISTORY LOG:
C   96-04-10  IREDELL
C   06-01-10  GAYNO    CALCULATIONS FOR ELLIPTICAL EARTH.
C
C USAGE:    CALL GDSWIZ05(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
C     &                   LROT,CROT,SROT)
C
C   INPUT ARGUMENT LIST:
C     KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
C     IOPT     - INTEGER OPTION FLAG
C                (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
C                (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
C     NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
C     FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
C                (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
C     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
C     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
C     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
C                (ACCEPTABLE RANGE: -360. TO 360.)
C     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
C                (ACCEPTABLE RANGE: -90. TO 90.)
C     LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
C
C   OUTPUT ARGUMENT LIST:
C     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
C     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
C     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
C     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
C     NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
C     CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
C     SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
C                (UGRID=CROT*UEARTH-SROT*VEARTH;
C                 VGRID=SROT*UEARTH+CROT*VEARTH)
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C
C$$$
      IMPLICIT NONE
      INTEGER         :: IM, JM, IROT, KGDS(200), IPROJ
      INTEGER         :: ISCAN, JSCAN, NSCAN, NRET
      INTEGER         :: ITER, NPTS, IOPT, N, LROT
      REAL            :: XPTS(NPTS),YPTS(NPTS),RLON(NPTS),RLAT(NPTS)
      REAL            :: CROT(NPTS),SROT(NPTS)
      REAL, PARAMETER :: RERTH=6.3712E6
      REAL, PARAMETER :: RERTH_WGS84=6.378137E6
      REAL, PARAMETER :: PI=3.14159265358979
      REAL, PARAMETER :: DPR=180./PI
      REAL, PARAMETER :: PI2=PI/2.0
      REAL, PARAMETER :: PI4=PI/4.0
      REAL, PARAMETER :: E2=.00669437999013  ! wgs84 datum
      REAL            :: E=0.0          ! initialized below w/ intrinsic
      REAL            :: E_OVER_2=0.0   ! initialized below
      REAL, PARAMETER :: SLAT=60.0  ! standard latitude according
                                    ! to grib standard
      REAL, PARAMETER :: SLATR=SLAT/DPR
      REAL            :: MC, RLAT1, RLON1, ORIENT, DX, DY, H, HI, HJ
      REAL            :: DXS, DYS, DE, DR, XP, YP, DE2, ALAT, ALONG
      REAL            :: T, RHO, TC, XMIN, XMAX, YMIN, YMAX, DI, DJ
      REAL            :: ALAT1, DR2, DIFF, FILL
      LOGICAL         :: ELLIPTICAL
C tms moved initialization here pgf90 won't allow use of intrinsic during
C     declaration
      E=SQRT(E2)      !E2 set in declariation
      E_OVER_2=E*0.5
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      IF(KGDS(1).EQ.005) THEN
        ELLIPTICAL=MOD(KGDS(6)/64,2).EQ.1
        IM=KGDS(2)
        JM=KGDS(3)
        RLAT1=KGDS(4)*1.E-3
        RLON1=KGDS(5)*1.E-3
        IROT=MOD(KGDS(6)/8,2)
        ORIENT=KGDS(7)*1.E-3
        DX=KGDS(8)
        DY=KGDS(9)
        IPROJ=MOD(KGDS(10)/128,2)
        ISCAN=MOD(KGDS(11)/128,2)
        JSCAN=MOD(KGDS(11)/64,2)
        NSCAN=MOD(KGDS(11)/32,2)
        H=(-1.)**IPROJ
        IF(H.EQ.-1)ORIENT=ORIENT+180.
        HI=(-1.)**ISCAN
        HJ=(-1.)**(1-JSCAN)
        DXS=DX*HI
        DYS=DY*HJ
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C FIND X/Y OF POLE
        IF (.NOT.ELLIPTICAL) THEN
          DE=(1.+SIN(SLAT/DPR))*RERTH
          DR=DE*COS(RLAT1/DPR)/(1+H*SIN(RLAT1/DPR))
          XP=1-H*SIN((RLON1-ORIENT)/DPR)*DR/DXS
          YP=1+COS((RLON1-ORIENT)/DPR)*DR/DYS
          DE2=DE**2
        ELSE
          ALAT=H*RLAT1/DPR
          ALONG = (RLON1-ORIENT)/DPR
          T=TAN(PI4-ALAT/2.)/((1.-E*SIN(ALAT))/
     &      (1.+E*SIN(ALAT)))**(E_OVER_2)
          TC=TAN(PI4-SLATR/2.)/((1.-E*SIN(SLATR))/
     &      (1.+E*SIN(SLATR)))**(E_OVER_2)
          MC=COS(SLATR)/SQRT(1.0-E2*(SIN(SLATR)**2))
          RHO=RERTH_WGS84*MC*T/TC
          YP = 1.0 + RHO*COS(H*ALONG)/DYS
          XP = 1.0 - RHO*SIN(H*ALONG)/DXS
        ENDIF ! ELLIPTICAL
        XMIN=0
        XMAX=IM+1
        YMIN=0
        YMAX=JM+1
        NRET=0
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  TRANSLATE GRID COORDINATES TO EARTH COORDINATES
        IF(IOPT.EQ.0.OR.IOPT.EQ.1) THEN
          IF(.NOT.ELLIPTICAL)THEN
            DO N=1,NPTS
              IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND.
     &           YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN
                DI=(XPTS(N)-XP)*DXS
                DJ=(YPTS(N)-YP)*DYS
                DR2=DI**2+DJ**2
                IF(DR2.LT.DE2*1.E-6) THEN
                  RLON(N)=0.
                  RLAT(N)=H*90.
                ELSE
                  RLON(N)=MOD(ORIENT+H*DPR*ATAN2(DI,-DJ)+3600,360.)
                  RLAT(N)=H*DPR*ASIN((DE2-DR2)/(DE2+DR2))
                ENDIF
                NRET=NRET+1
                IF(LROT.EQ.1) THEN
                  IF(IROT.EQ.1) THEN
                    CROT(N)=H*COS((RLON(N)-ORIENT)/DPR)
                    SROT(N)=SIN((RLON(N)-ORIENT)/DPR)
                  ELSE
                    CROT(N)=1
                    SROT(N)=0
                  ENDIF
                ENDIF
              ELSE
                RLON(N)=FILL
                RLAT(N)=FILL
              ENDIF
            ENDDO
          ELSE ! ELLIPTICAL
            DO N=1,NPTS
              IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND.
     &           YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN
                DI=(XPTS(N)-XP)*DXS
                DJ=(YPTS(N)-YP)*DYS
                RHO=SQRT(DI*DI+DJ*DJ)
                T=(RHO*TC)/(RERTH_WGS84*MC)
                IF(ABS(YPTS(N)-YP)<0.01)THEN
                  IF(DI>0.0) ALONG=ORIENT+H*90.0
                  IF(DI<=0.0) ALONG=ORIENT-H*90.0
                ELSE
                  ALONG=ORIENT+H*ATAN(DI/(-DJ))*DPR
                  IF(DJ>0) ALONG=ALONG+180.
                END IF
                ALAT1=PI2-2.0*ATAN(T)
                DO ITER=1,10
                  ALAT = PI2 - 2.0*ATAN(T*(((1.0-E*SIN(ALAT1))/
     &                  (1.0+E*SIN(ALAT1)))**(E_OVER_2)))
                  DIFF = ABS(ALAT-ALAT1)*DPR
                  IF (DIFF < 0.000001) EXIT
                  ALAT1=ALAT
                ENDDO
                RLAT(N)=H*ALAT*DPR
                RLON(N)=ALONG
                IF(RLON(N)<0.0) RLON(N)=RLON(N)+360.
                IF(RLON(N)>360.0) RLON(N)=RLON(N)-360.0
                NRET=NRET+1
                IF(LROT.EQ.1) THEN
                  IF(IROT.EQ.1) THEN
                    CROT(N)=H*COS((RLON(N)-ORIENT)/DPR)
                    SROT(N)=SIN((RLON(N)-ORIENT)/DPR)
                  ELSE
                    CROT(N)=1
                    SROT(N)=0
                  ENDIF
                ENDIF
              ELSE
                RLON(N)=FILL
                RLAT(N)=FILL
              ENDIF
            ENDDO
          ENDIF ! ELLIPTICAL
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  TRANSLATE EARTH COORDINATES TO GRID COORDINATES
        ELSEIF(IOPT.EQ.-1) THEN
          IF(.NOT.ELLIPTICAL)THEN
            DO N=1,NPTS
              IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90.AND.
     &                                   H*RLAT(N).NE.-90) THEN
                DR=DE*TAN((90-H*RLAT(N))/2/DPR)
                XPTS(N)=XP+H*SIN((RLON(N)-ORIENT)/DPR)*DR/DXS
                YPTS(N)=YP-COS((RLON(N)-ORIENT)/DPR)*DR/DYS
                IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND.
     &             YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN
                  NRET=NRET+1
                  IF(LROT.EQ.1) THEN
                    IF(IROT.EQ.1) THEN
                      CROT(N)=H*COS((RLON(N)-ORIENT)/DPR)
                      SROT(N)=SIN((RLON(N)-ORIENT)/DPR)
                    ELSE
                      CROT(N)=1
                      SROT(N)=0
                    ENDIF
                  ENDIF
                ELSE
                  XPTS(N)=FILL
                  YPTS(N)=FILL
                ENDIF
              ELSE
                XPTS(N)=FILL
                YPTS(N)=FILL
              ENDIF
            ENDDO
          ELSE ! ELLIPTICAL CALCS
            DO N=1,NPTS
              IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90.AND.
     &                                   H*RLAT(N).NE.-90) THEN
                ALAT = H*RLAT(N)/DPR
                ALONG = (RLON(N)-ORIENT)/DPR
                T=TAN(PI4-ALAT*0.5)/((1.-E*SIN(ALAT))/
     &            (1.+E*SIN(ALAT)))**(E_OVER_2)
                RHO=RERTH_WGS84*MC*T/TC
                XPTS(N)= XP + RHO*SIN(H*ALONG) / DXS
                YPTS(N)= YP - RHO*COS(H*ALONG) / DYS
                IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND.
     &             YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN
                  NRET=NRET+1
                  IF(LROT.EQ.1) THEN
                    IF(IROT.EQ.1) THEN
                      CROT(N)=H*COS((RLON(N)-ORIENT)/DPR)
                      SROT(N)=SIN((RLON(N)-ORIENT)/DPR)
                    ELSE
                      CROT(N)=1
                      SROT(N)=0
                    ENDIF
                  ENDIF
                ELSE
                  XPTS(N)=FILL
                  YPTS(N)=FILL
                ENDIF
              ELSE
                XPTS(N)=FILL
                YPTS(N)=FILL
              ENDIF
            ENDDO
          ENDIF
        ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  PROJECTION UNRECOGNIZED
      ELSE
        NRET=0
        IF(IOPT.GE.0) THEN
          DO N=1,NPTS
            RLON(N)=FILL
            RLAT(N)=FILL
          ENDDO
        ENDIF
        IF(IOPT.LE.0) THEN
          DO N=1,NPTS
            XPTS(N)=FILL
            YPTS(N)=FILL
          ENDDO
        ENDIF
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END
