# 月平均した物理量の緯度経度分布を描画するスクリプト # dcpam と NCEP を 2 枚並べて描画する用 # # * 2011/12/04 ファイル名, タイトルを変更 # * 2011/12/01 OLRA にも対応 # * 2011/11/29 viewport などの細かい設定を修正 # * 2011/11/28 drawxy_clim_mon.rb をもとに作成 require "numru/ggraph" include NumRu ###################################### def clim_mean( iyrs, iyre, imon, gphys, data ) if data == "dcpam" then dataNumOfDay = 4 daysOfMonth = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] daysofyear = 0 daysOfMonth.each do |i| daysofyear += i end time = gphys.coord('time').val for iyr in iyrs..iyre t = daysofyear*(iyr-1) for i in 1..imon-1 t = t + daysOfMonth[i-1] end ts = t te = t + daysOfMonth[imon-1] ts = ts * dataNumOfDay te = te * dataNumOfDay - 1 if iyr == iyrs gphystmpD = gphys.cut('time'=>time[ts]..time[te]).mean('time') else gphystmpD = gphystmpD + gphys.cut('time'=>time[ts]..time[te]).mean('time') end gphystmpD = gphystmpD / ( iyre - iyrs + 1 ) end return gphystmpD elsif data == "NCEP" then hoursOfDay = 24 daysofyear = 0 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 gphystmpN = gphys.cut('time'=>time[ts]) else gphystmpN = gphystmpN + gphys.cut('time'=>time[ts]) end end gphystmpN = gphystmpN / ( iyre - iyrs + 1 ) return gphystmpN end end ###################################### iyrsD = ARGV[4].to_i iyreD = ARGV[5].to_i iyrsN = ARGV[6].to_i iyreN = ARGV[7].to_i imon = ARGV[8].to_i title = ARGV[9] dirD = ARGV[0] dirN = ARGV[1] # 物理変数の設定と gphys の準備 ## dcpam 計算結果用 vnameD = ARGV[2] gphysD = GPhys::NetCDF_IO.open(dirD+'/'+vnameD+".nc", vnameD) tmD = gphysD.shape[2] # number of elements for time dimension timeD = gphysD.coord('time').val ## NCEP 用 if vnameD == "OLRA" || vnameD == "OSRA" # 放射フラックスの場合 vnameN1 = ARGV[3] vnameN1_longname = vnameN1 + "top" vnameN2 = 'uswrf' vnameN2_longname = vnameN2 + "top" gphysN1 = GPhys::NetCDF_IO.open(dirN+'/'+vnameN1_longname + ".mon.mean"+".nc", vnameN1) gphysN2 = GPhys::NetCDF_IO.open(dirN+'/'+vnameN2_longname + ".mon.mean"+".nc", vnameN2) if vnameD == "OSRA" gphysD = - gphysD gphysN = gphysN1 - gphysN2 else gphysD = gphysD gphysN = gphysN1 end else # 放射フラックス以外の場合 vnameN = ARGV[3] gphysN = GPhys::NetCDF_IO.open(dirN+'/'+vnameN+ ".mon.mean"+".nc", vnameN) end tmN = gphysN.shape[0] # number of elements for time dimension timeN = gphysN.coord('time').val #< DCLのオープンと設定 > DCL.gropn(2) DCL.sgpset('lcntl', false) # 制御文字を解釈しない DCL.sgpset('lfull',true) # 全画面表示 DCL.sgpset('lcorner', false) # コーナーマークを書かない DCL.uzfact(0.5) # 座標軸の文字列サイズを 0.5 倍 DCL.sgpset('lfprop',true) # プロポーショナルフォントを使う DCL.udpset('lmsg',false) # コンター間隔非表示 #< GGraph による 描画 > vpt = NArray[0.2, 0.7, 0.01, 0.21] vpt1 = (vpt + ([0.050]*2 + [0.32]*2) ).to_a vpt2 = (vpt + ([0.050]*2 + [0.10]*2) ).to_a GGraph.set_axes('xlabelint'=>90) # x 軸に値を打つ間隔 GGraph.set_axes('ylabelint'=>30) # y 軸に値を打つ間隔 GGraph.set_fig 'itr'=>10, 'window'=>[0,360,-90,90] GGraph.set_map 'coast_world'=>true, 'grid'=>false data = "dcpam" gphystmpD = clim_mean( iyrsD, iyreD, imon, gphysD, data ) data = "NCEP" gphystmpN = clim_mean( iyrsN, iyreN, imon, gphysN, data ) if vnameD == 'SurfTemp' GGraph.tone( gphystmp, true, 'lev'=>[210,220,230,240,250,260,270,280,290,300], # レベル＆パターンを陽に指定 'pat'=>[20999,30999,40999,45999,50999,55999,60999,70999,80999,85999,90999], # パタンの方が1つ多→±∞まで 'map_axes'=>true ) #GGraph.contour( gphystmp, false, 'lev'=>[200,210,220,230,240,250,260,270,280,290,300], 'index'=>3 ) # 参考まで elsif vnameD == 'Rain' # 1 枚目の描画の設定 GGraph.set_fig('viewport'=>vpt1) GGraph.set_axes('xunits'=>'', 'yunits'=>'', 'xtitle'=>'') DCL.uzpset('labelxb',false) # 描画 gphystmpD = gphystmpD / 2.5e6 GGraph.tone( gphystmpD, true, 'lev'=>[0,1e-5,2e-5,4e-5,8e-5,12e-5,16e-5,20e-5,24e-5,28e-5,32e-5], # レベル＆パターンを陽に指定 'pat'=>[1,20999,30999,35999,40999,50999,60999,65999,70999,75999,80999,85999], # パタンの方が1つ多→±∞まで 'map_axes'=>true, 'annot'=>false, 'titl'=>'' ) # 図の左上のタイトル DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','dcpam5',-1) ; DCL.uzpset('pad1',0.7) # 2 枚目の描画の設定 GGraph.set_fig('viewport'=>vpt2,'new_frame'=>false) GGraph.set_axes('xunits'=>nil, 'yunits'=>'', 'xtitle'=>'longitude','ytitle'=>'latitude') DCL.uzpset('labelxb',true) # 描画 GGraph.tone( gphystmpN, true, 'lev'=>[0,1e-5,2e-5,4e-5,8e-5,12e-5,16e-5,20e-5,24e-5,28e-5,32e-5], # レベル＆パターンを陽に指定 'pat'=>[1,20999,30999,35999,40999,50999,60999,65999,70999,75999,80999,85999], # パタンの方が1つ多→±∞まで 'map_axes'=>true, 'annot'=>false, 'titl'=>'' ) # 図の左上のタイトル DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','NCEP',-1) ; DCL.uzpset('pad1',0.7) #GGraph.contour( gphystmp, false, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], 'index'=>3 ) # 参考まで # カラーバーの設定 GGraph.color_bar( 'voff'=>0.05, 'vcent'=>0.3,'units'=>'(kg/m/s)' ) # 図のタイトルの設定 DCL::sgtxzv(0.5,vpt1[3]+0.06, 'Monthly mean precipitation '+'('+title+')', 1.15*DCL.uzpget('rsizec2'),0,0,3) elsif vnameD == 'SoilMoist' GGraph.tone( gphystmp, true, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], # レベル＆パターンを陽に指定 'pat'=>[15999,20999,25999,30999,35999,40999,45999,50999,55999,60999,65999,70999,75999,80999,85999,90999,95999], # パタンの方が1つ多→±∞まで 'map_axes'=>true ) #GGraph.contour( gphystmp, false, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], 'index'=>3 ) # 参考まで elsif vnameD == 'SurfSnow' GGraph.tone( gphystmp, true, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], # レベル＆パターンを陽に指定 'pat'=>[15999,20999,25999,30999,35999,40999,45999,50999,55999,60999,65999,70999,75999,80999,85999,90999,95999], # パタンの方が1つ多→±∞まで 'map_axes'=>true ) #GGraph.contour( gphystmp, false, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], 'index'=>3 ) # 参考まで elsif vnameD == 'OLRA' # 1 枚目の描画の設定 GGraph.set_fig('viewport'=>vpt1) GGraph.set_axes('xunits'=>'', 'yunits'=>'', 'xtitle'=>'') DCL.uzpset('labelxb',false) GGraph.tone( gphystmpD, true, 'lev'=>[100,120,140,160,180,200,220,240,260,280,300], # レベル＆パターンを陽に指定 'pat'=>[10999,15999,20999,25999,30999,35999,40999,50999,60999,70999,80999,90999], # パタンの方が1つ多→±∞まで 'map_axes'=>true, 'annot'=>false, 'titl'=>'') # 図の左上のタイトル DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','dcpam5',-1) ; DCL.uzpset('pad1',0.7) # 2 枚目の描画の設定 GGraph.set_fig('viewport'=>vpt2,'new_frame'=>false) GGraph.set_axes('xunits'=>nil, 'yunits'=>'', 'xtitle'=>'longitude','ytitle'=>'latitude') DCL.uzpset('labelxb',true) # 描画 GGraph.tone( gphystmpN, true, 'lev'=>[100,120,140,160,180,200,220,240,260,280,300], # レベル＆パターンを陽に指定 'pat'=>[10999,15999,20999,25999,30999,35999,40999,50999,60999,70999,80999,90999], # パタンの方が1つ多→±∞まで 'map_axes'=>true, 'annot'=>false, 'titl'=>'') # 図の左上のタイトル DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','NCEP',-1) ; DCL.uzpset('pad1',0.7) #GGraph.contour( gphystmp, false, 'lev'=>[0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150], 'index'=>3 ) # 参考まで # カラーバーの設定 GGraph.color_bar( 'voff'=>0.05, 'vcent'=>0.3,'units'=>'(W/m^2)' ) # 図のタイトルの設定 DCL::sgtxzv(0.5,vpt1[3]+0.06, 'Monthly OLR '+'('+title+')', 1.15*DCL.uzpget('rsizec2'),0,0,3) end DCL.grcls