#! /usr/bin/ruby # -*- coding: utf-8-with-signature -*- # # =begin # == 概要 # # 横軸に自転角速度, 縦軸に熱輸送量 [W] をとった図を作成する # 水, 大気別, かつ全球, 昼, 夜の合計 9 枚. # # 太陽定数変更実験用. # # == bug # 複数ファイルを生成しようとしているが seqv になる. # そのため, 現在では 1 枚描くごとに終了している. # 根気良く複数回実行すること. # # 現状では # $ ruby ./omega_heat-wat_atmos-surf_by_Watt_SolCon.rb # としないと実行できない. # # == 更新履歴 # # * 2015-02-25 石渡 太陽定数変更実験用に変更 # * 2014-09-22 石渡 複数実験シリーズ用に変更 # * 2013-11-15 石渡 地球放射計算用に変更 # * 2011/08/12 納多 哲史 Omega = 0.799 (0.8 で南北非対称が出た場合) を追加 # * 2011/04/27 納多 哲史 新規作成 # # =end require "numru/ggraph" require 'fileutils' include NumRu # SR_Omega_*_T42L26 シリーズの場合 #DATA_HOME = "/GFD_Dennou_Work10/momoko/SyncRotEarthRad-2" ##DIR_HEADER = ["SR_CLT0_Omega", "SR_CLT1500_Omega", "SR_CLT1500000_Omega"] ##DIR_HEADER = ["SR_CLT0_Omega", "SR_CLT1500_Omega"] #DIR_HEADER = ["SR_CLT0_Omega"] #DIR_FOOTER = "_T42L26" #DATA_HOME_LOCAL = DATA_HOME # SR_Omega_*_T42L26 シリーズの場合 DATA_HOME = "/GFD_Dennou_Work10/momoko/SyncRotEarthRad-2" DIR_HEADER = ["SR_CLT1500_Omega1.0_S"] DIR_FOOTER = "_T42L26" DATA_HOME_LOCAL = DATA_HOME nexp = DIR_HEADER.size LatentHeat = 2.5e6 TMPFILE = "dcl.ps" TIME_START = 731 TIME_END = 1095 #TIME_START = 620 #TIME_END = 920 VAL_MAX = 300 VAL_MIN = 0 WSN = 2 MARK_SIZE = 0.02 YNAME = 'Heat Flux' TITLE = '' 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'] = 9 #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 = [0.0, 0.1, 0.5, 1.0] na_omegas = NArray.to_na(omegas) nomega = omegas.size solcons = [1366, 1700, 1800] na_solcons = NArray.to_na(solcons) nsolcon = solcons.size flux = Hash.new flux_all = Hash.new flux_east = Hash.new flux_west = Hash.new hash_flux = Hash.new # 描画設定 DCL.uzfact(0.6) # 各変数で回す for iexp in 0..nexp-1 for v in vars p var = v['name'] # all = NArray.sfloat(nomega) # west = NArray.sfloat(nomega) # east = NArray.sfloat(nomega) all = NArray.sfloat(nsolcon) west = NArray.sfloat(nsolcon) east = NArray.sfloat(nsolcon) ## 各 omega で回す. # for i in 0..nomega-1 for i in 0..nsolcon-1 # omega = na_omegas[i] solcon = na_solcons[i] # p var.to_s + ' ' + omega.to_s ### 読み込み # dir = DIR_HEADER[iexp] + omega.to_s + DIR_FOOTER dir = DIR_HEADER[iexp] + solcon.to_s + DIR_FOOTER 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 end # end of loop of omega # p all # 領域平均値の計算 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 # end of loop of var temp = [flux_all, flux_east, flux_west] hash_flux[iexp] = Marshal.load(Marshal.dump(temp)) end # end of loop of exp flux_exps = Hash.new for iexp in 0..nexp-1 flux_all = hash_flux[iexp][0] flux_east = hash_flux[iexp][1] flux_west = hash_flux[iexp][2] # 昼から夜への輸送の計算 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 flux_exps[iexp] = Marshal.load(Marshal.dump(flux)) end ## gphys オブジェクトの生成 #na_omegas_4axis = na_omegas #tmp_long_name = DCL::csgi(151)+'|*"' #va_omega = VArray.new( na_omegas_4axis, # {"long_name"=>tmp_long_name}, # "omega" ) #axis_omega = Axis.new.set_pos(va_omega) na_solcons_4axis = na_solcons tmp_long_name = 'Solar Constant' va_solcon = VArray.new( na_solcons_4axis, {"long_name"=>tmp_long_name}, "SolCon" ) axis_solcon = Axis.new.set_pos(va_solcon) vars_draw = Hash.new p draw_types = ['atmos_heat'] # until 2014/10/20 #vars_draw['atmos_heat'] = ['OLRA', 'trans_lat' , 'trans_sens'] vars_draw['atmos_heat'] = ['OLRA', 'PRCP' , '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'] gp_exps = Hash.new ary = ['all', 'west', 'east'] for i in 0..ary.size-1 area = ary[i] for iexp in 0..nexp-1 gp = Hash.new for v in vars_wt_DtoN var = v['name'] long_name = YNAME units = 'W m|-2"' flux = flux_exps[iexp] tmp = flux[area][var] data = VArray.new( tmp, {"long_name"=>long_name, "units"=>units}, var ) # gp[var] = GPhys.new( Grid.new(axis_omega), data ) gp[var] = GPhys.new( Grid.new(axis_solcon), data ) end gp_exps[iexp] = Marshal.load(Marshal.dump(gp)) end # 描画準備 GGraph.set_fig('viewport'=>[0.3,0.6,0.3,0.5]) DCL.sgpset('lfull', true) # 全画面表示 DCL.sgpset('lcntl', true) # オプションのデフォルト値の設定 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 psfn = name + '.ps' epsfn = name + '.eps' vars_tmp = vars_draw[draw_type] # eps ファイルが存在する場合は終了 if File.exist?(epsfn) then puts("MESSAGE: #{epsfn} have already existed. Skipped.") next end DCL.gropn(WSN) p draw_type + '_' + area for iexp in 0..nexp-1 gp = gp_exps[iexp] 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 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'] = (iexp + 2) * 10 + 9 # 最後の 1 桁は太さ # 299 赤 # 499 青 # 509 緑 # 変数毎に色を変える場合 opt['index'] = index_type[v].to_i if iexp == 0 && i == 0 then GGraph.mark( gp_tmp, true, opt ) else GGraph.mark( gp_tmp, false, opt ) end end end # end of loop of vars end # end of loop of iexp DCL.grcls # ps ファイルの処理 # remove extra space of ps file if WSN == 2 then if File.exist?(TMPFILE) then # DCL.grcls FileUtils.mv(TMPFILE, psfn) puts("MESSAGE: after processing ... #{epsfn} will be generate.") `dclpsrmcm #{psfn} | dclpsrot | dclpsmargin > #{epsfn}` File.delete(psfn) exit end end end # end of loop of draw_type = 'atmos_wat' etc. end # end of loop of ary = ['all', 'west', 'east'] exit(0)