       SUBROUTINE COUV(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     * CALCULATE THE U,V AT THE TRACER LOCATION
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     * LATH = HALF CELL GRID POINT LATITUDES (U,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),LATS
      REAL LONF(NLON,NLAT),LONH(NLON,NLAT)
      REAL COF(NLON,NLAT),PI,RAD,DAP1,DAP2,DAM1,DAM2
      DATA PI/3.1415926535898D0/
C-----------------------------------------------------------------------
C
C     *
C     * CALCULATE THE COLOCATED U OR V VALUE AT THE TRACER POINTS
C     *  
C     *             I     
C     *  
C     *     U-------+-------U---- LATH(I,J)
C     *     +       +       +  
C     *     +  DAP1 +  DAP2 +  
C     *     +       +       +
C     *     +-------T-------+----- J --- LATF(I,J)
C     *     +       +       +
C     *     +  DAM1 +  DAM2  +  
C     *     +       +       +
C     *     U-------+-------U----- LATH(I,J-1)
C     *
C     *
C     *     DAP = FRACTIONAL AREA OF GRID CELL T FOR U(J)
C     *     DAM = FRACTIONAL AREA OF GRID CELL T FOR U(J-1)
C     *
C     *     E.G.:
C     *     DAP1 = DX1 * DY1
C     *     DX1  = R*(LONF(I,J)-LONH(I-1,J))
C     *     DY1  = R*COS((LATF(I,J)+LATH(I,J))/2)*(LATH(I,J)-LATF(I,J))
C     *
C     *     WE DO AN INTERPOLATION USING
C     *
C     *     COU(I,J)=(U(I,J)*DAM1+U(I-1,J)*DAM2+U(I,J-1)*DAP1+U(I-1,J-1)*DAP2)
C     *              ---------------------------------------------------------
C     *                              (DAP1+DAP2+DAM1+DAM2)
C     *
C
C     * DO COLOCATION
C
      RAD = PI / 180.0D0
      NLATM=NLAT-1      
C
      DO 310 I=2,NLON
C
C        * NO NEED TO USE THE EARTH'S RADIUS SINCE THEY CANCEL OUT, SAME THING WITH 
C        * THE TRANSFORMATION IN RADIAN OF THE LATITUDE SUBSTRACTION.
C        * SIMILARLY WITH THE LONGITUDE SUBSTRACTION.
C
         DO 300 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+INF(I-1,J)*DAM2)+
     1               (INF(I,J-1)*DAP1+INF(I-1,J-1)*DAP2))/
     2               ((DAP1+DAM1+DAP2+DAM2))
 300     CONTINUE
C        * SOUTH POLE
         LATS=-LATH(I,NLAT)
         ANGM=RAD*(LATS+LATF(I,1))/2.E0
         DLON1=(LONF(I,1)-LONH(I-1,1))
         DLON2=(LONH(I,1)-LONF(I,1))
         DAM1=COS(ANGM)*(LATF(I,1)-LATS)*DLON1
         DAM2=COS(ANGM)*(LATF(I,1)-LATS)*DLON2
         COF(I,1)=BETA(I,1)*(INF(I,1)*DAM1+INF(I-1,1)*DAM2)/(DAM1+DAM2)
C        * NORTH POLE
C        * THE U AND V VALUE AT THE NORTH POLE DO NOT MAKE SENSE
C        * SO WE TAKE THE VALUE AT NLAT-1 FOR THE COLOCATION
         LATS=-LATH(I,NLAT)
         ANGM=RAD*(LATF(I,NLAT)+LATH(I,NLATM))/2.E0
         DLON1=(LONF(I,NLAT)-LONH(I-1,NLATM))
         DLON2=(LONH(I,NLATM)-LONF(I,NLAT))
         DAM1=COS(ANGM)*(LATF(I,NLAT)-LATH(I,NLATM))*DLON1
         DAM2=COS(ANGM)*(LATF(I,NLAT)-LATH(I,NLATM))*DLON2
         COF(I,NLAT)=BETA(I,NLAT)*(INF(I,NLATM)*DAM1+
     1               INF(I-1,NLATM)*DAM2)/(DAM1+DAM2)
 310  CONTINUE
C     
C     * DO REPEATED LONGITUDE
C
      I=1
      DO 340 J=1,NLAT
         COF(1,J)=COF(NLON,J)
 340  CONTINUE
C     
      RETURN
      END
