#!/usr/bin/ruby # -*- coding: utf-8 -*- # 2013/??/?? noda # 2014/08/13 Ishiwatari modified require "numru/ggraph" include NumRu require "pp" wsn = 4 RPlanet = 6.371 * 10 ** 6 # [m] #Ps = 1e5 # [Pa] Grav = 9.8 # [m/s2] title = "" # vapor の描画範囲 val_min = 0 val_max = 300 val_interval = (val_max - val_min) / 30.0 #YUNIT_vapor = 'kg m-2' #epsfn = "qvap_area_column_mass.eps" # ベクトルの描画間隔 (単位は格子点数) #xintv = 6 #yintv = 2 xintv = 4 yintv = 3 #ITR=30 # 球 ITR=1 # 平面プロット LON_CNT=125 # 球のときに使用 #xfact1 = 0.01 #yfact1 = 0.01 xfact1 = 1.0 yfact1 = 1.0 # 凡例の単位ベクトルの長さ uxunit = 0.03 uyunit = 0.03 time_start = 400 time_end = 500 omega = 1.0 type = 1 # 1: mean, 2: disturbance dir = nil # 簡易 while ARGV.size > 0 case ARGV[0] when '--wsn' then ARGV.shift wsn = ARGV[0] ARGV.shift when '--xintv' then ARGV.shift xintv = ARGV[0] ARGV.shift when '--yintv' then ARGV.shift yintv = ARGV[0] ARGV.shift when '--xfact1' then ARGV.shift xfact1 = ARGV[0] ARGV.shift when '--yfact1' then ARGV.shift yfact1 = ARGV[0] ARGV.shift when '--uxunit' then ARGV.shift uxunit = ARGV[0] ARGV.shift when '--uyunit' then ARGV.shift uyunit = ARGV[0] ARGV.shift when '--type' then ARGV.shift type = ARGV[0].to_i ARGV.shift when '--range_time' then ARGV.shift tmp = ARGV[0].split(':') time_start = tmp[0].to_f time_end = tmp[1].to_f ARGV.shift when '--range_val' then ARGV.shift tmp = ARGV[0].split(':') val_min = tmp[0].to_f val_max = tmp[1].to_f ARGV.shift val_interval = (val_max - val_min) / 30.0 else raise "ERROR: invalid option: #{ARGV[0]}" end end if type == 1 then puts "all" elsif type == 2 then puts "disturbance" else raise "Type #{type} is not supported." end #if dir == nil then # dir = File.join(DATA_HOME, # "101220_diff-DryAir-Vapor_T21L16_sr_Omega" + omega.to_s + "_dt20m") #end dir='.' # input path = File.join(dir, 'U.nc') na_sig_weight = GPhys::IO.open(path, 'sig_weight').val nsig = na_sig_weight.size var = 'U' path = File.join(dir, var+'.nc') gp_u = GPhys::IO.open(path, var).cut('time'=>time_start..time_end) var = 'V' path = File.join(dir, var+'.nc') gp_v = GPhys::IO.open(path, var).cut('time'=>time_start..time_end) var = 'QVap' path = File.join(dir, var+'.nc') gp_q = GPhys::IO.open(path, var).cut('time'=>time_start..time_end) var = 'Ps' path = File.join(dir, var+'.nc') gp_ps = GPhys::IO.open(path, var).cut('time'=>time_start..time_end) gp_u_q = gp_u * gp_q gp_v_q = gp_v * gp_q # 時間平均 gp_u_q = gp_u_q.mean('time') gp_v_q = gp_v_q.mean('time') if type == 2 then gp_u = gp_u.mean('time') gp_v = gp_v.mean('time') gp_q = gp_q.mean('time') end # 擾乱成分 if type == 2 then gp_ud_qd = gp_u_q - gp_u * gp_q gp_vd_qd = gp_v_q - gp_v * gp_q end # 鉛直積分 if type == 1 then gp_u_q = (gp_u_q * na_sig_weight.reshape(1,1,nsig)).sum('sig') gp_v_q = (gp_v_q * na_sig_weight.reshape(1,1,nsig)).sum('sig') elsif type == 2 then gp_ud_qd = (gp_ud_qd * na_sig_weight.reshape(1,1,nsig)).sum('sig') gp_vd_qd = (gp_vd_qd * na_sig_weight.reshape(1,1,nsig)).sum('sig') end # for debug #pp gp_u_q # 水蒸気量の鉛直積分 p gp_vapor = (gp_q.mean('time') * na_sig_weight.reshape(1,1,nsig)).sum('sig') \ * gp_ps.mean('time') / Grav #DCL.sgscmn(9) # colormap: white-blue-black DCL.gropn(wsn) DCL.sgpset('lcntl', false) ; DCL.uzfact(0.7) GGraph.set_axes('xtickint'=>10, 'xlabelint'=>30) GGraph.set_axes('ytickint'=>10, 'ylabelint'=>30) # 球の時 (平面プロットの時は 下の 2 行をコメントアウトするだけ?) #GGraph.set_fig 'itr'=>ITR, 'viewport'=>[0.15,0.95,0.2,0.8], 'map_axis'=>[LON_CNT, 0, 0] #GGraph.set_fig 'viewport'=>[0.15,0.95,0.2,0.8] # same as gpview vx0,vx1,vy0,vy1 = 0.15,0.85,0.2,0.55 GGraph.set_fig 'itr'=>ITR, 'viewport'=>[vx0,vx1,vy0,vy1], 'map_axis'=>[LON_CNT, 0, 0] DCL.uscset('cyspos', 'B' ) # move unit y axis DCL.udlset('LMSG', false) DCL.ugpset('XFACT1', xfact1) DCL.ugpset('YFACT1', yfact1) DCL.ugpset('UXUNIT', uxunit) # size of unit vector by dimentional value DCL.ugpset('UYUNIT', uyunit) DCL::ugrset('VXUNIT', 0.1) # size of unit vector DCL::ugrset('VYUNIT', 0.1) DCL::uxmttl('T', '', 0.0) # title (large character) DCL.ugpset('LNRMAL', false) # scale vector automatically #DCL.ugpset('LMSG', false) # no message DCL.ugpset('LUMSG', false) # no message DCL.ugpset('INDEX', 9) #DCL.ugpset('INDEX', 959) opt_vector = { 'xintv'=>xintv, 'yintv'=>yintv, 'flow_vect'=>false, 'annotate'=>false, 'unit_vect'=>true } opt_contour = { 'min'=>val_min, 'max'=>val_max, 'interval'=>val_interval, 'title'=>title, 'annotate'=>false } if type == 1 then # mean GGraph.tone_and_contour( gp_vapor, true, opt_contour) # GGraph.color_bar('landscape'=>true,'top'=>true, 'vlength'=>0.5, 'voff'=>0.03) GGraph.color_bar('landscape'=>true) GGraph.vector( gp_u_q, gp_v_q, false, opt_vector) elsif type == 2 then # disturbance GGraph.tone_and_contour( gp_vapor, true, opt_contour) # GGraph.color_bar('landscape'=>true,'top'=>true, 'vlength'=>0.5, 'voff'=>0.03) GGraph.color_bar('landscape'=>true) GGraph.vector( gp_ud_qd, gp_vd_qd, false, opt_vector) else raise "Type #{type} is not supported." end # label of vector # Ref.: # * definition of unit_vect in dclext.rb # * http://ruby.gfd-dennou.org/products/ruby-dcl/ruby-dcl-doc/grph1/node62.html space = 0.0 vxuloc = vx1 + space vyuloc = vy0 + space rsizet = 0.017 index = 1 p xunit = (uxunit.to_f * 1e7).round / 1e7 p yunit = (uyunit.to_f * 1e7).round / 1e7 if type == 1 then # mean # sfxunit = [xunit.to_s, gp_u_q.data.units.to_s].join(' ') # sfyunit = [yunit.to_s, gp_v_q.data.units.to_s].join(' ') # sfxunit = [xunit.to_s, "m s-1 kg kg -1"].join(' ') # sfyunit = [yunit.to_s, "m s-1 kg kg -1"].join(' ') sfxunit = [xunit.to_s, "m s-1"].join(' ') sfyunit = [yunit.to_s, "m s-1"].join(' ') elsif type == 2 then # disturbance sfxunit = [xunit.to_s, "m s-1"].join(' ') sfyunit = [yunit.to_s, "m s-1"].join(' ') else raise "Type #{type} is not supported." end DCL.sgtxzv(vxuloc, vyuloc-1.2*rsizet, sfxunit, rsizet, 0, -1, index) #DCL.sgtxzv(vxuloc+1.2*rsizet, vyuloc+0.5*rsizet, # sfyunit, rsizet, 90, -1, index) DCL.sgtxzv(vxuloc+2.5*rsizet, vyuloc+0.5*rsizet, sfyunit, rsizet, 90, -1, index) DCL.grcls