Module: NumRu::GAnalysis::Met

consider_ice() click to toggle source

whether or not ice is considered in the water phase change

# File ../../lib/numru/ganalysis/met.rb, line 458
def consider_ice
  @@consider_ice
end

consider_ice=(t_or_f) click to toggle source

set whether or not ice is considered in the water phase change

# File ../../lib/numru/ganalysis/met.rb, line 463
def consider_ice=(t_or_f)
  @@consider_ice=t_or_f
end

convert_units2Pa(prs) click to toggle source

Convert units into Pa. To deal with old versions of NumRu::Units that do not support “millibar”.

ARGUMENT

prs [GPhys and UNumeric]

# File ../../lib/numru/ganalysis/met.rb, line 398
def convert_units2Pa(prs)
  pa = Units["Pa"]
  un = prs.units
  if ((un =~ pa) and !(un === pa))
    prs = convert_units(prs, pa)
  elsif un.to_s=="millibar"
    if UNumeric === prs
      ret = UNumeric[prs.val*100, "Pa"]
    else
      ret = prs*100
      ret.units = "Pa"
    end
    ret
  else
    prs
  end
end

df_dx_vialogscale(f, x, dim) click to toggle source

# File ../../lib/numru/ganalysis/met.rb, line 440
def df_dx_vialogscale(f, x, dim)
  z = Misc::EMath.log(x)
  if GPhys === f
    mdl = NumRu::GPhys::Derivative
    dfdz = mdl.threepoint_O2nd_deriv(f, dim, mdl::LINEAR_EXT, z)
  else
    mdl = NumRu::Derivative
    dfdz = mdl.threepoint_O2nd_deriv(f, z, dim, mdl::LINEAR_EXT)
  end
  dfdz / x
end

e2q(e,prs=nil) click to toggle source

water vapor pressure -> specific humidity

ARGUMENTS

e: water vapor pressure
prs: pressure

RETURN VALUE

q: specific humidity

# File ../../lib/numru/ganalysis/met.rb, line 567
def e2q(e,prs=nil)
  prs = get_prs(e) if !prs
  prs = convert_units(prs,e.units)
  rratio = R / Rv
  q = rratio * e / (prs-(1-rratio)*e) 
  q.name = "q"
  q.long_name = "specific humidity"
  q
end

e2r(e,prs=nil) click to toggle source

water vapor pressure -> mixing ratio

ARGUMENTS

e: water vapor pressure
prs: pressure

RETURN VALUE

r: mixing ratio

# File ../../lib/numru/ganalysis/met.rb, line 531
def e2r(e,prs=nil)
  prs = get_prs(e) if !prs
  prs = convert_units(prs,e.units)
  rratio = R / Rv
  r = rratio * e / (prs-e)
  r.name = "r"
  r.long_name = "mixing ratio"
  r
end

e_sat(temp) click to toggle source

Calculates saturation water vapor pressure using enhanced

ARGUMENTS

temp: temperature

RETURN VALUE

es: saturation water vapor pressure

# File ../../lib/numru/ganalysis/met.rb, line 695
def e_sat(temp)

  ice = @@consider_ice && ( temp.lt(@@ice_thres) )
  #ice = ice.to_na
  water = !@@consider_ice || ( (ice==true||ice==false) ? !ice : ice.not)

  if water
    es = e_sat_water(temp)
  end

  case ice
  when true
    es = e_sat_ice(temp)
  when NArray, NArrayMiss 
    es[ice] = e_sat_ice(temp,ice).val[ice]
  end
  es.name = "e_sat"
  es.long_name = "e_sat"
  es
end

e_sat_bolton(temp) click to toggle source

Bolton formula for saturation water vapor pressure against water

ARGUMENTS

temp: temperature

RETURN VALUE

es: saturation water vapor pressure [Pa]

# File ../../lib/numru/ganalysis/met.rb, line 601
def e_sat_bolton(temp)
  tempC = temp.convert_units("degC")
  es = UNumeric[6.112e2,"Pa"] * 
    Misc::EMath.exp( 17.67 * tempC / (tempC + UNumeric[243.5,"degC"] ) )
  es.name = "e_sat"
  es.long_name = "e_sat_water bolton"
  es
end

