#!/usr/bin/env ruby =begin = mkfig_25time_series.rb -- make time-series figure of prate intensity =end require "getopts" # for option_parse require "numru/ggraph" include NumRu include Misc::EMath def cal_correlation(x, y) # 相関係数を求める xa = x - x.mean(0) # x の偏差 ya = y - y.mean(0) # y の偏差 xstddev = x.stddev # x の標準偏差(分散の平方根) ystddev = y.stddev # y :: crc = (( xa*ya ).mean(0))/(xstddev*ystddev) # 相関係数 return crc end ################################################################ # make gphys ################################################################ L = UNumeric.new(2.5008e6, 'J.Kg-1') year = [1979, 2003] #year = [1979, 1980] month = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"] # 初期値 #month = ["02", "08"] # 初期値 prate_long_name = "Precipitation" prate_law_10_10_long_name = prate_long_name + " (N10 - S10)" prate_law_s5_15_long_name = prate_long_name + " (S5 - S15)" prate_law_n5_15_long_name = prate_long_name + " (N5 - N15)" prate_law_30_30_long_name = prate_long_name + " (N30 - S30)" prate_nh_midlat_long_name = prate_long_name + " (N30 - N60)" prate_sh_midlat_long_name = prate_long_name + " (S30 - S60)" ## make compisit-GPhys object # make array which set the gphys objects. ## ガウス重み path_gw = "../../../GAUSSIAN_WEIGHT/gaussian_weight.nc" gp_gw = GPhys::NetCDF_IO.open(path_gw, 'gw') axmonth = [] i = 0 prate_name = ""; prate_unit = "" prate_law_10_10 = [] prate_law_s5_15 = [] prate_law_n5_15 = [] prate_law_30_30 = [] prate_nh_mid = [] prate_sh_mid = [] year[0].upto(year[1]) do |y| month.each do |m| # NetCDFVar => GPhys path_prate = "../../../PRATE.NCEP/PRATE.#{y}.NCEP/PRATE_#{y}-#{m}_NCEP.nc" gp_prate = GPhys::NetCDF_IO.open(path_prate, 'prate').mean('lon')*L # 時系列 VArray 作り ## prate prate_law_10_10 << ( ( ( ( gp_prate ) * gp_gw ).cut( 'lat'=>-10..10 ) ).sum.val ) prate_law_30_30 << ( ( ( ( gp_prate ) * gp_gw ).cut( 'lat'=>30..-30 ) ).sum.val ) # if m == "02" prate_law_s5_15 << ( ( ( ( gp_prate ) * gp_gw ).cut( 'lat'=>-5..-15 ) ).sum.val ) # elsif m == "08" prate_law_n5_15 << ( ( ( ( gp_prate ) * gp_gw ).cut( 'lat'=>5..15 ) ).sum.val ) # end prate_nh_mid << ( ( ( ( gp_prate ) * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val ) prate_sh_mid << ( ( ( ( gp_prate ) * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val ) prate_name = gp_prate.data.name prate_unit = gp_prate.data.units.to_s axmonth << DCL::dateg3(year[0],1,1,y.to_i,m.to_i,1) i += 1 p y if i%12 == 0 end end ### for intensity # 日付軸作成 namonth = NArray.to_na(axmonth) # 日付軸の NArray vamonth = VArray.new(namonth).rename!("month") # VArray 化. 名前は "month" vamonth.set_att('long_name','month') # 属性設定. long_name と vamonth.set_att('units','month since 1979 JAN') # units を設定. grid = Grid.new( Axis.new(true).set_cell_guess_bounds(vamonth).set_pos_to_center ) # prate を GPhys 化 ## VArray -> GPhys gp_prate_nh_mid = GPhys.new(grid, VArray.new(NArray.to_na(prate_nh_mid)).rename!(prate_name)) / gp_gw.cut( 'lat'=> 60..30 ).sum.val gp_prate_sh_mid = GPhys.new(grid, VArray.new(NArray.to_na(prate_sh_mid)).rename!(prate_name)) / gp_gw.cut( 'lat'=>-30..-60 ).sum.val #NArray.to_na(prate_law_10_10)).rename!(prate_name)) / gp_gw.cut( 'lat'=> 10..-10 ).sum.val na_prate_law_s5_15 = NArray.to_na(prate_law_s5_15) / gp_gw.cut( 'lat'=> -5..-15 ).sum.val na_prate_law_n5_15 = NArray.to_na(prate_law_n5_15) / gp_gw.cut( 'lat'=> 5..15 ).sum.val na_prate_law_30_30 = NArray.to_na(prate_law_30_30) / gp_gw.cut( 'lat'=> -30..30 ).sum.val ############ stream function ######################### strm_nh_hadley_long_name = "hadley circulation intencity (NH)" strm_sh_hadley_long_name = "hadley circulation intencity (SH)" ## make compisit-GPhys object # make array which set the gphys objects. axmonth = [] i = 0 strm_name = "" strm_unit = "" strm_hadley_north = [] strm_hadley_south = [] year[0].upto(year[1]) do |y| month.each do |m| path = "../../../STRM.NCEP/STRM.#{y}.NCEP/STRM_#{y}-#{m}_NCEP.nc" # stream function gp_climat = GPhys::NetCDF_IO.open("../../25MEANS.NCEP/25MEANS.1979-2003.NCEP/100hPa_PHYS_1979-2003_#{m}_NCEP.nc", 'strm') gp_strm = GPhys::NetCDF_IO.open(path, 'strm').cut('level'=>100..1000) gp_strm = GPhys::NetCDF_IO.open(path, 'strm') # if m == "02" strm_hadley_north << gp_strm.cut( 'lat'=> 0..30 ).data.val.max - gp_climat.cut( 'lat'=> 0..30 ).data.val.max # elsif m == "08" strm_hadley_south << gp_strm.cut( 'lat'=> 0..-30 ).data.val.min - gp_climat.cut( 'lat'=> 0..-30 ).data.val.min # end strm_name = gp_strm.data.name strm_unit = gp_strm.data.get_att("units") end end ### for intensity na_strm_hadley_north = NArray.to_na(strm_hadley_north) na_strm_hadley_south = NArray.to_na(strm_hadley_south) ################################################################ # 描画ルーチン ################################################################ ## # 事前準備 DCL.uscset('cyspos', 'B' ) # y 軸の単位の位置を下方へ rsizel2 = DCL.uzrget('rsizel2') # 現在のラベルサイズを取得 DCL.uzrset('rsizel2', rsizel2*0.42 ) # ラベルサイズをデフォルトの 0.5 倍に ## # お絵描きメイン DCL.swpset('lsep', true) # ページ毎別々のファイルに落す DCL.gropn(2) DCL.sldiv('T', 2, 2) DCL.sgpset('lcntl', false) # DCL.sgpset('lfull',true) # フルスクリーン DCL.uzfact(0.55) # DCL.sgpset('lcorner',false) # コーナーを取っ払う DCL.sgpset('lfprop',true) # DCL.udpset('lmsg',false) # DCL.uscset('cyspos', 'B' ) # y 軸の単位の位置を下方へ ## y 軸範囲 axx = NArray[-5e10, 5e10] # x 軸範囲. -5e10 から 5e10 axy = NArray[100, 200] # y 軸範囲. 0 から 200 W/m^2 ### 1 枚目 : 北半球ハドレー循環強度 & 降水 (s5 - s15) DCL::grfrm # ページ確定 x = na_strm_hadley_north y = na_prate_law_s5_15 DCL::usspnt(axx, axy) DCL::uspfit DCL::grstrf DCL::ussttl('Mass stream function', strm_unit, 'precipitation (S5-S15) ', prate_unit ) DCL::usdaxs DCL::uxmttl('T', 'Hadley intensity (NH) & Precipitation ', 0.0) DCL::uusmkt(2) DCL::uumrk(x, y) print "北半球ハドレー循環強度と降水(N5-N15)", cal_correlation(x, y), "\n" ### 2 枚目 : 南半球ハドレー循環強度 & 降水 (n5 - n15) axy = NArray[150, 250] # y 軸範囲. 0 から 200 W/m^2 DCL::grfrm # ページ確定 x = na_strm_hadley_south y = na_prate_law_n5_15 DCL::usspnt(axx, axy) DCL::uspfit DCL::grstrf DCL::ussttl('Mass stream function', strm_unit, 'precipitation (N5-N15)', prate_unit ) DCL::usdaxs DCL::uxmttl('T', 'Hadley intensity (SH) & Precipitation (N5 - N15)', 0.0) DCL::uusmkt(2) DCL::uumrk(x, y) print "南半球ハドレー循環強度と降水(S5-S15)", cal_correlation(x, y), "\n" ### 3 枚目 : 北半球ハドレー循環強度 & 南北 ITCZ 降水差 axy = NArray[-150, 150] # y 軸範囲. -100 から 200 W/m^2 DCL::grfrm # ページ確定 x = na_strm_hadley_north y = na_prate_law_s5_15 - na_prate_law_n5_15 DCL::usspnt(axx, axy) DCL::uspfit DCL::grstrf DCL::ussttl('Mass stream function', strm_unit, 'precipitation ', prate_unit ) DCL::usdaxs DCL::uxmttl('T', 'Hadley intensity (NH) & Precipitation anomaly [(S5-S15) - (N5-N15)]', 0.0) DCL::uusmkt(2) DCL::uumrk(x, y) print "北半球ハドレー循環強度と南北降水差[(S5-S15) - (N5-N15)]", cal_correlation(x, y), "\n" ### 4 枚目 : 南半球ハドレー循環強度 & 南北 ITCZ 降水差 DCL::grfrm # ページ確定 x = na_strm_hadley_south y = na_prate_law_n5_15 - na_prate_law_s5_15 DCL::usspnt(axx, axy) DCL::uspfit DCL::grstrf DCL::ussttl('Mass stream function', strm_unit, 'precipitation ', prate_unit ) DCL::usdaxs DCL::uxmttl('T', 'Hadley intensity (SH) & Precipitation anomaly [(N5-N15) - (S5-S15)]', 0.0) DCL::uusmkt(2) DCL::uumrk(x, y) print "南半球ハドレー循環強度と南北降水差[(N5-N15) - (S5-S15)]", cal_correlation(x, y), "\n" =begin # マーカー情報記述 charsize = 1.0 * DCL.uzpget('rsizec1') dvx = 0.01 dvy = charsize*1.5 vxmin,vxmax,vymin,vymax = DCL.sgqvpt vx = 0.12 vy = 0.1 + charsize/2 DCL::sgpmu([vx], [vy]) DCL::sgtxzv(vx+dvx,vy,uwnd_nh_20_40_long_name, charsize, 0, -1, 1) =end DCL.grcls ################################################## # 画像ファイル名を変更 File.rename('dcl_001.ps', 'PRATE_1973-2003_TIME-SERIES_NCEP_5_15_sanpuzu.ps')