○ main

require "narray"
require "data"
require "runge-kutta-4"

include NMath

def calc_relative_volticity(u, v)
　u.volticity(v)
end

def check(time, pv)
　psi = pv.to_stream_func
　printf("%4.2f %12.10f %12.10f\n", time, psi.energy, psi.enstrophy)
end

omega = 1.0

mm = 85
jm = 128
im = 256

#
# 強制の設定
#

pv = DATA.new(mm, jm, im)
pv.omega = omega

pv0 = pv.dup

u0 = 180.0 * 86.4 / 40000.0
phi0 = 50.0 * PI / 180.0
b0 = 10.0 * PI / 180.0

u = DATA.new(mm, jm, im)
v = u.dup

u.grid = NArray.float(jm*im)
u.grid[0..jm-1] = u0 * cos(u.latitude) * sech(2.0 * (u.latitude - phi0) / b0)

for i in 0..im-1
　u.grid[(i*jm)..(i*jm+jm-1)] = u.grid[0..jm-1]
end

v.grid = NArray.float(jm*im)

pv0.spect = calc_relative_volticity(u, v)
pv0.spect = pv0.absolute_volticity

#
# 初期値の設定
#

pv.spect = pv0.spect + 0.001000000

time = 0.0
dt = 1.0 / 80.0

check(time, pv)

for i in 1..10
　for j in 1..8
　　time, pv = integral(time, dt, pv, pv0)
　end
　check(time, pv)
end

次