# 温度場, 渦度場, 圧力場のデータを用いて, 温位, 渦位の値を求め,
# その東西-鉛直断面図を描画するスクリプト.

# 履歴情報 (新しい更新履歴を上に追記していくこと)
# * 2011/03/22 (北野 太朗) 加筆/修正
# * 2011/03/19 (北野 太朗) 修正
# * 2011/02/01 (北野 太朗) 新規作成

require "numru/ggraph"

include NumRu

DCL.gropn(4)

file1 = 'Basic_Temp.nc' # 温度場のデータファイル
file2 = 'Basic_Vor.nc'  # 渦度場のデータファイル
file3 = 'Basic_Ps.nc'   # 圧力場のデータファイル
# 違うデータを利用して計算したい場合は,
# ここで指定するファイルを変更する

# file 1 - 3 中の各変数に対し,
# GPhys オブジェクトを構成し,
# それぞれ変数名を割り当てている
time = GPhys::IO.open( file1,'time').val
lon = GPhys::IO.open( file1,'lon')
lat = GPhys::IO.open( file1,'lat')
sig = GPhys::IO.open(file1,'sig').val

temp = GPhys::IO.open(file1, 'Temp')
vor = GPhys::IO.open(file2, 'Vor')
ps = GPhys::IO.open(file3, 'Ps')

# 着目する経度, 緯度を指定する
setlon1 = 100..200
setlat1 = 55..70

# 指定した経度, 緯度の範囲の各値を置き換える
sublon = lon.cut('lon'=>setlon1).val
sublat = lat.cut('lat'=>setlat1).val
subtemp = temp.cut('lon'=>setlon1,'lat'=>setlat1)
subvor = vor.cut('lon'=>setlon1,'lat'=>setlat1)
subps = ps.cut('lon'=>setlon1,'lat'=>setlat1)

# i,j,k,t の大きさと最大値, 最小値を指定する
# i は経度方向に繰り返し計算する作業変数
# j は緯度方向に繰り返し計算する作業変数
# k は鉛直方向に繰り返し計算する作業変数
# t は時間方向に繰り返し計算する作業変数
ni = sublon.size
imin,imax = 1, sublon.size
nj = sublat.size
jmin, jmax = 1, sublat.size
nk = sig.size
kmin, kmax = 1, sig.size
nt = time.size
tmin, tmax = 1, time.size

subtemp = subtemp[false]
subps = subps[false]
subvor = subvor[false].val

# sig と subps は下で掛け合わせるため,
# 配列の型をあわせている
sig = sig.reshape(1,1,nk,1)
subps = subps.val.reshape(ni,nj,1,nt)

lat = subtemp.coord("lat").val*PI/180
subtemp = subtemp.val

# 温位を計算するための準備
p0 =  100000.0
pr = sig*subps
pfact = (pr/p0)**(-2.0/7.0)

# 定数値を指定する
Grav = 9.8
Omega = 7.292*(10**(-5))

# 配列を示す
theta = NArray.float(ni,nj,nk,nt).indgen     # 温位
p_ka = NArray.float(ni,nj,nk,nt)             # 渦位
ab_vor1 = NArray.float(ni,nj,nk,nt).indgen   # 絶対渦度
dthetadp1 = NArray.float(ni,nj,nk,nt).indgen # 温位の圧力微分

############### 温位を求める ##################
for t in (tmin-1)..(tmax-1)
  for i in (imin-1)..(imax-1)
    for j in (jmin-1)..(jmax-1)
      for k in (kmin-1)..(kmax-1)
        theta[i,j,k,t] = subtemp[i,j,k,t] * pfact[i,j,k,t]
      end
    end
  end
end

