############################################################

=begin
=method for class NumRu::GPhys and  NumRu::GphysFunc in libgphys-n.rb

==Index
* ((<class NumRu::NArray (methods added)>))
  * ((<cshift(n)>))
    要素全体を n 個ずらした NArray オブジェクトを返す.

* ((<class NumRu::GPhys (methods added)>))
  * ((<enhanced_strm_function(vname)>))
    南北風のデータから質量流線関数を計算. 用いた式は...
    daily データしか扱えない.
  * ((<strm_function_ncep(vname)>))
    南北風のデータから質量流線関数を計算. 用いた式は...
    時間平均をとった NCEP データしか扱えない. (圧力軸の向きが NCEP 固有.)
  * ((<theta(vname)>))
    温度から温位を計算.
  * ((<p_bar(vname)>))
    温位から等温位面に沿った圧力の平均を計算.
  * ((<cd_dp(v, k, dp)>))
    一階微分を計算. 中央差分.
  * ((<fd_dp(v, k, dp)>))
    一階微分を計算. 前方差分.
  * ((<bd_dp(v, k, dp)>))
    一階微分を計算. 後方差分.
  * ((<ape(vname)>))
    温位から有効位置エネルギーを計算.    
  * ((<stspct_fft(vname)>))
    時空間スペクトルを計算.
  * ((<stspct_coh(vname)>))
    P-R コヒーレンスを計算.
  * ((<stspct_st_tr(vname)>))
    時空間スペクトルのうち, 定在成分(standing wave) と 移動成分(traveling wave) を計算.
  * ((<make_sfc_mask(p_gphys)>))
    引数に与えた地表面圧力の gphys オブジェクトから,self のgphys オブジェクトと同じ配列構造をもつ, 3 次元のマスク gphys オブジェクトを作成する.
    実際は地面下なのに, 外装してある部分を排除するためのマスクを作成.
  * ((<save(output, global_attr=nil, vname=nil, var_attr=nil, history=$history, new_vname=nil)>))
    GPhys オブジェクトを netCDF ファイルに保存.
  * ((<add_lost_axes( lost )>))
    データ, 座標軸操作記録をいじる
  * ((<set_lost_axes( lost )>))
    データ, 座標軸操作記録をいじる
  * ((<rename(name)>))
    ((<rename(name)>)) を gphys_class を返すよう再定義
  * ((<set_att(name,val)>))
    ((<set_att>)) を gphys_class を返すよう再定義

* other meshod
  * ((<global_attr(nc_file)>))
    * nc_file で指定した netCDF ファイルの大域属性を返す
  * ((<mean_gphys>))
    * 引数に指定した gphys オブジェクト群の平均をとる
  * ((<ask_overwrite>))
    * netCDF ファイルに上書き保存してよいか確認
=end
require "numru/ggraph"
require "numru/netcdf_miss"
require "date"
include NumRu

class NArray # NArray への追加メソッド
  def cshift(n) # Fortran90 の cshiftの代わり
    #if(n > x.size) ...
    y = self.dup
    y[n..-1] = self[0..-1-n]
    y[0..n-1] = self[-n..-1]
    return y
  end
end

class NArray # for deep_copy
  def self._load(o); to_na(*Marshal::load(o)).ntoh; end
  def _dump(limit); Marshal::dump([hton.to_s, typecode, *shape]); end
end

