#!/usr/bin/env ruby =begin = mkfig.rb -- make figure for zonal mean. (a)オイラー平均質量流線関数, (b)東西風, (c)温度, (d)比湿, (e)放射収支(黒線:大気上端正味短波フラックス, 赤線:大気上端上向き長波フラックス, 緑線:地表面正味短波フラックス, 青線:地表面正味長波フラックス), (f)熱収支(黒線:地表面正味潜熱フラックス, 赤線:地表面正味顕熱フラックス, 緑線:降水), (g)南北乾燥静的エネルギー輸送量 (黒線:合計, 赤線:帯状平均成分, 緑線:帯状擾乱成分), (h)南北潜熱輸送量 (黒線:合計, 赤線:帯状平均成分, 緑線:帯状擾乱成分), (i)EP フラックスとその発散 (矢羽: EP フラックス, トーン:発散), (j)残差循環の質量流線関数 を一枚にまとめた図. 以下の実行時オプションを指定することで, 各年平均, 25 年平均の絵がかける. --kikan : 年平均, 季節平均の絵かを指定. 以下の文字列の何れかを入力 "annual", "winter", "spring", "summer", "fall" --range : 平均する期間の指定. (開始年):(終了年) の形式で指定. とある一年間の平均の場合, 開始年と終了年は同一年を入力. == 実行例 1979 - 2003 年 の年平均図 $ ruby mkfig_25years.rb --kikan annual --range 1979:2003 1979 年 の冬平均図 $ ruby mkfig_25years.rb --kikan winter --range 1979:1979 ... =end require "getopts" # for option_parse require "numru/ggraph" require "libgphys-n" include NumRu include Misc::EMath ## enable to the option min, max module NumRu module GGraph module_function def vector_for_dcchart(fx, fy, newframe=true, options=nil) if ! defined?(@@vector_for_dcchart_options) @@vector_for_dcchart_options = Misc::KeywordOptAutoHelp.new( ['newfig', true, 'if false, do not cleared before figure setting.'], ['title', nil, 'Title of the figure(if nil, internally determined)'], ['annotate', true, 'if false, do not put texts on the right margin even when newframe==true'], ['transpose', false, 'if true, exchange x and y axes'], ['xintv', 1, '(Effective only if flow_vect) interval of data sampling in x'], ['yintv', 1, '(Effective only if flow_vect) interval of data sampling in y'], ['xfact', 1, 'scale factor for x'], ['yfact', 1, 'scale factor for y'], ['itr', 1, 'log or linear'], ['unit_vect', false, 'Show the unit vector'], ['max_unit_vect', false, '(Effective only if flow_vect && unit_vect) If true, use the maximum arrows to scale the unit vector; otherwise, normalize in V coordinate.'] ) end opts = @@vector_for_dcchart_options.interpret(options) fx = fx.first2D.copy fy = fy.first2D.copy sh = fx.shape if sh != fy.shape raise ArgumentError, "shapes of fx and fy do not agree with each other" end fx = fx.transpose(1,0) if opts['transpose'] fy = fy.transpose(1,0) if opts['transpose'] if ((xi=opts['xintv']) >= 2) idx = NArray.int(sh[0]/xi).indgen!*xi # [0,xi,2*xi,..] fx = fx[idx, true] fy = fy[idx, true] end if ((yi=opts['yintv']) >= 2) idx = NArray.int(sh[1]/yi).indgen!*yi # [0,yi,2*yi,..] fx = fx[true, idx] fy = fy[true, idx] end xax = fx.coord(0) yax = fy.coord(1) if newframe nextfig = @@next_fig.dup if ( @@next_fig != nil ) # backup next_fig fig(xax, yax, {'itr'=>itr}) axes(xax, yax) if opts['title'] ttl = opts['title'] else fxnm = fx.data.get_att('long_name') || fx.name fynm = fy.data.get_att('long_name') || fy.name ttl = '('+fxnm+','+fynm+')' end title( ttl ) annotate(fx.lost_axes) if opts['annotate'] @@next_fig = nextfig if ( @@next_fig != nil ) end fig(xax, yax, {'new_frame'=>false, 'itr'=>1, 'yreverse'=>false}) \ if (opts['newfig']) DCL.uwsgxa(xax.val) DCL.uwsgya(yax.val) xfact = opts["xfact"]; yfact = opts["yfact"] before2=DCLExt.ug_set_params( {'LNRMAL'=>false, 'LMSG'=>false, 'XFACT1'=>xfact, 'YFACT1'=>yfact} ) # before=DCLExt.ug_set_params({'lunit'=>true}) if opts['unit_vect'] DCL.ugvect(fx.val[true, -1..0], fy.val[true, -1..0]) DCLExt.unit_vect(xfact, yfact, nil, nil, "vxuoff"=>0.01, "vyuoff"=>-0.1) if opts['unit_vect'] # DCLExt.ug_set_params(before) if opts['unit_vect'] nil end end end def GGraph::line(gphys, newframe=true, options=nil) if newframe!=true && newframe!=false raise ArgumentError, "2nd arg (newframe) must be true or false" end if ! defined?(@@line_options) @@line_options = Misc::KeywordOptAutoHelp.new( ['title', nil, 'Title of the figure(if nil, internally determined)'], ['annotate', true, 'if false, do not put texts on the right margin even when newframe==true'], ['exchange', false, 'whether to exchange x and y axes'], ['index', 1, 'line/mark index'], ['type', 1, 'line type'], ['label', nil, 'if a String is given, it is shown as the label'], ['max', nil, 'maximam value for draw line'], ['min', nil, 'minimam value for draw line'] ) end opts = @@line_options.interpret(options) gp = gphys.first1D if !opts['exchange'] x = gp.coord(0) y = gp.data xax = x if opts['min'] ymin = opts['min'].to_f else ymin = gp.data.val.min end if opts['max'] ymax = opts['max'].to_f else ymax = gp.data.val.max end yax = VArray.new(NArray[ymin, ymax], gp.data, gp.data.name) else y = gp.coord(0) x = gp.data yax = y if opts['min'] xmin = opts['min'] else xmin = gp.data.val.min end if opts['max'] xmax = opts['max'] if opts['max'] else xmax = gp.data.val.max end xax = VArray.new(NArray[xmin, xmax], gp.data, gp.data.name) end if newframe fig(xax, yax) axes(xax, yax) title( (opts['title'] || gp.data.get_att('long_name')) ) annotate(gp.lost_axes) if opts['annotate'] end if opts['label'] lcharbk = DCL.sgpget('lchar') DCL.sgpset('lchar',true) DCL.sgsplc(opts['label']) end DCL.uulinz(x.val, y.val, opts['type'], opts['index']) DCL.sgpset('lchar',lcharbk) if opts['label'] nil end ##----------------------- # ラインメソッド再定義? def plot_line_main(gphys_array, line_opts={}) line_index_ary = Array.new; line_type_ary = Array.new; name_ary = Array.new default_index = 12 default_type = 1 line_hash = { "index"=> default_index }.update(line_opts) line_index_ary.push(default_index) line_type_ary.push(default_type) name_ary.push(gphys_array[0].data.get_att("long_name")) GGraph.line( gphys_array[0], true, line_hash ) gphys_array.size.times{ |num| unless num == 0 if num > 3 line_type = default_type + 2 index = ( (num - 3)*10 + default_index) else line_type = default_type index = (num*10 + default_index) end if type = gphys_array[num].data.get_att("line_type") line_type = type end line_hash = { "index"=>index, "type"=>line_type } line_index_ary.push(index); line_type_ary.push(line_type); name_ary.push(gphys_array[num].data.get_att("long_name")) GGraph.line( gphys_array[num], false, line_hash ) end } return name_ary, line_index_ary, line_type_ary end ##----------------------- # ラインインデックスの種類を表示 def __show_line_index(str_ary,line_index_ary, line_type_ary, x=0.15, y=0.1) charsize = 0.85 * DCL.uzpget('rsizec1') len = 0.05 dvx = 0.01 dvy = charsize*1.6 raise TypeError,"Array expected" if ! str_ary.is_a?(Array) vxmin,vxmax,vymin,vymax = DCL.sgqvpt vx = x vy = y + charsize/2 str_ary.size.times{|num| DCL::sgplzv([vx, vx + len],[vy, vy], line_type_ary[num], line_index_ary[num]) DCL::sgtxzv(vx + len + 0.02 ,vy,str_ary[num], charsize, 0, -1, line_index_ary[num]) vy -= dvy if num == 3 vx = 0.50 vy = 0.12 - charsize/2 end } nil end ##----------------------- ## 25 年平均のオブジェクトを作成 def make_mean_gphys(phys, src, season, year, var=phys) pathary = [] y1 = year[0]; y2 = year[1] path = "../../#{phys}.#{src}/" season.each do |m| y1.upto(y2) do |y| y = y-1 if m == "12" fn = path + "#{phys}.#{y}.#{src}/#{phys}_#{y}-#{m}_#{src}.nc" pathary << fn end end # gp = mean_gphys_date_weight(pathary, phys.downcase) gp = mean_gphys(pathary, var.downcase) return gp end ## GGraph の margin_info を参考に, 余白の上下につけるようにした. def margin_info(program=nil, data_source=nil, char_height=nil, date=nil, xl=0.0, xr=0.0, yb=nil, yt=nil) program = File.expand_path($0) if !program if date program = program + ' ' + Date.today.to_s else program = program if date.nil? program = program + ' ' + Date.today.to_s end end if !data_source data_source = Dir.pwd data_source = '' if data_source == program.sub(/\/[^\/]+$/,'') end sz = 1.0/( program.length + data_source.length ) char_height = [0.008, sz].min if !char_height yb = 2.0 * char_height if !yb yt = 2.0 * char_height if !yt DCL.slmgn(xl, xr, yb, yt) # 上につける DCL.slsttl(program, 't', -1.0, 1.0, char_height, 1) DCL.slsttl(data_source, 't', 1.0, 1.0, char_height, 2) # 下につける DCL.slsttl(program, 'b', -1.0, 1.0, char_height, 3) DCL.slsttl(data_source, 'b', 1.0, 1.0, char_height, 4) end ################################################################## unless getopts( "", "info", # 余白情報をつけるか否か "climat", # 気候値読み込み "deviation", # 年平均季節変動からの偏差 "kikan:", # "annual", "winter", "spring", "summer", "fall", 月を指定 "range:" # 平均する年を指定. とある一年の平均の場合 ) print "#{$0}:illegal options. please exec this program with -h/--help. \n" exit 1 end ## define Constants Infty = 1.0/0.0 L = UNumeric.new(2.5008e6, 'J.Kg-1') Rad = UNumeric.new(6.37E6, 'm') G = UNumeric.new(9.81 , 'm.s-2') # make gphys objects ## set kikan year = $OPT_range.split(/:/).collect{|y| y.to_i } if $OPT_range ############################################### if $OPT_climat ## set kikan winter = "DJF"; spring = "MAM" summer = "JJA"; fall = "SON" annual = "ANNUAL" case $OPT_kikan when "annual" kikan = annual when "winter" kikan = winter when "spring" kikan = spring when "summer" kikan = summer when "fall" kikan = fall else kikan = $OPT_kikan.to_s end ## 質量流線関数 strm = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "strm") ## 東西風速 uwnd = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "uwnd") ## 温度場 temp = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "t") ## 比湿 shum = GPhys::IO.open("300hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "shum") ## 放射 nswrf = GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan}_NCEP.nc", "uswrf") ulwrf = GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan}_NCEP.nc", "ulwrf") nlwrs = GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan}_NCEP.nc", "nlwrs") nswrs = GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan}_NCEP.nc", "nswrs") nsw = (-nswrf + nswrs).abs # 大気の吸収量(短波放射) nsw.data.set_att('long_name', 'RHEAT (long)') nlw = (-ulwrf + nlwrs).abs # 大気の吸収量(長波放射) nlw.data.set_att('long_name', 'RCOOL (short)') nrw = (nlw - nsw).abs # 大気の吸収量(トータル) nrw.data.set_att('long_name', 'RCOOL (total)') nrw.data.set_att('line_type', 3) radiation = [ nswrf.abs, ulwrf, nswrs.abs, nlwrs, nsw, nlw, nrw ] ## 熱 #### 潜熱フラックス lhtfl = GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan}_NCEP.nc", "lhtfl") #### 顕熱フラックス shtfl = GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan}_NCEP.nc", "shtfl") #### 降水(潜熱に換算したもの) lheat = GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan}_NCEP.nc", "prate") heat = [ lhtfl, shtfl, lheat, nrw ] ## 乾燥静的エネルギー dse = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "dsefl") dse_bar = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "dsefl_bar")*1 dse_dash = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "dsefl_dash")*1 dse_bar.data.set_att('long_name', 'dry static energy transport (mean)') dse_dash.data.set_att('long_name', 'dry static energy transport (eddy)') ## 潜熱輸送量 lhe = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "lhefl") lhe_bar = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "lhefl_bar")*1 lhe_dash = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "lhefl_dash")*1 lhe_bar.data.set_att('long_name', 'latent heat energy transport (mean)') lhe_dash.data.set_att('long_name', 'latent heat energy transport (eddy)') # 湿潤静的エネルギー mse = lhe + dse mse.data.set_att('long_name', 'moist static energy transport (total)') dry_ene = [ mse, dse, dse_bar, dse_dash ] lat_ene = [ mse, lhe, lhe_bar, lhe_dash ] ## EP フラックス epflx_phi = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "epflx_phi") epflx_p = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "epflx_p") epflx_div = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "epflx_div") ## 残差循環 strm_rmean = GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan}_NCEP.nc", "strm_rmean") mgn_str = kikan+" MEAN (1979-2003)" elsif $OPT_deviation winter = ["12", "01", "02"]; spring = ["03", "04", "05"] summer = ["06", "07", "08"]; fall = ["09", "10", "11"] annual = winter + spring + summer + fall case $OPT_kikan when "annual" kikan_ary = annual kikan_str = "ANNUAL" when "winter" kikan_ary = winter kikan_str = "DJF" when "spring" kikan_ary = spring kikan_str = "MAM" when "summer" kikan_ary = summer kikan_str = "JJA" when "fall" kikan_ary = fall kikan_str = "SON" else kikan_ary = $OPT_kikan kikan_str = $OPT_kikan.to_s end ## 質量流線関数 strm = make_mean_gphys('STRM', 'NCEP', kikan_ary, year).cut( 'level'=>1000..100 ) \ - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "strm") ## 東西風速 uwnd = make_mean_gphys('UWND', 'NCEP', kikan_ary, year).mean('lon').cut( 'level'=>1000..100 ) \ - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "uwnd") ## 温度場 temp = make_mean_gphys('T', 'NCEP', kikan_ary, year).mean('lon').cut( 'level'=>1000..100 ) \ - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "t") ## 比湿 shum = make_mean_gphys('SHUM', 'NCEP', kikan_ary, year).mean('lon').cut( 'level'=>1000..100 ) \ - GPhys::IO.open("300hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "shum") ## 放射 dswrf = make_mean_gphys('DSWRF', 'NCEP', kikan_ary, year).mean('lon') uswrf = make_mean_gphys('USWRF', 'NCEP', kikan_ary, year).mean('lon') nswrf = uswrf - dswrf # 正味の短波放射フラックス nswrf = nswrf - GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "uswrf") nswrf.data.set_att('long_name', 'NSR at TOP') ulwrf = make_mean_gphys('ULWRF', 'NCEP', kikan_ary, year).mean('lon') \ - GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "ulwrf") ulwrf.data.set_att('long_name', 'OLR at TOP') nlwrs = make_mean_gphys('NLWRS', 'NCEP', kikan_ary, year).mean('lon') \ - GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "nlwrs") nlwrs.data.set_att('long_name', 'NLR at SFC') nswrs = make_mean_gphys('NSWRS', 'NCEP', kikan_ary, year).mean('lon') \ - GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "nswrs") nswrs.data.set_att('long_name', 'NSR at SFC') nsw = (-nswrf + nswrs).abs # 大気の吸収量(短波放射) nsw.data.set_att('long_name', 'RHEAT (long)') nlw = (-ulwrf + nlwrs).abs # 大気の吸収量(長波放射) nlw.data.set_att('long_name', 'RCOOL (short)') nrw = (nlw - nsw).abs # 大気の吸収量(トータル) nrw.data.set_att('long_name', 'RCOOL (total)') nrw.data.set_att('line_type', 3) radiation = [ nswrf.abs, ulwrf, nswrs.abs, nlwrs, nsw, nlw, nrw ] ## 熱 #### 潜熱フラックス lhtfl = make_mean_gphys('LHTFL', 'NCEP', kikan_ary, year).mean('lon') \ - GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "lhtfl") lhtfl.data.set_att('long_name', 'net latent heat flux at surface') #### 顕熱フラックス shtfl = make_mean_gphys('SHTFL', 'NCEP', kikan_ary, year).mean('lon') \ - GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "shtfl") shtfl.data.set_att('long_name', 'net sensible heat flux at surface') #### 降水(潜熱に換算したもの) prate = make_mean_gphys('PRATE', 'NCEP', kikan_ary, year).mean('lon') lheat = prate * L lheat = lheat - GPhys::IO.open("GAUSSIAN_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "prate") lheat.data.set_att('long_name', 'latent heat with precipitation') heat = [ lhtfl, shtfl, lheat, nrw ] ## 乾燥静的エネルギー dsefl = make_mean_gphys('DSEFL', 'NCEP', kikan_ary, year, 'dsefl') dsefl_bar = make_mean_gphys('DSEFL', 'NCEP', kikan_ary, year, 'dsefl_bar') dsefl_dash = make_mean_gphys('DSEFL', 'NCEP', kikan_ary, year, 'dsefl_dash') dse_bar.data.set_att('long_name', 'dry static energy transport (mean)') dse_dash.data.set_att('long_name', 'dry static energy transport (eddy)') ## 潜熱輸送量 lhefl = make_mean_gphys('LHEFL', 'NCEP', kikan_ary, year, 'lhefl') lhefl_bar = make_mean_gphys('LHEFL', 'NCEP', kikan_ary, year, 'lhefl_bar') lhefl_dash = make_mean_gphys('LHEFL', 'NCEP', kikan_ary, year, 'lhefl_dash') lhe_bar.data.set_att('long_name', 'latent heat energy transport (mean)') lhe_dash.data.set_att('long_name', 'latent heat energy transport (eddy)') # get cos lat = lhefl.coord('lat') grid = Grid.new(lhefl.axis('lat')) cos_lat = GPhys.new(grid, cos(lat*2*PI/360.0)) dse = dsefl.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 dse_bar = dsefl_bar.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 dse_dash = dsefl_dash.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 dse.data.set_att('units', 'W') dse_bar.data.set_att('units', 'W') dse_bar.data.set_att('long_name', 'dry static energy transport (mean)') dse_dash.data.set_att('units', 'W') dse_bar.data.set_att('long_name', 'dry static energy transport (eddy)') dse = dse - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "dsefl") dse_bar = dse - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "dsefl_bar") dse_dash = dse - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "dsefl_dash") dse_bar.data.set_att('long_name', 'dry static energy transport (mean)') dse_dash.data.set_att('long_name', 'dry static energy transport (eddy)') lhe = lhefl.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 lhe_bar = lhefl_bar.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 lhe_dash = lhefl_dash.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 lhe.data.set_att('units', 'W') lhe.data.set_att('long_name', 'latent heat energy transport (total)') lhe_bar.data.set_att('units', 'W') lhe_bar.data.set_att('long_name', 'latent heat energy transport (mean)') lhe_dash.data.set_att('units', 'W') lhe_dash.data.set_att('long_name', 'latent heat energy transport (eddy)') lhe = lhe - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "lhefl") lhe_bar = lhe - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "lhefl_bar") lhe_dash = lhe - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "lhefl_dash") lhe_bar.data.set_att('long_name', 'latent heat energy transport (mean)') lhe_dash.data.set_att('long_name', 'latent heat energy transport (eddy)') # 湿潤静的エネルギー mse = lhe + dse mse.data.set_att('long_name', 'moist static energy transport (total)') dry_ene = [ mse, dse, dse_bar, dse_dash ] lat_ene = [ mse, lhe, lhe_bar, lhe_dash ] ## EP フラックス epflx_phi = make_mean_gphys('EPFLX', 'NCEP', kikan_ary, year, 'epflx_phi') \ - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "epflx_phi") epflx_p = make_mean_gphys('EPFLX', 'NCEP', kikan_ary, year, 'epflx_p') \ - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "epflx_p") epflx_div = make_mean_gphys('EPFLX', 'NCEP', kikan_ary, year, 'epflx_div') \ - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "epflx_div") ## 残差循環 strm_rmean = make_mean_gphys('STRM_RMEAN', 'NCEP', kikan_ary, year) \ - GPhys::IO.open("100hPa_PHYS_1979-2003_#{kikan_str}_NCEP.nc", "strm_rmean") ############################################################################# # マージンに載せる年月情報 mgn_kikan = { 'annual' => 'ANNUAL', 'winter' => 'DJF', 'spring' => 'MAM', 'summer' => 'JJA', 'fall' => 'SON', '01' => 'JAN', '02' => 'FEB', '03' => 'MAR', '04' => 'APR', '05' => 'MAY', '06' => 'JUN', '07' => 'JUL', '08' => 'AUG', '09' => 'SEP', '10' => 'OCT', '11' => 'NOV', '12' => 'DEC' } mgn_str = mgn_kikan[$OPT_kikan]+" MEAN (#{year.uniq})" else ############################################### winter = ["12", "01", "02"]; spring = ["03", "04", "05"] summer = ["06", "07", "08"]; fall = ["09", "10", "11"] annual = winter + spring + summer + fall case $OPT_kikan when "annual" kikan = annual when "winter" kikan = winter when "spring" kikan = spring when "summer" kikan = summer when "fall" kikan = fall else kikan = $OPT_kikan end ## 質量流線関数 strm = make_mean_gphys('STRM', 'NCEP', kikan, year).cut( 'level'=>1000..100 ) ## 東西風速 uwnd = make_mean_gphys('UWND', 'NCEP', kikan, year).mean('lon').cut( 'level'=>1000..100 ) ## 温度場 temp = make_mean_gphys('T', 'NCEP', kikan, year).mean('lon').cut( 'level'=>1000..100 ) ## 比湿 shum = make_mean_gphys('SHUM', 'NCEP', kikan, year).mean('lon').cut( 'level'=>1000..100 ) ## 放射 dswrf = make_mean_gphys('DSWRF', 'NCEP', kikan, year).mean('lon') uswrf = make_mean_gphys('USWRF', 'NCEP', kikan, year).mean('lon') nswrf = uswrf - dswrf # 正味の短波放射フラックス nswrf.data.set_att('long_name', 'NSR at TOP') ulwrf = make_mean_gphys('ULWRF', 'NCEP', kikan, year).mean('lon') ulwrf.data.set_att('long_name', 'OLR at TOP') nlwrs = make_mean_gphys('NLWRS', 'NCEP', kikan, year).mean('lon') # 正味の長波放射フラックス nlwrs.data.set_att('long_name', 'NLR at SFC') nswrs = make_mean_gphys('NSWRS', 'NCEP', kikan, year).mean('lon') # 地表面正味の短波放射フラックス nswrs.data.set_att('long_name', 'NSR at SFC') nsw = (-nswrf + nswrs).abs # 大気の吸収量(短波放射) nsw.data.set_att('long_name', 'RHEAT (long)') nlw = (-ulwrf + nlwrs).abs # 大気の吸収量(長波放射) nlw.data.set_att('long_name', 'RCOOL (short)') nrw = (nlw - nsw).abs # 大気の吸収量(トータル) nrw.data.set_att('long_name', 'RCOOL (total)') nrw.data.set_att('line_type', 3) radiation = [ nswrf.abs, ulwrf, nswrs.abs, nlwrs, nsw, nlw, nrw ] ## 熱 #### 潜熱フラックス lhtfl = make_mean_gphys('LHTFL', 'NCEP', kikan, year).mean('lon') lhtfl.data.set_att('long_name', 'net latent heat flux at surface') #### 顕熱フラックス shtfl = make_mean_gphys('SHTFL', 'NCEP', kikan, year).mean('lon') shtfl.data.set_att('long_name', 'net sensible heat flux at surface') #### 降水(潜熱に換算したもの) prate = make_mean_gphys('PRATE', 'NCEP', kikan, year).mean('lon') lheat = prate*L lheat.data.set_att('long_name', 'latent heat with precipitation') heat = [ lhtfl, shtfl, lheat, nrw ] ## 乾燥静的エネルギー dsefl = make_mean_gphys('DSEFL', 'NCEP', kikan, year, 'dsefl') dsefl_bar = make_mean_gphys('DSEFL', 'NCEP', kikan, year, 'dsefl_bar') dsefl_dash = make_mean_gphys('DSEFL', 'NCEP', kikan, year, 'dsefl_dash') ## 潜熱輸送量 lhefl = make_mean_gphys('LHEFL', 'NCEP', kikan, year, 'lhefl') lhefl_bar = make_mean_gphys('LHEFL', 'NCEP', kikan, year, 'lhefl_bar') lhefl_dash = make_mean_gphys('LHEFL', 'NCEP', kikan, year, 'lhefl_dash') # get cos lat = lhefl.coord('lat') grid = Grid.new(lhefl.axis('lat')) cos_lat = GPhys.new(grid, cos(lat*2*PI/360.0)) dse = dsefl.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 dse_bar = dsefl_bar.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 dse_dash = dsefl_dash.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 dse.data.set_att('units', 'W') dse_bar.data.set_att('units', 'W') dse_bar.data.set_att('long_name', 'dry static energy transport (mean)') dse_dash.data.set_att('units', 'W') dse_bar.data.set_att('long_name', 'dry static energy transport (eddy)') lhe = lhefl.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 lhe_bar = lhefl_bar.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 lhe_dash = lhefl_dash.integrate('level') /G * Rad * cos_lat * 2 * PI * 100 lhe.data.set_att('units', 'W') lhe.data.set_att('long_name', 'latent heat energy transport (total)') lhe_bar.data.set_att('units', 'W') lhe_bar.data.set_att('long_name', 'latent heat energy transport (mean)') lhe_dash.data.set_att('units', 'W') lhe_dash.data.set_att('long_name', 'latent heat energy transport (eddy)') # 湿潤静的エネルギー mse = lhe + dse mse.data.set_att('long_name', 'moist static energy transport (total)') dry_ene = [ mse, dse, dse_bar, dse_dash ] lat_ene = [ mse, lhe, lhe_bar, lhe_dash ] ## EP フラックス epflx_phi = make_mean_gphys('EPFLX', 'NCEP', kikan, year, 'epflx_phi') epflx_p = make_mean_gphys('EPFLX', 'NCEP', kikan, year, 'epflx_p') epflx_div = make_mean_gphys('EPFLX', 'NCEP', kikan, year, 'epflx_div') ## 残差循環 strm_rmean = make_mean_gphys('STRM_RMEAN', 'NCEP', kikan, year) ############################################################################# # マージンに載せる年月情報 mgn_kikan = { 'annual' => 'ANNUAL', 'winter' => 'DJF', 'spring' => 'MAM', 'summer' => 'JJA', 'fall' => 'SON', '01' => 'JAN', '02' => 'FEB', '03' => 'MAR', '04' => 'APR', '05' => 'MAY', '06' => 'JUN', '07' => 'JUL', '08' => 'AUG', '09' => 'SEP', '10' => 'OCT', '11' => 'NOV', '12' => 'DEC' } mgn_str = mgn_kikan[$OPT_kikan]+" MEAN (#{year.uniq})" end ############################################################################## # set User Path for dcldatabase DCL.glcset('DUPATH','/home/daktu32/.dcldir/') DCL.swpset('lsep', true) # ページ毎別々のファイルに落す DCL.gropn(-2) margin_info($0, mgn_str) if $OPT_info DCL.sgpset('lcntl', false) # DCL.uzfact(0.75) # DCL.sgpset('lcorner',false) # コーナーを取っ払う DCL.sgpset('lfprop',true) # # DCL.udpset('lmsg',false) # DCL.uscset('cyspos', 'B' ) # y 軸の単位の位置を下方へ DCL.sgpset('lfull',true) # フルスクリーン DCL.sldiv('y',2,5) ############################################################################## ##-------------------------------------- ## HEAT BALANCE (seasonal) # ビューポート設定 vpt00 = NArray[0.12, 0.88, 0.14, 0.5] vpt01 = NArray[0.11, 0.87, 0.14, 0.5] vpt10 = vpt00 vpt11 = vpt01 vpt20 = vpt00 vpt21 = vpt01 vpt30 = vpt00 vpt31 = vpt01 vpt40 = vpt00 vpt41 = vpt01 itr = 1 # 放射(左下の絵) ## draw GGraph.set_fig('viewport'=>vpt20, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('xunits'=>'', 'yunits'=>'', 'xtitle'=>'', 'ytitle'=>'') DCL.uzpset('labelxb',false) names, idxs, types = plot_line_main(radiation, 'annot'=>false, 'titl'=>'', 'min'=>-100, 'max'=>400) __show_line_index(names,idxs, types) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t',"(a) Radiation Budget(#{radiation[0].data.units})", -1) ; DCL.uzpset('pad1',0.7) # 熱(右下の絵) GGraph.set_fig('viewport'=>vpt21, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('xunits'=>'', 'yunits'=>'', 'xtitle'=>'', 'ytitle'=>'') names, idxs, types = plot_line_main(heat, 'annot'=>false, 'titl'=>'', 'min'=>-100, 'max'=>400) __show_line_index(names,idxs, types) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t',"(b) Heat Budget(#{heat[0].data.units})", -1) ; DCL.uzpset('pad1',0.7) # 乾燥静的エネルギーフラックス(左下の絵) GGraph.set_fig('viewport'=>vpt30, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('ytitle'=>'','xtitle'=>nil) names, idxs, types = plot_line_main(dry_ene, 'annot'=>false, 'titl'=>'', 'min'=>-8.0e15, 'max'=>+8.0e15) __show_line_index(names,idxs, types) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','(c) Dry Static Energy Transport(W)',-1) ; DCL.uzpset('pad1',0.7) # 南北エネルギーフラックス(左下の絵) GGraph.set_fig('viewport'=>vpt31, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('ytitle'=>'') names, idxs, types = plot_line_main(lat_ene, 'annot'=>false, 'titl'=>'', 'min'=>-8.0e15, 'max'=>+8.0e15) __show_line_index(names,idxs, types) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','(d) Latent Heat Transport(W)',-1) ; DCL.uzpset('pad1',0.7) itr = 2 # dcl に欠損値情報を教える before = DCLExt.gl_set_params('lmiss'=>true,'rmiss'=>temp.data.get_att("missing_value")[0]) # 温度場(左下の絵) ## set param GGraph.next_linear_contour_options( 'min'=>180, 'max'=>330, 'int'=>10) level = NArray.int(8).indgen!(0, 20)+180 level[0] = 150 pattern = NArray[25999, 30999, 40999, 55999, 70999, 75999, 85999] ## draw GGraph.set_fig('viewport'=>vpt10, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('ytitle'=>'','xtitle'=>nil) GGraph.tone( temp, true , 'annot'=>false, 'titl'=>'', 'pat'=>pattern, 'lev'=>level ) GGraph.contour( temp, false , 'annot'=>false, 'titl'=>'' ) GGraph::color_bar('constwidth'=>true) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t',"(e) Temperature(#{temp.data.units})", -1) ; DCL.uzpset('pad1',0.7) # 比湿(右下の絵) ## set param GGraph.next_linear_contour_options('int'=>0.001) level = NArray[0, 0.002, 0.004, 0.006, 0.008, 0.01, 0.012, 0.014, 0.016, 0.018, 0.02] pattern = NArray[15999, 25999, 30999, 35999, 40999, 50999, 70999, 75999, 80999, 85999] ## draw GGraph.set_fig('viewport'=>vpt11, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('ytitle'=>'') GGraph.tone( shum, true , 'annot'=>false, 'titl'=>'', 'pat'=>pattern, 'lev'=>level ) GGraph.contour( shum, false , 'annot'=>false, 'titl'=>'' ) GGraph::color_bar('constwidth'=>true) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t',"(f) Specific Humidity(#{shum.data.units})", -1) ; DCL.uzpset('pad1',0.7) # 東西風速(右上の絵) ## set param GGraph.next_linear_contour_options('min'=>-50,'int'=>2.5, 'max'=>70) level = [-15, -5, 0, 5, 15, 25] pattern = [30999, 35999, 40999, 55999, 70999, 75999, 85999] ## draw GGraph.set_fig('viewport'=>vpt01, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('ytitle'=>'') GGraph.tone( uwnd, true , 'annot'=>false, 'titl'=>'', 'pat'=>pattern, 'lev'=>level ) GGraph.contour( uwnd, false , 'annot'=>false, 'titl'=>'' ) GGraph::color_bar('constwidth'=>true) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t',"(g) Zonal Wind(#{uwnd.data.units})", -1) ; DCL.uzpset('pad1',0.7) # 質量流線関数(左上の絵) level = [-1.0e12, 0, +1.0e12]; pattern = [35999, 35999, 70999, 70999] GGraph.next_linear_contour_options('int'=>1.0e10) ## draw GGraph.set_fig('viewport'=>vpt00, 'itr'=>itr) GGraph.set_axes('xunits'=>'', 'yunits'=>'', 'xtitle'=>'', 'ytitle'=>'') DCL.uzpset('labelxb',false) GGraph.tone( strm, true , 'annot'=>false, 'titl'=>'', 'pat'=>pattern, 'lev'=>level ) GGraph.contour( strm, false , 'annot'=>false, 'titl'=>'' ) GGraph::color_bar('constwidth'=>true) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t',"(h) Mass Stream Function(#{strm.data.units})", -1) ; DCL.uzpset('pad1',0.7) # EP フラックス(左下の絵) ## settings level = NArray[-5.0e1, -1.0e1, 0, 1.0e1, 5.0e1]*1.0e-6 pattern = NArray[30999, 35999, 40999, 55999, 75999, 85999] vect_ratio = NArray[1.0, 100.0]*8e-12 vect_fact = 10.0 ## draw GGraph.set_fig('viewport'=>vpt40, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('ytitle'=>'','xtitle'=>nil) DCL.uzpset('labelxb',true) GGraph.tone( epflx_div, true , 'annot'=>false, 'titl'=>'', 'pat'=>pattern, 'lev'=>level ) GGraph::vector_for_dcchart(epflx_phi, - epflx_p, false, 'xfact'=>vect_ratio[0]*vect_fact, 'yfact'=>vect_ratio[-1]*vect_fact, 'unit'=>true, 'annot'=>false, 'itr'=>2, 'xintv'=>2 ) GGraph::color_bar('constwidth'=>true) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','(i) EP Flux(kg.s-2), Div F(m.s-2)',-1) ; DCL.uzpset('pad1',0.7) # 残差循環(左下の絵) ## settings level = [-1.0e12, 0, +1.0e12]; pattern = [35999, 35999, 70999, 70999] GGraph.next_linear_contour_options('int'=>1.0e10) ## draw GGraph.set_fig('viewport'=>vpt41, 'new_frame'=>true, 'itr'=>itr) GGraph.set_axes('ytitle'=>'') GGraph.tone( strm_rmean, true , 'annot'=>false, 'titl'=>'', 'pat'=>pattern, 'lev'=>level ) GGraph::color_bar('constwidth'=>true) GGraph.contour( strm_rmean, false , 'annot'=>false, 'titl'=>'' ) DCL.uzpset('pad1',0.2) ; DCL.uxsttl('t','(j) Residual Mass Stream Function(kg/s)',-1) ; DCL.uzpset('pad1',0.7) DCL.grcls #######################################################