#!/usr/bin/ruby # -*- coding: euc-jp -*- # =begin # == 概要 # # 横軸に自転角速度, 縦軸に熱輸送量 [W] をとった図を作成する # 水, 大気別, かつ全球, 昼, 夜の合計 9 枚. # # == bug # # 複数ファイルを生成しようとしているが seqv になる. # そのため, 現在では 1 枚描くごとに終了している. # 根気良く複数回実行すること. # # == 更新履歴 # # * 2011/08/12 納多 哲史 Omega = 0.799 (0.8 で南北非対称が出た場合) を追加 # * 2011/04/27 納多 哲史 新規作成 # # =end require "numru/ggraph" require 'fileutils' include NumRu require '../load_config.rb' # x 軸をログプロットする場合 #axis_x_log = true axis_x_log = false #DATA_HOME_LOCAL = DATA_HOME # config.yml のを使う DATA_HOME_ENS = "/home/noda/data/exps_101220_diff-DryAir-Vapor_T21L16_sr_ensemble" DIR_HEADER_ENS = "101220_diff-DryAir-Vapor_T21L16_sr_Omega" DIR_HOOTER_ENS = "E_rndTemp_" TMPFILE = "dcl.ps" nexp = 10 draw_mean = 'yes' TIME_START = 1000 TIME_END = 2000 VAL_MAX = 400 VAL_MIN = 0 WSN = 2 MARK_SIZE = 0.01 YNAME = 'Heat Flux' TITLE = '' RPlanet =6.371 * 10.0 ** 6 # [m] vars = [ {'name'=>'OLR'}, {'name'=>'Evap'}, {'name'=>'Rain'}, {'name'=>'Sens'}, {'name'=>'OSR'}, {'name'=>'SLR'} ] 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['OLR'] = 4 mark_type['Rain'] = 5 mark_type['Evap'] = 6 mark_type['Sens'] = 7 mark_type['OSR'] = 9 mark_type['SLR'] = 10 mark_type['trans_all'] = 11 mark_type['trans_sens'] = 14 mark_type['trans_lat'] = 16 omegas.delete(0.0) if axis_x_log na_omegas = NArray.to_na(omegas) nomega = omegas.size flux_std = Hash.new flux_all_std = Hash.new flux_east_std = Hash.new flux_west_std = Hash.new flux_ens = Hash.new flux_all_ens = Hash.new flux_east_ens = Hash.new flux_west_ens = Hash.new # 描画設定 #p i # if i == 0 then DCL.uzfact(0.6) # end ## TOP of ensemble data handling # Variable Loop for v in vars p var = v['name'] all_ens = NArray.sfloat(nomega, nexp) west_ens = NArray.sfloat(nomega, nexp) east_ens = NArray.sfloat(nomega, nexp) ## Omega Loop for i in 0..nomega-1 omega = na_omegas[i] (0..nexp-1).each{|iexp| num = "00" + iexp.to_s # i is up to 9 dir = File.join(DATA_HOME_ENS, DIR_HEADER_ENS + omega.to_s + DIR_HOOTER_ENS + num) p var.to_s + ' ' + omega.to_s file = var + '.nc' ### Read Ensemble Data path = File.join(dir, file) gphys = GPhys::IO.open(path, var) 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 # time average gphys = gphys.cut('time'=>TIME_START..TIME_END).mean('time') # Horizontal Mean gphys = gphys * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum all_ens[i,iexp] = gphys.mean('lon').sum('lat').val west_ens[i,iexp] = gphys.cut('lon'=>na_lon[0]..na_lon[nlon/2 -1]).mean('lon').sum('lat').val east_ens[i,iexp] = gphys.cut('lon'=>na_lon[nlon/2]..na_lon[nlon-1]).mean('lon').sum('lat').val } end ### End of Omega loop case var when 'OLR' flux_all_ens[var] = all_ens flux_east_ens[var] = east_ens flux_west_ens[var] = west_ens when 'OSR' flux_all_ens[var] = -1.0 * all_ens flux_east_ens[var] = -1.0 * east_ens flux_west_ens[var] = -1.0 * west_ens when 'Evap' flux_all_ens[var] = -1.0 * all_ens flux_east_ens[var] = -1.0 * east_ens flux_west_ens[var] = -1.0 * west_ens else flux_all_ens[var] = all_ens flux_east_ens[var] = east_ens flux_west_ens[var] = west_ens end end # End of variable loop # 昼から夜への輸送の計算 d2n = Hash.new d2n['trans_all'] = flux_east_ens['OLR'] d2n['trans_lat'] = flux_west_ens['Evap'] + flux_west_ens['Rain'] d2n['trans_sens'] = d2n['trans_all'] + d2n['trans_lat'] # 後で図を描くときの便利のため, 全て同一の配列に押し込める for var in ['trans_all', 'trans_lat', 'trans_sens'] flux_all_ens[var] = d2n[var] flux_east_ens[var] = d2n[var] flux_west_ens[var] = -1.0 * d2n[var] end # form of Hash flux_ens # flux_ens[region][var][iomega,iexp] flux_ens['all'] = flux_all_ens flux_ens['east'] = flux_east_ens flux_ens['west'] = flux_west_ens ## BOTTOM of ensemble data setting ## TOP of standard data handling # Variable Loop for v in vars p var = v['name'] all = NArray.sfloat(nomega) west = NArray.sfloat(nomega) east = NArray.sfloat(nomega) ## Omega Loop for i in 0..nomega-1 omega = na_omegas[i] # p var.to_s + ' ' + omega.to_s ### Read STANDARD Data file = var + '.nc' p dir = exp_dirs[omega] path = File.join(dir, file) gphys = GPhys::IO.open(path, var) 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 # time average gphys = gphys.cut('time'=>TIME_START..TIME_END).mean('time') # Horizontal Mean gphys = gphys * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum # all, west, east : NArray 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 Omega loop case var when 'OLR' flux_all_std[var] = all flux_east_std[var] = east flux_west_std[var] = west when 'OSR' flux_all_std[var] = -1.0 * all flux_east_std[var] = -1.0 * east flux_west_std[var] = -1.0 * west when 'Evap' flux_all_std[var] = -1.0 * all flux_east_std[var] = -1.0 * east flux_west_std[var] = -1.0 * west else flux_all_std[var] = all flux_east_std[var] = east flux_west_std[var] = west end end # End of variable loop # 昼から夜への輸送の計算 d2n = Hash.new d2n['trans_all'] = flux_east_std['OLR'] d2n['trans_lat'] = flux_west_std['Evap'] + flux_west_std['Rain'] d2n['trans_sens'] = d2n['trans_all'] + d2n['trans_lat'] # 後で図を描くときの便利のため, 全て同一の配列に押し込める for var in ['trans_all', 'trans_lat', 'trans_sens'] flux_all_std[var] = d2n[var] flux_east_std[var] = d2n[var] flux_west_std[var] = -1.0 * d2n[var] end flux_std['all'] = flux_all_std flux_std['east'] = flux_east_std flux_std['west'] = flux_west_std ## BOTTOM of std data setting ### ### Enasemble Mean ### ensemble_mean_rain_small = NArray.sfloat(nomega) ensemble_mean_rain_large = NArray.sfloat(nomega) ensemble_mean_trans_sens_small = NArray.sfloat(nomega) ensemble_mean_trans_sens_large = NArray.sfloat(nomega) threshold_rain = 100.0 #for iomega in 0..3 do # ensemble_mean_rain_small[iomega] = 999 # ensemble_mean_rain_large[iomega] = 999 #end #for iomega in 4..nomega-1 do for iomega in 0..nomega-1 do counter_small_branch = 0 counter_large_branch = 0 ensemble_mean_rain_small[iomega] = 0.0 ensemble_mean_rain_large[iomega] = 0.0 # Standard Experiment if flux_east_std['Rain'][iomega] < threshold_rain && iomega > 9 ensemble_mean_rain_large[iomega] = ensemble_mean_rain_large[iomega] \ + flux_east_std['Rain'][iomega] counter_large_branch = counter_large_branch + 1 else ensemble_mean_rain_small[iomega] = ensemble_mean_rain_small[iomega] \ + flux_east_std['Rain'][iomega] counter_small_branch = counter_small_branch + 1 end # Ensemble Experiment for iexp in 0..nexp-1 do if flux_ens['east']['Rain'][iomega,iexp] < threshold_rain && iomega > 9 ensemble_mean_rain_large[iomega] = ensemble_mean_rain_large[iomega] \ + flux_ens['east']['Rain'][iomega,iexp] counter_large_branch = counter_large_branch + 1 else ensemble_mean_rain_small[iomega] = ensemble_mean_rain_small[iomega] \ + flux_ens['east']['Rain'][iomega,iexp] counter_small_branch = counter_small_branch + 1 end end if counter_small_branch == 0 then ensemble_mean_rain_small[iomega] = 999 else ensemble_mean_rain_small[iomega] = ensemble_mean_rain_small[iomega] \ / counter_small_branch end if counter_large_branch == 0 then ensemble_mean_rain_large[iomega] = 999 else ensemble_mean_rain_large[iomega] = ensemble_mean_rain_large[iomega] \ / counter_large_branch end end # trans_sens threshold_trans_sens = 180.0 #for iomega in 0..3 do # ensemble_mean_trans_sens_small[iomega] = 999 # ensemble_mean_trans_sens_large[iomega] = 999 #end #for iomega in 4..nomega-1 do for iomega in 0..nomega-1 do counter_small_branch = 0 counter_large_branch = 0 ensemble_mean_trans_sens_small[iomega] = 0.0 ensemble_mean_trans_sens_large[iomega] = 0.0 if flux_east_std['trans_sens'][iomega] > threshold_trans_sens && iomega > 9 ensemble_mean_trans_sens_large[iomega] = ensemble_mean_trans_sens_large[iomega] \ + flux_east_std['trans_sens'][iomega] counter_large_branch = counter_large_branch + 1 else ensemble_mean_trans_sens_small[iomega] = ensemble_mean_trans_sens_small[iomega] \ + flux_east_std['trans_sens'][iomega] counter_small_branch = counter_small_branch + 1 end for iexp in 0..nexp-1 do if flux_ens['east']['trans_sens'][iomega,iexp] > threshold_trans_sens && iomega > 9 ensemble_mean_trans_sens_large[iomega] = ensemble_mean_trans_sens_large[iomega] \ + flux_ens['east']['trans_sens'][iomega,iexp] counter_large_branch = counter_large_branch + 1 else ensemble_mean_trans_sens_small[iomega] = ensemble_mean_trans_sens_small[iomega] \ + flux_ens['east']['trans_sens'][iomega,iexp] counter_small_branch = counter_small_branch + 1 end end if counter_small_branch == 0 then ensemble_mean_trans_sens_small[iomega] = 999 else ensemble_mean_trans_sens_small[iomega] = ensemble_mean_trans_sens_small[iomega] \ / counter_small_branch end if counter_large_branch == 0 then ensemble_mean_trans_sens_large[iomega] = 999 else ensemble_mean_trans_sens_large[iomega] = ensemble_mean_trans_sens_large[iomega] \ / counter_large_branch end end DCL.glpset('lmiss',true) # handle data missing DCL::sglset('LCNTL', true) ## Omega axis na_omegas_4axis = na_omegas na_omegas_4axis = NMath::log10( na_omegas ) if axis_x_log # ORIGINAL VER. #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'] #vars_draw['atmos_heat'] = ['OLR', 'Sens', 'SLR', 'Rain', 'trans_sens'] vars_draw['atmos_heat'] = ['OLR' , 'Rain' , 'trans_sens'] vars_draw['atmos_wat'] = ['Rain', 'Evap', 'trans_lat'] vars_draw['surf_heat'] = ['Evap', 'Sens', 'SLR', 'OSR'] #ary = ['all', 'west', 'east'] ary = ['east'] GGraph.set_fig('viewport'=>[0.3,0.6,0.3,0.5]) DCL.sgpset('lfull', true) # 全画面表示 DCL.sgpset('lcntl', false) DCL.gropn(WSN) opt = { 'max'=>VAL_MAX, 'min'=>VAL_MIN, 'index'=>5, 'size'=>MARK_SIZE, 'legend'=>false, 'annotate'=>false, 'title'=>TITLE } # BEGIN region loop for i in 0..ary.size-1 area = ary[i] ### exp LOOP # std gp_std = Hash.new for v in vars_wt_DtoN var = v['name'] long_name = YNAME # ORIGINAL Ver. # units = 'W m-2' units = 'W m|-2"' tmp_std = flux_std[area][var] data_std = VArray.new( tmp_std, {"long_name"=>long_name, "units"=>units}, var) gp_std[var] = GPhys.new( Grid.new(axis_omega), data_std ) end for draw_type in draw_types name = draw_type + '_' + area + '_unit-watt' name = name + '_xlog' if axis_x_log psfn = name + '.ps' epsfn = name + '.eps' vars_tmp = vars_draw[draw_type] if File.exist?(epsfn) then puts("MESSAGE: #{epsfn} have already existed. Skipped.") next end # variable loop for i in 0..vars_tmp.size-1 v = vars_tmp[i] gp_tmp = gp_std[v] opt['type'] = mark_type[v].to_i if i == 0 then GGraph.mark( gp_tmp, true, opt ) else GGraph.mark( gp_tmp, false, opt ) end end # Marks of Ensemble data (0..nexp-1).each{|iexp| for i in 0..vars_tmp.size-1 long_name = YNAME # units = 'W m-2' units = 'W m|-2"' v = vars_tmp[i] tmp_ens = flux_ens[area][v][true,iexp] data_ens = VArray.new( tmp_ens, {"long_name"=>long_name, "units"=>units}, v) gp_tmp = GPhys.new( Grid.new(axis_omega), data_ens) opt['type'] = mark_type[v].to_i GGraph.mark( gp_tmp, false, opt ) end } # Ensemble Mean of rain data_rain_large = VArray.new( ensemble_mean_rain_large[true], {"long_name"=>'Rain', "units"=>units}, 'Rain' ) gp_rain_large = GPhys.new( Grid.new(axis_omega), data_rain_large ) data_rain_small = VArray.new( ensemble_mean_rain_small[true], {"long_name"=>'Rain', "units"=>units}, 'Rain' ) gp_rain_small = GPhys.new( Grid.new(axis_omega), data_rain_small ) if draw_mean == 'yes' then GGraph.line( gp_rain_small, false, {'index' => 299} ) GGraph.line( gp_rain_large, false, {'index' => 899} ) end # Ensemble Mean of trans_sens data_trans_sens_large = VArray.new( ensemble_mean_trans_sens_large[true], {"long_name"=>'trans_sens', "units"=>units}, 'trans_sens' ) gp_trans_sens_large = GPhys.new( Grid.new(axis_omega), data_trans_sens_large ) data_trans_sens_small = VArray.new( ensemble_mean_trans_sens_small[true], {"long_name"=>'trans_sens', "units"=>units}, 'trans_sens' ) gp_trans_sens_small = GPhys.new( Grid.new(axis_omega), data_trans_sens_small ) if draw_mean == 'yes' then GGraph.line( gp_trans_sens_small, false, {'index' => 899, 'type'=> 3} ) GGraph.line( gp_trans_sens_large, false, {'index' => 299, 'type'=> 3} ) end # DCL::sglset('LCNTL', true) # DCL::sgtxv(0.6, 0.6, 'W m|-2"') # DCL::sgtxv(0.7, 0.7, DCL::csgi(151)+'_E"') # DCL::sgtxv(0.7, 0.7, DCL::csgi(151)+'|*"') DCL.grcls # remove extra space of ps file if WSN == 2 then if File.exist?(TMPFILE) then FileUtils.mv(TMPFILE, psfn) puts("MESSAGE: after processing ... #{epsfn} will be generate.") `dclpsrmcm #{psfn} | dclpsrot | dclpsmargin > #{epsfn}` File.delete(psfn) exit end end end # draw_type LOOP end # END region LOOP # debug # DCL.grcls exit(0)