=begin
module NumRu
  class NetCDFVar
     def get_with_miss_and_scaling2(*args) # make mask before scaling
      __interpret_missing_params if !defined?(@missval)
      packed_data = simple_get(*args)
      scaled_data = scaled_get(*args)
      sf = att('scale_factor')
      ao = att('add_offset')
      if @vmin || @vmax
	if sf && ao
	  csf = sf.get
	  cao = ao.get
	  vmin = (@vmin-cao)/csf if @vmin
	  vmax = (@vmax-cao)/csf if @vmax
	elsif
	  vmin = @vmin; vmax = @vmax
	end
	if vmin
	  mask = (packed_data >= vmin) 
	mask = mask.and(packed_data <= vmax) if vmax
	else
	  mask = (packed_data <= vmax)
	end
	data = NArrayMiss.to_nam(scaled_data, mask)
      elsif @missval	# only missing_value is present.
	eps = 1e-6
	missval = @missval[0].to_f
	vmin = missval - eps*missval
	vmax = missval + eps*missval
	mask = (packed_data <= vmin) 
	mask = mask.or(packed_data >= vmax)
	data = NArrayMiss.to_nam(scaled_data, mask)
      else
	data = scaled_data
      end
      data
    end

    def put_with_miss_after_scaling(data, *args)
      if data.is_a?( NArrayMiss )
	__interpret_missing_params if !defined?(@missval)
	if @missval
	  scaled_put_without_missval(data, *args)
	else
	  scaled_put(data.to_na, *args)
	end
      else
	scaled_put(data, *args)
      end
    end

    def scaled_put_without_missval(var,hash=nil)
      sf = att('scale_factor')
      ao = att('add_offset')
      if ( sf == nil && ao == nil ) then
	# no scaling --> just call put
	simple_put(var,hash)
      else
	if (sf != nil) 
	  csf = sf.get
	  if csf.is_a?(NArray) then  # --> should be a numeric
	    csf = csf[0]
	  elsif csf.is_a?(String)
	    raise TypeError, "scale_factor is not a numeric"
	  end
	  if(csf == 0) then; raise NetcdfError, "zero scale_factor"; end
	else
	  csf = 1.0      # assume 1 if not defined
	end
	if (ao != nil) 
	  cao = ao.get
	if cao.is_a?(NArray) then  # --> should be a numeric
	  cao = cao[0]
	elsif csf.is_a?(String)
	  raise NetcdfError, "add_offset is not a numeric"
	end
	else
	  cao = 0.0      # assume 0 if not defined
	end
	var = var.to_na( @missval[0]*csf + cao)
	simple_put( (var-cao)/csf, hash )
      end
    end  
  end
end
=end


#------------------------------------------------------------------------------#


class GPhys
  # gtstrm メソッド GPhys に追加
  def strm_function_ncep_with_time(vname="strm")

  end

  # gtstrm 用下請けメソッド

#------------------------------------------------------------------------------#


  # gtstrm メソッド GPhys に追加
  def strm_function_ncep(vname="strm")
    grid = self.grid
    
    lat = grid.axis(0).pos.val              # [90.0, 87.5, ... , -87.5, -90.0]
    ps  = grid.axis(1).pos.val*100          # [100000,  ... , 2000, 1000] milibar => N/m^2
    va_data = self.data
    strm = va_data.val.dup.fill!(0.0)
    mask = va_data.val.get_mask if va_data.val.is_a?(NArrayMiss)
    miss = self.data.get_att("missing_value")[0] if self.data.get_att("missing_value")

    axs = grid.shape[0]
    ays = grid.shape[1]

    pi = NMath::PI
    g  = UNumeric.new(9.81, "m.s-2")
    r  = UNumeric.new(6.37E6,"m") # radius of the planet6.370e6
    cos = Misc::EMath::cos(lat*pi*2.0/360.0)

    keisuu = 2*pi*r*cos/g

    # 台形公式による積分
    (ays-1).downto(0){|j|
      if j == (ays-1) then
	strm[*([true] + [j, false])] = 0.5*(va_data.val[*([true] + [j, false])]+0)*ps[ays-1]
      else
	dp = (ps[j] - ps[j+1])
	strm[*([true] + [j, false])] = 0.5*(va_data.val[*([true] + [j, false])]+va_data.val[*([true] + [j+1, false])])*dp + strm[*([true] + [j+1, false])]
      end
    }
    strm = GPhys.new(grid, VArray.new(strm, va_data, vname)*keisuu)
    strm.units = "kg.s-1"
    strm.set_att("long_name", "Zonal mean mass stream function")
    strm.del_att("title")
    strm.del_att("var_desc")

    return strm
  end

