[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[dennou-ruby:001833] Re: GPhys::EP_Flux



堀之内です

塚原君は、その後さらにチュートリアル用プログラム作りを進めた上で
帰っていきました。お疲れ様でした。

さて、塚原様、製品の改良に向けてです。

* 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
@xxxxxx@xxxxxx -105,10 +105,10 @xxxxxx@xxxxxx
       R   = 287.0   # gas constant per unit mass for dry air of the earth [J.K-1.kg-1]
 
       #<<< module variable >>>
-      @xxxxxx@xxxxxx = UNumeric.new(7000,  "m")   # log pressure scale height in [m]  
-      @xxxxxx@xxxxxx       = UNumeric.new(6.37E6,"m")   # radius in [m]
-      @xxxxxx@xxxxxx   = UNumeric.new(8.64E4,"s")   # rotation period in [s]
-      @xxxxxx@xxxxxx          = UNumeric.new(1.0E5,"Pa")   # surface pressure in [Pa]
+      @xxxxxx@xxxxxx = UNumeric.new(7000,  "m")   # log-pressure scale height
+      @xxxxxx@xxxxxx       = UNumeric.new(6.37E6,"m")   # radius of the planet
+      @xxxxxx@xxxxxx   = UNumeric.new(8.64E4,"s")   # rotation period of the plnt
+      @xxxxxx@xxxxxx          = UNumeric.new(1.0E5,"Pa")   # Reference surface pressure
 
       module_function
 
@xxxxxx@xxxxxx -131,7 +131,7 @xxxxxx@xxxxxx
       def rot_period=(rp)
 	@xxxxxx@xxxxxx = rp
       end
-      def p00  # note that surface pressure cannot changeed, only read.
+      def p00  # note that surface pressure cannot be changeed, only read.
 	@xxxxxx@xxxxxx
       end
       def set_constants(scale_height, radius, rot_period)
@xxxxxx@xxxxxx -153,47 +153,44 @xxxxxx@xxxxxx
         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 = @xxxxxx@xxxxxx * 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 / @xxxxxx@xxxxxx) * 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')
 
@xxxxxx@xxxxxx -203,52 +200,54 @xxxxxx@xxxxxx
 	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 = @xxxxxx@xxxxxx * 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
@xxxxxx@xxxxxx -258,32 +257,20 @xxxxxx@xxxxxx
 	}
 	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 = @xxxxxx@xxxxxx*exp(-z/@xxxxxx@xxxxxx)
-	  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*(-@xxxxxx@xxxxxx/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
@xxxxxx@xxxxxx -294,10 +281,7 @xxxxxx@xxxxxx
 	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/@xxxxxx@xxxxxx)
 	  gp_t.data.set_att('long_name', "Potential Temperature")
 	end
@xxxxxx@xxxxxx -305,22 +289,13 @xxxxxx@xxxxxx
       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 = @xxxxxx@xxxxxx(va_z.units)
-	  ln_va_z = -@xxxxxx@xxxxxx*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 = @xxxxxx@xxxxxx(gp_z.units)
+	  gp_z = -@xxxxxx@xxxxxx*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
@xxxxxx@xxxxxx -329,18 +304,13 @xxxxxx@xxxxxx
 
       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
@xxxxxx@xxxxxx -371,13 +341,13 @xxxxxx@xxxxxx
       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 >>>
      @xxxxxx@xxxxxx = UNumeric.new(7000,  "m")   # log-pressure scale height
      @xxxxxx@xxxxxx       = UNumeric.new(6.37E6,"m")   # radius of the planet
      @xxxxxx@xxxxxx   = UNumeric.new(8.64E4,"s")   # rotation period of the plnt
      @xxxxxx@xxxxxx          = UNumeric.new(1.0E5,"Pa")   # Reference surface pressure

      module_function

      #<<< access to constants method >>> -------------------------------------
      def scale_height 
	@xxxxxx@xxxxxx 
      end
      def scale_height=(h)
	@xxxxxx@xxxxxx = h 
      end
      def radius
	@xxxxxx@xxxxxx
      end
      def radius=(r)
	@xxxxxx@xxxxxx = r
      end
      def rot_period
	@xxxxxx@xxxxxx
      end
      def rot_period=(rp)
	@xxxxxx@xxxxxx = rp
      end
      def p00  # note that surface pressure cannot be changeed, only read.
	@xxxxxx@xxxxxx
      end
      def set_constants(scale_height, radius, rot_period)
	@xxxxxx@xxxxxx = scale_height; @xxxxxx@xxxxxx = radius; @xxxxxx@xxxxxx = rot_period
      end
      def get_constants
	return @xxxxxx@xxxxxx, @xxxxxx@xxxxxx, @xxxxxx@xxxxxx
      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 = @xxxxxx@xxxxxx * 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 / @xxxxxx@xxxxxx) * 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/@xxxxxx@xxxxxx)
	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 = @xxxxxx@xxxxxx * 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 = @xxxxxx@xxxxxx*exp(-z/@xxxxxx@xxxxxx)
	  gp_un = gp_units
	  pr = ( pr.data.units ).convert2(pr, gp_un*Units.new('s'))
	  pr.units = gp_un*Units.new('s')
	  gp = gp*(-@xxxxxx@xxxxxx/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/@xxxxxx@xxxxxx)
	  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 = @xxxxxx@xxxxxx(gp_z.units)
	  gp_z = -@xxxxxx@xxxxxx*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