[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[dennou-ruby:001833] Re: GPhys::EP_Flux
- To: dennou-ruby@xxxxxxxxxxx
- Subject: [dennou-ruby:001833] Re: GPhys::EP_Flux
- From: Takeshi Horinouchi <horinout@xxxxxxxxxxxxxxxxxx>
- Date: Tue, 10 Aug 2004 13:56:39 +0900
堀之内です
塚原君は、その後さらにチュートリアル用プログラム作りを進めた上で
帰っていきました。お疲れ様でした。
さて、塚原様、製品の改良に向けてです。
* lib/以下のファイルは全部インストールされますんで、テストプログラ
ムに注意してください。（口頭了解済）
* 適当に手を入れたのを添付します。
ep_flux.rb 私が手を入れたの
patch-ep_flux.rb 電脳においてるバージョンとの差分
(古いほうを ep_flux_old.rb として diff -u ep_flux_old.rb ep_flux.rb)
テストプログラム test.rb の結果が完全に一致しているのは確認し
ました。しかし、エンバグしてない保証はできないです。参考用につ
かって下さい。
* Units#convet2 では VArray の単位まで自動設定しませんでした。
よって、今のところ、後からの単位を書き込みは残しています。
これが無変換の場合に問題を生じうるのはご承知のとおり。
VArray#convert_units(to), GPhys#convert_units(to) を作り
対応します。
--- ep_flux_old.rb Mon Aug 9 18:53:33 2004
+++ ep_flux.rb Tue Aug 10 13:42:35 2004
@@ -105,10 +105,10 @@
R = 287.0 # gas constant per unit mass for dry air of the earth [J.K-1.kg-1]
#<<< module variable >>>
- @@scale_height = UNumeric.new(7000, "m") # log pressure scale height in [m]
- @@radius = UNumeric.new(6.37E6,"m") # radius in [m]
- @@rot_period = UNumeric.new(8.64E4,"s") # rotation period in [s]
- @@p00 = UNumeric.new(1.0E5,"Pa") # surface pressure in [Pa]
+ @@scale_height = UNumeric.new(7000, "m") # log-pressure scale height
+ @@radius = UNumeric.new(6.37E6,"m") # radius of the planet
+ @@rot_period = UNumeric.new(8.64E4,"s") # rotation period of the plnt
+ @@p00 = UNumeric.new(1.0E5,"Pa") # Reference surface pressure
module_function
@@ -131,7 +131,7 @@
def rot_period=(rp)
@@rot_period = rp
end
- def p00 # note that surface pressure cannot changeed, only read.
+ def p00 # note that surface pressure cannot be changeed, only read.
@@p00
end
def set_constants(scale_height, radius, rot_period)
@@ -153,47 +153,44 @@
lon_nm, lat_nm, z_nm = ax_lon.pos.name, ax_lat.pos.name, ax_z.pos.name
gp_lon, gp_lat, gp_z = make_gphys(ax_lon, ax_lat, ax_z)
- ## convert
- gp_z = to_z_if_pressure(gp_z) # convert P => z{-H*log(P/P00)} if gp_w is P
- gp_lon = to_rad_if_deg(gp_lon) # convert deg => rad if lon unit is deg
- gp_lat = to_rad_if_deg(gp_lat) # convert deg => rad if lat unit is deg
- gp_w = to_w_if_omega(gp_w, gp_z) # convert dz/dt => dP/dt if gp_w is dP/dt
+ ## convert axes
+ gp_z = to_z_if_pressure(gp_z) # P => z=-H*log(P/P00) (units-based)
+ gp_lon = to_rad_if_deg(gp_lon) # deg => rad (unit convesion)
+ gp_lat = to_rad_if_deg(gp_lat) # deg => rad (unit convesion)
+ gp_w = to_w_if_omega(gp_w, gp_z) # dP/dt => dz/dt (units-based)
gp_t = to_theta_if_temperature(gp_t, gp_z, flag_temp_or_theta)
- # convert temp => \theta depend on flag
- ## transform gphys
- grid = gp_u.grid_copy; old_grid = gp_u.grid_copy
- grid.axis(lon_nm).pos = gp_lon.axis(lon_nm).pos
- grid.axis(lat_nm).pos = gp_lat.axis(lat_nm).pos
- grid.axis(z_nm).pos = gp_z.axis(z_nm).pos
+ # temperature => potential temperature (if flag is true)
+
+ ## replace grid (without duplicating data)
+ grid = gp_u.grid_copy
+ old_grid = gp_u.grid_copy # saved to use in outputs
+ grid.axis(lon_nm).pos = gp_lon.data # in radian
+ grid.axis(lat_nm).pos = gp_lat.data # in radian
+ grid.axis(z_nm).pos = gp_z.data # log-p height
gp_u = GPhys.new(grid, gp_u.data)
gp_v = GPhys.new(grid, gp_v.data)
gp_w = GPhys.new(grid, gp_w.data)
gp_t = GPhys.new(grid, gp_t.data)
- ## only for test
- #v_mean = gp_v.mean(lon_nm)
- #w_mean = gp_w.mean(lon_nm)
- #t_mean = gp_t.mean(lon_nm)
## get each term
- # F_y and F_z
+ # needed in F_y and F_z
uv_dash, vt_dash, uw_dash = eddy_products(gp_u, gp_v, gp_w, gp_t, lon_nm)
theta_mean = gp_t.mean(lon_nm)
dtheta_dz = GPhys::Derivative::cderiv(theta_mean, z_nm)
cos_lat = cos(gp_lat)
a_cos_lat = @@radius * cos_lat
a_cos_lat.data.rename!('a_cos_lat')
- a_cos_lat.data.set_att('long_name', 'radius mult cos_lat')
+ a_cos_lat.data.set_att('long_name', 'radius * cos_lat')
remove_0_at_poles(a_cos_lat)
- # F_y only
+ # needed in F_y only
u_mean = gp_u.mean(lon_nm)
du_dz = GPhys::Derivative::cderiv(u_mean, z_nm)
- # F_z only
+ # needed in F_z only
f_cor = 2 * (2 * PI / @@rot_period) * sin(gp_lat)
f_cor.data.rename!('f_cor')
- f_cor.data.set_att('long_name', 'corioli parameter')
+ f_cor.data.set_att('long_name', 'Coriolis parameter')
ducos_dphi = GPhys::Derivative::cderiv( u_mean * cos_lat, lat_nm)
- avort = (-ducos_dphi/a_cos_lat) + f_cor # -- absolute vorticity
- sec = Units.new("s-1")
- avort.data.units = sec
+ avort = (-ducos_dphi/a_cos_lat) + f_cor # -- absolute vorticity
+ avort.data.units = "s-1"
avort.data.rename!('avort')
avort.data.set_att('long_name', 'zonal mean absolute vorticity')
@@ -203,52 +200,54 @@
epflx_z = ( - uw_dash + avort*vt_dash/dtheta_dz ) * cos_lat * sigma
epflx_y.data.name = "epflx_y"; epflx_z.data.name = "epflx_z"
- epflx_y.data.set_att("long_name", "Eliassen_Palm Flux y component")
- epflx_z.data.set_att("long_name", "Eliassen_Palm Flux z component")
+ epflx_y.data.set_att("long_name", "EP flux y component")
+ epflx_z.data.set_att("long_name", "EP flux z component")
## convert with past grid
gp_ary = [] # grid convertes gphyss into
grid_xmean = old_grid.delete_axes(lon_nm)
- [epflx_y, epflx_z, gp_lat, gp_z, u_mean, theta_mean, uv_dash, vt_dash, uw_dash, dtheta_dz].each {|gp|
+ [epflx_y, epflx_z, gp_lat, gp_z, u_mean, theta_mean, uv_dash, vt_dash,
+ uw_dash, dtheta_dz].each {|gp|
if grid_xmean.shape.size != gp.shape.size
gp_ary << gp
else
- gp_ary << GPhys.new(grid_xmean, gp.data)
+ gp_ary << GPhys.new(grid_xmean, gp.data) #back to the original grid
end
}
return gp_ary
end
- def div_sphere(gp_fphi, gp_fz, yzdims=[0,1])
+ def div_sphere(gp_fy, gp_fz, yzdims=[0,1])
## get axis and name
- ax_lat = gp_fphi.axis(yzdims[0]) # Axis of latitude
- ax_z = gp_fphi.axis(yzdims[1]) # Axis of vertical
+ ax_lat = gp_fy.axis(yzdims[0]) # Axis of latitude
+ ax_z = gp_fy.axis(yzdims[1]) # Axis of vertical
lat_nm, z_nm = ax_lat.pos.name, ax_z.pos.name
gp_lat, gp_z = make_gphys(ax_lat, ax_z)
## convert
- gp_z = to_z_if_pressure(gp_z) # convert P => z{-H*log(P/P00)} if gp_w is P
- gp_lat = to_rad_if_deg(gp_lat) # convert deg => rad if lat unit is deg
+ gp_z = to_z_if_pressure(gp_z) # P => z=-H*log(P/P00) (units-based)
+ gp_lat = to_rad_if_deg(gp_lat) # deg => rad (unit convesion)
- ## transform gphys
- grid = gp_fphi.grid_copy; cp_grid = gp_fphi.grid_copy
- grid.axis(lat_nm).pos = gp_lat.axis(lat_nm).pos
- grid.axis(z_nm).pos = gp_z.axis(z_nm).pos
- gp_fphi = GPhys.new(grid, gp_fphi.data)
+ ## replace grid (without duplicating data)
+ grid = gp_fy.grid_copy
+ cp_grid = gp_fy.grid_copy # saved to use in outputs
+ grid.axis(lat_nm).pos = gp_lat.data
+ grid.axis(z_nm).pos = gp_z.data
+ gp_fy = GPhys.new(grid, gp_fy.data)
gp_fz = GPhys.new(grid, gp_fz.data)
## d_F_phi_dz
a_cos_lat = @@radius * cos (gp_lat)
remove_0_at_poles(a_cos_lat)
- d_gp_phi_d_phi = GPhys::Derivative::cderiv(gp_fphi * cos(gp_lat), lat_nm)
+ d_gp_fy_d_phi = GPhys::Derivative::cderiv(gp_fy * cos(gp_lat), lat_nm)
## d_F_z_dz
- d_gp_z_d_z = GPhys::Derivative::cderiv(gp_fz, z_nm)
+ d_gp_fz_d_z = GPhys::Derivative::cderiv(gp_fz, z_nm)
- f_div = ( d_gp_phi_d_phi / a_cos_lat ) + d_gp_z_d_z
+ f_div = ( d_gp_fy_d_phi / a_cos_lat ) + d_gp_fz_d_z
## convert with past grid
return GPhys.new(cp_grid, f_div.data)
end
- def make_gphys(*ax_ary) # it will be lost when new geid.rb released : to use delete_ax
+ def make_gphys(*ax_ary) # it will be lost when new grid.rb released : to use delete_ax
gphys_ary = []
ax_ary.each{|ax|
ax_data = ax.pos
@@ -258,32 +257,20 @@
}
return gphys_ary
end
-=begin
- def make_gphys(ax_lon, ax_lat, ax_z) # for operate
- lon_data = ax_lon.pos
- lat_data = ax_lat.pos
- z_data = ax_z.pos
- lon_grid = Grid.new(ax_lon)
- lat_grid = Grid.new(ax_lat)
- z_grid = Grid.new(ax_z)
- gp_lon = GPhys.new(lon_grid, lon_data)
- gp_lat = GPhys.new(lat_grid, lat_data)
- gp_z = GPhys.new(z_grid, z_data)
- return gp_lon, gp_lat, gp_z
- end
-=end
def to_w_if_omega(gp, z) # it is only for z coordinate!!!
- if gp.data.units =~ Units.new("Pa/s")
+ gp_units = gp.data.units
+ if gp_units =~ Units.new("Pa/s")
pr = @@p00*exp(-z/@@scale_height)
- gp_un = gp.data.units
+ gp_un = gp_units
pr = ( pr.data.units ).convert2(pr, gp_un*Units.new('s'))
+ pr.units = gp_un*Units.new('s')
gp = gp*(-@@scale_height/pr)
gp.data.rename!("wwnd")
- gp.data.set_att('long_name', "W")
- elsif gp.data.units =~ Units.new("m/s")
- gp = (gp.data.units).convert2( gp, Units.new('m/s') )
- gp.data.units = Units["m/s"]
+ gp.data.set_att('long_name', "log-P vertical wind")
+ elsif gp_units =~ Units.new("m/s")
+ gp = (gp_units).convert2( gp, Units.new('m/s') )
+ gp.units = 'm/s'
else
raise ArgumentError,"units of gp.data (#{gp.data.units}) must be dimention of pressure/time or length/time."
end
@@ -294,10 +281,7 @@
if flag_temp_or_theta
gp_un = gp_t.data.units
gp_t = ( gp_t.data.units ).convert2(gp_t, Units.new("K"))
- #va_t = gp_t.data.copy
- #gr_t = gp_t.grid_copy
- #va_t.units = Units.new('K')
- #gp_t = GPhys.new(gr_t, va_t)
+ gp_t.units = 'K'
gp_t = gp_t*exp((R/C_p)*z/@@scale_height)
gp_t.data.set_att('long_name', "Potential Temperature")
end
@@ -305,22 +289,13 @@
end
def to_z_if_pressure(gp_z) # number in units is not considerd operater as log.
- if ( gp_z.data.units =~ Units.new('Pa') )&& ( gp_z.axis(0).pos.units =~ Units.new('Pa') )
- grid = gp_z.grid_copy
- va_z = gp_z.data.copy
- p00 = @@p00.convert(va_z.units)
- ln_va_z = -@@scale_height*log(va_z/p00)
- ln_va_z.set_att('long_name', "z")
- ln_va_z.rename!("z")
- grid.axis(0).pos=ln_va_z
- gp_z = GPhys.new(grid, ln_va_z)
- elsif ( gp_z.data.units =~ Units.new('m') )&& ( gp_z.axis(0).pos.units =~ Units.new('m') )
- grid = gp_z.grid_copy
- va_z = gp_z.data.copy
- va_z = (va_z.units).convert2( va_z, Units.new('m') )
- va_z.units = Units["m"]
- grid.axis(0).pos = va_z
- gp_z = GPhys.new(grid, va_z)
+ if ( gp_z.data.units =~ Units.new('Pa') )
+ p00 = @@p00.convert(gp_z.units)
+ gp_z = -@@scale_height*log(gp_z/p00)
+ gp_z.data.set_att('long_name', "z").rename!("z")
+ elsif ( gp_z.data.units =~ Units.new('m') )
+ gp_z = (gp_z.data.units).convert2( gp_z, Units.new('m') )
+ gp_z.units = 'm'
else
raise ArgumentError,"units of gp (#{gp.data.units}) must be dimention of pressure or length."
end
@@ -329,18 +304,13 @@
def to_rad_if_deg(gp)
if gp.data.units =~ Units.new("degrees")
- grid = gp.grid_copy
- va = gp.data.copy
- newva = (va.units).convert2(va, Units["rad"])
- newva.units = Units[""]
- grid.axis(0).pos=newva
- gp = GPhys.new(grid, newva)
+ gp = (gp.units).convert2(gp, Units["rad"])
+ gp.units = Units[""]
+ gp
elsif gp.data.units =~ Units.new('rad')
- grid = gp.grid_copy
- va = gp.data
- va.units = Units[""]
- grid.axis(0).pos = va
- gp = GPhys.new(grid, va)
+ gp.data = gp.data.copy
+ gp.data.units = Units[""]
+ gp
else
raise ArgumentError,"units of gp #{gp.data.units} must be equal to deg or radian."
end
@@ -371,13 +341,13 @@
end
def remove_0_at_poles(a_cos_lat)
- (a_cos_lat.eq 0).where.each{ |zero|
- if zero == 0
- a_cos_lat[zero] = (a_cos_lat[zero] + a_cos_lat[zero + 1])/2
- else
- a_cos_lat[zero] = (a_cos_lat[zero] + a_cos_lat[zero - 1])/2
- end
- }
+ a_cos_lat[0] = (a_cos_lat[0] + a_cos_lat[1])/2 if a_cos_lat[0] == 0
+ a_cos_lat[-1] = (a_cos_lat[-1] + a_cos_lat[-2])/2 if a_cos_lat[-1] == 0
+ if a_cos_lat.min.val <= 0
+ raise "Illegal cos(phi) data. phi must between -pi/2 and +pi/2 " +
+ "and aligned in increasing or decreasing order."
+ end
+ nil
end
end
require 'narray'
require 'numru/gphys'
require 'numru/gphys/derivative' # 'numru/derivative' loaded by GPhys::Derivative
############################################################
=begin
=module NumRu::GPhys::EP_Flux in ep_flux.rb
==todo
* 1 行を 80 文字以内に納める
* Recast document in English
* I must cheack every case
* I must care remove_0_at_poles => cos(PI/2).nq 0. it is resonsible for NArray?
==Index
* ((<module NumRu::GPhys::EP_Flux>))
* ((<ep_full_sphere>))
* 球面座標データから E-P_Flux を計算
* ((<div_sphere>))
* 球面座標の発散を計算
* ((<scale_height>))
* スケールハイトを参照
* ((<scale_height=(h)>))
* スケールハイトを与える
* ((<radius>))
* 半径を参照
* ((<radius=(r)>))
* 半径を与える
* ((<rot_period>))
* 回転角速度を与える
* ((<rot_period=(rp)>))
* スケールハイトを与える
* ((<get_constants>))
* スケールハイト, 半径, 自転角速度を参照
* ((<set_constants(scale_height, radius, rot_period)>))
* スケールハイト, 半径, 自転角速度を設定
* ((<make_gphys(ax_ary)>))
引数に与えられた配列に格納される Axis から GPhys オブジェクトを作成.
* ((<to_w_if_omega(gp)>))
* ((<gp>))が圧力速度の場合, 速度( ex. m/s ) に変換
* ((<to_z_if_pressure(gp)>))
* ((<gp>))が圧力の場合, z( ex. m ) に変換
* ((<to_theta_if_temperature(gp, flag)>))
* ((<flag>))が true の場合, gp を温位に変換
* ((<to_rad_if_deg(gp)>))
* ((<gp>))の単位が deg の場合, rad に変換
* ((<eddy_products(gp_u, gp_v, gp_w, gp_t, dimname)>))
* ((<dimname>)) に関する eddy flux を計算する.
* ((<remove_0_at_poles(cos_gp)>))
* ((<cos_gp>)) が極で 0 の場合に, 隣接する格子との平均値を代入し回避する.
=module NumRu::GPhys::EP_Flux
Module functions of EP_Flux Operater for NArray.
EP_Flux 計算モジュール.
---ep_flux(gp_u, gp_v, gp_w_or_omega, gp_temp_or_theta, flag1(temp_or_theta))
Calculate Eliassen-Palm Flux on the spherical coordinate.
memo:
各物理量のグリッドは全て同一であるもののみ扱う.
ARGUMENTS
*
RETURN VALUE
* F_y, F_z, rv1, ...
---div_sphere(data)#(u, v, w_or_omega, temp_or_theta, flag1(w_or_omega), flag2(temp_or_theta))
....description...
ARGUMENTS
* hoge...
RETURN VALUE
* GPhys
---scale_height
---scale_height=(h)
---radius
---radius=(r)
---rot_period
---rot_period=(rp)
---set_constants(scale_height, radius, rot_period)
---get_constants
=end
############################################################
module NumRu
class GPhys
module EP_Flux
include Misc::EMath
#<<< module constants >>>
C_p = 1004 # specific heat at constant pressure of the earth's atmosphere [J.K-1.kg-1]
R = 287.0 # gas constant per unit mass for dry air of the earth [J.K-1.kg-1]
#<<< module variable >>>
@@scale_height = UNumeric.new(7000, "m") # log-pressure scale height
@@radius = UNumeric.new(6.37E6,"m") # radius of the planet
@@rot_period = UNumeric.new(8.64E4,"s") # rotation period of the plnt
@@p00 = UNumeric.new(1.0E5,"Pa") # Reference surface pressure
module_function
#<<< access to constants method >>> -------------------------------------
def scale_height
@@scale_height
end
def scale_height=(h)
@@scale_height = h
end
def radius
@@radius
end
def radius=(r)
@@radius = r
end
def rot_period
@@rot_period
end
def rot_period=(rp)
@@rot_period = rp
end
def p00 # note that surface pressure cannot be changeed, only read.
@@p00
end
def set_constants(scale_height, radius, rot_period)
@@scale_height = scale_height; @@radius = radius; @@rot_period = rot_period
end
def get_constants
return @@scale_height, @@radius, @@rot_period
end
#<<< keisan method >>> ---------------------------------------------
def ep_full_sphere(gp_u, gp_v, gp_w, gp_t,
flag_temp_or_theta=true, xyzdims=[0,1,2])
## get axis and name
raise ArgumentError,"xyzdims's size (#{xyzdims.size}) must be 3." if xyzdims.size != 3
ax_lon = gp_u.axis(xyzdims[0]) # Axis of longitude
ax_lat = gp_u.axis(xyzdims[1]) # Axis of latitude
ax_z = gp_u.axis(xyzdims[2]) # Axis of vertical
lon_nm, lat_nm, z_nm = ax_lon.pos.name, ax_lat.pos.name, ax_z.pos.name
gp_lon, gp_lat, gp_z = make_gphys(ax_lon, ax_lat, ax_z)
## convert axes
gp_z = to_z_if_pressure(gp_z) # P => z=-H*log(P/P00) (units-based)
gp_lon = to_rad_if_deg(gp_lon) # deg => rad (unit convesion)
gp_lat = to_rad_if_deg(gp_lat) # deg => rad (unit convesion)
gp_w = to_w_if_omega(gp_w, gp_z) # dP/dt => dz/dt (units-based)
gp_t = to_theta_if_temperature(gp_t, gp_z, flag_temp_or_theta)
# temperature => potential temperature (if flag is true)
## replace grid (without duplicating data)
grid = gp_u.grid_copy
old_grid = gp_u.grid_copy # saved to use in outputs
grid.axis(lon_nm).pos = gp_lon.data # in radian
grid.axis(lat_nm).pos = gp_lat.data # in radian
grid.axis(z_nm).pos = gp_z.data # log-p height
gp_u = GPhys.new(grid, gp_u.data)
gp_v = GPhys.new(grid, gp_v.data)
gp_w = GPhys.new(grid, gp_w.data)
gp_t = GPhys.new(grid, gp_t.data)
## get each term
# needed in F_y and F_z
uv_dash, vt_dash, uw_dash = eddy_products(gp_u, gp_v, gp_w, gp_t, lon_nm)
theta_mean = gp_t.mean(lon_nm)
dtheta_dz = GPhys::Derivative::cderiv(theta_mean, z_nm)
cos_lat = cos(gp_lat)
a_cos_lat = @@radius * cos_lat
a_cos_lat.data.rename!('a_cos_lat')
a_cos_lat.data.set_att('long_name', 'radius * cos_lat')
remove_0_at_poles(a_cos_lat)
# needed in F_y only
u_mean = gp_u.mean(lon_nm)
du_dz = GPhys::Derivative::cderiv(u_mean, z_nm)
# needed in F_z only
f_cor = 2 * (2 * PI / @@rot_period) * sin(gp_lat)
f_cor.data.rename!('f_cor')
f_cor.data.set_att('long_name', 'Coriolis parameter')
ducos_dphi = GPhys::Derivative::cderiv( u_mean * cos_lat, lat_nm)
avort = (-ducos_dphi/a_cos_lat) + f_cor # -- absolute vorticity
avort.data.units = "s-1"
avort.data.rename!('avort')
avort.data.set_att('long_name', 'zonal mean absolute vorticity')
## F_y, F_z
sigma = exp(-gp_z/@@scale_height)
epflx_y = ( - uv_dash + du_dz*vt_dash/dtheta_dz ) * cos_lat * sigma
epflx_z = ( - uw_dash + avort*vt_dash/dtheta_dz ) * cos_lat * sigma
epflx_y.data.name = "epflx_y"; epflx_z.data.name = "epflx_z"
epflx_y.data.set_att("long_name", "EP flux y component")
epflx_z.data.set_att("long_name", "EP flux z component")
## convert with past grid
gp_ary = [] # grid convertes gphyss into
grid_xmean = old_grid.delete_axes(lon_nm)
[epflx_y, epflx_z, gp_lat, gp_z, u_mean, theta_mean, uv_dash, vt_dash,
uw_dash, dtheta_dz].each {|gp|
if grid_xmean.shape.size != gp.shape.size
gp_ary << gp
else
gp_ary << GPhys.new(grid_xmean, gp.data) #back to the original grid
end
}
return gp_ary
end
def div_sphere(gp_fy, gp_fz, yzdims=[0,1])
## get axis and name
ax_lat = gp_fy.axis(yzdims[0]) # Axis of latitude
ax_z = gp_fy.axis(yzdims[1]) # Axis of vertical
lat_nm, z_nm = ax_lat.pos.name, ax_z.pos.name
gp_lat, gp_z = make_gphys(ax_lat, ax_z)
## convert
gp_z = to_z_if_pressure(gp_z) # P => z=-H*log(P/P00) (units-based)
gp_lat = to_rad_if_deg(gp_lat) # deg => rad (unit convesion)
## replace grid (without duplicating data)
grid = gp_fy.grid_copy
cp_grid = gp_fy.grid_copy # saved to use in outputs
grid.axis(lat_nm).pos = gp_lat.data
grid.axis(z_nm).pos = gp_z.data
gp_fy = GPhys.new(grid, gp_fy.data)
gp_fz = GPhys.new(grid, gp_fz.data)
## d_F_phi_dz
a_cos_lat = @@radius * cos (gp_lat)
remove_0_at_poles(a_cos_lat)
d_gp_fy_d_phi = GPhys::Derivative::cderiv(gp_fy * cos(gp_lat), lat_nm)
## d_F_z_dz
d_gp_fz_d_z = GPhys::Derivative::cderiv(gp_fz, z_nm)
f_div = ( d_gp_fy_d_phi / a_cos_lat ) + d_gp_fz_d_z
## convert with past grid
return GPhys.new(cp_grid, f_div.data)
end
def make_gphys(*ax_ary) # it will be lost when new grid.rb released : to use delete_ax
gphys_ary = []
ax_ary.each{|ax|
ax_data = ax.pos
ax_grid = Grid.new(ax)
gp = GPhys.new(ax_grid, ax_data)
gphys_ary << gp
}
return gphys_ary
end
def to_w_if_omega(gp, z) # it is only for z coordinate!!!
gp_units = gp.data.units
if gp_units =~ Units.new("Pa/s")
pr = @@p00*exp(-z/@@scale_height)
gp_un = gp_units
pr = ( pr.data.units ).convert2(pr, gp_un*Units.new('s'))
pr.units = gp_un*Units.new('s')
gp = gp*(-@@scale_height/pr)
gp.data.rename!("wwnd")
gp.data.set_att('long_name', "log-P vertical wind")
elsif gp_units =~ Units.new("m/s")
gp = (gp_units).convert2( gp, Units.new('m/s') )
gp.units = 'm/s'
else
raise ArgumentError,"units of gp.data (#{gp.data.units}) must be dimention of pressure/time or length/time."
end
return gp
end
def to_theta_if_temperature(gp_t, z, flag_temp_or_theta) # it is only for z coordinate!!!
if flag_temp_or_theta
gp_un = gp_t.data.units
gp_t = ( gp_t.data.units ).convert2(gp_t, Units.new("K"))
gp_t.units = 'K'
gp_t = gp_t*exp((R/C_p)*z/@@scale_height)
gp_t.data.set_att('long_name', "Potential Temperature")
end
return gp_t
end
def to_z_if_pressure(gp_z) # number in units is not considerd operater as log.
if ( gp_z.data.units =~ Units.new('Pa') )
p00 = @@p00.convert(gp_z.units)
gp_z = -@@scale_height*log(gp_z/p00)
gp_z.data.set_att('long_name', "z").rename!("z")
elsif ( gp_z.data.units =~ Units.new('m') )
gp_z = (gp_z.data.units).convert2( gp_z, Units.new('m') )
gp_z.units = 'm'
else
raise ArgumentError,"units of gp (#{gp.data.units}) must be dimention of pressure or length."
end
return gp_z
end
def to_rad_if_deg(gp)
if gp.data.units =~ Units.new("degrees")
gp = (gp.units).convert2(gp, Units["rad"])
gp.units = Units[""]
gp
elsif gp.data.units =~ Units.new('rad')
gp.data = gp.data.copy
gp.data.units = Units[""]
gp
else
raise ArgumentError,"units of gp #{gp.data.units} must be equal to deg or radian."
end
return gp
end
def eddy_products(gp_u, gp_v, gp_w, gp_t, dimname)
# get zonal_eddy
u_dash = gp_u - gp_u.mean(dimname)
v_dash = gp_v - gp_v.mean(dimname)
w_dash = gp_w - gp_w.mean(dimname)
t_dash = gp_t - gp_t.mean(dimname)
# get eddy_product
uv_dash = u_dash*v_dash # u'v'
vt_dash = v_dash*t_dash # v't'
uw_dash = u_dash*w_dash # u'w'
# set attribute
uv_dash.data.set_att("long_name", "U'V'")
vt_dash.data.set_att("long_name", "V'T'")
uw_dash.data.set_att("long_name", "U'W'")
uv_dash.data.rename!("uv_dash")
vt_dash.data.rename!("vt_dash")
uw_dash.data.rename!("uw_dash")
return uv_dash.mean(dimname), vt_dash.mean(dimname), uw_dash.mean(dimname)
end
def remove_0_at_poles(a_cos_lat)
a_cos_lat[0] = (a_cos_lat[0] + a_cos_lat[1])/2 if a_cos_lat[0] == 0
a_cos_lat[-1] = (a_cos_lat[-1] + a_cos_lat[-2])/2 if a_cos_lat[-1] == 0
if a_cos_lat.min.val <= 0
raise "Illegal cos(phi) data. phi must between -pi/2 and +pi/2 " +
"and aligned in increasing or decreasing order."
end
nil
end
end
end
end