#------------------------------------------------------------------------------#

  # gttheta メソッド GPhys に追加
  def theta(vname)
    grid = self.grid
    
    lon = grid.axis(0).pos.val  # [0.0, 2.5, ... , -359.5, -360.0]
    lat = grid.axis(1).pos.val  # [90.0, 87.5, ... , -87.5, -90.0]
    ps  = grid.axis(2).pos.val  # [10, 20, ... , 1000]
    data = self.data.val        

    axs = grid.shape[0]
    ays = grid.shape[1]
    azs = grid.shape[2]

    # 温位 theta の初期化
#    theta = NArray.float(axs, ays, azs).fill!(0.0)
    theta = NArray.float(axs, ays, azs).fill!(0.0)

=begin
    # theta
    0.upto(axs-1) {|i|
      0.upto(ays-1){|j|
	0.upto(azs-1){|k|
	  t = data[i, j, k]
	  theta[i, j, k] = t*(1000.0/ps[k])**0.288
	}
      }
    }
=end

    # theta
    0.upto(azs-1){|k|
      t = data[true, true, k]
      theta[true, true, k] = t*(1000.0/ps[k])**0.288
    }

    gttheta = GPhys.new(grid, VArray.new(theta).rename(vname))
    return gttheta
  end

#------------------------------------------------------------------------------#

  def bvfreq(vname) # ブラントバイサラ振動数
   
    g = 9.8

    grid = self.grid
    
    lon = grid.axis(0).pos.val  # [0.0, 2.5, ... , -359.5, -360.0]
    lat = grid.axis(1).pos.val  # [90.0, 87.5, ... , -87.5, -90.0]
    ps  = grid.axis(2).pos.val  # [10, 20, ... , 1000]
    data = self.data.val        

    axs = grid.shape[0]
    ays = grid.shape[1]
    azs = grid.shape[2]

    bvfreq = NArray.float(axs, ays, azs).fill!(0.0)


    # brunt-vaisala
    0.upto(azs-1){|k|
    theta = data[true, true, k]
      if k == 0 then
	 bvfreq[true, true, k] = g/theta*(-data[true, true, k+1]+data[true, true, k])/(ps[k+1]-ps[k])
      elsif k == (azs-1) then
	 bvfreq[true, true, k] = g/theta*(-data[true, true, k]+data[true, true, k-1])/(ps[k]-ps[k-1])
      else
	 bvfreq[true, true, k] = g/theta*(-data[true, true, k+1]+data[true, true, k-1])/(ps[k+1]-ps[k-1])
      end
    }
    gtbrunt = GPhys.new(grid, VArray.new(bvfreq).rename(vname))
    return gtbrunt
  end

#------------------------------------------------------------------------------#

  # p_mean メソッド GPhys に追加
  def p_bar(vname)
    theta_bar = self.mean(0, 1).data.val
    grid = self.grid
    
    lon = grid.axis(0).pos.val  # [0.0, 2.5, ... , -359.5, -360.0]
    lat = grid.axis(1).pos.val  # [90.0, 87.5, ... , -87.5, -90.0]
    ps  = grid.axis(2).pos.val  # [10, 20, ... , 1000]
    data = self.data.val        

    axs = grid.shape[0]
    ays = grid.shape[1]
    azs = grid.shape[2]

    # p_bar の初期化
    p_bar = NArray.float(axs, ays, azs).fill!(0.0)

    # p_bar
    0.upto(axs-1) {|i|
      0.upto(ays-1){|j|
	theta = data[i, j, true]
	0.upto(azs-1){|k|
	  if k == 0 then
	    p_bar[i, j, k] = 1.0 - ((ps[k] - (theta[k] - theta_bar[k])/(fd_dp(theta, k, ps[k+1]-ps[k])))/ps[k])
	  elsif k == (azs-1) then
	    p_bar[i, j, k] = 1.0 - ((ps[k] - (theta[k] - theta_bar[k])/(bd_dp(theta, k, ps[k]-ps[k-1])))/ps[k])
	  else
	    p_bar[i, j, k] = 1.0 - ((ps[k] - (theta[k] - theta_bar[k])/(cd_dp(theta, k, ps[k+1]-ps[k-1])))/ps[k])
	  end
	}
      }
    }

    gtpbar = GPhys.new(grid, VArray.new(p_bar).rename(vname))
    return gtpbar
  end

