FUNCTION eos, t, s
;
; potential density (p=0)
;
   sr=sqrt(abs(s))
   r1=((((6.536332E-9*t-1.120083E-6)*t+1.001685E-4)*t $
        -9.095290E-3)*t+6.793952E-2)*t+999.842594
   r2=(((5.3875E-9*t-8.2467E-7)*t+7.6438E-5)*t-4.0899E-3)*t+8.24493E-1
   r3=(-1.6546E-6*t+1.0227E-4)*t-5.72466E-3
   rhopn = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1)  

   return,  rhopn
END