# -*- coding: utf-8 -*- # # NCEP データの月平均・東西平均分布をまとめて描画する # # * 2012/02/07 カラーコンターを修正 # * 2012/01/25 カラーコンターを修正 # * 2012/01/24 カラーコンターを修正 # * 2011/12/06 drawzm_clim_mon.rb を基に作成 # require "numru/ggraph" include NumRu include NMath # cos を認識するために必要 =begin # NetCDF の ファクター等を認識しない, simple_get の使用 module NumRu class NetCDFVar alias :get :simple_get end end =end ###################################### def clim_mean( iyrs, iyre, imon, gphys ) hoursOfDay = 24 time = gphys.coord('time').val for iyr in iyrs..iyre t = (iyr - 1948) * 12 t = t + (imon - 1) ts = t # p iyrs, imon, ts if iyr == iyrs gphystmp = gphys.cut('time'=>time[ts]) else gphystmp = gphystmp + gphys.cut('time'=>time[ts]) end end gphystmp = gphystmp / ( iyre - iyrs + 1 ) return gphystmp end ############################################################## def calcmsf_notzm( gphys ) grav = 9.8 # 重力加速度 im = gphys.shape[0] # 経度の数 jm = gphys.shape[1] # 緯度の数 km = gphys.shape[2] # level の数 # p im,jm,km press = gphys.coord(2).val # level の値 # p press msf = gphys.copy # p msf.coord(2).name # 質量南北フラックス: V * (\delta p / g) の鉛直積分. 大気上端から順に求める i = km - 1 # NCEP は鉛直 17 層であり, 1 層目の値は (0) に格納されているのでループは km -1 まで. while i >= 0 # 1 層目以上について以下を行う # p i if i == km - 1 # 大気上端 (0 Pa) では南北風 0 m/s とする. (10 hPa から 0 hPa までに南北風の鉛直シアーがあると考える) msf[true,true,i,true] = gphys[true,true,i,true] * 0.5 * press[i] / grav else # 気圧の変化量は, 正の値になる様に i 番目 (下層) の値から i+1 番目 (上層) の値を引く. msf[true,true,i,true] = msf[true,true,i+1,true] + (gphys[true,true,i+1,true] + gphys[true,true,i,true] ) * 0.5 * (press[i] - press[i+1]) / grav end i = i - 1 end # 緯度依存のファクター: 2 \pi a cos(\phi) をかける pradi = 6378e3 # 惑星半径 lat = gphys.coord(1).val # 緯度の値 # /im は, 平均 = 割り算 -> 足し算 としていたために行っていた操作. factor = (cos( lat * PI / 180.0 ) * pradi * PI * 2 ).reshape!(1,jm,1) # p factor msf = msf * factor * 1e-10 # 名前と単位をつける msf.name = 'mass stream function' msf.long_name = 'mass stream function' msf.units = '1e10 kg/s' return msf end ############################################################## dir = ARGV[0] iyrs = ARGV[1].to_i iyre = ARGV[2].to_i imon = ARGV[3].to_i title = ARGV[4] vname0 = 'pres' vname1 = 'air' vname2 = 'uwnd' vname3 = 'vwnd' vname4 = 'shum' gphysPs = GPhys::NetCDF_IO.open(dir+'/'+vname0+".mon.mean.nc", vname0) gphys1 = GPhys::NetCDF_IO.open(dir+'/'+vname1+".mon.mean.nc", vname1) gphys2 = GPhys::NetCDF_IO.open(dir+'/'+vname2+".mon.mean.nc", vname2) gphys3 = GPhys::NetCDF_IO.open(dir+'/'+vname3+".mon.mean.nc", vname3) gphys4 = GPhys::NetCDF_IO.open(dir+'/'+vname4+".mon.mean.nc", vname4) #tm = gphysPs.shape[0] # number of elements for time dimension time = gphysPs.coord('time').val #< DCLのオープンと設定 > DCL.gropn(2) DCL.sldiv('y',2,2) # 2x2に画面分割, 'y'=yoko: 左上→右上→左下... DCL.sgpset('lcntl', false) # 制御文字を解釈しない DCL.sgpset('lfull',true) # 全画面表示 DCL.uzfact(0.75) # 座標軸の文字列サイズを 0.75 倍 DCL.sgpset('lfprop',true) # プロポーショナルフォントを使う #< GGraph による 描画 > GGraph.set_axes('xlabelint'=>30) # x 軸に値を打つ間隔 GGraph.set_axes('ylabelint'=>1) # y 軸に値を打つ間隔 GGraph.set_fig 'window'=>[-90,90,1e5,1e3] GGraph.set_fig 'itr'=>2, 'viewport'=>[0.10,0.80,0.15,0.55], 'yrev'=>'units:Pa' ## 鉛直座標軸の単位の変換 z = gphys1.axis("level").pos.convert_units( Units["Pa"] ) z.long_name = 'pressure' gphys1.axis("level").set_pos(z) z = gphys2.axis("level").pos.convert_units( Units["Pa"] ) z.long_name = 'pressure' gphys2.axis("level").set_pos(z) z = gphys3.axis("level").pos.convert_units( Units["Pa"] ) z.long_name = 'pressure' gphys3.axis("level").set_pos(z) z = gphys4.axis("level").pos.convert_units( Units["Pa"] ) z.long_name = 'pressure' gphys4.axis("level").set_pos(z) ## 東西・月平均値の計算 gphystmp1 = clim_mean( iyrs, iyre, imon, gphys1 ).mean('lon') gphystmp2 = clim_mean( iyrs, iyre, imon, gphys2 ).mean('lon') gphystmp3 = calcmsf_notzm( gphys3 ) gphystmp3 = clim_mean( iyrs, iyre, imon, gphystmp3 ).mean('lon') gphystmp4 = clim_mean( iyrs, iyre, imon, gphys4 ).mean('lon') # first panel (温度) gphystmp1 = gphystmp1 + 273.15 # 絶対温度への変換 gphystmp1.units = 'K' GGraph.tone( gphystmp1, true, 'lev'=>[170,180,190,200,210,220,230,240,250,260,270,280,290,300], # レベル＆パターンを陽に指定 'pat'=>[20999,30999,35999,40999,45999,50999,60999,65999,70999,75999,80999,85999,90999,95999,99999], 'annot'=>false,'titl'=>'temperature' ) # パタンの方が1つ多→まで GGraph.color_bar('vcent'=>0.35,'units'=>'('+gphystmp1.units.to_s+')') DCL::uxmttl('T', ' ', 1.0) DCL::uxmttl('T', ' ', 1.0) DCL::uxmttl('T', title, -1.0) # second panel (東西風) GGraph.tone( gphystmp2, true, 'lev'=>[-110,-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0,10,20,30,40,50,60,70], # レベル＆パターンを陽に指定 'pat'=>[20999,25999,30999,34999,38999,40999,43999,46999,50999,55999,60999,65999,73999,76999,79999,82999,88999,90999,95999,99999], 'annot'=>false,'titl'=>'eastward wind' ) # パタンの方が1つ多→まで GGraph.color_bar('vcent'=>0.35,'units'=>'('+gphystmp2.units.to_s+')') GGraph.contour( gphystmp2, false, 'lev'=>[0], # 'index'=>[0], # 'line_type'=>true, 'label'=>true ) # third panel (質量流線関数) GGraph.tone( gphystmp3, true, 'lev'=>[-24e0,-22e0,-20e0,-18e0,-16e0,-14e0,-12e0,-10e0,-8e0,-6e0,-4e0,-2e0,0,2e0,4e0,6e0,8e0,10e0,12e0,14e0,16e0,18e0,20e0,22e0,24e0], # レベル＆パターンを陽に指定 'pat'=>[20999,25999,30999,34999,36999,38999,40999,45999,50999,55999,60999,63999,65999,72999,74999,76999,78999,80999,83999,86999,88999,90999,92999,94999,96999,98999], 'annot'=>false,'titl'=>'mass stream function' ) # パタンの方が1つ多→まで GGraph.contour( gphystmp3, false, 'lev'=>[0], # 'index'=>[0], # 'line_type'=>true, 'label'=>true ) GGraph.color_bar('vcent'=>0.35,'units'=>'('+gphystmp3.units.to_s+')') # forth panel (比湿) gphystmp4 = gphystmp4 / 1000.0 # kg/kg への変換 gphystmp4.units = '1' GGraph.tone( gphystmp4, true, 'lev'=>[1e-4,2e-3,4e-3,6e-3,8e-3,10e-3,12e-3,14e-3,16e-3,18e-3,20e-3], # レベル＆パターンを陽に指定 'pat'=>[1,20999,30999,35999,40999,50999,60999,65999,70999,80999,90999,95999], 'annot'=>false,'titl'=>'specific humidity' ) GGraph.color_bar('vcent'=>0.35,'units'=> '(' +gphystmp4.units.to_s + ')') DCL.grcls