#------------------------------------------------------------------------------#

  def cd_dp(v, k, dp) # center 
    (v[k+1] - v[k-1])/(dp)
  end

#------------------------------------------------------------------------------#

  def fd_dp(v, k, dp) # forward
    (v[k+1] - v[k])/(dp) 
  end

#------------------------------------------------------------------------------#

  def bd_dp(v, k, dp) # backward
    (v[k] - v[k-1])/(dp) 
  end

#------------------------------------------------------------------------------#

  def ape(vname)
    grid = self.grid
    
    kappa = 0.288

    lon = grid.axis(0).pos.val  # [0.0, 2.5, ... , -359.5, -360.0]
    lat = grid.axis(1).pos.val  # [90.0, 87.5, ... , -87.5, -90.0]
    ps  = grid.axis(2).pos.val  # [10, 20, ... , 1000]
    date= grid.axis(3).pos.val  # [10, 20, ... , 1000]
    time  = grid.axis(4).pos.val  # [10, 20, ... , 1000]
    theta = self.data.val        

    axs = grid.shape[0]
    ays = grid.shape[1]
    azs = grid.shape[2]

    # theta の平均値を 3 次元に拡張
    theta_bar1 = self.mean(0, 1).data.val
    theta_bar = NArray.float(axs, ays, azs).fill!(0.0)
    for i in 0..axs-1
      for j in 0..ays-1
	theta_bar[i, j, true] = theta_bar1
      end
    end

    # ps を 3 次元に拡張
    psn = NArray.sint(axs, ays, azs).fill!(0)
    for i in 0..axs-1
      for j in 0..ays-1
	psn[i, j, true] = ps
      end
    end

    # ape
    ape = theta_bar**2*psn**(kappa-1)*((theta-theta_bar)/theta_bar)**2

    gtape = GPhys.new(grid, VArray.new(ape).rename(vname))
    ngtape = gtape.integrate(0)
    ngtape = ngtape.integrate(0)
    ngtape = ngtape.integrate(0)
    return ngtape
  end

