SUBROUTINE RHO(REFPRS,PRSDB,TEMP,SALTY,RHOSTP)
C
C     TITLE:
C     *****
C
C       RHO  -- CALCULATE DENSITY USING INTERNATIONAL EQUATION OF
C                 STATE
C
C                 FROM TEXT FURNISHED BY J. GIESKES
C
C     PARAMETERS:
C     **********
C
C       REFPRS -> REFERENCE PRESSURE IN DECIBARS
C       PRSDB  -> INSITU PRESSURE IN DECIBARS
C       TEMP   -> TEMPERATURE IN CELSIUS DEGREES
C       SALTY  -> SALINITY PSS 78
C       RHOSTP <- DENSITY IN Kg/M**3 : RHO(S,T,P)
C
        REAL REFPRS,PRSDB,TEMP,SALTY,RHOSTP
C
C     VARIABLES:
C     *********
C
        REAL BARS,TERMA,TERMB,KST0,KW,RHOW,KSTP
C
C     CODE:
C     ****
C
      BARS = REFPRS*.10
C       /* REFERENCE PRESSURE IN BARS */
C
C     /* CALCULATE DENSITY AT GIVEN SALINITY, TEMP AND PRESSURE=0.0 */
C
      RHOW = -.00909529+TEMP*(1.001685E-4+TEMP*(-1.120083E-6+
     1       TEMP*6.536332E-9))
      RHOW = 999.842594+TEMP*(.06793952+TEMP*RHOW)
C       /* RHOW = DENSITY OF PURE WATER Kg/M**3 */
C
      KW  = TEMP*(-.0040899+TEMP*(7.6438E-5+TEMP*(-8.2467E-7
     1       +TEMP*5.3875E-9)))
C           /* PURE WATER SECANT BULK MODULUS */
C
      KW  = (KW+.824493)*SALTY
C
      KST0 = (-.00572466+TEMP*(1.0227E-4+TEMP*(-1.6546E-6)))
     1       * ABS(SALTY)**1.5
C       /* K(S,T,0) */
C
      RHOSTP = RHOW + KW + KST0 + 4.8314E-4*SALTY*SALTY
C       /* RHO(S,T,P)  Kg/M**3 */
C
C     SELECT  /* IF REFERENCE PRESSURE=0.0, DONE */
C       (REFERENCE PRESSURE = 0.0):
          IF(REFPRS.EQ.0.0) GO TO 999
C       (OTHERWISE):  /* CALCULATE PRESSURE EFFECT */
C
C         /* RHO(S,T,0)/(1.0-P/K(S,T,P)) */
C
         KW = TEMP*(148.4206+TEMP*(-2.327105+TEMP
     1        *(.01360477+TEMP*(-5.155288E-5))))
         KW = KW + 19652.21
C
         KST0 = (54.6746+TEMP*(-.603459+TEMP*(.0109987+TEMP
     1          *(-6.167E-5))))*SALTY
C
         KST0 = KST0 + KW+(.07944+TEMP*(.016483+TEMP*(-5.3009E-4)))
     1          *ABS(SALTY)**1.5
C
C        /* CALCULATE PRESSURE TERMS */
C
         TERMA = 3.239908+TEMP*(.00143713+TEMP*(1.16092E-4+TEMP
     1           *(-5.77905E-7)))
         TERMA = TERMA + (.0022838+TEMP*(-1.0981E-5+TEMP
     1                 *(-1.6078E-6)))*SALTY
         TERMA = TERMA + 1.91075E-4*ABS(SALTY)**1.5
C
         TERMB = 8.50935E-5+TEMP*(-6.12293E-6+TEMP*5.2787E-8)
         TERMB = TERMB + (-9.9348E-7+TEMP*(2.0816E-8+TEMP
     1                 * 9.1697E-10))*SALTY
C
         KSTP  = KST0+BARS*(TERMA + BARS*TERMB)
C          /* SECANT BULK MODULUS K(S,T,P) */
C
         RHOSTP = RHOSTP/(1.0-BARS/KSTP)
C
C     END SELECT
C
  999 CONTINUE
      RETURN
C
C     END RHO
C
      END