はじめに
以前の記事では、Python の windspharm パッケージを用いた流線関数・速度ポテンシャルの計算方法を紹介した（なお、流線関数や速度ポテンシャル自体についての解説はこの記事では行わないので、リンク先の記事を参照してほしい）。
windspharm は Intel Mac では conda コマンドでインストールすることができるが、Apple Silicon向けにはインストールできない（最終更新が2018年とメンテナンスがなされていないようである）。
そこでCDOを用いた計算方法を調べたのでメモしておく。
なお、以下の記事を大きく参考にした。
実行環境は以下のとおりである
- Macbook Air (M1, 2020)
- macOS Monterey 12.4
- CDO version 2.0.5 (Homebrewでインストール)
計算方法
手元に東西風(ps_u.nc)と南北風(ps_v.nc)が含まれているファイルがあるとする。ファイルに1時刻、複数気圧面が含まれているとする。水平1度格子で、出力の流線関数と速度ポテンシャルも同じ格子で得たいとする。
サンプルコードを以下に示す。CDOのパイプラインの機能を使っており一つのコマンドが長くなるので適宜改行している。また、後で参照するときにわかりやすいように各行にコメントで番号を振っている。番号はパイプラインで処理される順番（右から）に対応している。
cdo -L \ # (0)
sp2gp \ # (8)
-dv2ps \ # (7)
-uv2dv \ # (6)
-remapbil,F90 \ # (5)
-chname,ps_u,u,ps_v,v \ # (4)
-setmisstoc,0 \ # (3)
-sellevel,850 \ # (2)
-merge ps_u.nc ps_v.nc \ # (1)
tmp.nc # (9)
cdo -L remapbil,outgrid tmp.nc ps_sfvp200.nc # (10)
(0) -Lは Lock I/O (sequential access) の意味。マルチスレッド環境で非スレッドセーフなNetCDF4/HDF5ライブラリを使うとエラーが起きうるとの警告が出たのでつけている。
cdo sp2gp (Warning): Using a non-thread-safe NetCDF4/HDF5 library in a multi-threaded environment may lead to erroneous results!
cdo sp2gp (Warning): Use a thread-safe NetCDF4/HDF5 library or the CDO option -L to avoid such errors.
(1) 東西風と南北風のファイルを一つに結合する。気圧面の切り出しをしてから結合したほうが効率的だが、その場合はパイプラインで繋ぐことができないため、東西風と南北風の一時ファイルを作成することになる。
(2) 複数気圧面を一度に処理できないため、一つの気圧面データを切り出す。ファイル内に複数時刻が含まれている場合は、-seltimestep,1のように特定時刻を切り出す必要があるかもしれない(未確認)。
(3) uv2dvをする際に、欠損値があると計算ができないため、欠損値を埋める。気圧面が地表面よりも下にある場合に欠損が生じるため、埋める値を0とすることで実用上問題ないだろう。
(4) uv2dvをする際に、変数名がu,vにする必要があるため、変数名の変更を行う。
(5) uv2dvをする際に、入力データはガウス格子にする必要があるため、等緯度経度格子からバイリニア補間で格子変換する。F90のFはFull regular Gaussian gridの意味。90は極から赤道までの格子数で、この場合は東西360x南北180格子となる。
(6) ガウス格子上の東西風と南北風から、発散と渦度の球面調和係数を計算する。
(7) 発散と渦度の球面調和係数から、流線関数と速度ポテンシャルの球面調和係数を計算する。
(8) 流線関数と速度ポテンシャルの球面調和係数から、ガウス格子上での流線関数と速度ポテンシャルの値を計算する。
(9) 中間ファイルに出力する。後続のガウス格子から等緯度傾度格子へと変換する部分もパイプラインにできるとよかったのだが、以下のようなエラーが出てしまう。
cdo remapbil (Abort): Unsupported spectral coordinates (Variable: )!
(10) ガウス格子から等緯度傾度格子へと変換する。outgridは出力格子で、remapbil,r360x180とすれば、経度0度、南緯89.5度開始の1度格子となる。東経0.5度開始にしたい場合は、outgridという名前のファイルに以下のように記述すればよい。
gridtype = latlon
xsize = 360
ysize = 180
xunits = 'degree'
yunits = 'degree'
xfirst = 0.5
xinc = 1.0
yfirst = -89.5
yinc = 1.0