      SUBROUTINE COVISO(INF,BETA,NLON,NLAT,LATF,LONF,LATH,LONH,COF)
C
C     * MAY 14/02 - B. MIVILLE - REVISED FOR NEW DATA DESCRIPTION FORMAT
C     * APR 22/02 - B. MIVILLE - NEW VALUE OF PI
C     * JUN 08/01 - B. MIVILLE
C
C     * VSIO COLOCATION
C
C     * CALCULATE THE COLOCATION VALUE OF THE V ISO VELOCITY
C
C     * INPUT PARAMETERS...
C
C     * INF  =  INPUT FIELD
C     * BETA = HEAVISIDE FUNCTION
C     * NLON = NUMBER OF LONGITUDES INCLUDING CYCLIC
C     * NLAT = NUMBER OF LATITUDES
C     * LATF = TRACER GRID POINT LATITUDES
C     * LONF = TRACER GRID POINT LONGITUDES
C     * LATH = HALF CELL GRID POINT LATITUDES (V)
C     * LONH = HALF CELL GRID POINT LONGITUDES (V)
C
C     * OUTPUT PARAMETERS...
C
C     * COF = OUTPUT COLOCATED FIELD
C
C--------------------------------------------------------------------------
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      INTEGER NLON,NLAT,NLATM
      REAL INF(NLON,NLAT),BETA(NLON,NLAT)
      REAL LATF(NLON,NLAT),LATH(NLON,NLAT)
      REAL LONF(NLON,NLAT),LONH(NLON,NLAT)
      REAL COF(NLON,NLAT),DAP1,DAP2,DAM1,DAM2,PI,RAD
      DATA PI/3.1415926535898D0/
C-----------------------------------------------------------------------
C
C     *  
C     *             I     
C     *  
C     *     +-------VI------+---- LATH(I,J)
C     *     +       +       +  
C     *     +  DAP1 +  DAP2 +  
C     *     +       +       +
C     *     +-------T-------+----- J --- LATF(I,J)
C     *     +       +       +
C     *     +  DAM1 +  DAM2 +  
C     *     +       +       +
C     *     +-------VI------+----- LATH(I,J-1)
C     *
C     *
C     *     DAP1 +DAP2 = FRACTIONAL AREA OF GRID CELL T FOR VI(J)
C     *     DAM1 +DAM2 = FRACTIONAL AREA OF GRID CELL T FOR VI(J-1)
C     *
C     *     WE DO AN INTERPOLATION USING
C     *
C     *     COVI(I,J)=(VI(I,J)* (DAM1+DAM2) + VI(I,J-1)* (DAP1+DAP2))
C     *               -----------------------------------------------
C     *                         (DAM1 + DAM2 + DAP1 + DAP2)
C     *
C
C     * DO COLOCATION
C
      RAD = PI / 180.0D0
      NLATM=NLAT-1
C
      DO 510 I=1,NLON
         DO 500 J=2,NLATM
            ANGP=RAD*(LATH(I,J)+LATF(I,J))/2.E0
            ANGM=RAD*(LATH(I,J-1)+LATF(I,J))/2.E0
            DLON1=(LONF(I,J)-LONH(I-1,J))
            DLON2=(LONH(I,J)-LONF(I,J))
            DAP1=COS(ANGP)*(LATH(I,J)-LATF(I,J))*DLON1
            DAP2=COS(ANGP)*(LATH(I,J)-LATF(I,J))*DLON2
            DAM1=COS(ANGM)*(LATF(I,J)-LATH(I,J-1))*DLON1
            DAM2=COS(ANGM)*(LATF(I,J)-LATH(I,J-1))*DLON2
            COF(I,J)=BETA(I,J)*((INF(I,J)*(DAM1+DAM2))+
     1               (INF(I,J-1)*(DAP1+DAP2)))/
     2               (DAM1+DAM2+DAP1+DAP2)
 500     CONTINUE
C        * SOUTH POLE
         COF(I,1)=BETA(I,1)*INF(I,1)
C        * NORTH POLE
         COF(I,NLAT)=BETA(I,NLATM)*INF(I,NLATM)
 510  CONTINUE
C
      RETURN
      END
