# # 相当温位をだす # ########################################################################## require "numru/gphys" # Gphys を使う include NumRu include NMath # NArray で計算できるように #-- 軸を読み込む/たぶんGphysオブジェクトを維持 #z = gphys.coord(2).val #-- 定数 Grav = 9.8 #CpDry = 1039.6423436145739742642 CpDry = 1004.0 # 標準的 CpDry T_0 = 273.16 #[K], 水の融点 P_0 = 1010.0 #deepconv 計算中での基準気圧 hPa R_d = 287.0 # J/kg 乾燥気体定数 R_v = 461.0 # J/kg 水蒸気の気体定数 GasRDry = 288.17502069917014750899 # deepconv のRd #--領域サイズ #----大文字変数は定数(ruby) #-- 幅設定 DelX = 500 DelZ =250 DelTime = 150 NT = (21600/150).to_i # どれだけ時間ステップか(NetCDFの軸を生成するために) NX = 256 NZ = 120 #-- NArray.を用意 zTemp = NArray.sfloat(NX) p_sat = NArray.sfloat(NX) #x_fNl = NArray.sfloat(NX) #x_fAl = NArray.sfloat(NX) #-- NArray オブジェクト #-- indgen(0,500)<- 0 から始めて 500 ずつ値を入れる x_a = VArray.new(NArray.sfloat(NX).indgen(0,500), {"long_name"=>"x-coodinate", "units"=>"m"}, # ハッシュ(ruby) "x") #---- GPhys の軸のオブジェクトに格上げ x = Axis.new.set_pos(x_a) z_a = VArray.new(NArray.sfloat(NZ).indgen(0,250), {"long_name"=>"z-coodinate", "units"=>"m"}, # ハッシュ(ruby) "z") #---- GPhys の軸のオブジェクトに格上げ z = Axis.new.set_pos(z_a) #-- 時間のオブジェクト t_a = VArray.new(NArray.sfloat(NT+1).indgen(0,150), {"long_name"=>"Time", "units"=>"sec"}, # ハッシュ(ruby) "t" ) #---- 軸のオブジェクトに格上げ t = Axis.new.set_pos(t_a) #-- 物理場オブジェクト theta_e_a = VArray.new(NArray.sfloat(NX,NZ,NT+1), {"long_name"=>"equivalent potential temperature", "units"=>"(1)"}, "theta_e") #---- 軸のオブジェクトに格上げ theta_e_gphys = GPhys.new(Grid.new(x,z,t), theta_e_a) #-- netCDF づくり ncFile = NetCDF.create("EPT.nc") GPhys::NetCDF_IO.write_grid(ncFile, Grid.new(x,z,t)) ###--- 必要な netCDF ファイルを open する ---### gphysPTemp = GPhys::IO.open('./thermal_PTempAll.nc','PTempAll') gphysExner = GPhys::IO.open('./thermal_ExnerAll.nc','ExnerAll') gphysH2Og = GPhys::IO.open('./thermal_H2O-gAll.nc','H2O-gAll') #-- 軸を読み込む/たぶんGphysオブジェクトを維持 tarray = gphysPTemp.coord(3).val zarray = gphysPTemp.coord(2).val ### y gphysPTemp = gphysPTemp.mean( 'y' ) gphysExner = gphysExner.mean( 'y' ) gphysH2Og = gphysH2Og.mean( 'y' ) for t in tarray #t=0 p t ### 時間でカット tPTemp = gphysPTemp.cut( 't'=>t ) tExner = gphysExner.cut( 't'=>t ) tH2Og = gphysH2Og.cut( 't'=>t ) for z in zarray ### zでカット zPTemp = tPTemp.cut( 'z'=>z ) zExner = tExner.cut( 'z'=>z ) zH2Og = tH2Og.cut( 'z'=>z ) zPTemp.units = '1' ### NArray オブジェクトにする zPTemp = zPTemp.val zExner = zExner.val zH2Og = zH2Og.val #-- 温度を計算 zTemp = zPTemp*zExner #-- 圧力を計算 p = P_0*(zExner**(CpDry/GasRDry)) #-- 飽和水蒸気圧 e_sat を計算 #-- e_sat = 6.112 exp[{17.67(T-273.15)}/(T-29.65)] [hPa] a = (17.67*(zTemp-273.15))/(zTemp-29.65) e_sat = 6.112*exp(a) #-- 飽和混合比 q_sat を計算 #-- q_sat = (R_d/R_v)(e_sat/p) q_sat = (R_d/R_v)*(e_sat/p) #-- 飽和温度 temp_sat を計算 #-- temp_sat = [1/(T-55)-ln(q_v/q_sat)/2840]^{-1}+55 [k] temp_sat = (1.0/(1.0/(zTemp-55.0)-(log(zH2Og/q_sat))/2840.0)) + 55.0 #-- 相当温位 theta_e を計算 #-- theta_e = T(Pref/p)**{0.2854(1-0.28 q_v)} *exp[q_v(1+0.81q_v)((3376/T_sat)-2.54)] b = (zH2Og*(1.0+(0.81*zH2Og))*((3376.0/temp_sat)-2.54)) theta_e = zTemp*(P_0/p)**(0.2854*(1-(0.28*zH2Og)))*exp(b) #-- 結果を最終変数に代入 k = ((z-125)/DelZ).to_i n = (t/DelTime).to_i theta_e_gphys[0..NX-1,k,n] = theta_e[0..NX-1] end end theta_e_gphys.data.set_att('long_name', "Equivalent potential temerature") theta_e_gphys.units = 'K' GPhys::NetCDF_IO.write(ncFile, theta_e_gphys.rename('theta_e')) # netCDFファイルにライト #rename で名前を変える ncFile.close