Source code for MITgcmutils.jmd95

#
# created by mlosch on 2002-08-09
# converted to python by jahn on 2010-04-29

import sys
import numpy as np

__doc__ = """
Density of Sea Water using the Jackett and McDougall 1995 (JAOT 12) polynomial
"""

# coefficients nonlinear equation of state in pressure coordinates for
# 1. density of fresh water at p = 0
eosJMDCFw = [ 999.842594,
              6.793952e-02,
           -  9.095290e-03,
              1.001685e-04,
           -  1.120083e-06,
              6.536332e-09,
            ]
# 2. density of sea water at p = 0
eosJMDCSw = [     8.244930e-01,
               -  4.089900e-03,
                  7.643800e-05,
               -  8.246700e-07,
                  5.387500e-09,
               -  5.724660e-03,
                  1.022700e-04,
               -  1.654600e-06,
                  4.831400e-04,
            ]
# coefficients in pressure coordinates for
# 3. secant bulk modulus K of fresh water at p = 0
eosJMDCKFw = [    1.965933e+04,
                  1.444304e+02,
                - 1.706103e+00,
                  9.648704e-03,
                - 4.190253e-05,
             ]
# 4. secant bulk modulus K of sea water at p = 0
eosJMDCKSw = [    5.284855e+01,
                - 3.101089e-01,
                  6.283263e-03,
                - 5.084188e-05,
                  3.886640e-01,
                  9.085835e-03,
                - 4.619924e-04,
             ]
# 5. secant bulk modulus K of sea water at p
eosJMDCKP = [     3.186519e+00,
                  2.212276e-02,
                - 2.984642e-04,
                  1.956415e-06,
                  6.704388e-03,
                - 1.847318e-04,
                  2.059331e-07,
                  1.480266e-04,
                  2.102898e-04,
                - 1.202016e-05,
                  1.394680e-07,
                - 2.040237e-06,
                  6.128773e-08,
                  6.207323e-10,
            ]

[docs]def densjmd95(s,theta,p): """ Computes in-situ density of sea water Density of Sea Water using Jackett and McDougall 1995 (JAOT 12) polynomial (modified UNESCO polynomial). Parameters ---------- s : array_like salinity [psu (PSS-78)] theta : array_like potential temperature [degree C (IPTS-68)]; same shape as s p : array_like pressure [dbar]; broadcastable to shape of s Returns ------- dens : array density [kg/m^3] Example ------- >>> densjmd95(35.5, 3., 3000.) 1041.83267 Notes ----- AUTHOR: Martin Losch 2002-08-09 ([email protected]) Jackett and McDougall, 1995, JAOT 12(4), pp. 381-388 """ # make sure arguments are floating point s = np.asfarray(s) t = np.asfarray(theta) p = np.asfarray(p) # convert pressure to bar p = .1*p t2 = t*t t3 = t2*t t4 = t3*t if np.any(s<0): sys.stderr.write('negative salinity values! setting to nan\n') # the sqrt will take care of this # if s.ndim > 0: # s[s<0] = np.nan # else: # s = np.nan s3o2 = s*np.sqrt(s) # density of freshwater at the surface rho = ( eosJMDCFw[0] + eosJMDCFw[1]*t + eosJMDCFw[2]*t2 + eosJMDCFw[3]*t3 + eosJMDCFw[4]*t4 + eosJMDCFw[5]*t4*t ) # density of sea water at the surface rho = ( rho + s*( eosJMDCSw[0] + eosJMDCSw[1]*t + eosJMDCSw[2]*t2 + eosJMDCSw[3]*t3 + eosJMDCSw[4]*t4 ) + s3o2*( eosJMDCSw[5] + eosJMDCSw[6]*t + eosJMDCSw[7]*t2 ) + eosJMDCSw[8]*s*s ) rho = rho / (1. - p/bulkmodjmd95(s,t,p)) return rho
[docs]def bulkmodjmd95(s,theta,p): """ Compute bulk modulus """ # make sure arguments are floating point s = np.asfarray(s) t = np.asfarray(theta) p = np.asfarray(p) t2 = t*t t3 = t2*t t4 = t3*t # if np.any(s<0): # sys.stderr.write('negative salinity values! setting to nan\n') # the sqrt will take care of this # if s.ndim > 0: # s[s<0] = np.nan # else: # s = np.nan s3o2 = s*np.sqrt(s) #p = pressure(i,j,k,bi,bj)*SItoBar p2 = p*p # secant bulk modulus of fresh water at the surface bulkmod = ( eosJMDCKFw[0] + eosJMDCKFw[1]*t + eosJMDCKFw[2]*t2 + eosJMDCKFw[3]*t3 + eosJMDCKFw[4]*t4 ) # secant bulk modulus of sea water at the surface bulkmod = ( bulkmod + s*( eosJMDCKSw[0] + eosJMDCKSw[1]*t + eosJMDCKSw[2]*t2 + eosJMDCKSw[3]*t3 ) + s3o2*( eosJMDCKSw[4] + eosJMDCKSw[5]*t + eosJMDCKSw[6]*t2 ) ) # secant bulk modulus of sea water at pressure p bulkmod = ( bulkmod + p*( eosJMDCKP[0] + eosJMDCKP[1]*t + eosJMDCKP[2]*t2 + eosJMDCKP[3]*t3 ) + p*s*( eosJMDCKP[4] + eosJMDCKP[5]*t + eosJMDCKP[6]*t2 ) + p*s3o2*eosJMDCKP[7] + p2*( eosJMDCKP[8] + eosJMDCKP[9]*t + eosJMDCKP[10]*t2 ) + p2*s*( eosJMDCKP[11] + eosJMDCKP[12]*t + eosJMDCKP[13]*t2 ) ) return bulkmod
# aliases dens = densjmd95