e_sat_emanuel_ice(temp, mask=nil) click to toggle source

Saturation water vapor pressure against ice.

Emanuel (1994) eq.(4.4.15)

ARGUMENTS

temp: temperature

RETURN VALUE

es: saturation water vapor pressure [Pa]

# File ../../lib/numru/ganalysis/met.rb, line 646
def e_sat_emanuel_ice(temp, mask=nil)
  es = temp.copy
  tempK = temp.convert_units("K").val   # units removed
  tempK = tempK[mask] if mask
  e = 23.33086 - 6111.72784/tempK               + 0.15215 * Misc::EMath.log(tempK)
  if !mask
    es.replace_val( Misc::EMath.exp(e) * 100.0 )
  else
    es[false] = 0
    es[mask] = Misc::EMath.exp(e) * 100.0
  end
  es.units = "Pa"
  es.name = "e_sat_ice"
  es.long_name = "e_sat_ice emanuel"
  es
end

Also aliased as: e_sat_ice

e_sat_emanuel_water(temp, mask=nil) click to toggle source

saturation water vapor pressure against ice.

Emanuel (1994) eq.(4.4.15)

ARGUMENTS

temp: temperature

RETURN VALUE

es: saturation water vapor pressure [Pa]

# File ../../lib/numru/ganalysis/met.rb, line 619
def e_sat_emanuel_water(temp, mask=nil)
  es = temp.copy
  tempK = temp.convert_units("K").val   # units removed
  tempK = tempK[mask] if mask
  e = 53.67957 - 6743.769/tempK               - 4.8451 * Misc::EMath.log(tempK)
  if !mask
    es.replace_val( Misc::EMath.exp(e) * 100.0 )
  else
    es[false] = 0
    es[mask] = Misc::EMath.exp(e) * 100.0
  end
  es.units = "Pa"
  es.name = "e_sat_ice"
  es.long_name = "e_sat_ice emanuel"
  es
end

find_prs_d(gphys, error=nil) click to toggle source

def frontogenesis_eulerian(theta, u, v, w=nil, fstodr=true )
  thx, thy = GAnalysis::Planet.grad_s( theta )
  uthxx,  uthxy =  GAnalysis::Planet.grad_s( u*thx )
  vthyx,  vthyy =  GAnalysis::Planet.grad_s( v*thy )
  va = GAnalysis::Planet.weight_tanphi( v, 1, -1 )
  frgf = - (uthxx+vthyx)*thx - (uthxy+vthyy)*thy
  frgf.name = "thgrd_tend"
  frgf.long_name = "Eluerian grad-theta tendency"
  if w 
    zdim = 2
    if (wun=w.units) !~ (ztun = theta.coord(zdim).units * Units["s-1"])
      raise "w in #{wun} is inconsistent with the vertical coordinate of theta in #{ztun}"
    else
      w = w.convert_units(ztun)    # For example, Pa/s -> hPa/s
    end
    z = theta.axis(zdim).to_gphys
    if z.units =~ Units["Pa"]
      thz = df_dx_vialogscale(theta, z, zdim)
    else
      thz = df_dx(theta, z, zdim)
    end
    wthzx,  wthzy =  GAnalysis::Planet.grad_s( w*thz ) 
    frgf -= wthzx*thx + wthzy*thy
  end
  if fstodr
    frgf /= (thx**2 + thy**2).sqrt 
  else
    frgf *= 2
  end
  frgf
end

# Find a pressure coordinate in a GPhys object # # ARGUMENT # * gphys [GPhys] # * error [nil/false or true] change the behavior if a # pressure coordinate is not found. Default: returns nil; # if error is true, an exception is raised. # RETURN VALUE # * Integer to indicate the dimension of the pressure coordinate, # or nil if not found by default (see above)

# File ../../lib/numru/ganalysis/met.rb, line 352
def find_prs_d(gphys, error=nil)
  pa = Units.new("Pa")
  (gphys.rank-1).downto(0) do |d|
    un = gphys.axis(d).pos.units
    if ( un =~ pa or un.to_s=="millibar" )
      return(d)
    end
  end
  if error
    raise("Could not find a pressure coordinate.")
  else
    nil
  end
end

frontogenesis_func(theta, u, v, w=nil, fstodr=true, full_adv = true) click to toggle source

