Module: NumRu::GAnalysis::MetZ

Defined in:: lib/numru/ganalysis/met_z.rb

Overview

Meterological analysis regarding vertical section, integration, etc.

Class Method Summary collapse

.mass_strm_any(v, ps, w, wcoord, vs = nil, ws = nil) ⇒ Object

mass stream function on any vertical coordinate.
.mass_strm_p(v, ps, pcoord = nil, vs = nil) ⇒ Object

Derive the mass stream function in the pressure coordinate.

Class Method Details

.mass_strm_any(v, ps, w, wcoord, vs = nil, ws = nil) ⇒ `Object`

mass stream function on any vertical coordinate

Similar to mass_strm_p, but it supports representation to have an arbitrary physical quantity, such as potential temperature, as the vertical coordinate (instead of pressure).

Applicable both to pressure- and sigma-coordinate input data

ARGUMENTS

v [GPhys] : meridional wind with a vertical dimension (p or sigma) It must have a latitudinal dimension too. Longitudinal and time dimensions are optional. If it has a longitudinal dimension, zonal mean is taken. The order of the dimensions is not restricted.
ps [GPhys] : surface pressure. Its must have the same grid as v but for the vertical dimension (ps.rank must be v.rank-1)
w [GPhys] : Grid-point values (at the same points as v) of the quantity used to represent the vertical coordinate. Its shape must be the same as that of v, as a matter of course.
wcoord [1D VArray] : Output vertical coordinate. It must have the same units as w.
vs [nil(default) or GPhys]: vs is not needed (neglected) when v has a sigma coordinate. It is an optional parameter to specify the surface values of v, when it is in the pressure coordinate. vs can be omitted (nil), even when v has a pressure coordinate; in that case, vs is set by interpolating v if ps is within the p range of v (e.g. when ps<=1000hPa), or it is naively extended (using the bottom values of v) if ps is out of the range (e.g. when ps>1000hPa). In other words, the current implementation assumes that v is available below the surface, as is customary for reanalysis data.
ws [nil(default) or GPhys]: same as vs but for the surface value of w.

Raises:

(ArgumentError)

# File 'lib/numru/ganalysis/met_z.rb', line 182

def mass_strm_any(v, ps, w, wcoord, vs=nil, ws=nil)

  pascal = Units["Pa"]
  grav = Met.g.to_f

  #< check >

  raise(ArgumentError,"v.shape != w.shape")  if v.shape != w.shape
  raise(ArgumentError,"ps.rank != v.rank-1")  if ps.rank != v.rank-1
  raise(ArgumentError,"w.units !~wcoord.units") if w.units !~ wcoord.units

  #< preprare data >

  if zdim = Met.find_prs_d(v)   # substitution, not comparison
    # has a pressure coordinate
    pcv = v.coord(zdim)   # v's p coord
    pcv_val = pcv.val
    v_val = v.val             # should be NArray or NArrayMiss
    v_val = v_val.to_na if v_val.is_a?(NArrayMiss)
    w_val = w.val             # should be NArray or NArrayMiss
    w_val = w_val.to_na if w_val.is_a?(NArrayMiss)
    if pcv_val[0] > pcv_val[-1]
      # reverse the p coordinate to the increasing order
      pcv_val = pcv_val[-1..0]
      v_val = v_val[ *([true]*zdim + [-1..0,false]) ]
      w_val = w_val[ *([true]*zdim + [-1..0,false]) ]
    end

    pcv_val = pcv.units.convert2(pcv_val, pascal) if pcv.units!=pascal
    pcv_over_g = pcv_val / grav

    ps_val = ps.val
    ps_val = ps_val.to_na if ps_val.is_a?(NArrayMiss)
    ps_val = ps.units.convert2(ps_val, pascal) if ps.units!=pascal
    ps_over_g = ps_val / grav

    vs_val = vs && vs.val   # nil (default) or vs.val (if vs is given)
    vs_val = vs_val.to_na if vs_val.is_a?(NArrayMiss)

    ws_val = ws && ws.val   # nil (default) or ws.val (if ws is given)
    ws_val = ws_val.to_na if ws_val.is_a?(NArrayMiss)

    v_val, p_over_g, nzbound = GPhys.c_cap_by_boundary(v_val, zdim, 
                                     pcv_over_g, true, ps_over_g, vs_val)

    w_val, p_over_g, nzbound = GPhys.c_cap_by_boundary(w_val, zdim, 
                                     pcv_over_g, true, ps_over_g, ws_val)

  elsif zdim = SigmaCoord.find_sigma_d(v)  # substitution, not comparison
    # has a sigma coordnate
    sig = v.coord(zdim)
    nz = sig.length
    nzbound = nil
    ps = ps.convert_units(pascal) if ps.units != pascal
    sig_val = sig.val
    v_val = v.val    # should be NArray, not NArrayMiss (coz sigma)
    w_val = w.val
    p_over_g = SigmaCoord.sig_ps2p(ps.val/grav, sig_val, zdim)
  else
    raise ArgumentError, "v does not have a p or sigma coordinate."
  end

  #< cumulative vertical integration >
 
  wc_val = wcoord.val
  if wc_val[0] > wc_val[-1]
    # change it to the increasing order
    wc_val = wc_val[-1..0]
    wcoord = wcoord.copy.replace_val(wc_val)
  end

  rho_v_cum = GPhys.c_cum_integ_irreg(v_val, p_over_g, zdim, nzbound,
                                   wc_val, w_val)

  #< zonal mean & latitudinal factor >

  lam, phi, lond, latd = Planet.get_lambda_phi(v, false)

  if latd.nil?
    raise(ArgumentError, "v appears not having a latitudinal dimension")
  end
  if lond
    rho_v_cum = rho_v_cum.mean(lond)
    latd -= 1 if lond<latd
  end

  a_cos = NMath.cos(phi.val) * ( 2 * Math::PI * Planet.radius.to_f )
  latd.times{a_cos.newdim!(0)}
  (rho_v_cum.rank - latd -1).times{a_cos.newdim!(-1)}

  mstrm_val = rho_v_cum * a_cos

  #< make a GPhys >

  axes = Array.new
  for d in 0...v.rank
    case d
    when lond
      # lost by zonal mean
    when zdim
      wax = Axis.new().set_pos(wcoord)
      axes.push(wax)
    else
      axes.push(v.axis(d).copy)   # kept
    end
  end
  grid = Grid.new( *axes )

  units = Units["kg.m-1"]    # p/g*a : Pa / (m.s-2) * m = kg.m-1
  units *= v.units
  mstrm_va = VArray.new(mstrm_val, {"long_name"=>"mass stream function",
                          "units"=>units.to_s}, "mstrm")
  mstrm = GPhys.new(grid, mstrm_va)
  mstrm
end

.mass_strm_p(v, ps, pcoord = nil, vs = nil) ⇒ `Object`

Derive the mass stream function in the pressure coordinate

Applicable both to pressure- and sigma-coordinate input data (the output is always on the pressure coordinate).