Module: NumRu::GAnalysis::QG_common

Included in:
QG, QG_sphere, QG_sphere_div
Defined in:
lib/numru/ganalysis/qg.rb

Overview

QG_common: collection of common methods for QG, QG_sphere, and QG_sphere_div.

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.cut_bottom(z) ⇒ Object 



154
155
156
157
158
159
160
161
162
 
# File 'lib/numru/ganalysis/qg.rb', line 154

def cut_bottom(z)
  pdim = Met.find_prs_d(z)
  plev = z.coord(pdim).val
  if plev[0] - plev[1] > 0
    z[ *([true]*pdim + [1..-1,false]) ]
  else
    z[ *([true]*pdim + [0..-2,false]) ]
  end
end

.extend_bottom(z, val_extended = nil) ⇒ Object

Extend the bottom pressure level by the lowest thickness (a hypothetical “Under-ground” level is created) If value of the extended bottom level is set to val_extended (Numeric or NArray etc), if it is specified (non nil). If nil, the value at the original bottom level is simply copied.



127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
 
# File 'lib/numru/ganalysis/qg.rb', line 127

def extend_bottom(z, val_extended=nil)
  pdim = Met.find_prs_d(z)
  plev = z.coord(pdim).val
  raise("Only one pressure level is found; 2 or more needed") if (plev==1)
  bottom_first = ( plev[0] - plev[1] > 0 )
  np = z.shape[pdim]
  idx = (0...np).collect{|i| i}
  if bottom_first       # The first level is the bottom one
    idx.unshift(0)      # idx => [0,0,1,2,...,np-1]
    ihb = 0             # index of the extended bottom level
    dp = plev[0] - plev[1]
    phb = plev[0] + dp  # pressure of the extended bottom level
  else                # The last level is the bottom one
    idx.push(np-1)     # idx => [0,1,2,...,np-1,np-1]
    ihb = np            # index of the extended bottom level
    dp = plev[-1] - plev[-2]
    phb = plev[-1] + dp  # pressure of the extended bottom level
  end
  ze = z[ *([true]*pdim + [idx,false]) ].copy  # add one level below
  ze.coord(pdim)[ihb] = phb
  if val_extended
    ze[ *([true]*pdim + [ihb,false]) ] = val_extended
  end
  ze
end

.gp2gpref(gp) ⇒ Object

geopotential (multi-D) -> reference geopotential profile (1D)



24
25
26
27
28
29
30
31
32
33
34
35
36
 
# File 'lib/numru/ganalysis/qg.rb', line 24

def gp2gpref(gp)
  gpref = Planet::ave_s(gp)  # horizontal ave (spherical)
  if gpref.rank >= 2
    # likely a time sequence. need to reduce more.
    pdim = Met.find_prs_d(gpref)
    idxs = (0...gpref.rank).collect{|i| i}
    idxs.delete(pdim)
    gpref = gpref.mean(*idxs)
  end
  gpref.name = "gpref"
  gpref.long_name = "Reference geopotential"
  gpref
end

.gpd2qzz(gp, b) ⇒ Object

[ (p/b) gp_z ]_z /p



65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
 
# File 'lib/numru/ganalysis/qg.rb', line 65

def gpd2qzz(gp, b)
  bunits = Units["s-2"]   # this is assumed!
  pdim = Met.find_prs_d(gp)
  p = gp.coord(pdim)
  z = LogP.p2z(p)
  zunits = z.units
  g = Derivative::b_expand_linear_ext( gp.val, pdim )
  z = Derivative::b_expand_linear_ext( z.val, 0 )
  p = Derivative::b_expand_linear_ext( p.val, 0 )
  b = b.val
  b = b.to_na if b.respond_to?(:to_na)  # likely a NArrayMiss
  b = Derivative::b_expand_linear_ext( b, 0 )
  pb = p/b

  pbm = (pb[0..-2] + pb[1..-1]) / 2.0  # pb_{i+1/2} (for i=0..-2)
  pbm01 = pbm[0..-2]                   # pb_{i-1/2} (for i=1..-2)
  pbm12 = pbm[1..-1]                   # pb_{i+1/2} (for i=1..-2)
  dz20 = z[2..-1] - z[0..-3]           # z_{i+1} - z_{i-1} (for i=1..-2)
  dz21 = z[2..-1] - z[1..-2]           # z_{i+1} - z_{i}   (for i=1..-2)
  dz10 = z[1..-2] - z[0..-3]           # z_{i} - x_{i-1}   (for i=1..-2)
  pc = p[1..-2]                        # p_{i}   (for i=1..-2)

  a2 = 2*pbm12/(dz21*dz20)/pc
  a0 = 2*pbm01/(dz10*dz20)/pc
  a1 = -a2 - a0

  to_rankD = [1]*pdim + [true] + [1]*(gp.rank-pdim-1)
  a2 = a2.reshape(*to_rankD)
  a1 = a1.reshape(*to_rankD)
  a0 = a0.reshape(*to_rankD)

  vqzz =   g[ *([true]*pdim+[2..-1,false]) ] * a2 \
         + g[ *([true]*pdim+[1..-2,false]) ] * a1 \
         + g[ *([true]*pdim+[0..-3,false]) ] * a0

  qzz = gp.copy
  qzz.data.replace_val(vqzz)
  qzz.name = "qzz"
  qzz.long_name = "z-deriv term in QG PV"
  qzz.units = qzz.units / zunits**2 / bunits
  qzz