Adiabatic frontogenesis function over the sphere. – D/Dt(|gradH theta|) or D/Dt(|gradH theta|^2), where gradH express the horizontal component of gradient.

if full_adv is true (default),

D/Dt(|gradH theta|) = 
  [ -(ux-va)*thetax^2 - (vx+uy)*thetax*thetay - vy*thetay^2
    - (wx*thetax + wy*thetay)*theta_z ]
  / |gradH theta|

or else,

(\del/\del t + u gradx + v grady )(|gradH theta|) = 
  [ -(ux-va)*thetax^2 - (vx+uy)*thetax*thetay - vy*thetay^2
    - (w*theta_z)_x*thetax - (w*theta_z)_y*thetay ]
  / |gradH theta|

Here, the 2nd line (vertical advection) is optional;

vx, vy: gradH v; [thetax, thetay] = gradH theta;
ux, uy: cos_phi * gradH (u/cos_phi)

va = v*tan_phi/a (a=radius). z and w is the vertical coordinate and the lagrangian “velocity” in that coordinate — Typically they are p and omega, or log-p height and log-p w.

This formulation is adiabatic; the diabatic heating effect can be easily included if needed.

ARGUMENTS

theta [GPhys] : potential temperature
u [GPhys] : zonal wind
v [GPhys] : meridional wind
w [nil (default) or GPhys] : (optional) “vertical wind”, which must be dimensionally consistent with the vertical coordiante (e.g., omega for the pressure coordinate). If w is given, the vertical cooridnate is assumed to be the 3rd one (dim=2).
fstodr [true (default) or false] (optional) if true D/Dt(|NablaH theta|) returned; if false D/Dt(|NablaH theta|^2) is returned.
full_adv [true (default) or false] : whether to calculate full lagrangian tendency or lagragian tendency only in horizontal direction

RETURN VALUE

frontogenesis function [GPhys]

# File ../../lib/numru/ganalysis/met.rb, line 268
def frontogenesis_func(theta, u, v, w=nil, fstodr=true, full_adv = true)
  thx, thy = GAnalysis::Planet.grad_s( theta )
  ux = GAnalysis::Planet.grad_sx( u )
  uy = GAnalysis::Planet.grad_sy_cosphifact( u, -1 )
  vx, vy =  GAnalysis::Planet.grad_s( v )
  va = GAnalysis::Planet.weight_tanphi( v, 1, -1 )
  frgf = - (ux-va)*thx*thx - (vx+uy)*thx*thy - vy*thy*thy
  frgf.name = "frgen"
  frgf.long_name = "frontogenesis function"
  if w 
    zdim = 2
    if (wun=w.units) !~ (ztun = theta.coord(zdim).units * Units["s-1"])
      raise "w in #{wun} is inconsistent with the vertical coordinate of theta in #{ztun}"
    else
      w = w.convert_units(ztun)    # For example, Pa/s -> hPa/s
    end
    z = theta.axis(zdim).to_gphys
    if z.units =~ Units["Pa"]
      thz = df_dx_vialogscale(theta, z, zdim)
    else
      thz = df_dx(theta, z, zdim)
    end
    if full_adv
      # full lagragian tendency of theta-gradient strength
      wx,  wy =  GAnalysis::Planet.grad_s( w ) 
      frgf -= (wx*thx + wy*thy)*thz
    else
      # lagragian tendency only in horizontal direction
      wthzx,  wthzy =  GAnalysis::Planet.grad_s( w*thz ) 
      frgf -= wthzx*thx + wthzy*thy
    end
  end
  if fstodr
    frgf /= (thx**2 + thy**2).sqrt 
  else
    frgf *= 2
  end
  frgf
end

g() click to toggle source

Returns gravity acceleration in the module (default: 9.8 ms-1).

# File ../../lib/numru/ganalysis/met.rb, line 60
def g
  @@g.dup
end

get_prs(gphys) click to toggle source

Find and return a pressure coordinate in a GPhys object

ARGUMENT

gphys [GPhys]

RETURN VALUE

pressure in a 1D GPhys object created from the pressure axis in the GPhys object or a UNumeric object if the pressure coordinated has been deleted by subsetting.

