# -*- coding: euc-jp -*- require "numru/gphys" include NumRu # # = 概要 # # * 圧力方程式における熱膨張項を計算する # 計算結果は zzz_test-thermal-moist_NetsuBoutyouTerm.nc に出力される # * 今のところは, 非断熱加熱の term だけ # #gphys = GPhys::NetCDF_IO.open('T.jan.nc', 'T') gphysPTempBZ = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'PTempBZ') gphysDensBZ = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'DensBZ') gphysVelSoudBZ = GPhys::NetCDF_IO.open('thermal-moist_restart.nc', 'VelSoundBZ') gphysExnerAll = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll.nc', 'ExnerAll') gphysPTempCond = GPhys::NetCDF_IO.open('thermal-moist_PTempCond.nc', 'PTempCond') gphysPTempDisp = GPhys::NetCDF_IO.open('thermal-moist_PTempDisp.nc', 'PTempDisp') gphysPTempRad = GPhys::NetCDF_IO.open('thermal-moist_PTempRad.nc', 'PTempRad') gphysPTempTurb = GPhys::NetCDF_IO.open('thermal-moist_PTempTurb.nc', 'PTempTurb') gphysSfcHeatFlux = GPhys::NetCDF_IO.open('thermal-moist_SfcHeatFlux.nc', 'SfcHeatFlux') #gphysExnerAll = GPhys::NetCDF_IO.open('thermal-moist_ExnerAll_IntValue30day.nc', 'ExnerAll') #gphysPTempCond = GPhys::NetCDF_IO.open('thermal-moist_PTempCond_IntValue30day.nc', 'PTempCond') #gphysPTempDis = GPhys::NetCDF_IO.open('thermal-moist_PTempDisp_IntValue30day.nc', 'PTempDisp') #gphysPTempRad = GPhys::NetCDF_IO.open('thermal-moist_PTempRad_IntValue30day.nc', 'PTempRad') #gphysPTempTurb = GPhys::NetCDF_IO.open('thermal-moist_PTempTurb_IntValue30day.nc', 'PTempTurb') #gphysSfcHeatFlux = GPhys::NetCDF_IO.open('thermal-moist_SfcHeatFlux_IntValue30day.nc', 'SfcHeatFlux') #-- 必要な定数を設定 --# xmin = 0 xmax = 256000 dx = 1000 xn = xmax / dx zmin = 0 #zmax = 24000 #zmax = 30000 #zmax = 36000 zmax = 48000 dz = 300 zn = zmax / dz tmin = 0 tmax = 86400 * 50 dt = 1800 tn = tmax / dt CpDry = 1039.642343614574 #-- cut --# gphysPTempBZ = gphysPTempBZ.mean( 'y' ) gphysPTempBZ = gphysPTempBZ.cut( 'x'=>xmin..xmax ) #gphysPTempBZ = gphysPTempBZ.mean( 'x' ) gphysPTempBZ = gphysPTempBZ.cut( 'z'=>zmin..zmax ) #gphysPTempBZ = gphysPTempBZ.cut( 'z'=>150 ) print "\n Basic field of Potential Temperature :\n" p gphysPTempBZ.val gphysDensBZ = gphysDensBZ.mean( 'y' ) gphysDensBZ = gphysDensBZ.cut( 'x'=>xmin..xmax ) #gphysDensBZ = gphysDensBZ.mean( 'x' ) gphysDensBZ = gphysDensBZ.cut( 'z'=>zmin..zmax ) #gphysDensBZ = gphysDensBZ.cut( 'z'=>150 ) print "\n Basic field of Density :\n" p gphysDensBZ.val gphysVelSoudBZ = gphysVelSoudBZ.mean( 'y' ) gphysVelSoudBZ = gphysVelSoudBZ.cut( 'x'=>xmin..xmax ) #gphysVelSoudBZ = gphysVelSoudBZ.mean( 'x' ) gphysVelSoudBZ = gphysVelSoudBZ.cut( 'z'=>zmin..zmax ) #gphysVelSoudBZ = gphysVelSoudBZ.cut( 'z'=>150 ) print "\n Basic field of Sound Velocity :\n" p gphysVelSoudBZ.val gphysExnerAll = gphysExnerAll.mean( 'y' ) #gphysExnerAll = gphysExnerAll.cut( 'z'=>150 ) print "\n Exner function :\n" p gphysExnerAll.val gphysPTempCond = gphysPTempCond.mean( 'y' ) #gphysPTempCond = gphysPTempCond.cut( 'z'=>150 ) print "\n Condensation term of Potential Temperature :\n" p gphysPTempCond.val gphysPTempDisp = gphysPTempDisp.mean( 'y' ) #gphysPTempDisp = gphysPTempDisp.cut( 'z'=>150 ) print "\n Dispation term of Potential Temperature :\n" p gphysPTempDisp.val gphysPTempRad = gphysPTempRad.mean( 'y' ) #gphysPTempRad = gphysPTempRad.cut( 'z'=>150 ) print "\n Radiation term of Potential Temperature :\n" p gphysPTempRad.val gphysPTempTurb = gphysPTempTurb.mean( 'y' ) #gphysPTempTurb = gphysPTempTurb.cut( 'z'=>150 ) print "\n Turbulence term of Potential Temperature :\n" p gphysPTempTurb.val gphysSfcHeatFlux = gphysSfcHeatFlux.mean( 'y' ) #---------------------------- xzt_f = gphysSfcHeatFlux.val p xzt_f #exit #-- 定数の設定 --# xmin = 500 zmin = 150 #-- 軸の設定 --# x_a = gphysPTempCond.coord(0) x = Axis.new.set_pos(x_a) z_a = gphysPTempCond.coord(1) z = Axis.new.set_pos(z_a) t_a = gphysPTempCond.coord(2) t = Axis.new.set_pos(t_a) data = VArray.new( NArray.sfloat(xn,zn,tn+1), {"long_name"=>"surface heat flux", "units"=>"W.m-2"}, "SfcHeatFlux" ) xzt_f_gphys = GPhys.new( Grid.new(x,z,t), data ) p xzt_f_gphys #-- 値を配列に入れる --# #xzt_f_gphys[0..xn-1,0..zn-1,0..tn] = xzt_f[0..xn-1,0..zn-1,0..tn] xzt_f_gphys[0..xn-1,0,0..tn] = xzt_f[0..xn-1,0..tn] p xzt_f_gphys #-- netCDF ファイルの生成とファイルへの書き出し --# #outfile = NetCDF.create("thermal-moist_ExnerThermExp.nc") outfile = NetCDF.create("thermal-moist_2D-SfcHeatFlux.nc") GPhys::NetCDF_IO.write(outfile,xzt_f_gphys) # outfile を閉じる outfile.close #---------------------------- #exit gphysSfcHeatFlux = gphysSfcHeatFlux / (CpDry * gphysDensBZ * gphysExnerAll * dz ) print "\n Surface Potential Temperature Flux :\n" p gphysSfcHeatFlux.val #-- 非断熱加熱の項を計算 (密度) --# #gphysdtheta = ( gphysDensBZ / gphysPTempBZ ) * ( gphysPTempCond + gphysPTempDisp + gphysPTempRad + gphysPTempTurb ) * dt gphysdtheta = ( gphysDensBZ / gphysPTempBZ ) * ( gphysPTempCond + gphysPTempDisp + gphysPTempRad + gphysPTempTurb + gphysSfcHeatFlux ) * dt #print "\n Sumention of diabatic heat term in 30 days :\n" print "\n Diabatic heat term in 50 days :\n" p gphysdtheta.val gphysdrho = -gphysdtheta #-- 非断熱加熱の項を計算 (圧力) --# #gphysdpidt = ( gphysVelSoudBZ**2 / (CpDry * gphysPTempBZ) ) * ( 1 / gphysPTempBZ ) * ( gphysPTempCond + gphysPTempDisp + gphysPTempRad + gphysPTempTurb ) * dt gphysdpidt = ( gphysVelSoudBZ**2 / (CpDry * gphysPTempBZ) ) * ( 1 / gphysPTempBZ ) * ( gphysPTempCond + gphysPTempDisp + gphysPTempRad + gphysPTempTurb + gphysSfcHeatFlux ) * dt print "\n Sumention of diabatic heat term in Pressure Eq. :\n" p gphysdpidt.val p gphysdpidt #-- NetCDF ファイルへの書き込み --# #-- 定数の設定 --# xmin = 500 zmin = 150 #-- 値を配列に入れておく --# #xzt_f = gphysdrho.val #-- 軸の設定 --# x_a = gphysPTempCond.coord(0) x = Axis.new.set_pos(x_a) z_a = gphysPTempCond.coord(1) z = Axis.new.set_pos(z_a) t_a = gphysPTempCond.coord(2) #t_a = gphysPTempCond.coord(1) t = Axis.new.set_pos(t_a) #-- drho 用 --# # #-- データ用の配列を準備 --# #data = VArray.new( NArray.sfloat(xn,zn,tn+1), # {"long_name"=>"Contribution of Thermal Expansion Term in Density", "units"=>"Kg.m-3"}, # "DensThermExp" ) #xzt_f_gphys = GPhys.new( Grid.new(x,z,t), data ) #p xzt_f_gphys # # #-- 値を配列に入れる --# #xzt_f_gphys[0..xn-1,0..zn-1,0..tn] = xzt_f[0..xn-1,0..zn-1,0..tn] #p xzt_f_gphys # # #-- netCDF ファイルの生成とファイルへの書き出し --# ##outfile = NetCDF.create("thermal-moist_DensThermExp.nc") #outfile = NetCDF.create("zzz-thermal-moist_DensThermExp.nc") #GPhys::NetCDF_IO.write(outfile,xzt_f_gphys) # ## outfile を閉じる #outfile.close #-- dpidt 用 --# #-- 値を配列に入れておく --# xzt_f = gphysdpidt.val p xzt_f #-- データ用の配列を準備 --# data = VArray.new( NArray.sfloat(xn,zn,tn+1), {"long_name"=>"Contribution of Thermal Expansion Term in Pressure Eq.", "units"=>"1"}, "ExnerThermExpTerm" ) xzt_f_gphys = GPhys.new( Grid.new(x,z,t), data ) #data = VArray.new( NArray.sfloat(xn,tn+1), # {"long_name"=>"Contribution of Thermal Expansion Term in Pressure Eq.", "units"=>"1"}, # "ExnerThermExp" ) #xzt_f_gphys = GPhys.new( Grid.new(x,t), data ) p xzt_f_gphys #-- 値を配列に入れる --# xzt_f_gphys[0..xn-1,0..zn-1,0..tn] = xzt_f[0..xn-1,0..zn-1,0..tn] #xzt_f_gphys[0..xn-1,0..tn] = xzt_f[0..xn-1,0..tn] p xzt_f_gphys #-- netCDF ファイルの生成とファイルへの書き出し --# #outfile = NetCDF.create("thermal-moist_ExnerThermExp.nc") outfile = NetCDF.create("thermal-moist_Exner-ThermExpTerm-DthetaDt.nc") GPhys::NetCDF_IO.write(outfile,xzt_f_gphys) # outfile を閉じる outfile.close