#------------------------------------------------------------------------------#

  def stspct_fft(vname)
  #------------------------------------------------------
  # 準備                                                 
  
    # グリッド情報取得   
    grid = self.grid    
    space = grid.axis(0).pos.val  # 空間グリッド
    time  = grid.axis(1).pos.val  # 時間グリッド
    axs = grid.shape[0] # 空間方向のデータ数
    ays = grid.shape[1] # 時間方向のデータ数

  #------------------------------------------------------
  # スペクトル計算

    # 複素空間フーリエ係数 f_k(t) を計算
    w = self.data.val                         # 元データ w(x,t)
    f_k = NArray.complex(axs, ays).fill!(0.0)  # f_k の箱を用意
    for i in 0..(axs-1)
      f_k[i,true] = FFTW.fftw(w[i, true], 1)/(ays/2)
    end

    # f_k(t) の実部 c_k, 虚部 s_k を取得
    c_k = f_k.real
    s_k = f_k.imag

    # c_k,s_k のパワースペクトル pc, ps を求める.
    pwc = NArray.complex(axs, ays).fill!(0.0)  # P_w(c) の箱を用意
    pws = NArray.complex(axs, ays).fill!(0.0)  # P_w(s) の箱を用意
    for i in 0..(ays-1)
      pwc[true,i] = FFTW.fftw(c_k[true, i], 1)/(axs/2) # forward FT
      pws[true,i] = FFTW.fftw(s_k[true, i], 1)/(axs/2) # forward FT
    end

    pc = (pwc).abs**2/2.0
    ps = (pws).abs**2/2.0


    # クアドラチャースペクトルを求める.
    qcs = (pwc.conj*pws).imag/2.0

    # パワースペクトル    
    p_ww = (pc + ps - 2*qcs)/4.0     # 西進波のパワースペクトル
    p_we = (pc + ps + 2*qcs)/4.0     # 東進波のパワースペクトル
    p_ww = p_ww[-1..0, true]         # 符号を逆向きに

  #------------------------------------------------------
  # 軸の設定

    if ays/2 == 0
      nays = (ays+1)/2
    else
      nays = ays/2
    end
    if axs/2 == 0
      naxs = (axs+1)/2
    else
      naxs = axs/2
    end

    axs = 2*naxs - 1

    stspct = NArray.float(axs, nays)
    stspct[0..naxs-2, true] = p_ww[(naxs)..-2, 0..(nays-1)]
    stspct[naxs-1, true] = p_ww[-1, 0..(nays-1)] + p_we[0, 0..(nays-1)] # 波数 0
    stspct[naxs..-1, true] = p_we[1..naxs-1, 0..(nays-1)]
    stspct[true, 0] = 0 # 波数 0, 周波数 0
    stspct[0, true] = 0 # 波数 0, 周波数 0
    data = VArray.new(stspct).rename(vname)

    # 軸の生成
    aks = VArray.new(NArray.sfloat(axs).indgen!(0)-(naxs-1)).rename("wvn")
    aws = VArray.new(NArray.sfloat(nays).indgen!(0)/ays).rename("freq")

    xax = Axis.new().set_pos(aks)
    yax = Axis.new(true).set_cell_guess_bounds(aws).set_pos_to_center

    grid = Grid.new(xax, yax)

    # GPhys オブジェクトを生成
    gtspw = GPhys.new(grid, data)
  
    gtspw.coord(0).set_att('units', '1')
    gtspw.coord(1).set_att('units', '1/DAY')
   
    return gtspw
  end

