温度・塩分から海水密度の計算(主にpython)

海水の密度を計算する場面というのは多々あると思います。 私はユネスコの状態方程式を計算するfortranサブルーチンをこれまで利用してきましたが、 pythonでの解析に適したものが必要になったので、その内容をまとめました。

---

• ユネスコの状態方程式 GitHub - bjornaa/seawater: Python functions for physical properties of sea water numpyさえいらず、サーバー上でのセットアップ等が楽そうなので導入。 しかし、"this seawater package is becoming obsolete."

• TEOS-10(Thermodynamic Equation of SeaWater TEOS-10)によるもの 上記ページにてMatlab, C, Fortranのコードをダウンロード可能。 Pythonは以下のGitHubリポジトリから利用可能。 推奨版?コードGitHub - TEOS-10/GSW-Python: Python implementation of TEOS-10 GSW based on ufunc wrappers of GSW-C：コア部分がcで書かれているらしい（ゆえに速い？） 非推奨版?コードGitHub - TEOS-10/python-gsw: Python implementation of the Thermodynamic Equation Of Seawater - 2010 (TEOS-10)：pythonだけで書かれている?

• TEOS-10での変更点 What every oceanographer needs to know about TEOS-10

研究で使う場合は以下引用(詳細:Thermodynamic Equation of SeaWater TEOS-10)

If you use the GSW Oceanographic Toolbox we ask that you include a reference to McDougall and Barker (2011), whose full citation is: McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127, ISBN 978-0-646-55621-5.

---

手持ちのモデルアウトプットで使ってみる。t, s（温度、塩分）から密度等を計算してみます。

• 関数

"""
gsw.rho(SA,CT,p)

INPUT:

SA  =  Absolute Salinity                                        [ g/kg ]
CT  =  Conservative Temperature                                [ deg C ]
p   =  sea pressure                                             [ dbar ]
       (i.e. absolute pressure - 10.1325 dbar)
SA & CT need to have the same dimensions.
p may have dimensions 1x1 or Mx1 or 1xN or MxN, where SA & CT are MxN.
OUTPUT:

rho  =  in-situ density                                      [ kg m^-3 ]
"""
# EXAMPLE:
>>> SA = [34.7118, 34.8915, 35.0256, 34.8472, 34.7366, 34.7324]
>>> CT = [28.8099, 28.4392, 22.7862, 10.2262,  6.8272,  4.3236]
>>> p =  [     10,      50,     125,     250,     600,    1000]
>>> rho = gsw.rho(SA,CT,p)
>>> print(rho)
[1021.83993574 1022.26245797 1024.42719541 1027.79015276 1029.837779
 1032.00245322]

gsw.cp_t_exact(SA,t,p)

INPUT:
SA   =   Absolute Salinity                                      [ g/kg ]
t    =   in-situ temperature (ITS-90)                          [ deg C ]
p    =   sea pressure                                           [ dbar ]
        ( ie. absolute pressure - 10.1325 dbar )
SA & t need to have the same dimensions.
p may have dimensions 1x1 or Mx1 or 1xN or MxN, where SA & t are MxN.
OUTPUT:
cp_t_exact   =   heat capacity of seawater                     [ J/(kg*K) ]

---

• やってみた

"""
t: 温度
s: 塩分
p: 圧力(dbar)
"""
import time
import numpy as np
import gsw

""
t,s,pを読み込む部分(省略)
"""

print(t.shape)

rho = np.zeros(t.shape)
cp  = np.zeros(t.shape)
start = time.time()
for k in range(len(p)):
  rho[k,:,:] = gsw.rho(s[k,:,:],t[k,:,:],p[k])
  cp[k,:,:]  = gsw.cp_t_exact(cp[k,:,:],t[k,:,:],p[k])
print("elapsed_time:{0}".format(time.time()-start) + "[sec]")

結果

(63, 1280, 1440) # 配列のサイズ
elapsed_time:171.92075514793396[sec] 

まあこんなものでしょうか。例えば、CMIPモデルはもっと解像度低いですし、高解像モデルの場合でも必要な領域だけ切り出したりすれば良いでしょう。

何より、地衡流流線関数や浮力振動数計算などT,Sから計算できるものは、だいたい関数が用意されているのがありがたいです。