
C EOS80 FCN ******* APRIL 6 1981 *****
      FUNCTION rrho(P,T,S)
c      FUNCTION EOS80(P,T,S)
C ************************************
C EQUATION OF STATE FOR SEAWATER PROPOSED BY JPOTS 1980
C REFERENCES
C MILLERO ET AL 1980, DEEP-SEA RES.,27A,255-264
C JPOTS NINTH REPORT 1978,TENTH REPORT 1980
C UNITS:
CC      PRESSURE        P        DECIBARS
C       TEMPERATURE     T        DEG CELSIUS (IPTS-68)
C       SALINITY        S        NSU (IPSS-78)
CC      DENSITY         RHO      G/CM**3
CC      SPEC. VOL.      EOS80    CM**3/G
C CHECK VALUE: EOS80 = 9.435561E-1 CM**3/G FOR S = 40 NSU,
C T = 40 DEG C, P = 10000 DBARS.
C
CC UNITS CHANGED TO DECIBARS AND CGS FOR USE AT MIT
C N FOFONOFF REVISED OCT 7 1980
C *****************************
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*4 P,T,S,rRHO,EOS80
      REAL*8 K,K0,KW
C EQUIV
      EQUIVALENCE (E,D,B1,R4),(BW,B,R3),(C,A1,R2)
      EQUIVALENCE (AW,A,R1,R0),(KW,K0,K)
CC CHANGE PRESSURE FROM INPUT UNITS OF DECIBARS TO BARS
      DP=DBLE(P)
      DT=DBLE(T)
      DS=DBLE(S)
      DP=DP/10.D0
C  SQUARE ROOT SALINITY.
      SR = DSQRT(DABS(DS))
C COMPUTE DENSITY PURE WATER AT ATM PRESSURE
      R1 = ((((6.536332D-9*DT-1.120083D-6)*DT+1.001685D-4)*DT
     X-9.095290D-3)*DT+6.793952D-2)*DT+999.842594D0
C SEAWATER DENSITY ATM PRESS.
      R2 = (((5.3875D-9*DT-8.2467D-7)*DT+7.6438D-5)*DT-4.0899D-3)*DT
     X+8.24493D-1
      R3 = (-1.6546D-6*DT+1.0227D-4)*DT-5.72466D-3
      R4 = 4.8314D-4
      rRHO = SNGL((R4*DS + R3*SR + R2)*DS + R1)
C SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
      EOS80 = 1.0/rRHO
      IF(ABS(P).GE..00001)GO TO 1
CC CHANGE EOS80 TO CM**3/G,rRHO TO G/CM**3
      EOS80=EOS80*1000.
      rRHO=rRHO*.001
      RETURN
C COMPUTE COMPRESSION TERMS
 1    E = (9.1697D-10*DT+2.0816D-8)*DT-9.9348D-7
      BW = (5.2787D-8*DT-6.12293D-6)*DT+8.50935D-5
      B = BW + E*DS
C
      D = 1.91075D-4
      C = (-1.6078D-6*DT-1.0981D-5)*DT+2.2838D-3
      AW = ((-5.77905D-7*DT+1.16092D-4)*DT+1.43713D-3)*DT
     X+3.239908D0
      A = (D*SR + C)*DS + AW
C
      B1 = (-5.3009D-4*DT+1.6483D-2)*DT+7.944D-2
      A1 = ((-6.1670D-5*DT+1.09987D-2)*DT-0.603459D0)*DT+54.6746D0
      KW = (((-5.155288D-5*DT+1.360477D-2)*DT-2.327105D0)*DT
     X+148.4206D0)*DT+19652.21D0
      K0 = (B1*SR + A1)*DS + KW
C EVALUATE PRESSURE POLYNOMIAL AND RETURN
      K = (B*DP + A)*DP + K0
      EOS80 = EOS80*SNGL(1.D0 - DP/K)
CC  CHANGE EOS80 TO CM**3/G,rRHO TO G/CM**3
      EOS80=EOS80*1000.
      rRHO=rRHO*.001
      RETURN
      END