end

.gph2gpd_gpref(gph) ⇒ Object

geopotential height to geopotential deviation from the global&time mean



40
41
42
43
44
45
46
47
 
# File 'lib/numru/ganalysis/qg.rb', line 40

def gph2gpd_gpref(gph)
  gp = gph * Met::g
  gpref = gp2gpref(gp)
  gpd = gp - gpref
  gpd.name = "gpd"
  gpd.long_name = "Geopotential deviation"
  [gpd, gpref]
end

.gph2gpref(gph) ⇒ Object

geopotential height (multi-D) -> reference geopotential profile (1D)



18
19
20
 
# File 'lib/numru/ganalysis/qg.rb', line 18

def gph2gpref(gph)
  gp2gpref(gph) * Met::g
end

.gpref2n2(gpref) ⇒ Object

reference geopotential -> buoyancy frequency squared



51
52
53
54
55
56
57
58
59
60
61
 
# File 'lib/numru/ganalysis/qg.rb', line 51

def gpref2n2(gpref)
  gp_z = LogP.pcdata_dz( gpref )
  gp_zz = LogP.pcdata_dz2( gpref )
  gp_zz[0] = gp_zz[1]   # At boundary, it's safer to extend lapse rate
  gp_zz[-1] = gp_zz[-2] # At boundary, it's safer to extend lapse rate
  n2 = gp_zz + gp_z * (Met::Kappa / LogP.h)
  n2.name = "N2"
  n2.long_name = "N**2 (log-p)"
  #p "@@@@@ N2 @@@@",n2.coord(0).val.to_a, n2.val.sqrt.to_a
  n2
end

Instance Method Details

#div_waf(fx, fy, fz, bottom_treatment = true) ⇒ Object

div of WAF

(p cos_phi)^-1 div(p waf) = (cos_phi)^-1 ( div_h(fx,fy) + p^-1 d_z (p fz) )

  • fx, fy, fz (GPhys) : the x, y and z components of waf

  • bottom_treatment (true (==default) or false) : If true, the lowest level vertical divergence is computed by assuming that fz is zero at the extended “underground” level. The thickness assumed (=p-p) is consistent with the ((<extend_bottom>)) method.



185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
 
# File 'lib/numru/ganalysis/qg.rb', line 185

def div_waf(fx, fy, fz, bottom_treatment=true)
  cosphi = cos_phi(fx)
  p = Met.get_prs(fx)

  fz_z = LogP.pcdata_dz( fz*p ) / p

  #>>>>>> the lowest layer treatment consistent with qb, in which
  # geopotential (or stream function) is extended by extend_bottom. 
  # Assumption: the first level is the lowest (bottom) one
  if bottom_treatment
    # using the relation p^{-1} d/dz = -H^{-1} d/dp 
    # and assuming fz=0 below the bottom (the "underground" level),
    # p^{-1} d/dz (p fz) = -H^{-1} d/dp (p fz),
    # which is H^{-1} p fz / delta_p, at the lowest level with a
    # "thickness" of delta_p. 
    w = p[0..1].val
    dp = w[1] - w[0]
    p0 = w[0]
    pdim = Met.find_prs_d(fz)
    sel0 = [true]*pdim + [0,false]  # to specify the first level
    fz_z[*sel0] = fz[*sel0]*p0 / (LogP.h*dp)
  end
  #<<<<<<

  divh = ( div_h(fx, fy) + fz_z ) / cosphi
         # ^ div_h is defined in QG, QG_sphere,..., but not in QG_common
  divh.name = "divwaf"
  divh.long_name = "div of waf (#{fx.name},..)"
  divh
