#! /usr/bin/ruby # -*- coding: utf-8-with-signature -*- # # =begin # == 概要 # # 横軸に自転角速度, 縦軸に熱輸送量 [W] をとった図を作成する # 水, 大気別, かつ全球, 昼, 夜の合計 9 枚. # # == bug # watt 単位に変更. # # 複数ファイルを生成しようとしているが seqv になる. # そのため, 現在では 1 枚描くごとに終了している. # 根気良く複数回実行すること. # # 現状では # $ ruby ./omega_heat-wat_atmos-surf_by_Watt_v2.rb # としないと実行できない. # # == 更新履歴 # # * 2013-11-15 石渡 地球放射計算用に変更 # * 2011/08/12 納多 哲史 Omega = 0.799 (0.8 で南北非対称が出た場合) を追加 # * 2011/04/27 納多 哲史 新規作成 # # =end require "numru/ggraph" require 'fileutils' include NumRu require './load_config.rb' LatentHeat = 2.5e6 # x 軸をログプロットする場合 #axis_x_log = true axis_x_log = false #read_meant_data = false DATA_HOME_LOCAL = DATA_HOME # config.yml で設定されているディレクトリ TMPFILE = "dcl.ps" TIME_START = 730 TIME_END = 1095 VAL_MAX = 300 VAL_MIN = 0 WSN = 2 # もともとの設定 #MARK_SIZE = 0.01 # 2014 年秋気象学会の予稿用に変更 MARK_SIZE = 0.02 YNAME = 'Heat Flux' #TITLE = 'Heat Budget' TITLE = '' RPlanet =6.371 * 10.0 ** 6 # [m] vars = [ {'name'=>'OLRA'}, {'name'=>'EvapU'}, {'name'=>'PRCP'}, {'name'=>'SensA'}, {'name'=>'OSRA'}, {'name'=>'SLRA'} ] vars_wt_DtoN = vars.dup # 単なる代入だと参照渡しになるため vars_wt_DtoN.push({'name'=>'trans_all'}) vars_wt_DtoN.push({'name'=>'trans_sens'}) vars_wt_DtoN.push({'name'=>'trans_lat'}) #1: 点 #2: 十字 #3: アスタリスク #4: 白抜き丸 #5: バツ #6: 白抜き正方形 #7: 白抜き上向き三角 #8: 白抜き菱形 #9: 白抜き星 #10: 黒塗り丸 #11: 黒塗り正方形 #12: 黒塗り上向き三角 #13: 黒塗り左向き三角 #14: 黒塗り下向き三角 #15: 黒塗り右向き三角 #16: 黒塗り星 #17: 黒塗り右向き旗 mark_type = Hash.new mark_type['OLRA'] = 4 mark_type['PRCP'] = 5 mark_type['EvapU'] = 6 mark_type['SensA'] = 7 mark_type['OSRA'] = 9 mark_type['SLRA'] = 10 mark_type['trans_all'] = 11 mark_type['trans_sens'] = 14 mark_type['trans_lat'] = 16 # 色をつける場合 index_type = Hash.new index_type['OLRA'] = 33 index_type['PRCP'] = 43 index_type['EvapU'] = 13 index_type['SensA'] = 13 index_type['OSRA'] = 13 index_type['SLRA'] = 13 index_type['trans_all'] = 3 index_type['trans_sens'] = 23 index_type['trans_lat'] = 43 omegas.delete(0.0) if axis_x_log #omegas = [0.0, 0.01, 0.05, 0.15, 0.5, 0.6, 0.75, 1.0] #omegas = [0.0, 0.01, 0.05, 0.15, 0.5, 0.75, 1.0] omegas = [0.0, 0.01, 0.05, 0.15, 0.5, 0.75, 1.0] #omegas = [0.0, 0.01, 0.05, 0.15, 0.5, 1.0] #omegas = [0.0, 0.01, 0.05, 0.5, 1.0] na_omegas = NArray.to_na(omegas) nomega = omegas.size flux = Hash.new flux_all = Hash.new flux_east = Hash.new flux_west = Hash.new # 描画設定 #p i # if i == 0 then DCL.uzfact(0.6) # end # 各変数で回す for v in vars p var = v['name'] all = NArray.sfloat(nomega) west = NArray.sfloat(nomega) east = NArray.sfloat(nomega) ## 各 omega で回す. for i in 0..nomega-1 omega = na_omegas[i] p var.to_s + ' ' + omega.to_s ### 読み込み dir = DIR_HEADER + omega.to_s + DIR_FOOTER # dir = 'SR_NoCloud_Omega' + omega.to_s # dir = 'SR_CLT1500_Omega' + omega.to_s file = var + '.nc' path = File.join(DATA_HOME_LOCAL, dir, file) p path gphys = GPhys::IO.open(path, var) if var == 'PRCP' then gphys = gphys*LatentHeat end 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 ### 平均の計算 #### 時間平均 gphys = gphys.cut('time'=>TIME_START..TIME_END).mean('time') #### 空間平均 ##### 経度方向の重み付け. gphys = gphys * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum all[i] = gphys.mean('lon').sum('lat').val west[i] = gphys.cut('lon'=>na_lon[0]..na_lon[nlon/2 -1]).mean('lon').sum('lat').val east[i] = gphys.cut('lon'=>na_lon[nlon/2]..na_lon[nlon-1]).mean('lon').sum('lat').val ## check # p all - (west + east)/2.0 end # 領域平均値の計算 case var # when 'OLRA' # flux_all[var] = -1.0 * all # flux_east[var] = -1.0 * east # flux_west[var] = -1.0 * west when 'OSRA' flux_all[var] = -1.0 * all flux_east[var] = -1.0 * east flux_west[var] = -1.0 * west when 'EvapU' flux_all[var] = -1.0 * all flux_east[var] = -1.0 * east flux_west[var] = -1.0 * west else flux_all[var] = all flux_east[var] = east flux_west[var] = west end end # 昼から夜への輸送の計算 d2n = Hash.new #d2n['trans_all'] = -1.0 * flux_east['OLRA'] d2n['trans_all'] = flux_east['OLRA'] # 潜熱, 夜半球で正になるように一時的に書き換えた. # 昼半球の図はコンシステントになってない. d2n['trans_lat'] = -1.0 * (flux_west['EvapU'] + flux_west['PRCP']) d2n['trans_sens'] = d2n['trans_all'] - d2n['trans_lat'] # 後で図を描くときの便利のため, 全て同一の配列に押し込める for var in ['trans_all', 'trans_lat', 'trans_sens'] flux_all[var] = d2n[var] flux_east[var] = d2n[var] flux_west[var] = -1.0 * d2n[var] end flux['all'] = flux_all flux['east'] = flux_east flux['west'] = flux_west ## gphys オブジェクトの生成 na_omegas_4axis = na_omegas na_omegas_4axis = NMath::log10( na_omegas ) if axis_x_log # (2013-6-14 石渡) 本当は, omega 軸の long_name は Omega/Omega_E # にしたかったのだけど, なぜかうまくいかない. # しかたなく, 今は Omega/OmegaE にしてある. #tmp_long_name = "Omega/OmegaE" tmp_long_name = DCL::csgi(151)+'|*"' tmp_long_name = 'log10(' + tmp_long_name + ')' if axis_x_log va_omega = VArray.new( na_omegas_4axis, {"long_name"=>tmp_long_name}, "omega" ) axis_omega = Axis.new.set_pos(va_omega) # testing vars_draw = Hash.new #p draw_types = ['atmos_heat', 'atmos_wat', 'surf_heat', 'trans'] p draw_types = ['atmos_heat'] #vars_draw['atmos_heat'] = ['OLRA', 'SensA', 'SLRA', 'PRCP', 'trans_sens'] vars_draw['atmos_heat'] = ['OLRA', 'trans_lat' , 'trans_sens'] vars_draw['atmos_wat'] = ['PRCP', 'EvapU', 'trans_lat'] vars_draw['surf_heat'] = ['EvapU', 'SensA', 'SLRA', 'OSRA'] vars_draw['trans'] = ['OLRA', 'trans_sens', 'trans_lat'] ary = ['all', 'west', 'east'] for i in 0..ary.size-1 area = ary[i] # epsfn = 'atmos_wat_' + area + '.eps' # if File.exist?(epsfn) then # puts(epsfn + " exists. process skipped.") # next # end gp = Hash.new for v in vars_wt_DtoN var = v['name'] long_name = YNAME # units = 'W m-2' units = 'W m|-2"' # units = '10^15 W' tmp = flux[area][var] # tmp = 2.0 * Math::PI * (RPlanet **2) * tmp / (10.0 ** 15) # if area == 'all' then # tmp = 2.0 * tmp # end data = VArray.new( tmp, {"long_name"=>long_name, "units"=>units}, var ) gp[var] = GPhys.new( Grid.new(axis_omega), data ) end ## 描画 # 描画準備 # testing # DCL.gropn(WSN) GGraph.set_fig('viewport'=>[0.3,0.6,0.3,0.5]) DCL.sgpset('lfull', true) # 全画面表示 DCL.sgpset('lcntl', true) # 初期化? # 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 ) # オプションのデフォルト値の設定 opt = { 'max'=>VAL_MAX, 'min'=>VAL_MIN, 'index'=>5, 'size'=>MARK_SIZE, 'index'=>13, 'legend'=>false, 'annotate'=>false, 'title'=>TITLE } for draw_type in draw_types p draw_type name = draw_type + '_' + area name = name + '_xlog' if axis_x_log psfn = name + '.ps' epsfn = name + '.eps' vars_tmp = vars_draw[draw_type] # if area == 'all' then #p vars_tmp.delete('trans_all') #p vars_tmp.delete('trans_sens') #p vars_tmp.delete('trans_lat') # end # eps ファイルが存在する場合は終了 if File.exist?(epsfn) then puts("MESSAGE: #{epsfn} have already existed. Skipped.") next end # testing DCL.gropn(WSN) p draw_type + '_' + area for i in 0..vars_tmp.size-1 v = vars_tmp[i] gp_tmp = gp[v] if (draw_type == 'surf_heat') and ((v == 'SensA') or (v == 'SLRA')) then gp_tmp.val = -1.0 * gp_tmp.val elsif (draw_type == 'atmos_wat') then gp_tmp.val = -1.0 * gp_tmp.val elsif (draw_type == 'atmos_wat') then gp_tmp.val = -1.0 * gp_tmp.val elsif (draw_type == 'trans') and (v == 'trans_sens') then gp_tmp.val = -1.0 * gp_tmp.val end # (2013-11-16 石渡) # 次は必要ないはずだが.... # if (draw_type == 'atmos_heat') and (v == 'trans_sens') then # gp_tmp.val = -1.0 * gp_tmp.val # end p gp_tmp.val if not area == 'all' or not v.include?('trans') then p area + " " + v opt['type'] = mark_type[v].to_i # opt['index'] = index_type[v].to_i opt['index'] = 1 if i == 0 then GGraph.mark( gp_tmp, true, opt ) else GGraph.mark( gp_tmp, false, opt ) end end end # ps ファイルの処理 # remove extra space of ps file if WSN == 2 then if File.exist?(TMPFILE) then # debug DCL.grcls # DCL.grfrm # DCL.grfig FileUtils.mv(TMPFILE, psfn) puts("MESSAGE: after processing ... #{epsfn} will be generate.") `dclpsrmcm #{psfn} | dclpsrot | dclpsmargin > #{epsfn}` File.delete(psfn) # testing exit end end end end # debug # DCL.grcls exit(0)