#!/usr/bin/ruby # -*- coding: utf-8 -*- # =begin # = Nakajima et al. (1992) の Tg vs OLR の(ような)図を描く # # == 概要 # # == 更新履歴 # # * 2011/04/25 納多 哲史 新規作成 # # =end require "numru/ggraph" require "pp" include NumRu SCATTER_MAX = 3 SCATTER_MIN = 0 TMPFILE = "dcl.ps" epsfn = "OLR-viewpoint.eps" RPlanet = 6.371 * 10 ** 6 # [m] TITLE = 'OLR' YNAME = 'OLR' #YUNIT = '10^15 W' longname_rh='relative humidity' name_rh= 'RH' units_rh= '1' P_0 = 1.4 * 1.0e11 # NHA92 で用いる基準気圧 [Pa] Grav = 9.8 # 重力加速度 [m s-2] GasRUniv = 8.314 # 普遍気体定数 [J K-1 mol-1] MolWtWet = 18.01528e-3 # 凝結成分の平均分子量 [kg mol-1] MolWtDry = 28.964e-3 # 乾燥成分の平均分子量 [kg mol-1] R_v = GasRUniv / MolWtWet # 凝結成分の気体定数, EpsV = MolWtWet / MolWtDry # 凝結成分と大気の分子量比 TIME_START = 1000 TIME_END = 2000 XMIN = 170 XMAX = 340 YMIN = 90 YMAX = 450 WSN = 2 #MARK_SIZE = 0.01 vars = [ # {'name'=>'Evap', 'min'=>0, 'max'=>600}, # {'name'=>'Rain', 'min'=>0, 'max'=>500}, {'name'=>'OLR', 'min'=>270, 'max'=>400}, # {'name'=>'Sens', 'min'=>-20, 'max'=>100}, # {'name'=>'SLR', 'min'=>0, 'max'=>200} ] #omegas = [0.0, 0.01, 0.02, 0.03, 0.05, 0.06, 0.07, # 0.1, 0.15, 0.2, 0.25, 0.33, 0.4, 0.5, 0.67, 0.8, 1.0] omegas = [0.0, 0.15, 0.5, 1.0] omegas = NArray.to_na(omegas) nomega = omegas.size # 2 次元以上のデータを 1 次元配列に落とす (GGraph.scatter 対策) def flatten_gphys(gp) na_dat = gp.val.flatten n = na_dat.size axis = dummy_axis(n) data = VArray.new( na_dat, {"long_name"=>gp.data.long_name, "units"=>gp.data.units.to_s}, gp.data.name ) gp = GPhys.new( Grid.new(axis), data ) end # ダミーの gphys を作成 (GGraph.scatter での図の範囲指定のため) def dummy_gphys(na) n = na.size axis = dummy_axis(n) data = VArray.new( na, {"long_name"=>"dummy", "units"=>"1"}, "dummy" ) gp = GPhys.new( Grid.new(axis), data ) end def dummy_axis(n) axis_a = VArray.new( NArray.sfloat(n).indgen, {"long_name"=>"dummy","units"=>"1"}, "dummy" ) axis = Axis.new.set_pos(axis_a) end # 領域平均. ただし lon のみ cut したものを使っても構わない def ave_area(gp, na_lat_weight) gp = gp.mean('lon') * na_lat_weight gp = 0.5 * gp.sum('lat') gp.val = [gp.val] end hash_gp_Tg = Hash.new hash_gp_OLR = Hash.new ### OLR 読み込み var = 'OLR' file = var + '.nc' gp_OLR = GPhys::IO.open(file, var) lon = GPhys::IO.open(file, 'lon') sig = GPhys::IO.open(file, 'sig') sig_weight = GPhys::IO.open(file, 'sig_weight') na_lon = lon.val na_lat = GPhys::IO.open(file, 'lat').val na_lat_weight = GPhys::IO.open(file, 'lat_weight').val nlon = na_lon.size nlat = na_lat.size nsig = sig.val.size # OLR 時間平均, 緯度平均 gp_OLR = gp_OLR.cut('time'=>TIME_START..TIME_END).mean('time') gp_OLR = gp_OLR * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum gp_OLR = gp_OLR.sum('lat') # まずはダミー gp_OLR_viewpoints = gp_OLR.copy gp_OLR_viewpoints.val = 0 for i in 0..nlon-1 if lon.val[i] < 90 then # 西のはじっこ # 幅は i + 90 gp_part1 = gp_OLR.cut('lon'=>0..lon.val[i]+90).mean('lon') # 90-lon だけはみ出す よって part2 の開始点は 360 - (90-i) = 270 + i # part2 の幅は 90-i gp_part2 = gp_OLR.cut('lon'=>270+lon.val[i]..360).mean('lon') # p (gp_part1 * (i+90) + gp_part2 * (90-i))/180.0 gp_OLR_viewpoints[i].val=(gp_part1.val*(lon.val[i]+90)+gp_part2.val * (90-lon.val[i]))/180.0 # p lon.val[i] p gp_OLR_viewpoints[i].val elsif lon.val[i] > 270 then # 東のはじっこ # 東側では i+90 - 360 = i-270 だけはみ出す. #よって 西側 part1 は 0 から i-270 # part1 の幅は i-270 gp_part1 = gp_OLR.cut('lon'=>0..lon.val[i]-270).mean('lon') # 東側は i -90 から 360 # 幅は 450-i gp_part2 = gp_OLR.cut('lon'=>lon.val[i]-90..360).mean('lon') gp_OLR_viewpoints[i].val=(gp_part1.val*(lon.val[i]-270)+gp_part2.val*(450-lon.val[i]))/180.0 # p lon.val[i] p gp_OLR_viewpoints[i].val else gp_OLR_viewpoints[i].val = gp_OLR.cut('lon'=>lon.val[i]-90..lon.val[i]+90).mean('lon').val # p lon.val[i] p gp_OLR_viewpoints[i].val end end wsn=2 # 描画準備 DCL.gropn(wsn) GGraph.set_fig('viewport'=>[0.3,0.6,0.3,0.5]) DCL.sgpset('lfull', true) # 全画面表示 DCL.sgpset('lcntl', false) DCL.uzfact(0.6) GGraph.line( gp_OLR_viewpoints, true, {'index'=>899, 'annotate'=>false} ) DCL.grcls exit exit # 範囲指定のためのダミーデータ na_xrange = NArray.to_na([XMIN, XMAX]) na_yrange = NArray.to_na([YMIN, YMAX]) gp_xrange = dummy_gphys(na_xrange) gp_yrange = dummy_gphys(na_yrange) gp_xrange.data.long_name = gp_Tg.data.long_name gp_yrange.data.long_name = gp_OLR.data.long_name gp_xrange.data.units = gp_Tg.data.units gp_yrange.data.units = gp_OLR.data.units ## 描画 # 描画準備 DCL.gropn(WSN) #GGraph.set_fig('viewport'=>[0.15,0.80,0.15,0.6]) GGraph.set_fig('viewport'=>[0.3,0.6,0.3,0.5]) GGraph.set_axes('ytitle'=>YNAME) DCL.sgpset('lfull', true) # 全画面表示 DCL.sgpset('lcntl', false) # 初期化? # DCL.grfrm #GGraph.set_axes('xtickint'=>xlabel_int) #GGraph.set_axes('xlabelint'=>xlabel_int) #GGraph.set_axes('ytickint'=>ylabel_int) #GGraph.set_axes('ylabelint'=>ylabel_int) # 図の範囲を指定する方法が分からないのでダミーの点を打つ GGraph.scatter( gp_xrange, gp_yrange, true, 'type'=>1) opt = { # 'max'=>VAL_MAX, 'min'=>VAL_MIN, \ 'annotate'=>false, 'title'=>TITLE, 'index'=>1 } # draw opt['type'] = 1 # 4 # GGraph.scatter( flatten_gphys( gp_Tg.cut('lon'=>0..179.99) ), # flatten_gphys( gp_OLR.cut('lon'=>0..179.99) ), # false, opt) # opt['min'] = RH_SCATTER_MIN # opt['max'] = RH_SCATTER_MAX opt['min'] = SCATTER_MIN opt['max'] = SCATTER_MAX for i in 0..omegas.size-1 omega = omegas[i] gp_Tg = flatten_gphys( hash_gp_Tg[omega].cut('lon'=>0..179.99) ) gp_OLR = flatten_gphys( hash_gp_OLR[omega].cut('lon'=>0..179.99) ) gp_color = gp_Tg.copy gp_color.val = gp_color.val.fill(i.to_f) GGraph.color_scatter( gp_Tg, gp_OLR, gp_color, false, opt) end # opt['type'] = 2 # 5 # GGraph.scatter( flatten_gphys( gp_Tg.cut('lon'=>180.01..360) ), # flatten_gphys( gp_OLR.cut('lon'=>180.01..360) ), # false, opt) # # opt['type'] = 4 # circle # GGraph.scatter( gp_Tg.cut('lon'=>90..90,'lat'=>0), # gp_OLR.cut('lon'=>90..90,'lat'=>0), # false, opt) # opt['type'] = 6 # cross # GGraph.scatter( gp_Tg.cut('lon'=>270..270,'lat'=>0), # gp_OLR.cut('lon'=>270..270,'lat'=>0), # false, opt) # opt['type'] = 7 # GGraph.scatter( ave_area( gp_Tg.cut('lon'=>0..179.99), na_lat_weight), # ave_area( gp_OLR.cut('lon'=>0..179.99), na_lat_weight), # false, opt) # opt['type'] = 8 # GGraph.scatter( ave_area( gp_Tg.cut('lon'=>180.01..360), na_lat_weight), # ave_area( gp_OLR.cut('lon'=>180.01..360), na_lat_weight), # false, opt) ## line # opt['index'] = 7 # opt['type'] = 1; GGraph.line( gp_all, true, opt ) ## mark # opt['index'] = 5; opt['size'] = MARK_SIZE # opt['type'] = 4 # circle # GGraph.mark( gp_all, true, opt ) # opt['type'] = 5 # cross # GGraph.mark( gp_west, false, opt ) # opt['type'] = 6 # square # GGraph.mark( gp_east, false, opt ) DCL.grcls # remove extra space of ps file if WSN == 2 then if File.exist?(TMPFILE) then puts("MESSAGE: after processing ... #{epsfn} will be generate.") `dclpsrmcm #{TMPFILE} | dclpsrot | dclpsmargin > #{epsfn}` end end `rm -f #{TMPFILE}` exit(0)