#!/usr/bin/env ruby load "/work11/ape/yukiko/lib/ape-view.rb" # -------------------------------------------------------------------- END{ $rezol = ["T39L48_eml", "T39L24_eml", "T39L96_eml", "T79L48_eml" ] $t= ape_new "/work11/ape/NetCDF/#{$rezol[0]}/AGUforAPE-03a_TR_control.nc" # ダンプする場合 DCL.swpset('LDUMP', true) DCL.swpset('LWAIT',false) ; DCL.swpset('LWAIT0',false) DCL.swpset('LWAIT1',false) DCL.swpset('IPOSX', 50) ; DCL.swpset('IPOSY',50) `rm dcl_* tmp.pnm` $t.gropn(1) =begin # ps にする場合 $t.gropn(2) 'rm dcl*.ps' =end $rezol.each{ |rezol| $t= ape_new "/work11/ape/NetCDF/#{rezol}/AGUforAPE-03a_TR_control.nc" g = $t.go("tr_tppn") gs = g.cut(true,0,true)[true,-401..-1].stspct_fft("#{g.name}_spct"). set_att("ape_name","precipitation flux S-T power spectrum"). set_att("units","(kg m-2 s-1 day-1)2") dim1 = gs.coord(0).shape.to_s.to_i dim2 = gs.coord(1).shape.to_s.to_i $file_label= "#{rezol}_control" p $file_label $t.plot_spct_bunsan(gs[((dim1+1)/2-31)..((dim1+1)/2+29),0..80]) file_name = "#{$file_label}_#{gs.name}_wk1999.gif" =begin # ps にする場合 `dclpsrmcm dcl*.ps | dclpsrot > tmp.ps` `eps2eps tmp.ps #{$fig_path}#{file_name}` `rm tmp.ps dcl*.ps` =end # ダンプする場合 GGraph.tone( gs[((dim1+1)/2-31)..((dim1+1)/2+29),0..80]) `for file in *.xwd ; do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm; done` `ppmtogif tmp.pnm > #{$fig_path}#{file_name}` `rm dcl_* tmp.pnm` } $rezol = ["T159L48_eml","T319L48_eml"] $rezol.each{ |rezol| $gr2ncfile_path ="/work11/ape/yukiko/data/#{rezol}_expID01/" file_name = Array.new Dir.foreach($gr2ncfile_path) { |item| file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/ } $t = Ape.new(file_name) g = $t.go("tr_tppn") gs = g.cut(true,0,true)[true,-401..-1].stspct_fft("#{g.name}_spct"). set_att("ape_name","precipitation flux S-T power spectrum"). set_att("units","(kg m-2 s-1 day-1)2") dim1 = gs.coord(0).shape.to_s.to_i dim2 = gs.coord(1).shape.to_s.to_i $file_label= "#{rezol}_control" p $file_label $t.plot_spct_bunsan(gs[((dim1+1)/2-31)..((dim1+1)/2+29),0..80]) file_name = "#{$file_label}_#{gs.name}_wk1999.gif" =begin # ps にする場合 `dclpsrmcm dcl*.ps | dclpsrot > tmp.ps` `eps2eps tmp.ps #{$fig_path}#{file_name}` `rm tmp.ps dcl*.ps` =end # ダンプする場合 GGraph.tone( gs[((dim1+1)/2-31)..((dim1+1)/2+29),0..80]) `for file in *.xwd ; do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm; done` `ppmtogif tmp.pnm > #{$fig_path}#{file_name}` `rm dcl_* tmp.pnm` } $t.grcls `rm tmp.pnm dcl_*` } # -------------------------------------------------------------------- class Ape def plot_spct_bunsan(gphys) gropn(1) if $gropn == nil plot_set(gphys) # attribute の long_name を消す (タイトル描画位置の都合) if gphys.data.get_att("long_name") gphys = gphys.set_att("long_name","") end GGraph.tone( gphys, true , $tone_hash ) # 分散曲線 $x_size = gphys.grid_copy.axis(0).pos.val h_array = NArray.to_na([ 12, 25, 50, 100, 200 ]) omega = 2*3.1415/24/60/60 beta = 2*(omega)/(4e+7) * ( 60*60*24*(6370e+3) ) # (2*Omega/a)*cos(0) name_ary = ["h = 12 ", "h = 25 ","h = 50 ","h =100", "h =200 "] # line_index_ary = [4,24,104,64,724] line_index_ary = [1,21,101,61,721] for i in (0..4) h = h_array[i] bunsan_calc(h) GGraph.line( $kelvin, false, "index" => line_index_ary[i] ) GGraph.line( $rossby, false, "index" => line_index_ary[i] ) GGraph.line( $grav1, false, "index" => line_index_ary[i] ) end GGraph.line( $grav1*0.0, false, "index" => line_index_ary[0] ) __line_index(name_ary,line_index_ary) # タイトル tani = "(#{gphys.data.get_att("units")})" DCL::uxsttl("t", tani , -1.0 ) DCL::uxmttl("t", gphys.data.get_att("ape_name").gsub("_", " "), -1.0 ) # nc ファイル名 if $gropn == 2 charsize = 0.79 - $file_label.size * 0.0115 DCL::sgtxzv(charsize,0.62,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1) else charsize = 0.81 - $file_label.size * 0.0095 DCL::sgtxzv(charsize,0.6,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1) end # トーンバー unless gphys.axnames.size == 1 DCL::Util::color_bar($cbar_conf) unless $tone_flag == false end end def __line_index(str_ary,line_index_ary) charsize = 0.7 * DCL.uzpget('rsizec1') dvx = 0.01 dvy = charsize*1.5 raise TypeError,"Array expected" if ! str_ary.is_a?(Array) vxmin,vxmax,vymin,vymax = DCL.sgqvpt vx = 0.15 vy = 0.12 - charsize/2 str_ary.size.times{|num| DCL::sgtxzv(vx,vy,"--- #{str_ary[num]}", charsize, 0, -1,line_index_ary[num]) vy -= dvy if num == 2 vx = 0.30 vy = 0.12 - charsize/2 end } nil end def bunsan_calc(h=200) omega = 2*3.1415/24/60/60 beta = 2*(omega)/(4e+7) * ( 60*60*24*(6370e+3) ) # (2*Omega/a)*cos(0) c_g = ((h*9.8)**0.5)*60*60*24/(4e+7) # 軸 knum = VArray.new($x_size).rename("knum").put_att("units","1") knum = Axis.new().set_pos(knum) grid = Grid.new(knum) ## 虚数を含む軸に変換 x_size = NArray.complex($x_size.size).fill!(1.0) * $x_size # モード num = 1 ## 3 次方程式解法: カルダノの公式 ## x^3 + 3px + q ## t^2 + qt - p =0, t = -q/2 +- (q**2/4 + p)**0.5 p = - ( (c_g*x_size)**2 + c_g*beta*(2*num+1)) / 3 q = - c_g**2*beta*x_size t1 = -q/2 + (q**2/4 + p**3)**0.5 t2 = -q/2 - (q**2/4 + p**3)**0.5 t1 = t1**(1.0/3) t2 = t2**(1.0/3) # 重力波 $grav1 = (t1 + t2).real $grav1 = VArray.new($grav1).rename("grav1"). put_att("units","1").put_att("ape_name","gravity") $grav1 = GPhys.new( grid, $grav1) # 反対重力波 (図示しない) $grav2 = (t1 * Complex.polar(1, Math::PI*2/3) + t2 * Complex.polar(1, Math::PI*4/3)).real $grav2 = VArray.new($grav2).rename("grav2"). put_att("units","1").put_att("ape_name","gravity") $grav2 = GPhys.new( grid, $grav2) # ロスビー波 $rossby = (t2 * Complex.polar(1, Math::PI*2/3) + t1 * Complex.polar(1, Math::PI*4/3) ).real $rossby2 = $rossby $rossby = ($rossby > 0) * $rossby $rossby = VArray.new($rossby).rename("rossby"). put_att("units","1").put_att("ape_name","rossby") $rossby = GPhys.new( grid, $rossby) $rossby2 = VArray.new($rossby2).rename("rossby"). put_att("units","1").put_att("ape_name","rossby") $rossby2 = GPhys.new( grid, $rossby2) # ケルビン波 $kelvin = (c_g*x_size).real $kelvin = ($kelvin > 0) * $kelvin $kelvin = VArray.new($kelvin).rename("kelvin"). put_att("units","1").put_att("ape_name","kelvin") $kelvin = GPhys.new( grid, $kelvin) end end