#------------------------------------------------------------------------------#

  def stspct_coh(vname)
    # グリッド情報取得   
    grid = self.grid    
    space = grid.axis(0).pos.val  # 空間グリッド
    time  = grid.axis(1).pos.val  # 時間グリッド
    axs = grid.shape[0] # 空間方向のデータ数
    ays = grid.shape[1] # 時間方向のデータ数

    # 複素空間フーリエ係数 f_k(t) を計算
    w = self.data.val                         # 元データ w(x,t)
    f_k = NArray.complex(axs, ays).fill!(0.0)  # f_k の箱を用意
    for i in 0..(axs-1)
      f_k[i,true] = FFTW.fftw(w[i, true], 1)/(ays/2)
    end

    # f_k(t) の実部 c_k, 虚部 s_k を取得
    c_k = f_k.real
    s_k = f_k.imag

    # c_k,s_k のパワースペクトル pc, ps を求める.
    pwc = NArray.complex(axs, ays).fill!(0.0)  # P_w(c) の箱を用意
    pws = NArray.complex(axs, ays).fill!(0.0)  # P_w(s) の箱を用意
    for i in 0..(ays-1)
      pwc[true,i] = FFTW.fftw(c_k[true, i], 1)/(axs/2) # forward FT
      pws[true,i] = FFTW.fftw(s_k[true, i], 1)/(axs/2) # forward FT
    end

    pc = (pwc).abs**2/2.0            # C_k のパワースペクトル
    ps = (pws).abs**2/2.0            # S_k のパワースペクトル

    qcs = (pwc.conj*pws).imag/2.0    # クアドラチャースペクトル
    kcs = (pwc.conj*pws).real/2.0    # コスペクトル

    # パワースペクトル    
    p_ww = (pc + ps - 2*qcs)/4.0     # 西進波のパワースペクトル
    p_we = (pc + ps + 2*qcs)/4.0     # 東進波のパワースペクトル
    p_ww = p_ww[-1..0, true]         # 符号逆向き

    # P-R コヒーレンス
    coh_cs = (kcs**2 + qcs**2)/pc/ps # C_k, C_s のコヒーレンス
    coh_pr = NMath::sqrt(1-(1-coh_cs**2)/4*pc*ps/p_ww/p_we)
                                     # P-R コヒーレンス


    # 以下, 空間軸が偶数の場合と奇数の場合で分けなければならない.

    if ays/2 == 0
      nays = (ays+1)/2
    else
      nays = ays/2
    end
    if axs/2 == 0
      naxs = (axs+1)/2
    else
      naxs = axs/2
    end

    axs = 2*naxs - 1

    stspct = NArray.float(axs, nays)
    stspct[0..naxs-2, true] = p_ww[(naxs)..-2, 0..(nays-1)]
    stspct[naxs-1, true] = p_ww[-1, 0..(nays-1)] + p_we[0, 0..(nays-1)] # 波数 0
    stspct[naxs..-1, true] = p_we[1..naxs-1, 0..(nays-1)]
    stspct[true, 0] = 0 # 波数 0, 周波数 0
    stspct[0, true] = 0 # 波数 0, 周波数 0
    data = VArray.new(stspct).rename(vname)

    # 軸の生成
    aks = VArray.new(NArray.sfloat(axs).indgen!(0)-(naxs-1)).rename("wvn")
    aws = VArray.new(NArray.sfloat(nays).indgen!(0)/ays).rename("freq")

    xax = Axis.new().set_pos(aks)
    yax = Axis.new(true).set_cell_guess_bounds(aws).set_pos_to_center

    grid = Grid.new(xax, yax)

    # GPhys オブジェクトを生成
    gtspw = GPhys.new(grid, data)
    p gtspw.max
    return gtspw
  end

#------------------------------------------------------------------------------#

  def stspct_st_tr(vname)
    #!!注意 !!#
    # 上記の点に関して修正すべし.

    # グリッド情報取得   
    grid = self.grid    
    space = grid.axis(0).pos.val  # 空間グリッド
    time  = grid.axis(1).pos.val  # 時間グリッド
    axs = grid.shape[0] # 空間方向のデータ数
    ays = grid.shape[1] # 時間方向のデータ数

    # 複素空間フーリエ係数 f_k(t) を計算
    w = self.data.val                         # 元データ w(x,t)
    f_k = NArray.complex(axs, ays).fill!(0.0)  # f_k の箱を用意
    for i in 0..(axs-1)
      f_k[i,true] = FFTW.fftw(w[i, true], 1)/(ays/2)
    end

    # f_k(t) の実部 c_k, 虚部 s_k を取得
    c_k = f_k.real
    s_k = f_k.imag

