#!/usr/bin/env ruby =begin = mkfig.rb -- make figure for energy balance エネルギーバランス関連図. 以下の物理量の緯度平均値をプロット. * RADIATION_BALANCE * DSWRF : 大気上端下向き短波放射フラックス * USWRF : 大気上端上向き短波放射フラックス * ULWRF : 大気上端上向き長波放射フラックス * NSWRS : 地表面正味短波放射フラックス * NLWRS : 地表面正味長波放射フラックス * HEAT_BALANCE * PRATE : 降水強度(潜熱の形に書き直す) * LHTFL : 地表面における潜熱フラックス * SHTFL : 地表面における顕熱フラックス =end require "getopts" # for option_parse require "numru/netcdf_miss" require "numru/ggraph" require "colorbar" require "libgphys-n" include NumRu include Misc::EMath ## enable to the option min, max 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; name_ary = Array.new default_index = 12 line_hash = { "index"=> default_index }.update(line_opts) line_index_ary.push(default_index) 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 line_hash = { "index"=> (num*10 + default_index)} line_index_ary.push(num*10 +default_index) 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 end ##----------------------- # ラインインデックスの種類を表示 def __show_line_index(str_ary,line_index_ary, x=0.15, y=0.12) charsize = 1.0 * 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 = x vy = y + 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 == 6 vx = 0.30 vy = 0.12 - charsize/2 end } nil end ##----------------------- ## 25 年平均のオブジェクトを作成 def make_mean_gphys(phys, src, season,year) pathary = [] y1 = year[0]; y2 = year[1] path = "../../#{phys}.#{src}/" season.each do |m| y1.upto(y2) do |y| 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, phys.downcase) return gp.average('lon') end ##----------------------- # 全球平均値求めてテキストに書き込む (TeX 形式). # 各配列は先頭から OLR_top(トップ上向き放射), NLR_sfc(地表面ネット長波), LH(潜熱), SH(顕熱), PR(降水) def write_tablar(file, ann, djf, mam, jja, son) ## -- 全球平均値を求める -- ann = cal_global_mean(ann) djf = cal_global_mean(djf) mam = cal_global_mean(mam) jja = cal_global_mean(jja) son = cal_global_mean(son) ## -- 結果を書き込む -- fmt = "& %10.1f & %10.1f & %10.1f & %10.1f & %10.1f & %10.1f & %10.1f" file.puts <winter, 'spring'=>spring, 'summer'=>summer, 'fall' =>fall } season_gp = {# 各季節平均の GPhys オブジェクトを持つハッシュ 'winter'=>djf, 'spring'=>mam, 'summer'=>jja, 'fall' =>son } ['spring', 'summer', 'fall', 'winter'].each {|ss| ## 大気上端放射フラックス dswrf = make_mean_gphys('DSWRF', 'NCEP', season_kikan[ss], year) uswrf = make_mean_gphys('USWRF', 'NCEP', season_kikan[ss], year) nswrf = uswrf - dswrf # 正味の短波放射フラックス nswrf.data.set_att('long_name', 'net short wave radiation flux at nominal top of atmosphere') ulwrf = make_mean_gphys('ULWRF', 'NCEP', season_kikan[ss], year) # 正味の長波放射フラックス ## 地表面放射フラックス nlwrs = make_mean_gphys('NLWRS', 'NCEP', season_kikan[ss], year) # 地表面正味の長波放射フラックス nswrs = make_mean_gphys('NSWRS', 'NCEP', season_kikan[ss], year) # 地表面正味の短波放射フラックス ## 潜熱フラックス lhtfl = make_mean_gphys('LHTFL', 'NCEP', season_kikan[ss], year) lhtfl.data.set_att('long_name', 'net latent heat flux at surface') ## 顕熱フラックス shtfl = make_mean_gphys('SHTFL', 'NCEP', season_kikan[ss], year) shtfl.data.set_att('long_name', 'net sensible heat flux at surface') ## 降水(潜熱に換算したもの) prate = make_mean_gphys('PRATE', 'NCEP', season_kikan[ss], year) lheat = prate*L lheat.data.set_att('long_name', 'latent heat with precipitation') season_gp[ss] << ulwrf season_gp[ss] << nlwrs season_gp[ss] << nswrf season_gp[ss] << nswrs season_gp[ss] << lhtfl season_gp[ss] << shtfl season_gp[ss] << lheat } ann = (NArray.to_na(djf) + NArray.to_na(mam) + NArray.to_na(jja) + NArray.to_na(son))/4 ################################################################### # set User Path for dcldatabase DCL.glcset('DUPATH','/home/daktu32/.dcldir/') itr = 1 DCL.swpset('lsep', true) # ページ毎別々のファイルに落す DCL.gropn(2) DCL.sgpset('lcntl', false) # DCL.sgpset('lfull',true) # フルスクリーン # DCL.uzfact(0.35) # DCL.uzfact(0.55) # DCL.sgpset('lcorner',false) # コーナーを取っ払う DCL.sgpset('lfprop',true) # DCL.udpset('lmsg',false) # DCL.uscset('cyspos', 'B' ) # y 軸の単位の位置を下方へ ############################################################################## ##--------------------------------------------------- ## HEAT BALANCE (annual) # ビューポート設定 vpt = NArray[0.15,0.85,0.23,0.58] # dcl に欠損値情報を教える before = DCLExt.gl_set_params('lmiss'=>true,'rmiss'=>ann[0].data.get_att("missing_value")[0]) # annual mean GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) GGraph.set_axes('ytitle'=>'') DCL.uzpset('labelyl',true) names, idxs = plot_line_main(ann[-3..-1], 'annot'=>false, 'titl'=>'25 years annual mean Heat Balance (1979-2003)', 'min'=>min, 'max'=>max) __show_line_index(names,idxs) ## -- extra information -- # DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years annual mean Heat Balance (1979-2003)", # 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{ann[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ##-------------------------------------- ## HEAT BALANCE (seasonal) ## ------- DJF -------- GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) names, idxs = plot_line_main(djf[-3..-1], 'annot'=>false, 'titl'=>'25 years DJF mean Heat Balance (1979-2003)', 'min'=>min, 'max'=>max) __show_line_index(names,idxs) # DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years DJF mean Heat Balance (1979-2003)", # 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{djf[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ## ------- MAM -------- GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) names, idxs = plot_line_main(mam[-3..-1], 'annot'=>false, 'titl'=>'25 years MAM mean Heat Balance (1979-2003)', 'min'=>min, 'max'=>max) __show_line_index(names,idxs) # DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years MAM mean Heat Balance (1979-2003)", # 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{mam[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ## ------- JJA -------- GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) names, idxs = plot_line_main(jja[-3..-1], 'annot'=>false, 'titl'=>'', 'min'=>min, 'max'=>max) __show_line_index(names,idxs) DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years JJA mean Heat Balance (1979-2003)", 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{jja[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ## ------- SON -------- GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) names, idxs = plot_line_main(son[-3..-1], 'annot'=>false, 'titl'=>'', 'min'=>min, 'max'=>max) __show_line_index(names,idxs) DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years SON mean Heat Balance (1979-2003)", 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{son[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ############################################################################## ##--------------------------------------------------- ## RADIATION BALANCE (annual) # ビューポート設定 vpt = NArray[0.15,0.85,0.23,0.58] # dcl に欠損値情報を教える before = DCLExt.gl_set_params('lmiss'=>true,'rmiss'=>ann[0].data.get_att("missing_value")[0]) # annual mean GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) GGraph.set_axes('ytitle'=>'') DCL.uzpset('labelyl',true) names, idxs = plot_line_main(ann[0..3], 'annot'=>false, 'titl'=>'25 years annual mean Radiation Balance (1979-2003)', 'min'=>rad_min, 'max'=>rad_max) __show_line_index(names,idxs) ## -- extra information -- # DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years annual mean Radiation Balance (1979-2003)", # 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{ann[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ##-------------------------------------- ## RADIATION BALANCE (seasonal) ## ------- DJF -------- GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) names, idxs = plot_line_main(djf[0..3], 'annot'=>false, 'titl'=>'25 years DJF mean Radiation Balance (1979-2003)', 'min'=>rad_min, 'max'=>rad_max) __show_line_index(names,idxs) # DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years DJF mean Radiation Balance (1979-2003)", # 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{djf[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ## ------- MAM -------- GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) names, idxs = plot_line_main(mam[0..3], 'annot'=>false, 'titl'=>'25 years MAM mean Radiation Balance (1979-2003)', 'min'=>rad_min, 'max'=>rad_max) __show_line_index(names,idxs) # DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years MAM mean Radiation Balance (1979-2003)", # 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{mam[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ## ------- JJA -------- GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) names, idxs = plot_line_main(jja[0..3], 'annot'=>false, 'titl'=>'25 years JJA mean Radiation Balance (1979-2003)', 'min'=>rad_min, 'max'=>rad_max) __show_line_index(names,idxs) # DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years JJA mean Radiation Balance (1979-2003)", # 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{jja[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) ## ------- SON -------- GGraph.set_fig('viewport'=>vpt, 'itr'=>itr, 'new_frame'=>true) names, idxs = plot_line_main(son[0..3], 'annot'=>false, 'titl'=>'25 years SON mean Radiation Balance (1979-2003)', 'min'=>rad_min, 'max'=>rad_max) __show_line_index(names,idxs) # DCL::sgtxzv(0.5,vpt[3]+0.028,"25 years SON mean Radiation Balance (1979-2003)", # 1.15*DCL.uzpget('rsizec2'),0,0,3) DCL::sgtxzv(0.78,vpt[3]+0.01,"(units:#{son[0].data.units.to_s})", 0.8*DCL.uzpget('rsizec2'),0,0,3) DCL.grcls ####################################################### ################################################## # 全球平均値をファイルへ書き込む file = File.open('ENERGY_BALANCE_1973-2003_GLOBAL-MEAN_NCEP.txt', 'w') file = write_tablar(file, ann, djf, mam, jja, son) file.close ################################################## # 画像ファイル名を変更 File.rename('dcl_001.ps', 'HEAT_BALANCE_1973-2003-ANNUAL_NCEP.ps') File.rename('dcl_002.ps', 'HEAT_BALANCE_1973-2003-DJF_NCEP.ps') File.rename('dcl_003.ps', 'HEAT_BALANCE_1973-2003-MAM_NCEP.ps') File.rename('dcl_004.ps', 'HEAT_BALANCE_1973-2003-JJA_NCEP.ps') File.rename('dcl_005.ps', 'HEAT_BALANCE_1973-2003-SON_NCEP.ps') File.rename('dcl_006.ps', 'RADIATION_BALANCE_1973-2003-ANNUAL_NCEP.ps') File.rename('dcl_007.ps', 'RADIATION_BALANCE_1973-2003-DJF_NCEP.ps') File.rename('dcl_008.ps', 'RADIATION_BALANCE_1973-2003-MAM_NCEP.ps') File.rename('dcl_009.ps', 'RADIATION_BALANCE_1973-2003-JJA_NCEP.ps') File.rename('dcl_010.ps', 'RADIATION_BALANCE_1973-2003-SON_NCEP.ps') exit