      SUBROUTINE BVINTRP(FIELDA,NLONA,NLATA,ELONA,ELATA,INDI,INDJ,
     1                   NLONO,NLATO,LONO,LATO,DAO,BETAO,BETAA,
     2                   DAA,FIELDI,EFIELDA)
C
C     * JAN 20/04 - B. MIVILLE - FIXED EFIELDA NOT BEING INITIALIZED
C                                AT THE CORNERS OF THE GRID
C     * APR 04/02 - B. MIVILLE - REVISED FOR NEW DATA DESCRIPTION FORMAT
C     * JUN 08/01 - B. MIVILLE
C
C     * DOES A BIVARIATE INTERPOLATION TO FIND A FIRST APPROXIMATION FOR
C     * THE OCEAN POINTS
C
C     * INPUT PARAMETERS...
C
C     * FIELDA= INPUT ATMOSPHERE FIELD
C     * NLONA = NUMBER OF ATMOSPHERE LONGITUDES INCLUDING CYCLIC
C     * NLATA = NUMBER OF ATMOSPHERE LATITUDES 
C     * ELONA = ATMOSPHERE LONGITUDES (PLUS ONE EXTRA TO THE WEST OF THE 
C     *         FIRST LONGITUDE
C     * ELATA = ATMOSPHERE LATITUDES (PLUS ONE EXTRA TO THE NORTH AND 
C     *         SOUTH)
C     * INDI  = INDEX THE ATMOSPHERE II POINT TO OCEAN I POINT FOR 
C     *         INTERPOLATION
C     * INDJ  = INDEX THE ATMOSPHERE JJ POINT TO OCEAN J POINT FOR 
C     *         INTERPOLATION
C     * NLONO = NUMBER OF OCEAN LONGITUDES INCLUDING CYCLIC
C     * NLATO = NUMBER OF OCEAN LATITUDES
C     * LONO  = OCEAN LONGITUDES
C     * LATO  = OCEAN LATITUDES
C     * DAO   = AREA OF OCEAN GRID CELL
C     * BETAO = HEAVISIDE FUNCTION FOR OCEAN GRID CELL
C     * BETAA = HEAVISIDE FUNCTION FOR ATMOSPHERE GRID CELL
C     * DAA   = AREA OF ATMOSPHERE GRID CELL
C
C     * OUTPUT PARAMETERS...
C
C     * FIELDI  = INTERPOLATED OCEAN FIELD
C     * EFIELDA = EXPANDED ATMOSPHERE FIELD
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      INTEGER NLONA,NLATA,NLONO,NLATO
      INTEGER INDI(NLONO,NLATO),INDJ(NLONO,NLATO)
      REAL FIELDA(NLONA,NLATA),EFIELDA(0:NLONA,0:NLATA+1)
      REAL FIELDI(NLONO,NLATO)
      REAL DAO(NLONO,NLATO),BETAO(NLONO,NLATO)
      REAL LONO(NLONO,NLATO),LATO(NLONO,NLATO)
      REAL ELONA(0:NLONA,0:NLATA+1),ELATA(0:NLONA,0:NLATA+1)
      REAL BETAA(0:NLONA,0:NLATA+1),DAA(0:NLONA,0:NLATA+1)
      REAL XF,SUMA,SUMF
      REAL A00,A01,A10,A11
      REAL INTER,P,Q
C-----------------------------------------------------------------------
C
C     * THE WEST BOUNDARY (I=0) FOR THE ATMOSPHERE
C
      DO 222 J=1,NLATA
         DO 220 I=1,NLONA
            EFIELDA(I,J)=FIELDA(I,J)
 220     CONTINUE
         EFIELDA(0,J)=FIELDA(NLONA-1,J)
 222  CONTINUE
C
C     * THE NORTH AND SOUTH BOUNDARY FOR THE ATMOSPHERE
C
      DO 224 I=1,NLONA
         EFIELDA(I,0)=FIELDA(I,1)
         EFIELDA(I,NLATA+1)=FIELDA(I,NLATA)
 224  CONTINUE
C
C     * INITIALIZE CORNERS
C
      EFIELDA(0,0)=FIELDA(NLONA-1,1)
      EFIELDA(0,NLATA+1)=FIELDA(NLONA-1,NLATA)
C     
C     * FILL IN ATMOSPHERE DATA OVER LAND DATA WITH AVERAGE WEIGHTED OF 
C     * 8 ATMOSPHERE POINTS AROUND IT
C     
      DO 232 J=1,NLATA
         DO 230 I=1,NLONA-1
            SUMF=0.E0
            SUMA=0.E0
            DO 228 JJ=J-1,J+1
               DO 226 II=I-1,I+1
                  SUMA=SUMA+BETAA(II,JJ)*DAA(II,JJ)
                  SUMF=SUMF+BETAA(II,JJ)*EFIELDA(II,JJ)*DAA(II,JJ)
 226           CONTINUE
 228        CONTINUE
            IF(SUMA.EQ.0.E0) THEN
               XF=0.E0
            ELSE
               XF=SUMF/SUMA
            ENDIF
            EFIELDA(I,J)=(1.0E0-BETAA(I,J))*XF+BETAA(I,J)*EFIELDA(I,J)
 230     CONTINUE
         EFIELDA(0,J)=EFIELDA(NLONA-1,J)
         EFIELDA(NLONA,J)=EFIELDA(1,J)
 232  CONTINUE
C     
C     * REDO THE NORTH AND SOUTH BOUNDARY FOR THE ATMOSPHERE
C     
      DO 234 I=0,NLONA
         EFIELDA(I,0)=EFIELDA(I,1)
         EFIELDA(I,NLATA+1)=EFIELDA(I,NLATA)
 234  CONTINUE
C     
C     * INTERPOLATION
C     *
C     *   
C     *                  
C     *  
C     *     A-------+-------A---- LATA(II,JJ+1)
C     *     +       +       +  
C     *     +       +       +  
C     *     +       +       +
C     *     +-------+-------+---- 
C     *     +       +       +
C     *     +       +   F   +---- LATO(I,J)  
C     *     +       +       +
C     *     A-------+-------A---- LATA(II,JJ)
C     *     |           |   |
C     *     |           |   |
C     *  LONA(II,JJ)    |   LONA(II+1,JJ)     
C     *                 |
C     *                LONO(I,J)
C     *
C     *
C     *  WE USE BIVARIATE INTERPOLATION
C     *
C     *  IF WE HAVE 4 POINTS (ALL 4 ATMOSPHERE BETAS EQUAL TO ONE)
C     *  THEN WE USE:
C     *  
C     *  FIELDI(I,J)=(1-P)(1-Q)*FIELDA(II,JJ)+P(1-Q)*FIELDA(II+1,JJ)+
C     *               Q*(1-P)*FIELDA(II,JJ+1)+PQ*FIELDA(II+1,JJ+1)
C     *
C     *  WHERE
C     *        P=(LONO(I,J)-LONA(II,JJ))/(LONA(II+1,JJ)-LONA(II,JJ))
C     *        Q=(LATO(I,J)-LATA(II,JJ))/(LATA(II,JJ+1)-LATA(II,JJ))
C     *
C     * WE ASSUME THE FOLLOWING:
C     *
C     *        INTERPOLATION LONGITUDE AND LATITUDE FALL INSIDE 
C     *        THE RANGE OF EXPANDED ATMOSPHERE POINTS. THAT IS
C     *        THE ORIGINAL ATMOSPHERE RANGE PLUS ONE MORE POINT TO 
C     *        THE EAST,WEST,NORTH AND SOUTH.
C     *
C
C     * P AND Q ARE THE RELATIVE DISTANCE BETWEEN OCEAN POINT AND 
C     * ATMOSPHERE POINT FROM THE LOWER LEFT CORNER
C
      DO 238 J=1,NLATO
         DO 236 I=1,NLONO-1
            II=INDI(I,J)
            JJ=INDJ(I,J)
C           * NORTH AND SOUTH POLE WE MAKE Q=0.0
            IF((JJ.EQ.0).OR.(JJ.EQ.NLATA))THEN
               Q=0.0E0
            ELSE
               Q=(LATO(I,J)-ELATA(II,JJ))/(ELATA(II,JJ+1)-ELATA(II,JJ))
            ENDIF
C
            P=(LONO(I,J)-ELONA(II,JJ))/(ELONA(II+1,JJ)-ELONA(II,JJ))
C           * ATMOSPHERE FIELD
            A00=EFIELDA(II,JJ)
            A01=EFIELDA(II,JJ+1)
            A10=EFIELDA(II+1,JJ)
            A11=EFIELDA(II+1,JJ+1)
C     
C           * FORMULA FOR INTERPOLATION
C
            INTER=(((1.0E0-P)*(1.0E0-Q)*A00)+(P*(1.0E0-Q)*A10)+
     1           (Q*(1.0E0-P)*A01)+(P*Q*A11))
C     
            FIELDI(I,J)=BETAO(I,J)*INTER
C     
 236     CONTINUE
         FIELDI(NLONO,J)=FIELDI(1,J)   
 238  CONTINUE
C     
      RETURN
      END
