/* ** inteqnstate.c 1.0 Solaris 2.3 Unix 940906 SIO/ODF fmd ** ** Various calculations using the ** International Equation of State for seawater: ** ** double IESRho(double s, double t, double rp) ** calculate the potential density of seawater. ** double IESSpVolAn(double s, double t, double p, double *sigma) ** calculate specific volume anomaly and potential density anomaly. ** ** Plus, soon-to-be obsolete Fortran interface routines. ** ** References: ** ========== ** ** Millero, F. J., Chen, C. T., Bradshaw, A. and Schleicher, K., 1980. ** A New High Pressure Equation of State for Seawater. Deep-Sea ** Research 27A, 255-264. ** ** Fofonoff, N. P. and Millard, R. C., 1983. Algorithms for Computation ** of Fundamental Properties of Seawater. UNESCO Report No. 44, 15-24. */ #include #include #ifndef FORTRAN /* ** Calculate the potential density of seawater. */ double /* <- density (kg/m**3) : rho(s,t,p) */ IESRho(s, t, rp) double s; /* -> salinity (psu) */ double t; /* -> potential temperature (degrees C) */ double rp; /* -> reference pressure (decibars) */ { double bars, rhow, kstp, terma, termb, termc, kw, kst0, rho; bars = rp * 0.1; /* reference pressure in bars */ /* ** Calculate Rho(s,t,0.0) */ rhow = 999.842594 +t* ( 0.06793952 +t* (-0.00909529 +t* ( 1.001685e-4+t* (-1.120083e-6+t* 6.536332e-9)))); /* rhow = density of pure water kg/m**3 */ kw = ((t*(-0.0040899+ t*( 7.6438e-5+ t*(-8.2467e-7+ t* 5.3875e-9))))+ 0.824493)*s; /* pure water secant bulk modulus */ termc = fabs(s) * sqrt((fabs(s))); kst0 = (-0.00572466 + t*( 1.0227e-4 + t*(-1.6546e-6))) * termc; /* k(s,t,0) */ rho = rhow + kw + kst0 + 4.8314e-4 *s*s; /* rho(s,t,0.0) kg/m**3 */ /* ** Calculate pressure effects. */ if (bars > 0.0) { /* ** Rho(s,t,0) ** Rho(s,t,p) = ------------- ** p ** 1.0 - ------- ** k(s,t,p) */ kw = t*(148.4206 + t*( -2.327105 + t*( 0.01360477 + t*( -5.155288e-5)))) + 1.965221e4; kst0 = (54.6746 + t*(-0.603459 + t*( 0.0109987 + t*(-6.167e-5))))*s + kw + ( 0.07944 + t*( 0.016483 + t*(-5.3009e-4)))*termc; /* ** Calculate pressure terms. */ terma = 3.239908 + t*( 0.00143713 + t*( 1.16092e-4 + t*(-5.77905e-7)))+ ( 0.0022838 + t*(-1.0981e-5 + t*(-1.6078e-6)))*s + 1.91075e-4*termc; termb = 8.50935e-5 + t*(-6.12293e-6 + t* 5.2787e-8) + (-9.9348e-7 + t*( 2.0816e-8 + t* 9.1697e-10))*s; kstp = kst0 + bars*(terma + bars*termb); /* secant bulk modulus k(s,t,p) */ rho = rho/(1.0-bars/kstp); } return (rho); } /* IESRho() */ /* ** Calculate specific volume anomaly (steric anomaly) and ** potential density anomaly for seawater. */ double /* <- specific volume anomaly (m**3/kg) */ IESSpVolAn(s, t, p, sigma) double s; /* -> PSS-78 salinity (psu) */ double t; /* -> potential temperature (degrees C) */ double p; /* -> reference pressure (decibars) */ double *sigma; /* <- potential density anomaly (kg/m**3) */ { double bars, sigmaW, kstp, termA, termB, termC, kw, sr, sigmaRP, kst0, spVan; #define dr350 28.106331 #define r3500 (1000.0+dr350) #define v350p (1.0/r3500) /* ** check values: ** ===== ====== ** ** for s=40, t=40, p=10000: ** ** sva = 9.8130210e-6, sigma = 59.82037 ** */ bars = p * 0.1; /* reference pressure in bars */ /* ** Calculate the density anomaly at in-situ ** salinity, temperature and pressure=0.0. */ sigmaW = -28.26373741481084 + t*( 0.06793952 + t*(-0.00909529 + t*( 1.001685e-4+ t*(-1.120083e-6+ t*( 6.536332e-9))))); /* sigmaW = density anomaly of pure water kg/m**3 */ kw = (t*(-0.0040899 + t*( 7.6438e-5 + t*(-8.2467e-7 + t*( 5.3875e-9)))))+.824493; /* pure water secant bulk modulus */ sr = (s<=0.0) ? 0.0 : sqrt(s); kst0 = (-0.00572466+ t*( 1.0227e-4 + t*(-1.6546e-6)))*sr; /* k(s,t,0) */ sigmaRP= sigmaW + (kw + kst0 + s * 4.8314e-4) * s; /* sigma(s,t,0.0) kg/m**3 */ spVan = -sigmaRP*v350p/(r3500+sigmaRP); sigmaRP+=dr350; /* ** Calculate pressure effects. */ if (bars > 0.0) { /* ** sigma(s,t,0.0) ** sigma(s,t,p) = -------------- ** p ** 1.0 - -------- ** k(s,t,p) */ kw = t*(148.4206 + t*( -2.327105 + t*( 0.01360477+ t*( -5.155288e-5)))) - 1930.06; kst0 = (54.6746 + t*(-0.603459 + t*( 0.0109987 + t*(-6.167e-5))) + ( 0.07944 + t*( 0.016483 + t*(-5.3009e-4)))*sr)*s + kw; /* ** Calculate pressure terms. */ termA = -0.1194975 + t*( 0.00143713 + t*( 1.16092e-4 + t*(-5.77905e-7))) + ( 0.0022838 + t*(-1.0981e-5 + t*(-1.6078e-6)) + 1.91075e-4*sr)*s; termB = 3.47718e-5 + t*(-6.12293e-6 + t*( 5.2787e-8)) + (-9.9348e-7 + t*( 2.0816e-8 + t*( 9.1697e-10)))*s; kstp = kst0 + bars*(termA + bars*termB); /* secant bulk modulus k(s,t,p) */ termA = 2.158227e4 + bars*(3.359406 + bars* 5.03217e-5); termB = bars/termA; termC = 1.0 - termB; spVan = spVan*termC+(v350p+spVan)*bars*kstp/(termA*(termA+kstp)); termA = v350p*termC; termB = termB/termA; termC = spVan/(termA*(termA+spVan)); sigmaRP= dr350+termB-termC; } if (sigma) *sigma = sigmaRP; return (spVan); } /* IESSpVolAn() */ #else /* FORTRAN */ double rho_(s, t, rp) float *s, *t, *rp; { return (IESRho(*s, *t, *rp)); } void alpha_(prsdb, temp, salt, alph, delta) float *prsdb, *temp, *salt, *alph, *delta; { double densAn; *delta = IESSpVolAn(*salt, *temp, *prsdb, &densAn); *alph = 1.0 / (densAn + 1000.0); } /* alpha_() */ void dalpha_(prsdb, temp, salt, alph, delta) double *prsdb, *temp, *salt, *alph, *delta; { double densAn; *delta = IESSpVolAn(*salt, *temp, *prsdb, &densAn); *alph = 1.0 / (densAn + 1000.0); } /* dalpha_() */ void sigmp_(refprs, prsdb, temp, salt, sigma) float *refprs, *prsdb, *temp, *salt, *sigma; { double densAn; (void)IESSpVolAn(*salt, PotentialTemp(*salt, *temp, *prsdb, *refprs), *refprs, &densAn); *sigma = densAn; } /* sigmp_() */ void svan_(refprs, temp, salt, sva, sigma) float *refprs, *temp, *salt, *sva, *sigma; { double densAn; *sva = IESSpVolAn(*salt, *temp, *refprs, &densAn); *sigma = densAn; } void dsvan_(refprs, temp, salt, sva, sigma) double *refprs, *temp, *salt, *sva, *sigma; { double densAn; *sva = IESSpVolAn(*salt, *temp, *refprs, &densAn); *sigma = densAn; } #endif /* FORTRAN */