p c_k.abs.max
p s_k.abs.max
   
    # c_k,s_k のパワースペクトル pc, ps を求める.
    pwc = NArray.complex(axs, ays).fill!(0.0)  # P_w(c) の箱を用意
    pws = NArray.complex(axs, ays).fill!(0.0)  # P_w(s) の箱を用意
    for i in 0..(ays-1)
      pwc[true,i] = FFTW.fftw(c_k[true, i], 1)/(axs/2) # forward FT
      pws[true,i] = FFTW.fftw(s_k[true, i], 1)/(axs/2) # forward FT
    end

    pc = (pwc).abs**2/2.0
    ps = (pws).abs**2/2.0

    qcs = (pwc.conj*pws).imag/2.0    # クアドラチャースペクトル
    kcs = (pwc.conj*pws).real/2.0    # コスペクトル

    # パワースペクトル    
    p_ww = (pc + ps - 2*qcs)/4.0     # 西進波のパワースペクトル
    p_we = (pc + ps + 2*qcs)/4.0     # 東進波のパワースペクトル

    # 定在波, 移動波のパワースペクトル
    var = ((pc-ps)**2)/4 + kcs**2
    p_st = NMath::sqrt(var)          # 定在波のパワースペクトル
    p_tr = (p_ww + p_we) - p_st      # 移動波のパワースペクトル
    p_tr = p_tr[-1..0, true]

    # 以下, 空間軸が偶数の場合と奇数の場合で分けなければならない.

    if ays/2 == 0
      nays = (ays+1)/2
    else
      nays = ays/2
    end
    if axs/2 == 0
      naxs = (axs+1)/2
    else
      naxs = axs/2
    end

    axs = 2*naxs - 1

    stspct = NArray.float(axs, nays)
    stspct[0..naxs-2, true] = p_tr[(naxs)..-2, 0..(nays-1)]
    stspct[naxs-1, true] = p_tr[-1, 0..(nays-1)] + p_st[0, 0..(nays-1)] # 波数 0
    stspct[naxs..-1, true] = p_st[1..naxs-1, 0..(nays-1)]
    stspct[true, 0] = 0 # 波数 0, 周波数 0
    stspct[0, true] = 0 # 波数 0, 周波数 0
    data = VArray.new(stspct).rename(vname)

    # 軸の生成
    aks = VArray.new(NArray.sfloat(axs).indgen!(0)-(naxs-1)).rename("wvn")
    aws = VArray.new(NArray.sfloat(nays).indgen!(0)/ays).rename("freq")

    xax = Axis.new().set_pos(aks)
    yax = Axis.new(true).set_cell_guess_bounds(aws).set_pos_to_center

    grid = Grid.new(xax, yax)

    # GPhys オブジェクトを生成
    gtspw = GPhys.new(grid, data)
    p gtspw.max
    return gtspw
  end

#------------------------------------------------------------------------------#

  def make_sfc_mask(p_gphys)

    raise ArgumentError,"Self GPhys rank must have dimentions, longitude, latitude, time ." if ! p_gphys.is_a?(GPhys) 
    raise ArgumentError,"Arg is not GPhys." if ! p_gphys.is_a?(GPhys) 
    raise ArgumentError,"Arg's dimention is not 2nd order" if ! p_gphys.rank==2

    missval = 6.0e24 # missing_value

    p_data = p_gphys.data.val
    data = self.data.val.dup
    grid = self.grid

    ps  = grid.axis(2).pos.val  # pressure
    axs = grid.shape[0]
    ays = grid.shape[1]
    azs = grid.shape[2]

    if self.rank == 3
      0.upto(azs-1) do |z|
	data[true,true,z] = ps[z]
	data[true,true,z] = (data[true,true,z].ge p_data)
      end
    else
      0.upto(azs-1) do |z|
	data[true,true,z,true] = ps[z]
	data[true,true,z,true] = (data[true,true,z,true].ge p_data/100)
      end
    end
    mask = data*missval
    gtmask = GPhys.new(grid, VArray.new(mask,
					{"long_name"=>"surface mask", 
					  "missing_value"=>missval, 
					  "valid_max"=>1.0e24, 
					  "valid_min"=>-6.0e36},
					"sfc_mask")
		       )
    return gtmask
  end