end

#waf_tn2001(gph, gphbs, nx: nil, ny: nil, ubs: nil, vbs: nil, zonal: false, merid: false, n2: nil, derive_M0: false) ⇒ Object

Wave activity flux by Takaya and Nakamura (2001) (TN2001)

This method works by mix-in in QG (Cartesian) or QG_sphere (spherical). (so the call it either as QG::waf_tn2001(gph, gphbs, ..) or QG_sphere::waf_tn2001(gph, gphbs, ..) ).

DEFINITIONS

  • if called through QG: a modified version of Eq (32) of TN2001

  • if called through QG_sphere: a modified version of Eq (38) of TN2001

The modification of Eq (32) and (38) is as follows:

  • The multiplication by p is omitted. This is because

    it is more convenient to do it afterwords if needed
(for example, to integrate with p would require not to have
done this multiplication).

  • The C_u M term is not included as in convention. (but the multiplication by cos(phi) is done in QG_sphere). Note that this term vanishes for stationary waves, because C_u=0.

  • The normal vector (n1, n2) in TN2001 can be set flexibly. This vector is called (nx, ny) in this method, since n2 is used for the buoyancy frequency. In eqs. (32) and (38), (nx, ny) is prescribed as (nx, ny) = (U,V)/|U,V| where U is ubs and V is vbs. However, this relation is not necessarily valid, since (nx, ny) is originally attributed to the QGPV gradient of the basic flow. Therefore, in this method, (nx, ny) can be set arbitrary. Also, one can set (nx, ny) = (1, 0) or (0, 1) by setting the optional argument zonal or merid to true.

ARGUMENTS

  • gph: geopotential height

  • gphbs: geopotential height (basic state). Basically, it needs to have the same spacial grid as gph, but it can be a meridional section (such as zonal mean), if the normal vector is specified by (ubs, vbs) or (nx, ny) or zonal or merid.

  • nx, ny: normalized vector to designate the direction of wave propagation. By default, set to the direction of (ubs, vbs), or set by zonal/merid (see below).

  • ubs, vbs: Basic state zonal and meridional winds to compute the normal vector (nx,ny) as (ubs,vbs)/sqrt(ubs**2 + vbs**2). Neglected if nx and ny are specified. By default ubs and vbs are derived geostrophically from gphbs

  • zonal: derive the flux in the zonal direction assuming (nx,ny)=(1,0)

  • merid: derive the flux in the meridional dir assuming (nx,ny)=(0,1)

  • n2: Brunt-Vaisala (buoyancy) frequency squared. By default, derived from gphbs

  • derive_M0: derive M (the wave-activity pseudomomentum of TN20001) when the wave is stationary (cp=0). You need to give ubs and vbs (or have them derived from gphbs by not using nx, ny, zonal or merid) to make this option available.

RETURN VALUE

  • wx, wy, wz, ret_hash

    where [wx, wy, wz] are the WAF vector

    and ret_hash provides extra derived quantities such as the perturbation stream function (:pp), pseudomomentum M0 (:M0), etc.

REFERENCE

  • Takaya and Namaura (2001): A formulation of a phase-independent wave-activity flux for stationary and migratory quasigeostrophic eddies on a zonally varying basic flow. J. Atmos. Sci, 58, 608-627.



276
277
278
279
280
281
282
283
284
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
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
 
# File 'lib/numru/ganalysis/qg.rb', line 276

