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"]
@@consider_ice =

moist thermodynamics ###

true
@@ice_thres =

Celcius

T0

Class Method Summary collapse

Class Method Details

.consider_ice ⇒ Object

whether or not ice is considered in the water phase change



475
476
477
 
# File 'lib/numru/ganalysis/met.rb', line 475

def consider_ice
  @@consider_ice
end

.consider_ice=(t_or_f) ⇒ Object

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



480
481
482
 
# File 'lib/numru/ganalysis/met.rb', line 480

def consider_ice=(t_or_f)
  @@consider_ice=t_or_f
end

.convert_units(obj, un) ⇒ Object

Unit conversion good to both GPhys and UNumeric



436
437
438
439
440
441
442
 
# File 'lib/numru/ganalysis/met.rb', line 436

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]



415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
 
# File 'lib/numru/ganalysis/met.rb', line 415

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

.del_ngp(*gps) ⇒ Object

-g del/del p



185
186
187
188
189
190
191
 
# File 'lib/numru/ganalysis/met.rb', line 185

def del_ngp(*gps)  # gps: array of GPhys objects having the same grid
  d = find_prs_d(gps[0])
  prs = gps[0].coord(d).convert_units(Units["Pa"])
  dngps = gps.collect{|gp|
    gp.threepoint_O2nd_deriv(d,Derivative::LINEAR_EXT,prs) * (-@@g)
  }
end

.df_dx(f, x, dim) ⇒ Object

numerical derivative good to both GPhys and UNumeric



446
447
448
449
450
451
452
453
454
 
# File 'lib/numru/ganalysis/met.rb', line 446

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 



457
458
459
460
461
462
463
464
465
466
467
 
# File 'lib/numru/ganalysis/met.rb', line 457

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

water vapor pressure -> specific humidity

ARGUMENTS

  • e: water vapor pressure

  • prs: pressure

RETURN VALUE

  • q: specific humidity



584
585
586
587
588
589
590
591
592
 
# File 'lib/numru/ganalysis/met.rb', line 584

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

water vapor pressure -> mixing ratio

ARGUMENTS

  • e: water vapor pressure

  • prs: pressure

RETURN VALUE

  • r: mixing ratio



548
549
550
551
552
553
554
555
556
 
# File 'lib/numru/ganalysis/met.rb', line 548

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

Calculates saturation water vapor pressure using enhanced

ARGUMENTS

  • temp: temperature

RETURN VALUE

  • es: saturation water vapor pressure



712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
 
# File 'lib/numru/ganalysis/met.rb', line 712

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

Bolton formula for saturation water vapor pressure against water

ARGUMENTS

  • temp: temperature

RETURN VALUE

  • es: saturation water vapor pressure [Pa]



618
619
620
621
622
623
624
625
 
# File 'lib/numru/ganalysis/met.rb', line 618

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) ⇒ Object Also known as: e_sat_ice

Saturation water vapor pressure against ice.

Emanuel (1994) eq.(4.4.15)

ARGUMENTS

  • temp: temperature

RETURN VALUE

  • es: saturation water vapor pressure [Pa]



663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
 
# File 'lib/numru/ganalysis/met.rb', line 663

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

.e_sat_emanuel_water(temp, mask = nil) ⇒ Object

saturation water vapor pressure against ice.

Emanuel (1994) eq.(4.4.15)

ARGUMENTS

  • temp: temperature

RETURN VALUE

  • es: saturation water vapor pressure [Pa]



636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
 
# File 'lib/numru/ganalysis/met.rb', line 636

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

.e_sat_water ⇒ Object

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



369
370
371
372
373
374
375
376
377
378
379
380
381
382
 
# File 'lib/numru/ganalysis/met.rb', line 369

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]



285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
 
# File 'lib/numru/ganalysis/met.rb', line 285

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



60
61
62
 
# File 'lib/numru/ganalysis/met.rb', line 60

def g
  @@g.dup
end

.get_f(gp) ⇒ Object

Get the Coriolis parameter



231
232
233
234
235
236
237
238
 
# File 'lib/numru/ganalysis/met.rb', line 231

def get_f(gp)
  lam, phi, = GAnalysis::Planet.get_lambda_phi(gp)
  sin_phi = phi.sin
  f = 2 * GAnalysis::Planet.omega * sin_phi
  f.name = "f"
  f.long_name = "Coriolis parameter"
  f
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.



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

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