# File ../../lib/numru/ganalysis/met.rb, line 376
def get_prs(gphys)
  begin
    prs = gphys.axis( find_prs_d(gphys, true) ).to_gphys
  rescue
    regexp = /([\d\.]+) *(millibar|hPa)/
    cand = gphys.lost_axes.grep(regexp)
    if (str=cand[0])   # substitution, not ==
      regexp =~ str
      hpa = $1.to_f
      prs = UNumeric[hpa*100,"Pa"]
    else
      raise "No pressure axis was found"
    end
  end
  prs
end

ice_thres() click to toggle source

the threshold temperature for liquid/ice-phase treatment

# File ../../lib/numru/ganalysis/met.rb, line 468
def ice_thres
  @@ice_thres
end

ice_thres=(temp) click to toggle source

set the threshold temperature for liquid/ice-phase treatment (default: O degC)

# File ../../lib/numru/ganalysis/met.rb, line 473
def ice_thres=(temp)
  @@ice_thres=temp
end

interpolate_onto_theta(gphys, theta, theta_levs) click to toggle source

Interpolate onto the potential temperature coordinate

ARGUMENTS

gphys [GPhys] a gphys object that have a pressure dimension
theta [GPhys] potential temperature defined on the same grid as gphys
theta_vals : 1D NArray or Array

# File ../../lib/numru/ganalysis/met.rb, line 125
def interpolate_onto_theta(gphys, theta, theta_levs)
  theta_levs = NArray[*theta_levs].to_f if Array === theta_levs
  th_crd = VArray.new( theta_levs, 
         {"units"=>"K", "long_name"=>"potential temperature"}, "theta" )
  gphys.set_assoc_coords([theta])
  pdnm = gphys.coord(find_prs_d(gphys)).name
  gphys.interpolate(pdnm=>th_crd)
end

lat(temp) click to toggle source

temperature –> latent heat [J.kg-1]

good for -100<T<50

ARGUMENTS

temp: temperature

RETURN VALUE

lat: latent heat

# File ../../lib/numru/ganalysis/met.rb, line 586
def lat(temp)
  tempK = temp.convert_units("K")
  lat = Lat0*(T0/tempK)**(0.167+tempK.val*3.67E-4) 
  lat.name = "L"
  lat.long_name = "Latent heat"
  lat
end

pv_on_p(u, v, theta) click to toggle source

Derive Ertel’s potential vorticity on the pressure coordinate

ARGUMENTS

u [GPhys] : zonal wind on pressure coordinate (u, v, and theta must share same coordinates)
v [GPhys] : meridional wind on pressure coordinate
theta [GPhys] : potential temperature on pressure coordinate

RETURN VALUE

potential vorticity [GPhys] on the same grid as the inputs

# File ../../lib/numru/ganalysis/met.rb, line 173
def pv_on_p(u, v, theta)
  up,vp,thp = del_ngp(u,v,theta)    # -g del/del p
  pv = Planet.absvor_s(u, v) * thp               - vp * Planet.grad_sx(theta)               + up * Planet.grad_sy(theta)
  pv.long_name = "potential vorticity"
  pv.name = "PV"
  pv.units = pv.units.reduce5    # express in the MKS fundamental units
  pv
end

pv_on_theta(u, v, theta, theta_levs) click to toggle source

Derive Ertel’s potential vorticity on the theta (isentropic) coordinate

ARGUMENTS

u [GPhys] : zonal wind on pressure coordinate (u, v, and theta must share same coordinates)
v [GPhys] : meridional wind on pressure coordinate
theta [GPhys] : potential temperature on pressure coordinate
theta_levs [NArray or Array] : one dimensional array of potential temperature values (Kelvin) on which PV is derived.

RETURN VALUE

potential voticity [GPhys] on a theta coordinate, where levels are set to theta_levs

# File ../../lib/numru/ganalysis/met.rb, line 148
def pv_on_theta(u, v, theta, theta_levs)
  sigi = GAnalysis::Met.sigma_inv(theta)
  uth = GAnalysis::Met.interpolate_onto_theta(u, theta, theta_levs)
  vth = GAnalysis::Met.interpolate_onto_theta(v, theta, theta_levs)
  sigith = GAnalysis::Met.interpolate_onto_theta(sigi, theta, theta_levs)
  avorth = GAnalysis::Planet.absvor_s(uth,vth)
  pv = avorth*sigith
  pv.long_name = "potential vorticity"
  pv.name = "PV"
  pv