def waf_tn2001(gph, gphbs, nx: nil, ny: nil, ubs: nil, vbs: nil,
               zonal: false, merid: false, n2: nil,
               derive_M0: false)

  #< prepration >
  if self==QG
    f = f0
    wgt = nil
  elsif self==QG_sphere
    f = f_mask0(gph)
    wgt = cos_phi(gph)
  else
    raise("unsupported module to mix-in #{self}")
  end
  ret_hash = Hash.new

  #< perturbation stream function >
  gf = Met::g / f
  pp = (gph - gphbs)*gf   # psi' (perturbation stream function)
  if /gpm/i =~ (us=pp.units.to_s)
    pp.units = us.sub(/gpm/i,'m')
  end
  ret_hash[:pp] = pp

  #< derive nx, ny if not given >
  ul = psibs = nil
  if zonal
    nx = 1.0
    ny = 0.0
  elsif merid
    nx = 0.0
    ny = 1.0
  elsif nx.nil? && ny.nil?
    if !ubs || !vbs
      psibs, gpref = gph2psi_gpref(gphbs)
      if /gpm/i =~ (us=psibs.units.to_s)
        psibs.units = us.sub(/gpm/i,'m')
      end
      ubs, vbs = psi2ug_vg(psibs)
      ret_hash[:ubs] = ubs
      ret_hash[:vbs] = vbs
      ret_hash[:gpref] = gpref
    end
    ul = (ubs**2+vbs**2).sqrt
    nx = ubs/ul
    ny = vbs/ul
  end

  #< derive aspect ratio >
  unless n2
    n2 = LogP::gph2n2(gphbs)
    ret_hash[:n2] = n2
  end
  eps = f**2/n2

  #< cross terms --> WAF >
  px,py = grad_h(pp)
  pxx,pxy = grad_h(px)
  pyy = grad_hy(py) unless ny==0.0
  pz = LogP::pcdata_dz(pp)
  pxz = LogP::pcdata_dz(px) unless nx==0.0
  pyz = LogP::pcdata_dz(py) unless ny==0.0
  ret_hash[:px] = px
  ret_hash[:py] = py

  wx = nil
  unless nx==0.0
    wx = px*px - pp*pxx
    wy = px*py - pp*pxy
    wz = ( px*pz - pp*pxz ) * eps
    nx2 = nx/2.0
    wx = wx * nx2
    wy = wy * nx2
    wz = wz * nx2
  end

  unless ny==0.0
    wx_ = px*py - pp*pxy
    wy_ = py*py - pp*pyy
    wz_ = (py*pz - pp*pyz) * eps
    ny2 = ny/2.0
    wx_ = wx_ * ny2
    wy_ = wy_ * ny2
    wz_ = wz_ * ny2
    if wx
      wx += wx_
      wy += wy_
      wz += wz_
    else
      wx = wx_
      wy = wy_
      wz = wz_
    end
  end

  if wgt
    wx = wx * wgt
    wy = wy * wgt
    wz = wz * wgt
  end
  wx.name = "Wx"
  wy.name = "Wy"
  wz.name = "Wz"
  wx.long_name = "WAFx"
  wy.long_name = "WAFy"
  wz.long_name = "WAFz"

  #< derive M0 if required >
  if derive_M0
    raise("ubs and vbs are needed to derive M0") if !ubs || !vbs
    ul = (ubs**2+vbs**2).sqrt unless ul
    unless psibs
      psibs, gpref = gph2psi_gpref(gphbs)
      if /gpm/i =~ (us=psibs.units.to_s)
        psibs.units = us.sub(/gpm/i,'m')
      end
      ret_hash[:gpref] = gpref
    end
    n2p = Planet::ave_s(n2)
    if n2p.rank >= 2
      pdim = Met.find_prs_d(n2p)
      ds = Array.new
      (0...n2p.rank).each do |d| 
        ds.push(d) unless d == pdim
      end
      n2p = n2p.mean(*ds)
    end
    qp, = psi2q(pp,n2p)
    qbs, = psi2q(psibs,n2p)
    qbx, qby = grad_h(qbs)
    e2 = px**2 + py**2 + pz**2 * eps  # 2*e
    qbgrad = (qbx**2+qby**2).sqrt.running_mean(0,3).running_mean(1,3)
    a = qp**2 / qbgrad / 4.0   # A of TN2001
    e0 = e2 / ul / 4.0         # E of TN2001 when cp=0
    m0 = a + e0
    a.name = "A"
    a.long_name = "A"
    e0.name = "E0"
    e0.long_name = "stationary E"
    m0.name = "M0"
    m0.long_name = "stationary TN2001 pseudomomentum"
    ret_hash[:qp] = qp
    ret_hash[:qbs] = qbs
    ret_hash[:A] = a
    ret_hash[:E0] = e0
    ret_hash[:M0] = m0
  end

  #< return >
  [wx, wy, wz, ret_hash]
end