FORTRAN 77 風 Ruby 版モデル(test12d.rb)

require 'narray'
require 'ispack2.so' # ラッパー

include NMath # NArray 配列まるごと sin, exp などの演算ができるモジュール。

def check(mm,time,pv) # エネルギー・エンストロフィーをチェック。
　omg = 1.0
　psi = Ispack2.ndca2p(mm,omg,pv)
　ene = Ispack2.ndgeea(mm,psi)
　ens = Ispack2.ndgena(mm,psi)
　printf("%4.2f %12.10f %12.10f\n",time, ene, ens)
end

def pv0init(mm,jm,im,q,r,it,t,y) # 強制の設定
　omg = 1.0

　u0 = 180.0 * 86.4 / 40000.0 # m/s -> a Omega
　phi0 = 50.0 * PI / 180.0　　# deg -> rad
　b0　= 10.0 * PI / 180.0　　 # deg -> rad

　u = NArray.float(jm*im)
　u[0..jm-1] = u0*cos(y)*sech( 2.0*(y-phi0)/b0 ) # sech型ジェット

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

　v = NArray.float(jm*im).fill(0.0)

　pv0 = Ispack2.stvrsa(mm,jm,im,u,v,q,r,it,t) # u,vを相対渦度に変換
　pv0 = Ispack2.ndtv2a(mm,omg,pv0)　　　　　　# 相対渦度を絶対渦度に変換

　return pv0
end

def derivative(pv,pv0,mm,jm,im,q,r,it,t) # dq/dt = -J(ψ,q) -α(q-q0) +νΔ(q-q0)

　rank = 1
　nu = 1.0/(mm*(mm+1)-2.0)**rank
　alpha = 0.1

　omg = 1.0

　psi = Ispack2.ndca2p(mm,omg,pv)
　jacobian = Ispack2.stajba(mm,jm,im,psi,pv,q,r,it,t) * (-2*PI)
　relax = - alpha * (pv - pv0)
　viscous = nu * Ispack2.stclfa(mm, pv - pv0)

　dpvdt = jacobian + relax + viscous

　return dpvdt
end

def integral(time,pv,mm,jm,im,dt,q,r,it,t,pv0) # ルンゲクッタで積分

　d1 = derivative(pv,pv0,mm,jm,im,q,r,it,t)
　w = pv + d1*dt/2.0
　d2 = derivative(w,pv0,mm,jm,im,q,r,it,t)
　w = pv + d2*dt/2.0
　d3 = derivative(w,pv0,mm,jm,im,q,r,it,t)
　w = pv + d3*dt
　d4 = derivative(w,pv0,mm,jm,im,q,r,it,t)

　pv = pv + (d1 + d2*2.0 + d3*2.0 + d4)/6.0*dt
　time　= time + dt

　return time,pv
end

mm = 85
jm = 128
im = 256

time = 0.0
dt = 1.0/80.0

q,r,it,t = Ispack2.stinit(mm,jm,im) # 球面調和関数の初期化
y,x = Ispack2.stogrd(jm,im,q) # グリッドの座標

pv0 = pv0init(mm,jm,im,q,r,it,t,y)
pv = pv0 + 0.001000000 # 初期値の設定

check(mm,time,pv)

for i in 1..10
　for j in 1..8
　　time,pv = integral(time,pv,mm,jm,im,dt,q,r,it,t,pv0)
　end
　check(mm,time,pv)
end

次