end

q2e(q,prs=nil) click to toggle source

specific humidity -> water vapor pressure

ARGUMENTS

q: specific humidity
prs: pressure

RETURN VALUE

e: water vapor pressure

# File ../../lib/numru/ganalysis/met.rb, line 549
def q2e(q,prs=nil)
  prs = get_prs(q) if !prs
  prs = convert_units2Pa(prs)
  rratio = R / Rv
  e = prs*q/(rratio+(1-rratio)*q)     # water vapor pertial pressure
  e.name = "e"
  e.long_name = "water vapor pressure"
  e
end

q2r(q) click to toggle source

specific humidity -> mixing ratio

ARGUMENTS

q: specific humidty

RETURN VALUE

r: mixing ratio

# File ../../lib/numru/ganalysis/met.rb, line 484
def q2r(q)
  r = q/(1.0-q)
  r.name = "r"
  r.long_name = "mixing ratio"
  r
end

r2e(r,prs=nil) click to toggle source

water vapor mixing ratio -> water vapor pressure

ARGUMENTS

r: water vapor mixing ratio
prs: pressure

RETURN VALUE

e: water vapor pressure

# File ../../lib/numru/ganalysis/met.rb, line 513
def r2e(r,prs=nil)
  prs = get_prs(r) if !prs
  prs = convert_units2Pa(prs)
  rratio = R / Rv
  e = prs*r/(rratio+r)     # water vapor pertial pressure
  e.name = "e"
  e.long_name = "water vapor pressure"
  e
end

r2q(r) click to toggle source

mixing ratio -> specific humidity

ARGUMENTS

r: mixing ratio

RETURN VALUE

q: specific humidty

# File ../../lib/numru/ganalysis/met.rb, line 498
def r2q(r)
  q = r/(1.0+r)
  q.name = "q"
  q.long_name = "specific humidity"
  q
end

rh2e(rh,temp) click to toggle source

relative humidity -> water vapor pressure

ARGUMENTS

rh: relative humidity
temp: temperature

RETURN VALUE

e: water vapor pressure

# File ../../lib/numru/ganalysis/met.rb, line 724
def rh2e(rh,temp)
  es = e_sat(temp)
  rh = rh.convert_units("1")
  e = es * rh
  e.name = "e"
  e.long_name = "water vapor pressure"
  e
end

set_e_sat_water(formula=nil) click to toggle source

Selector of the formulat to compute saturation water vapor pressure against water (default: Bolton)

ARGUMENTS

formula: nil(default), “bolton”, “emanuel”

RETURN VALUE

nil

# File ../../lib/numru/ganalysis/met.rb, line 671
def set_e_sat_water(formula=nil)
  case formula
  when nil,"bolton"
    alias :e_sat_water :e_sat_bolton
  when "emanuel"
    alias :e_sat_water :e_sat_emanuel_water
    module_function :e_sat_water
  end
  nil
end

set_g(g) click to toggle source

Sets gravity acceleration in the module (default: 9.8 ms-1).

ARGUMENT

g [UNumeric]

EXAMPLE

Met::set_g( UNumeric[9.7,"m.s-2"] )

RETURN VALUE

The argument g

# File ../../lib/numru/ganalysis/met.rb, line 54
def set_g(g)
  @@g = g   # e.g., un[9.8,"m.s-2"]
end

sigma_inv(theta, dim=nil, prs=nil) click to toggle source

Inverse of the “sigma” density in the theta coordinate:

ARGUMENTS

theta [GPhys or VArray] potential temperature
dim [Integer] : the pressure dimension. If theta is a GPhys, this argument can be omitted, in which case a pressure dimension is searched internally by find_prs_d.
prs [GPhys or VArray] : the pressure values. Internally searched if omitted.

RETURN VALUE

1 / sigma = -g dtheta/dp [GPhys]

