      SUBROUTINE COUISO(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     * USIO COLOCATION
C
C     * CALCULATE THE COLOCATION VALUE OF THE U ISO VELOCITY
C
C     * INPUT PARAMETERS...
C
C     * INF  =  INPUT FIELD
C     * BETA = HEAVISIDE FUNCTION
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
      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),PI,RAD,DAP1,DAP2,DAM1,DAM2
      DATA PI/3.1415926535898D0/
C-----------------------------------------------------------------------
C
C     *             I     
C     *  
C     *     +-------+-------+---- LATH(I,J)
C     *     +       +       +  
C     *     +  DAP1 +  DAP2 +  
C     *     +       +       +
C     *     UI------T------UI--- J --- LATF(I,J)
C     *     +       +       +
C     *     +  DAM1 +  DAM2 +  
C     *     +       +       +
C     *     +-------+-------+----- LATH(I,J-1)
C     *
C     *
C     *     DAP + DAM = FRACTIONAL AREA OF GRID CELL T FOR UI
C     *
C     *     WE DO AN INTERPOLATION USING
C     *
C     *     COUI(I,J)=(UI(I,J)*(DAP1 + DAM1) + UI(I-1,J)* (DAP2 + DAM2) 
C     *               ------------------------------------------------
C     *                        (DAP1 + DAM1 + DAP2 + DAM2)
C     *
C
C     * DO COLOCATION
C      
      RAD = PI / 180.0D0
C
      DO 410 J=1,NLAT
         DO 400 I=2,NLON
            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)*(DAP1+DAM1)+
     1               INF(I-1,J)*(DAP2+DAM2))/
     2               (DAP1+DAM1+DAP2+DAM2)
 400     CONTINUE
         COF(1,J)=COF(NLON,J)
 410  CONTINUE
C
      RETURN
      END