#------------------------------------------------------------------------------#

  # gphys オブジェクトを nc ファイルに保存するメソッド
  def save(output, global_attr=nil, vname=nil, var_attr=nil, history=$history, new_vname=nil)

    if !vname
      vname = self.data.name
    end

    if File.exist?(output)
      nc = NetCDF.open(output, 'a+')
    else
      nc = NetCDF.open(output, 'w')
    end
    GPhys::IO.write(nc, self)
    nc.close

    nc = NetCDF.open(output, "a+")
    nc.redef                      # define mode に移行

    # 大域属性付加
    global_attr.to_a.each {|attr| 
      nc.put_att(attr[0], attr[1]) if attr[1]
    } if global_attr
    # 変数名変更
    if new_vname && vname
      nc.var(vname).name=new_vname 
      vname = new_vname
    end
    # 変数属性付加
    if var_attr
      var_attr.to_a.each {|attr|
	nc.var(vname).put_att(attr[0], attr[1])
      }   
    end  
    # history 付加
    now_hist = global_attr["history"]
    if now_hist then
      hist = now_hist + "\n" + history
      nc.put_att("history", hist)
    else
      nc.put_att("history", history)
    end

    nc.close

    return "saved as #{output}"

  end

  # 座標, データ操作ログ記録関数の定義 (Grid.lost_axes 参照.) <やまだ由さんのルーチンを拝借>
  def add_lost_axes( lost )
    GPhys.new( @grid.copy.add_lost_axes( lost.to_a ), @data)
  end
  def set_lost_axes( lost )
    GPhys.new( @grid.copy.set_lost_axes( lost.to_a ), @data)
  end

  # rename 再定義 (gphys class を返す)
  def rename(name)
    GPhys.new(@grid,@data.copy.rename(name))
  end

  # attribute 再定義 (gphys class 返す)
  def set_att(name,val)
    GPhys.new(@grid,@data.copy.set_att(name,val))
  end

end

def global_attr(ncfile)
  global_attr = Hash.new

  nc = NetCDF.open(ncfile, "r")
  nc.att_names.each do |attr| 
    global_attr[attr] = nc.att(attr).get
  end
  return global_attr
end



# 同じ次元構造をもつ複数の nc ファイルを指定した際に
# 平均したオブジェクトを返すメソッド
# ファイル名に 1990-10 のような文字が入ってた場合, 
# 月平均データを元に季節平均を計算. (日にちの重みを考慮)
def mean_gphys(file, var)

  gphys_file = Array.new()
  dateary =    Array.new()
  
  leap_year =      [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
  non_leap_year  = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] 

  file.each_index{|i|
    temp = GPhys::NetCDF_IO.open(file[i], var)
    ym = file[i].scan(/(\d\d\d\d)-(\d\d)/)
    if ym[0]
      year = ym[0][0].to_i
      month = ym[0][-1].to_i
      leap_flag = Date::new(year).leap?
      if leap_flag
	date = leap_year[month - 1]
      else
	date = non_leap_year[month - 1]
      end
      gphys_file << temp*date
      dateary << date
    else
      temp = GPhys::NetCDF_IO.open(f, var)
      gphys_file << temp
    end
  }

  sum = gphys_file[0]

  1.upto(gphys_file.size-1) {|f|
    sum += gphys_file[f] 
  }

  if dateary
    p sumdate = NArray.to_na(dateary).sum.to_f
    gphys = sum/sumdate
  else
    p gphys = sum/file.size.to_f
  end
  return gphys
end

#  旧 mean_gphys 足した数で平滑するだけ.
#
#def mean_gphys(file, var)
#
#  gphys_file = Array.new()
#
#  file.each {|f|
#  
#    temp = GPhys::NetCDF_IO.open(f, var)
#    gphys_file << temp
#
#  }
#
#  sum = gphys_file[0]#
#
#  1.upto(gphys_file.size-1) {|f|
#    sum += gphys_file[f] 
#  }
#
#  p gphys = sum/file.size.to_f
#  return gphys
#end

def ask_overwrite(output)
  
  if File.exist?(output) then
    
    print "\n"
    print %Q{the file name "#{output}" is already exists current directory.\n}
    print "Do you sure over write? (y/n)    "

    loop{
      flag = STDIN.gets.chomp

      if flag == "y" then
	break
      elsif flag == "n" then
	p "Bye."
	exit
      else
	print "please answer yes or no. (y/n)     "
      end

    }
    
  end
end