# File ../../lib/numru/ganalysis/met.rb, line 104
def sigma_inv(theta, dim=nil, prs=nil)
  dim = find_prs_d(theta) if !dim
  prs = get_prs(theta) if !prs
  prs = convert_units2Pa(prs)
  #dtheta_dp = df_dx(theta, prs, dim)
  dtheta_dp = df_dx_vialogscale(theta, prs, dim)
  sig_inv = dtheta_dp * (-@@g)
  if sig_inv.respond_to?(:long_name)  # VArray or GPhys
    sig_inv.name = "sig_inv"
    sig_inv.long_name = "1/sigma"
  end
  sig_inv
end

temp2theta(temp, prs=nil) click to toggle source

Convert temperature into potential temperature

ARGUMENTS

temp [UNumeric or GPhys or VArray, which supports method “units”] : temperature
prs [UNumeric or GPhys or VArray, which supports method “units”] : pressure

RETURN VALUE

potential temperature [UNumeric or GPhys or VArray,…] in Kelvin

# File ../../lib/numru/ganalysis/met.rb, line 76
def temp2theta(temp, prs=nil)
  prs = get_prs(temp) if !prs
  prs = convert_units2Pa(prs)
  if ( !(temp.units === Units["K"]) )
    temp = convert_units(temp,"K")
  end
  theta = temp * (prs/P00)**(-Kappa)
  if theta.respond_to?(:long_name)  # VArray or GPhys
    theta.name = "theta"
    theta.long_name = "potential temperature"
  end
  theta
end

theta_e(temp,q,prs=nil) click to toggle source

Derive equivalent potential temperature

ARGUMENTS

temp [GPhys] : temperature (ok whether degC or K)
q [GPhys] : specific humidity
prs [GPhys or VArray] : the pressure values. If nil, searched from coordinates (for data on the pressure coordinate)

RETURN VALUE

theta_e: equivalent potential temperature

# File ../../lib/numru/ganalysis/met.rb, line 744
def theta_e(temp,q,prs=nil)
  tempK = temp.convert_units("K")
  theta = temp2theta(tempK, prs)
  theta_e = theta * Misc::EMath.exp( lat(tempK)*q/(Cp*tempK) )
  theta_e.name = "theta_e"
  theta_e.long_name = "equivalent potential temperature"
  theta_e
end

theta_es(temp,prs=nil) click to toggle source

Derive the saturation equivalent potential temperature

ARGUMENTS

temp [GPhys] : temperature (ok whether degC or K)
prs [GPhys or VArray] : the pressure values. If nil, searched from coordinates (for data on the pressure coordinate)

RETURN VALUE

theta_es: saturation equivalent potential temperature

# File ../../lib/numru/ganalysis/met.rb, line 763
def theta_es(temp,prs=nil)
  tempK = temp.convert_units("K")
  theta = temp2theta(tempK, prs)
  q = e_sat(temp).e2q(prs)
  theta_e = theta * Misc::EMath.exp( lat(tempK)*q/(Cp*tempK) )
  theta_e.name = "theta_es"
  theta_e.long_name = "theta_e sat"
  theta_e
end

z2geostrophic_wind(z, f=nil) click to toggle source

Derive geostrophic wind from geopotential hight (spherical but fixed f)

ARGUMENTS

z [GPhys] : geopotential height on the pressure (or log-pressure) coordinate
f [nil or Numeric of UNumeric] : the constant f value (Coriolis parameter). If nil, the value at the north pole is assumed. If Numeric, units are assumed to be “s-1”.

# File ../../lib/numru/ganalysis/met.rb, line 204
def z2geostrophic_wind(z, f=nil)
  if f.nil?
    f = 2 * GAnalysis::Planet.omega
  elsif f.is_a?(Numeric)
    f = UNumeric[f,"s-1"]
  end
  z = z.convert_units("m")
  gx, gy = GAnalysis::Planet.grad_s( z * (g/f) )
  u = -gy
  v = gx
  u.name = "u"
  u.long_name = "geostrophic U"
  u.put_att("assumed_f",f.to_s)
  v.name = "v"
  v.long_name = "geostrophic V"
  v.put_att("assumed_f",f.to_s)
  [u, v]
end

Home Classes Methods

In Files

Methods

Class/Module Index

NumRu::GAnalysis::Met

Constants

Public Class Methods

Public Instance Methods

gradH v; [thetax, thetay] = gradH theta;

cos_phi * gradH (u/cos_phi)