結果、慣性周期の数倍ぐらいで落ち着いたものの、定常解ではなく、その周りで振動する解になってしまいました(流速の大きさも振動する)。
ちなみに設定は以下の通り。コリオリ係数 f = 10-4 s-1、鉛直粘性係数 K = 5 × 10-3 m2 s-1(なので慣性周期は 6 × 104 s、エクマン層の深さは 10 m)。Δz = 1 m とし、拡散の CFL 条件を考えて Δt = 102 s。あと、密度 ρ = 103 kg m-3。鉛直には100層とって、海面の境界条件は τx = 0.0, τy = 0.1 N m-2、いちばん底では ∂u/∂z = ∂v/∂z = 0。時間方向には前方差分、鉛直方向には中央差分で計算。
《追記》
慣性周期は 2π をかけなきゃダメだった、と気付いてよく考えたら一桁勘違いしていて……ちょうど慣性周期でした。あぁ(頭悪いなぁ……)。結局、エクマン流と慣性振動の和になるみたいです。