      SUBROUTINE GGPNL3(U,V,UG,VG,ILG,ILG1,ILAT,DLON,DLAT,
     1                  SLON,SLAT,USP,VSP,UNP,VNP,IFLD)
C
C     * MAR 25/13 - S.KHARIN
C
C     * BI-LINEAR INTERPOLATION FROM A SCALAR FIELD UG(ILG1,ILAT) OR
C     * A VECTOR FIELD (UG(ILG1,ILAT), VG(ILG1,ILAT)) TO A SINGLE
C     * POINT U OR (U,V) WITH LATITUDE DLAT AND LONGITUDE DLON.
C     * DLON IS MEASURED IN DEGREES
C     * DLAT IS MEASURED IN DEGREES FROM THE EQUATOR.
C     * SLAT AND SLON ARE INPUT GRID LATITUDES AND LONGITUDES.
C     * USP,VSP,UNP,VNP ARE SCALAR OR VECTOR ESTIMATES AT THE POLES.
C     * IFLD=0 USE SCALAR FIELD UG.
C     * IFLD=1 USE A VECTOR COMPONENT UG.
C     * IFLD=2 USE BOTH VECTOR COMPONENTS (UG,VG).
C     * RETURN SPVAL IF ANY OF NEIGHBORS HAS SPVAL.
C
      IMPLICIT NONE
      REAL U,V,DLAT,DLON,USP,VSP,UNP,VNP
      INTEGER ILG,ILG1,ILAT,IFLD,I,J
      REAL UG(ILG1,ILAT),VG(ILG1,ILAT),SLON(ILG1),SLAT(ILAT)
      REAL G(2,2),W(2,2),X,Y
      INTEGER IX(2),IY(2)
      REAL SPVAL,SPVALT
      DATA SPVAL/1.E+38/,SPVALT/1.E+32/
C-----------------------------------------------------------------------
      IF(IFLD.LT.0.OR.IFLD.GT.2) CALL              XIT('GGPNL3',-1)
      IF(ILG.LT.2.OR.ILAT.LT.2) CALL               XIT('GGPNL3',-2)
C
C     * FIND 2 CLOSEST INPUT LONGITUDES (TAKE INTO ACCOUNT PERIODICITY)
C
      IF(DLON.LT.SLON(1))THEN
        IX(1)=ILG
        IX(2)=1
        X=(DLON+360.E0-SLON(IX(1)))/(SLON(IX(2))+360.E0-SLON(IX(1)))
      ELSE IF(DLON.GE.SLON(ILG))THEN
        IX(1)=ILG
        IX(2)=1
        X=(DLON-SLON(IX(1)))/(SLON(IX(2))+360.E0-SLON(IX(1)))
      ELSE
        DO I=1,ILG1
          IF(DLON.LT.SLON(I)) GO TO 110
        ENDDO
 110    CONTINUE
        IX(1)=I-1
        IX(2)=I
        X=(DLON-SLON(IX(1)))/(SLON(IX(2))-SLON(IX(1)))
      ENDIF
C
C     * FIND 2 CLOSEST INPUT LATITUDES (RESTRICT TO 0...ILAT+1)
C
      IF(DLAT.LT.SLAT(1)) THEN
        IY(1)=0
        IY(2)=1
        Y=(DLAT+90.E0)/(SLAT(IY(2))+90.E0)
      ELSE IF(DLAT.GE.SLAT(ILAT)) THEN
        IY(1)=ILAT
        IY(2)=ILAT+1
        IF(SLAT(IY(1)).EQ.90.E0)THEN
          Y=0.E0
        ELSE
          Y=(DLAT-SLAT(IY(1)))/(90.E0-SLAT(IY(1)))
        ENDIF
      ELSE
        DO J=1,ILAT
          IF(DLAT.LT.SLAT(J)) GO TO 120
        ENDDO
 120    CONTINUE
        IY(1)=J-1
        IY(2)=J
        Y=(DLAT-SLAT(IY(1)))/(SLAT(IY(2))-SLAT(IY(1)))
      ENDIF
C
C     * SELECT 2X2 ARRAY FROM UG FOR BI-LINEAR INTERPOLATION
C     * (USE SOUTH OR NORTH POLE APPROXIMATIONS FOR J=0,ILAT+1)
C
      DO J=1,2
        DO I=1,2
          G(I,J)=SPVAL
          IF(IY(J).EQ.0)THEN
C
C           * SOUTH POLE
C
            IF(IFLD.EQ.0)THEN
              G(I,J)=USP
            ELSE
              IF(USP.NE.SPVAL)
     1             G(I,J)=USP*COSD(DLON)+VSP*SIND(DLON)
            ENDIF
          ELSE IF(IY(J).EQ.ILAT+1)THEN
C
C           * NORTH POLE
C
            IF(IFLD.EQ.0)THEN
              G(I,J)=UNP
            ELSE
              IF(UNP.NE.SPVAL)
     1             G(I,J)=UNP*COSD(DLON)+VNP*SIND(DLON)
            ENDIF
          ELSE
C
C           * AWAY FROM THE POLES
C
            G(I,J)=UG(IX(I),IY(J))
          ENDIF
        ENDDO
      ENDDO
C
C     * WEIGHT MATRIX FOR BI-LINEAR INTERPOLATION
C
      W(1,1)=(1.E0-X)*(1.E0-Y)
      W(1,2)=(1.E0-X)*Y
      W(2,1)=X*(1.E0-Y)
      W(2,2)=X*Y
C
C     * BI-LINEAR INTERPOLATION OF U
C
      U=0.
      DO J=1,2
        DO I=1,2
          IF(ABS(G(I,J)-SPVAL).LT.SPVALT.AND.W(I,J).NE.0.E0)THEN
            U=SPVAL
            RETURN
          ENDIF
          U=U+W(I,J)*G(I,J)
        ENDDO
      ENDDO
      IF(IFLD.LE.1)RETURN
C
C     * SELECT 2X2 ARRAY FROM VG FOR BI-LINEAR INTERPOLATION
C     * (USE SOUTH OR NORTH POLE APPROXIMATIONS FOR J=0, ILAT+1)
C
      DO J=1,2
        DO I=1,2
          IF(IY(J).EQ.0)THEN
C
C           * SOUTH POLE
C
            IF(USP.NE.SPVAL)
     1           G(I,J)=USP*SIND(DLON)-VSP*COSD(DLON)
          ELSE IF(IY(J).EQ.ILAT+1)THEN
C
C           * NORTH POLE
C
            IF(UNP.NE.SPVAL)
     1           G(I,J)=VNP*COSD(DLON)-UNP*SIND(DLON)
          ELSE
C
C           * AWAY FROM THE POLES
C
            G(I,J)=VG(IX(I),IY(J))
          ENDIF
        ENDDO
      ENDDO
C
C     * BI-LINEAR INTERPOLATION OF V
C
      V=0.
      DO J=1,2
        DO I=1,2
          IF(ABS(G(I,J)-SPVAL).LT.SPVALT.AND.W(I,J).NE.0.E0)THEN
            V=SPVAL
            RETURN
          ENDIF
          V=V+W(I,J)*G(I,J)
        ENDDO
      ENDDO
      RETURN
      END
