      REAL FUNCTION UNODEN(T,S,P0)
C
C     * JUN 07/2001 - B. MIVILLE
C
C     * CALCULATES THE OCEAN DENSITY USING THE UNESCO EQUATION OF STATE
C
C     * INPUT PARAMETERS...
C
C     * T  = TEMPERATURE (DEGREE CELSIUS)
C     * S  = SALINITY    
C     * P0 = PRESSURE AT CURRENT DEPTH (DECIBARS)
C
C     * OUTPUT PARAMETERS...
C
C     * RHO = OCEAN DENSITY (KG/M**3)
C
C-----------------------------------------------------------------------
C
C     * CHECK VALUE: RHO=999.96675 AT
C     *          T = 5 DEG C , S = 0, P0 = 0 DECIBARS.
C     * CHECK VALUE: RHO=1027.67547 AT
C     *          T = 5 DEG C , S = 35, P0 = 0 DECIBARS.
C     * CHECK VALUE: RHO=1062.53817 AT
C     *          T = 25 DEG C , S = 35, P0 = 10000 DECIBARS.
C     * (THESE ARE FROM GILL'S BOOK)
C
C-----------------------------------------------------------------------
C
C     * REFERENCES:
C
C     * MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629
C     *   (THERE WERE ERRORS IN THIS PAPER WHICH WERE FIXED IN THE
C     *   FOLLOWING UNESCO REPORT)
C     * UNESCO 1981. UNESCO TECHNICAL PAPERS IN MAR. SCI., NO. 38
C     * GILL, A., ATMOSPHERE-OCEAN DYNAMICS: INTERNATIONAL GEOPHYSICAL
C     *   SERIES NO. 30. ACADEMIC PRESS, LONDON, 1982, pp 599-600.
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      REAL T,S,P,P0
      REAL EXPO,RW,RSTO,XKW,XKSTO,XKSTP,RHO
C
C-----------------------------------------------------------------------
C
      EXPO = 1.5E0
C
C     * CONVERT PRESSURE FROM DECIBARS TO BARS
C     
      P = P0  / 10.E0
C
C     * CALCULATE DENSITY
C
C     * SAME CODE AS IN THE NCOM (VER 1.3) DENSITY.F PROGRAM
C     
C
      RW =     9.99842594E2 + 6.793952E-2*T - 9.095290E-3*T**2
     1        + 1.001685E-4*T**3 - 1.120083E-6*T**4 + 6.536332E-9*T**5
C
      RSTO =   RW + (8.24493E-1 - 4.0899E-3*T + 7.6438E-5*T**2
     1        - 8.2467E-7*T**3 + 5.3875E-9*T**4) * S
     2       + (-5.72466E-3 + 1.0227E-4*T - 1.6546E-6*T**2) * S**EXPO
     3       + 4.8314E-4 * S**2
C
      XKW =     1.965221E4 + 1.484206E2*T - 2.327105E0*T**2 +
     1         1.360477E-2*T**3 - 5.155288E-5*T**4
C
      XKSTO =   XKW + (5.46746E1 - 6.03459E-1*T + 1.09987E-2*T**2
     1        - 6.1670E-5*T**3) * S
     2       + (7.944E-2 + 1.6483E-2*T - 5.3009E-4*T**2) * S**EXPO
C
      XKSTP =   XKSTO + (3.239908E0 + 1.43713E-3*T + 1.16092E-4*T**2
     1        - 5.77905E-7*T**3) * P
     2       + (2.2838E-3 - 1.0981E-5*T - 1.6078E-6*T**2) * P * S
     3       + 1.91075E-4 * P * S**EXPO
     4       + (8.50935E-5 - 6.12293E-6*T + 5.2787E-8*T**2) * P**2
     5       + (-9.9348E-7 + 2.0816E-8*T + 9.1697E-10*T**2) * P**2 * S
C
      RHO =    RSTO / (1.0E0 - P/XKSTP)
C
      UNODEN=RHO
C
      RETURN
      END
