C-----------------------------------------------------------------------
      SUBROUTINE POLATES5(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
     &                    NO,RLAT,RLON,IBO,LO,GO,IRET)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  POLATES5   INTERPOLATE SCALAR FIELDS (FILTERED)
C   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS FILTERED INTERPOLATION
C           FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
C           IT REQUIRES THAT NO INPUT FIELDS HAVE BITMAPS (IBI=0).
C           OPTIONS ALLOW CHOICES BETWEEN STRAIGHT BICUBIC (IPOPT(1)=0)
C           AND CONSTRAINED BICUBIC (IPOPT(1)=1) WHERE THE VALUE IS
C           CONFINED WITHIN THE RANGE OF THE SURROUNDING 4 POINTS.
C           BILINEAR USED WITHIN ONE GRID LENGTH OF BOUNDARIES.
C           ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
C           THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS
C           (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63).
C           THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS:
C             (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
C             (KGDS(1)=001) MERCATOR CYLINDRICAL
C             (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
C             (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
C             (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
C             (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)
C           WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO.
C           AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
C           AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
C           ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
C           IF KGDSO(1)<0, IN WHICH CASE THE NUMBER OF POINTS
C           AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT.
C           OUTPUT BITMAPS WILL ONLY BE CREATED WHEN THE OUTPUT GRID
C           EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
C           THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
C        
C PROGRAM HISTORY LOG:
C   96-04-10  IREDELL
C
C USAGE:    CALL POLATES5(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
C    &                    NO,RLAT,RLON,IBO,LO,GO,IRET)
C
C   INPUT ARGUMENT LIST:
C     IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
C                IPOPT(1) IS FILTER TYPE.
C                  IPOPT(1)=0 FOR BUTTERWORTH FILTER;
C                  IPOPT(1)=1 FOR RATIONAL GAUSSIAN FILTER.
C                IPOPT(2) IS ORDER OF FILTER.
C                  AS ORDER APPROACHES INFINITY,
C                  BUTTERWORTH APPROACHES PURE SPECTRAL TRUNCATION AND
C                  RATIONAL GAUSSIAN APPROACHES PURE GAUSSIAN.)
C                IPOPT(3) IS GEOMETRY TYPE.
C                  IPOPT(3)=0 FOR SAME GEOMETRY AS INPUT GRID;
C                  IPOPT(3)=1 FOR SPHERICAL;
C                    (INPUT GRID MUST BE CYLINDRICAL)
C                  IPOPT(3)=2 FOR SAME GEOMETRY AS OUTPUT GRID.
C                    (INPUT GRID MUST BE CYLINDRICAL AND OUTPUT GRID
C                     MUST BE MERCATOR, STEREOGRAPHIC OR LAMBERT)
C                IPOPT(4) IS CHARACTERISTIC FILTER SCALE IN METERS.
C                  IPOPT(4)=0 FOR SAME SCALE AS OUTPUT GRID.
C                IPOPT(5) IS RESERVED FOR VECTOR INTERPOLATION OPTIONS.
C                IPOPT(6) IS INTERPOLATION TYPE.
C                  IPOPT(6)=0 FOR BILINEAR;
C                  IPOPT(6)=1 FOR BICUBIC.
C                IPOPT(7) IS MONOTONIC CONSTRAINT.
C                  IPOPT(7)=0 FOR NO CONSTRAINT;
C                  IPOPT(7)=1 FOR WITHIN RANGE OF SURROUNDING 4 POINTS;
C                  IPOPT(7)=2 FOR WITHIN TOTAL RANGE OF INPUT FIELD.
C     KGDSI    - INTEGER (22) INPUT GDS PARAMETERS AS DECODED BY W3FI63
C     KGDSO    - INTEGER (22) OUTPUT GDS PARAMETERS AS DECODED BY W3FI63
C                (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
C     MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
C                OR DIMENSION OF INPUT GRID FIELDS IF KM=1
C     MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
C                OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
C     KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
C     IBI      - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
C     LI       - LOGICAL (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
C     GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
C     NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
C     RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
C     RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
C
C   OUTPUT ARGUMENT LIST:
C     NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
C     RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
C     RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
C     IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
C     LO       - LOGICAL (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
C     GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
C     IRET     - INTEGER RETURN CODE
C                0    SUCCESSFUL INTERPOLATION
C                2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
C                3    UNRECOGNIZED OUTPUT GRID
C                51   INVALID GRIDS
C                52   INVALID INTERPOLATION OPTIONS
C
C SUBPROGRAMS CALLED:
C   GDSWIZ       GRID DESCRIPTION SECTION WIZARD
C   (IJKGDS)     RETURN FIELD POSITION FOR A GIVEN GRID POINT
C   POLFIXS      MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C
C$$$
CFPP$ EXPAND(IJKGDS)
      INTEGER IPOPT(20)
      INTEGER KGDSI(22),KGDSO(22)
      INTEGER IBI(KM),IBO(KM)
      LOGICAL LI(MI,KM),LO(MO,KM)
      REAL GI(MI,KM),GO(MO,KM)
      REAL RLAT(MO),RLON(MO)
      REAL XPTS(MO),YPTS(MO)
      INTEGER N11(MO),N21(MO),N12(MO),N22(MO)
      REAL W11(MO),W21(MO),W12(MO),W22(MO)
      REAL XPTI(MI),YPTI(MI),RLAI(MI),RLOI(MI)
      REAL YF(MI),GF(MI,KM)
      PARAMETER(FILL=-9999.)
      PARAMETER(RERTH=6.3712E6)
      PARAMETER(PI=3.14159265358979,DPR=180./PI)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
      IRET=0
      IF(KGDSO(1).GE.0) THEN
        CALL GDSWIZ(KGDSO, 0,MO,FILL,XPTS,YPTS,RLON,RLAT,NO,0,DUM,DUM)
        IF(NO.EQ.0) IRET=3
      ENDIF
      DO K=1,KM
        IF(IBI(K).NE.0) IRET=51
      ENDDO
      IDRTO=KGDSO(1)
      IDRTI=KGDSI(1)
      NX=KGDSI(2)
      NY=KGDSI(3)
      ISCAN=MOD(KGDSI(11)/128,2)
      JSCAN=MOD(KGDSI(11)/64,2)
      NSCAN=MOD(KGDSI(11)/32,2)
      IF(IDRTI.NE.0.AND.IDRTI.NE.1.AND.IDRTI.NE.3.AND.IDRTI.NE.4.AND.
     &   IDRTI.NE.5) IRET=51
      IF(NSCAN.EQ.1) IRET=51
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      IF(IRET.EQ.0) THEN
        CALL GDSWIZ(KGDSI, 0,MI,FILL,XPTI,YPTI,RLOI,RLAI,NI,0,DUM,DUM)
        MODEX=0
        IF(IJKGDS(-1,1,KGDSI).GT.0) MODEX=1
        MODEE=0
        IF(IPOPT(3).EQ.0) THEN
          MODEC1=0
          MODEC2=0
          MODEB1=0
          MODEB2=0
          MODEG=0
          IF(IDRTI.EQ.0.OR.IDRTI.EQ.1.OR.IDRTI.EQ.4) THEN
            RLAT1=KGDSI(4)*1.E-3
            RLON1=KGDSI(5)*1.E-3
            RLAT2=KGDSI(7)*1.E-3
            RLON2=KGDSI(8)*1.E-3
            IF(ISCAN.EQ.0) THEN
              XB=ABS(MOD(RLON2-RLON1-1+3600,360.)+1)/DPR
            ELSE
              XB=ABS(MOD(RLON1-RLON2-1+3600,360.)+1)/DPR
            ENDIF
            YB1=RLAT1/DPR
            YB2=RLAT2/DPR
          ELSE
            DX=KGDSI(8)
            DY=KGDSI(9)
            XB=(NX-1)*DX/RERTH
            YB1=-(NY-1)*DY/RERTH/2
            YB2=+(NY-1)*DY/RERTH/2
          ENDIF
        ELSEIF(IPOPT(3).EQ.1) THEN
          IF(IDRTI.EQ.0.OR.IDRTI.EQ.1.OR.IDRTI.EQ.4) THEN
            RLAT1=KGDSI(4)*1.E-3
            RLON1=KGDSI(5)*1.E-3
            RLAT2=KGDSI(7)*1.E-3
            RLON2=KGDSI(8)*1.E-3
            MODEC1=1
            MODEC2=1
            MODEB1=0
            IF(ABS(RLAT1).GE.89.9995) MODEB1=1
            MODEB2=0
            IF(ABS(RLAT2).GE.89.9995) MODEB2=1
            MODEG=1
            IF(ISCAN.EQ.0) THEN
              XB=ABS(MOD(RLON2-RLON1-1+3600,360.)+1)/DPR
            ELSE
              XB=ABS(MOD(RLON1-RLON2-1+3600,360.)+1)/DPR
            ENDIF
            YB1=RLAT1/DPR
            YB2=RLAT2/DPR
            IY=0
            DO N=1,NI
              IF(XPTI(N).EQ.1) THEN
                IY=IY+1
                YF(IY)=YPTI(NO)/DPR
              ENDIF
            ENDDO
          ELSE
            IRET=52
          ENDIF
        ELSEIF(IPOPT(3).EQ.2) THEN
          IF((IDRTI.EQ.0.OR.IDRTI.EQ.1.OR.IDRTI.EQ.4).AND.
     &       (IDRTO.EQ.1.OR.IDRTO.EQ.3.OR.IDRTO.EQ.5)) THEN
            RLAT1=KGDSI(4)*1.E-3
            RLON1=KGDSI(5)*1.E-3
            RLAT2=KGDSI(7)*1.E-3
            RLON2=KGDSI(8)*1.E-3
            MODEB1=0
            IF(ABS(RLAT1).GE.89.9995) MODEB1=1
            MODEB2=0
            IF(ABS(RLAT2).GE.89.9995) MODEB2=1
            MODEG=1
            IF(ISCAN.EQ.0) THEN
              XB=ABS(MOD(RLON2-RLON1-1+3600,360.)+1)/DPR
            ELSE
              XB=ABS(MOD(RLON1-RLON2-1+3600,360.)+1)/DPR
            ENDIF
            YB1=RLAT1/DPR
            YB2=RLAT2/DPR
            IY=0
            DO N=1,NI
              IF(XPTI(N).EQ.1) THEN
                IY=IY+1
                YF(IY)=YPTI(NO)/DPR
              ENDIF
            ENDDO
            IF(IDRTO.EQ.1) THEN
              MODEC1=0
              MODEC2=0
              DO IY=1,NY
                IF(YF(IY).GE.89.9995/DPR) THEN
                  YF(IY)=1.E20
                ELSEIF(YF(IY).LE.-89.9995/DPR) THEN
                  YF(IY)=-1.E20
                ELSE
                  YF(IY)=LOG(TAN(PI/4-YF(IY)/2))
                ENDIF
              ENDDO
            ELSEIF(IDRTO.EQ.3) THEN
              IPROJ=MOD(KGDSO(10)/128,2)
              RLATI1=KGDSO(12)*1.E-3
              RLATI2=KGDSO(13)*1.E-3
              H=(-1.)**IPROJ
              IF(H*(RLAT1-RLAT2).GT.0) THEN
                MODEC1=1
                MODEC2=0
              ELSE
                MODEC1=0
                MODEC2=1
              ENDIF
              IF(RLATI1.EQ.RLATI2) THEN
                AN=SIN(H*RLATI1/DPR)
              ELSE
                AN=LOG(COS(RLATI1/DPR)/COS(RLATI2/DPR))/
     &             LOG(TAN((H*RLATI1+90)/2/DPR)/
     &                 TAN((H*RLATI2+90)/2/DPR))
              ENDIF
              DO IY=1,NY
                IF(YF(IY).GE.89.9995/DPR) THEN
                  YF(IY)=1.E20*IPROJ
                ELSEIF(YF(IY).LE.-89.9995/DPR) THEN
                  YF(IY)=1.E20*(1-IPROJ)
                ELSE
                  YF(IY)=TAN(PI/4-H*YF(IY)/2)**AN
                ENDIF
              ENDDO
            ELSEIF(IDRTO.EQ.5) THEN
              IPROJ=MOD(KGDSO(10)/128,2)
              H=(-1.)**IPROJ
              IF(H*(RLAT1-RLAT2).GT.0) THEN
                MODEC1=1
                MODEC2=0
              ELSE
                MODEC1=0
                MODEC2=1
              ENDIF
              DO IY=1,NY
                IF(YF(IY).GE.89.9995/DPR) THEN
                  YF(IY)=1.E20*IPROJ
                ELSEIF(YF(IY).LE.-89.9995/DPR) THEN
                  YF(IY)=1.E20*(1-IPROJ)
                ELSE
                  YF(IY)=TAN(PI/4-H*YF(IY)/2)
                ENDIF
              ENDDO
            ENDIF
          ELSE
            IRET=52
          ENDIF
        ENDIF
        MODES=1
        MODEF=IPOPT(1)
        MODEP=0
        MODER=0
        NDEG=IPOPT(2)
        NORDH=0
        NYA=NY
        HWL=IPOPT(4)/RERTH
        IF(IPOPT(4).EQ.0) THEN
          IF(IDRTI.EQ.0.OR.IDRTI.EQ.1.OR.IDRTI.EQ.4) THEN
            HWL=RERTH*ABS(KGDSO(4)-KGDSO(7))*1.E-3/DPR/KGDSO(3)
          ELSE
            HWL=KGDSO(9)
          ENDIF
        ENDIF
        CALL QFPAR(MODEX,MODEE,MODEC1,MODEC2,MODEB1,MODEB2,MODEG,MODES,
     &             MODEF,MODEP,MODER,NX,NY,NY,NDEG,NORDH,XB,YB1,YB2,
     &             NQ,NQ3,NNNNQ,N,NW,NWAB,NJAB,NWORK,NB,NV)
        NTM1=MAX(NWORK,NB,NV)
        GF=GI
        CALL IPS51(MODEX,MODEE,MODEC1,MODEC2,MODEB1,MODEB2,MODEG,MODES,
     &             MODEF,MODEP,MODER,NX,NY,NY,NDEG,NORDH,XB,YB1,YB2,
     &             NQ,NQ3,NNNNQ,N,NW,NWAB,NJAB,NWORK,NB,NV,
     &             NTM1,YF,HWL,MI,KM,GF)
C INTERPOLATE GF TO GO HERE
      ENDIF
      END
      SUBROUTINE IPS51(MODEX,MODEE,MODEC1,MODEC2,
     &                 MODEB1,MODEB2,MODEG,MODES,
     &                 MODEF,MODEP,MODER,NX,NY,NYA,NDEG,NORDH,
     &                 XB,YB1,YB2,
     &                 NQ,NQ3,NNNNQ,N,NW,NWAB,NJAB,NWORK,NB,NV,
     &                 NTM1,Y,HWL,JM,KM,G)
      REAL Y(NY),G(JM,KM)
      REAL RTAU(NY),RTAUI(NY),FF(NNNNQ),W(NW),WAB(NWAB),TM1(NTM1)
      INTEGER NQC(NQ),JUMBLE(N),JAB(NJAB)
      CALL INFPAR(MODEX,MODEE,MODEC1,MODEC2,MODEB1,MODEB2,MODEG,MODES,
     &            MODEF,MODEP,MODER,NX,NY,NYA,NDEG,NORDH,XB,YB1,YB2,Y,
     &            HWL,RTAU,RTAUI,FF,NQC,NQ,NXU,NH,W,JUMBLE,WAB,JAB,TM1)
      DO K=1,KM
        CALL BRF(G,RTAU,RTAUI,FF,
     &           NX,NXU,NH,NY,NQC,NQ,WAB,JAB,W,JUMBLE,TM1)
      ENDDO
      END
