#!/usr/bin/ruby # -*- coding: euc-jp -*- # =begin # = Nakajima et al. (1992) の Tg vs OLR の(ような)図を描く # # == 概要 # # * 折れ線 # * Nakajima et al. (1992) の Tg vs OLR の図念頭に, # # * 散布図 # * 平均: 全球, 昼半球, 夜半球 # * 各格子点 (少なくとも昼半球と夜半球の区別, # 恒星直下点と対蹠点がどこかは分かるようにする) # # == 更新履歴 # # * 2011/04/25 納多 哲史 新規作成 # # =end require "numru/ggraph" require "pp" include NumRu #DATA_HOME = "/mnt/nas0/collatz/noda/work/dcpam/dcpam5-20100929_mod/practice" #DATA_HOME = "/home/noda/work/analysis/one_command_visualization2" #DATA_HOME = "/home/noda/data/exps_101220_diff_T21L16_sr_Omega_sweep" DATA_HOME = "/home/noda/data" prefix_line = "./dat/onishi2012" #ary_rh = [100,90,80,70,60,50,40] ary_rh = [100,90,80,70,60,50] #RH_SCATTER_MAX = 1.0 #RH_SCATTER_MIN = 0.4 SCATTER_MAX = 3 SCATTER_MIN = 0 TMPFILE = "dcl.ps" epsfn = "Tg_OLR_120x120region.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] #MolWtDry = MolWtWet # 乾燥成分の平均分子量 [kg mol-1] R_v = GasRUniv / MolWtWet # 凝結成分の気体定数, EpsV = MolWtWet / MolWtDry # 凝結成分と大気の分子量比 #LatentHeat = 2.5e+6 # 凝結成分の潜熱 [J kg-1] LatentHeat = 43655.0 # 凝結成分の潜熱 [J mol] TIME_START = 1000 TIME_END = 2000 # 昼半球全部のプロットの場合 #XMIN = 170 #XMAX = 340 #YMIN = 90 #YMAX = 450 XMIN = 250 XMAX = 320 YMIN = 100 YMAX = 330 WSN = 2 #MARK_SIZE = 0.01 vars = [ {'name'=>'SurfTemp', 'min'=>240, 'max'=>320}, # {'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 # 描画準備 DCL.uzfact(0.6) hash_gp_Tg = Hash.new hash_gp_OLR = Hash.new ## 各 omega で回す. for i in 0..nomega-1 omega = omegas[i] # epsfn = 'Omega' + omega.to_s + 'E_Tg_OLR.eps' if File.exist?(epsfn) puts("MESSAGE: #{epsfn} exist. skipped.") next end dir = '101220_diff-DryAir-Vapor_T21L16_sr_Omega' + omega.to_s + '_dt20m' #p s = 1450 # dir = '100929_noVF_NewRay_T21L48_S' + s.to_s + '_sr_Omega' + omega.to_s + '_dt00.5m' ### 読み込み # OLR and basic value var = 'OLR' file = var + '.nc' path = File.join(DATA_HOME, dir, file) gp_OLR = GPhys::IO.open(path, var) sig = GPhys::IO.open(path, 'sig') sig_weight = GPhys::IO.open(path, 'sig_weight') na_lon = GPhys::IO.open(path, 'lon').val na_lat = GPhys::IO.open(path, 'lat').val na_lat_weight = GPhys::IO.open(path, 'lat_weight').val nlon = na_lon.size nlat = na_lat.size nsig = sig.val.size # SurfTemp var = 'SurfTemp' file = var + '.nc' path = File.join(DATA_HOME, dir, file) gp_Tg = GPhys::IO.open(path, var) ## Temp # var = 'Temp' # file = var + '.nc' # path = File.join(DATA_HOME, dir, file) # gp_t = GPhys::IO.open(path, var) # ## QVap # var = 'QVap' # file = var + '.nc' # path = File.join(DATA_HOME, dir, file) # gp_q = GPhys::IO.open(path, var) # ## Ps # var = 'Ps' # file = var + '.nc' # path = File.join(DATA_HOME, dir, file) # gp_ps = GPhys::IO.open(path, var) ### 時間平均 # 昼半球全部 (オリジナル) # hash_gp_OLR[omega] = gp_OLR.cut('time'=>TIME_START..TIME_END).mean('time') # hash_gp_Tg[omega] = gp_Tg.cut('time'=>TIME_START..TIME_END).mean('time') # 直下点を中心とする 90度×90度領域 # hash_gp_OLR[omega] = gp_OLR.cut('lon'=>45..135, 'lat'=>-45..45).mean('time') # hash_gp_Tg[omega] = gp_Tg.cut('lon'=>45..135, 'lat'=>-45..45).mean('time') # 直下点を中心とする 120度×120度領域 hash_gp_OLR[omega] = gp_OLR.cut('lon'=>30..150, 'lat'=>-60..60).mean('time') hash_gp_Tg[omega] = gp_Tg.cut('lon'=>30..150, 'lat'=>-60..60).mean('time') end # omega のループ終わり # 1 次元モデルの結果 hash_gp_Tg_1d = Hash.new for rh in ary_rh file = "RH" + sprintf("%#03i", rh) + ".dat" path = File.join(prefix_line, file) arys_tmp = File.read(path).split("\n").map{|s| s.split(" ")} ary_tg = Array.new ary_olr = Array.new for ary in arys_tmp tg = ary[0].to_f olr = ary[1].to_f if tg <= XMAX then ary_tg.push(tg) ary_olr.push(olr) else ary_tg.push(XMAX) x0 = ary_tg[-1] x1 = tg y0 = ary_olr[-1] y1 = olr tmp = y0 + (y1-y0) / (x1-x0) * (XMAX - x0) ary_olr.push(tmp) break end end na_tg = NArray.to_na(ary_tg) na_olr = NArray.to_na(ary_olr) axis_a = VArray.new( na_tg, {"long_name"=>"surface temperature","units"=>"K"}, "SurfTemp" ) axis = Axis.new.set_pos(axis_a) data = VArray.new( na_olr, {"long_name"=>"Outgoing Longwave Flux", "units"=>"W m-2"}, "OLR" ) hash_gp_Tg_1d[rh] = GPhys.new( Grid.new(axis), data ) end ## relative humidity ## 対流圏は sigma=0.2-1.0 で決め打ち # ## 気圧の計算 # press = gp_ps.data.val.reshape(nlon,nlat,1) * sig.val.reshape(1,1,nsig) # ## 飽和水蒸気圧の計算 # sat_wat_vap_press = LatentHeat / (GasRUniv * gp_t) # sat_wat_vap_press *= -1.0 # sat_wat_vap_press = P_0 * sat_wat_vap_press.exp # ## 相対湿度の計算 # sat_spe_humidity = EpsV * sat_wat_vap_press / press # # sig_weight2 = sig_weight.cut('sig'=>0.2..1.0) # nsig2 = sig_weight2.val.size # # gp_q2 = gp_q.cut('sig'=>0.2..1.0) * sig_weight2.val.reshape(1,1,nsig2) # gp_q2 = gp_q2.sum('sig') # # sat_spe_humidity2 = sat_spe_humidity.cut('sig'=>0.2..1.0) * sig_weight2.val.reshape(1,1,nsig2) # sat_spe_humidity2 = sat_spe_humidity2.sum('sig') # # column_relative_humidity = gp_q2 / sat_spe_humidity2 # gp_rh = column_relative_humidity # #p gp_rh.max.val #p gp_rh.min.val # # 範囲指定のためのダミーデータ 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_axes('xtickint'=>10.0,'xlabelint'=>50.0) #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.line( 'max'=>val_max, 'min'=>val_min ) # 図の範囲を指定する方法が分からないのでダミーの点を打つ GGraph.scatter( gp_xrange, gp_yrange, true, 'type'=>1) opt = { # 'max'=>VAL_MAX, 'min'=>VAL_MIN, \ 'annotate'=>false, 'title'=>TITLE, 'index'=>1 } # draw for rh in ary_rh GGraph.line( hash_gp_Tg_1d[rh], false, 'index'=>1 ) end 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 # Omega0 だけ最後に書く for i in 1..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 omega = omegas[0] 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(0) GGraph.color_scatter( gp_Tg, gp_OLR, gp_color, false, opt) # 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)