############### 渦位を求める ##################
for t in (tmin-1)..(tmax-1)
  for i in (imin-1)..(imax-1)
    for j in (jmin-1)..(jmax-1)
      for k in (kmin)..(kmax-3)
	  ab_vor1[i,j,k,t] = subvor[i,j,k,t]+2*Omega*Math.sin(lat[j])
	  dthetadp1[i,j,k,t] = - (theta[i,j,k+1,t]-theta[i,j,k-1,t])/(pr[i,j,k+1,t]- pr[i,j,k-1,t])
	  p_ka[i,j,k,t] = Grav * ab_vor1[i,j,k,t] * dthetadp1[i,j,k,t] 
      end
    end
  end
end

############### 温位, 渦位の配列を Gphys オブジェクトにする ###############

p lon_a = VArray.new( sublon,
                    {"long_name"=>"longitude", "units"=>"degrees"},
                    "lon" )
p lon = Axis.new.set_pos(lon_a)
 
p lat_a = VArray.new( sublat,
                    {"long_name"=>"latitude","units"=>"degrees"},
                    "lat" )
lat = Axis.new.set_pos(lat_a)

 sig_a = VArray.new( GPhys::IO.open( file1,'sig').val,
                    {"long_name"=>"sigma","units"=>"hPa"},
                    "sig" )
 sig = Axis.new.set_pos(sig_a)

time_a = VArray.new( GPhys::IO.open(file1,'time').val,
                    {"long_name"=>"time","units"=>"day"},
                    "time" )
p time = Axis.new.set_pos(time_a)
#####################################################################

 pv_k = VArray.new( p_ka,
                   {"long_name"=>"potential vorticity", "units"=>"PVU"},
                   "PV" )
 PV_k = GPhys.new( Grid.new(lon,lat,sig,time), pv_k )

ptheta = VArray.new( theta,
                   {"long_name"=>"PotVor and PotTemp (x-z)", "units"=>"K"},
                   "PotTemp" )
pottheta = GPhys.new( Grid.new(lon,lat,sig,time), ptheta )
############################################################################

time = time.pos.val        # time を GPhys オブジェクトから配列に変更

DCL.sgpset('lcntl', false) #アンダースコアー等が特殊文字として解釈されないようにする
DCL.swiset('iwidth',600)   # 横方向の描画領域を広げる
DCL.swiset('iheight',600)  # 縦方向の描画領域を広げる
GGraph.set_fig('viewport'=>[0.1,0.7,0.3,0.7]) 
                           # viewport を指定
for it in [48]             # 12 日目を計算した場合
GGraph.set_axes('xtickint'=>10,'xlabelint'=>10)
                           # 横軸は経度 10 度ごとに目盛りをつける
gphys1 =pottheta.cut('lon'=>setlon1,'lat'=>65,'sig'=>1.0..0.22,'time'=>time[it])
                           # 温位の GPhys オブジェクトで描画したい箇所を切り取る
                           # ここでは, 緯度 65 度で切った断面の, 
                           # 経度 100 - 200 度, 鉛直方向に $\sigma$ = 1.0 - 0.2 
                           # の範囲を切り取っている 
gphys2 =PV_k.cut('lon'=>setlon1,'lat'=>65,'sig'=>1.0..0.22,'time'=>time[it])
                           # 渦位の GPhys オブジェクトで描画したい箇所を切り取る
GGraph.contour(gphys1, true, 'interval'=>2.5, 'color'=>true, 'title'=>'PotVor and PotTemp (x-z)')
                           # gphys1 に示される温位を描画
                           # 等値線は色つき
                           # 等値線の間隔は 2.5 K 
                           # タイトルをつけている
mj = DCL.udpget('indxmj')
mn = DCL.udpget('indxmn')
GGraph.contour(gphys2, false, 'levels'=>[1e-7,2e-7,3e-7,4e-7,5e-7,6e-7,7e-7,8e-7,9e-7,1e-6,1.5e-6,2e-6,4e-6,6e-6,8e-6],'line_type'=>[3],'index'=>[mn,mn,mn,mn,mj,mn,mn,mn,mn,mj,mj,mj,mj,mj,mj])
                           # gphys2 に示される渦位をコンターにより描画
                           # 等値線のの間隔, ラインタイプ, インデックス
                           # を個別に指定
end

DCL.grcls