the threshold temperature for liquid/ice-phase treatment



485
486
487
 
# File 'lib/numru/ganalysis/met.rb', line 485

def ice_thres
  @@ice_thres
end

.ice_thres=(temp) ⇒ Object

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



490
491
492
 
# File 'lib/numru/ganalysis/met.rb', line 490

def ice_thres=(temp)
  @@ice_thres=temp
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



125
126
127
128
129
130
131
132
 
# 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) ⇒ Object

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

good for -100<T<50

ARGUMENTS

  • temp: temperature

RETURN VALUE

  • lat: latent heat



603
604
605
606
607
608
609
 
# File 'lib/numru/ganalysis/met.rb', line 603

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

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



173
174
175
176
177
178
179
180
181
182
 
# 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) ⇒ 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 voticity [GPhys] on a theta coordinate, where levels are set to theta_levs



148
149
150
151
152
153
154
155
156
157
158
 
# 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) ⇒ Object

specific humidity -> water vapor pressure

ARGUMENTS

RETURN VALUE

  • e: water vapor pressure



566
567
568
569
570
571
572
573
574
 
# File 'lib/numru/ganalysis/met.rb', line 566

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

specific humidity -> mixing ratio

ARGUMENTS

  • q: specific humidty

RETURN VALUE

  • r: mixing ratio



501
502
503
504
505
506
 
# File 'lib/numru/ganalysis/met.rb', line 501

def q2r(q)
  r = q/(1.0-q)
  r.name = "r"
  r.long_name = "mixing ratio"
  r
end

.r2e(r, prs = nil) ⇒ Object

water vapor mixing ratio -> water vapor pressure

ARGUMENTS

  • r: water vapor mixing ratio

  • prs: pressure

RETURN VALUE

  • e: water vapor pressure



530
531
532
533
534
535
536
537
538
 
# File 'lib/numru/ganalysis/met.rb', line 530

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

mixing ratio -> specific humidity

ARGUMENTS

  • r: mixing ratio

RETURN VALUE

  • q: specific humidty



515
516
517
518
519
520
 
# File 'lib/numru/ganalysis/met.rb', line 515

def r2q(r)
  q = r/(1.0+r)
  q.name = "q"
  q.long_name = "specific humidity"
  q
end

.rh2e(rh, temp) ⇒ Object

relative humidity -> water vapor pressure

ARGUMENTS

  • rh: relative humidity

  • temp: temperature

RETURN VALUE

  • e: water vapor pressure



741
742
743
744
745
746
747
748
 
# File 'lib/numru/ganalysis/met.rb', line 741

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

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

ARGUMENTS

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

RETURN VALUE

  • nil



688
689
690
691
692
693
694
695
696
697
 
# File 'lib/numru/ganalysis/met.rb', line 688

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



54
55
56
 
# 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) ⇒ 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]



104
105
106
107
108
109
110
111
112
113
114
115
116
 
# 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) ⇒ 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



76
77
78
79
80
81
82
83
84
85
86
87
88
 
# 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) ⇒ Object

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



761
762
763
764
765
766
767
768
 
# File 'lib/numru/ganalysis/met.rb', line 761

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

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



780
781
782
783
784
785
786
787
788
 
# File 'lib/numru/ganalysis/met.rb', line 780

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

Derive geostrophic wind from geopotential hight (spherical)

ARGUMENTS

  • z [GPhys] : geopotential height on the pressure (or log-pressure) coordinate

  • f [nil or GPhys or Numeric or UNumeric] : Coriolis parameter If nil, local Corioris parameter is used. If Numeric or UNumeric, f is a constant as a matter of course (if Numeric, units are assumed to be “s-1”).



204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
 
# File 'lib/numru/ganalysis/met.rb', line 204

def z2geostrophic_wind(z, f=nil)
  if f.nil?
    f = get_f(z)
  elsif f.is_a?(Numeric)
    f = UNumeric[f,"s-1"]
  end
  #z = z.convert_units("m")
  gx, gy = GAnalysis::Planet.grad_s( z * g )
  u = -gy/f
  v = gx/f
  if f.respond_to? :val   # supposedly a GPhys
    if0 = f.val.eq(0).where
    if if0.length > 0
      u[true,if0,false] = 0   # replace Inf with 0
      v[true,if0,false] = 0   # replace Inf with 0
    end
  end
  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