Module: NumRu::GAnalysis::Met

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

Overview

Meteorological analysis

USEFUL METHODS

  • temp2theta

  • pv_on_theta

  • interpolate_onto_theta

Constant Summary collapse

R =

< Themodynamic constants for the Earth’s atmosphere >

UNumeric[287.04,"J.kg-1.K-1"]
Rv =

Gas constant of water vapor

UNumeric[461.50,"J.kg-1.K-1"]
Cp =

heat capacity at constant pressure for dry air

UNumeric[1004.6,"J.kg-1.K-1"]
Cpv =

heat capacity at constant pressure for water vapor

UNumeric[1870.0,"J.kg-1.K-1"]
Kappa =

2.0/7.0

(R/Cp).to_f
T0 =

K - Celcius

UNumeric[273.15,"K"]
Lat0 =

Latent heat of vaporization at 0 degC

UNumeric[2.500780e6,"J.kg-1"]
P00 = 
UNumeric[1e5,"Pa"]
@@g =

< Earth’s parameters >

UNumeric[9.8,"m.s-2"]

Class Method Summary collapse

Class Method Details

.convert_units(obj, un) ⇒ Object

Unit conversion good to both GPhys and UNumeric



375
376
377
378
379
380
381
 
# File 'lib/numru/ganalysis/met.rb', line 375

def convert_units(obj, un)
  begin
    obj.convert_units(un)   # for GPhys etc
  rescue
    obj.convert2(un)        # for UNumeric etc
  end
end

.convert_units2Pa(prs) ⇒ Object

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

ARGUMENT

  • prs [GPhys and UNumeric]



354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
 
# File 'lib/numru/ganalysis/met.rb', line 354

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(f, x, dim) ⇒ Object

numerical derivative good to both GPhys and UNumeric



385
386
387
388
389
390
391
392
393
 
# File 'lib/numru/ganalysis/met.rb', line 385

def df_dx(f, x, dim)
  if GPhys === f
    mdl = NumRu::GPhys::Derivative
    mdl.threepoint_O2nd_deriv(f, dim, mdl::LINEAR_EXT, x)
  else
    mdl = NumRu::Derivative
    mdl.threepoint_O2nd_deriv(f, x, dim, mdl::LINEAR_EXT)
  end
end

.df_dx_vialogscale(f, x, dim) ⇒ Object 



396
397
398
399
400
401
402
403
404
405
406
 
# File 'lib/numru/ganalysis/met.rb', line 396

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

.find_prs_d(gphys, error = nil) ⇒ Object

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)



308
309
310
311
312
313
314
315
316
317
318
319
320
321
 
# File 'lib/numru/ganalysis/met.rb', line 308

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) ⇒ Object

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]



224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
 
# File 'lib/numru/ganalysis/met.rb', line 224

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 ⇒ Object

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



50
51
52
 
# File 'lib/numru/ganalysis/met.rb', line 50

def g
  @@g.dup
end

.get_prs(gphys) ⇒ Object

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.



332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
 
# File 'lib/numru/ganalysis/met.rb', line 332

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

.interpolate_onto_theta(gphys, theta, theta_levs) ⇒ Object

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



115
116
117
118
119
120
121
122
 
# File 'lib/numru/ganalysis/met.rb', line 115

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

.pv_on_theta(u, v, theta, theta_levs) ⇒ Object

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 temperature [GPhys] on a theta coordinate, where levels are set to theta_levs



138
139
140
141
142
143
144
145
146
147
148
 
# File 'lib/numru/ganalysis/met.rb', line 138

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

.set_g(g) ⇒ Object

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



44
45
46
 
# File 'lib/numru/ganalysis/met.rb', line 44

def set_g(g)
  @@g = g   # e.g., un[9.8,"m.s-2"]
end

.sigma_inv(theta, dim = nil, prs = nil) ⇒ Object

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]



94
95
96
97
98
99
100
101
102
103
104
105
106
 
# File 'lib/numru/ganalysis/met.rb', line 94

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) ⇒ Object

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



66
67
68
69
70
71
72
73
74
75
76
77
78
 
# File 'lib/numru/ganalysis/met.rb', line 66

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

.z2geostrophic_wind(z, f = nil) ⇒ Object

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”.



160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
 
# File 'lib/numru/ganalysis/met.rb', line 160

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