Module: NumRu::GAnalysis::Planet

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

Overview

Library for spherical planets (thin spherical shell; default: Earth)

ASSUMPTIONS

  • longitude is assumed to increase in the eastward direction.

  • latitude is assumed to increase in the northward direction, and it is zero at the equator.

Constant Summary collapse

Earth =

< Pre-defined planets >

"Earth"
@@lonbc =

< Class variables regarding lon & lat >

GPhys::Derivative::CYCLIC_OR_LINEAR
@@lam =

@@latbc = GPhys::Derivative::LINEAR_EXT # this should be always fine

nil
@@phi =

lambda (lon in radian) obtaiend lately (see get_lambda_phi)

nil

Class Method Summary collapse

Class Method Details

.absvor_s(u, v) ⇒ Object 



181
182
183
184
185
186
 
# File 'lib/numru/ganalysis/planet.rb', line 181

def absvor_s(u,v)
  avor = vor_s(u,v) + @@phi.sin * (2*@@Ome)
  avor.long_name = "Absolute Vorticity"
  avor.name = "avor"
  avor
end

.ave_s(s) ⇒ Object

horizontal averaging considering the spherical geometry



66
67
68
69
70
71
72
73
74
 
# File 'lib/numru/ganalysis/planet.rb', line 66

def ave_s(s)
  lam, phi, lond, latd  = get_lambda_phi(s)
  xmean = s.mean(lond)
  cos_phi = phi.cos

  lond,latd = find_lon_lat_dims(xmean)  # find latd again
  wgt = cos_phi / cos_phi.sum
  (xmean * wgt).sum(latd)
end

.div_s(u, v) ⇒ Object 



87
88
89
90
91
92
93
94
95
96
 
# File 'lib/numru/ganalysis/planet.rb', line 87

def div_s(u,v)
  lam, phi, lond, latd  = get_lambda_phi(u)
  cos_phi = phi.cos
  du_dlam = u.cderiv(lond,@@lonbc,lam)
  dvc_dphi  = (v*cos_phi).cderiv(latd,latbc(phi),phi)
  rot = (du_dlam + dvc_dphi) / (@@R*cos_phi)
  rot.long_name = "div(#{u.name},#{v.name})"
  rot.name = "div"
  rot
end

.div_s_box(u, v, only_u: false, only_v: false, xs0: false, xe0: false, ys0: false, ye0: false) ⇒ Object

spherical divergence for entire rectagular regional data



99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
 
# File 'lib/numru/ganalysis/planet.rb', line 99

def div_s_box(u,v, only_u: false, only_v: false,
              xs0: false, xe0: false, ys0: false, ye0: false)
  lam, phi, xd, yd  = get_lambda_phi(u)
  lon = u.coord(xd)
  lat = u.coord(yd)
  cos_phi = phi.cos
  nx = lam.length
  ny = phi.length
  if nx>2
    ixends = {0..-1=>(nx-1)}
  else
    ixends = true
  end
  if ny>2
    iyends = {0..-1=>(ny-1)}
  else
    iyends = true
  end
  ub = u[*([true]*xd + [ixends, false])].copy  # ends along lon
  vb = v[*([true]*yd + [iyends, false])].copy  # ends along lat
  ub.axis(xd).set_pos(lam[ixends].data)
  ub.axis(yd).set_pos(phi.data)
  vb.axis(xd).set_pos(lam.data)
  vb.axis(yd).set_pos(phi[iyends].data)
  ubi = ub.integrate(yd)   # trapezoidal integration by default
  vbi = vb.integrate(xd) * cos_phi[iyends]
  dlam = lam[-1].val - lam[0].val
  sin_phib = phi[iyends].sin.val
  dsin_phi = sin_phib[-1] - sin_phib[0]
  if xd<yd
    xdi = xd
    ydi = yd-1
  else
    xdi = xd-1
    ydi = yd
  end
  scalar=false
  ubi.axis(xdi).set_pos(lon[ixends])  # to lon before the cut is recoreded
  vbi.axis(ydi).set_pos(lat[iyends])  # to lat before the cut is recoreded
  if ubi.rank>1
    xe = ubi[*([true]*xdi+[-1,false])]
    xs = ubi[*([true]*xdi+[0,false])]
    ye = vbi[*([true]*ydi+[-1,false])]
    ys = vbi[*([true]*ydi+[0,false])]
  else
    # to defer the result becomes a scalar
    scalar=true
    xe = ubi[*([true]*xdi+[-1..-1,false])]
    xs = ubi[*([true]*xdi+[0..0,false])]
    ye = vbi[*([true]*ydi+[-1..-1,false])]
    ys = vbi[*([true]*ydi+[0..0,false])]
  end
  xe *= 0.0 if xe0
  xs *= 0.0 if xs0
  ye *= 0.0 if ye0
  ys *= 0.0 if ys0
  dub = xe - xs
  dvb = ye - ys
  if only_u
    div = dub/(@@R*dlam*dsin_phi.abs) 
  elsif only_v
    div = dvb/(@@R*dlam.abs*dsin_phi)
  else
    div = dub/(@@R*dlam*dsin_phi.abs) +
          dvb/(@@R*dlam.abs*dsin_phi)
        # (dlam*dsin_phi).abs is the area for the unit radius
  end
  if scalar
    div = div[0]
  end
  div.long_name = "div(#{u.name},#{v.name})"
  div.name = "div"
  div
end

.find_lat_dim(gp, err_raise = false) ⇒ Object 



371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
 
# File 'lib/numru/ganalysis/planet.rb', line 371

def find_lat_dim(gp, err_raise=false)
  latd = nil
  (0...gp.rank).each do |dim|
    crd = gp.coord(dim)
    if /^degrees?_north$/i =~ crd.get_att("units")
      latd = dim
      break
    elsif ( ( /latitude/i =~ crd.long_name || 
              /^lat/i =~ crd.name ||
              (nm=crd.get_att("standard_name") && /latitude/i =~ nm ) &&
              (crd.units =~ Units["degrees_north"]) ) )
      latd = dim
      break
    end
  end
  if err_raise
    raise("Latitude dimension was not found") if !latd
  end
  latd
end

.find_lon_dim(gp, err_raise = false) ⇒ Object 



350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
 
# File 'lib/numru/ganalysis/planet.rb', line 350

def find_lon_dim(gp, err_raise=false)
  lond = nil
  (0...gp.rank).each do |dim|
    crd = gp.coord(dim)
    if /^degrees?_east$/i =~ crd.get_att("units")
      lond = dim
      break
    elsif ( ( /longitude/i =~ crd.long_name || 
              /^lon/i =~ crd.name ||
              (nm=crd.get_att("standard_name") && /longitude/i =~ nm ) &&
              (crd.units =~ Units["degrees_east"]) ) )
      lond = dim
      break
    end
  end
  if err_raise
    raise("Longitude dimension was not found") if !lond
  end
  lond
end

.find_lon_lat_dims(gp, err_raise = false) ⇒ Object

Find longitude and latitude coordinates.

ARGUMENTS

  • gp : GPhys to inspect

  • err_raise (OPTIONAL; default:false) : if true, exception is raised if longitude or latitude coordinate is not found.

SEARCH CRITERIA (1) Find coord having units “degrees_east” (lon) or

"degrees_north" (lat)

(2) Investigate coordinate name matches (to find a lonitude coord,

/longitude/i for long_name or standard_name, or /^lon/ for name)
and match units compatible with "degrees".

RETURN VALUE

  • lond,latd

    • lond: dimension number of longitude (0,1,..) if found / nil if not found

    • latd: dimension number of latitude (0,1,..) if found / nil if not found



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
 
# File 'lib/numru/ganalysis/planet.rb', line 313

def find_lon_lat_dims(gp, err_raise=false)
  lond = nil
  latd = nil
  (0...gp.rank).each do |dim|
    crd = gp.coord(dim)
    if /^degrees?_east$/i =~ crd.get_att("units")
      lond = dim
      break
    elsif ( ( /longitude/i =~ crd.long_name || 
              /^lon/i =~ crd.name ||
              (nm=crd.get_att("standard_name") && /longitude/i =~ nm ) &&
              (crd.units =~ Units["degrees_east"]) ) )
      lond = dim
      break
    end
  end
  (0...gp.rank).each do |dim|
    next if dim == lond
    crd = gp.coord(dim)
    if /^degrees?_north$/i =~ crd.get_att("units")
      latd = dim
      break
    elsif ( ( /latitude/i =~ crd.long_name || 
              /^lat/i =~ crd.name ||
              (nm=crd.get_att("standard_name") && /latitude/i =~ nm ) &&
              (crd.units =~ Units["degrees_north"]) ) )
      latd = dim
      break
    end
  end
  if err_raise
    raise("Longitude dimension was not found") if !lond
    raise("Latitude dimension was not found") if !latd
  end
  [lond,latd]
end

.get_lambda(gp, err_raise = true) ⇒ Object 



275
276
277
278
279
280
281
 
# File 'lib/numru/ganalysis/planet.rb', line 275

def get_lambda(gp, err_raise=true)
  lond = find_lon_dim(gp, err_raise)
  lam = lond && gp.axis(lond).to_gphys.convert_units("rad") # lon in rad
  lam.units = "" if lond                            # treat as non-dim
  @@lam = lam
  [lam, lond]
end

.get_lambda_phi(gp, err_raise = true) ⇒ Object

Find longitude and latitude coordinates and convert into radian.

RETURN VALUE

  • lam, phi, lond, latd

    • lam (GPhys): longitude in radian (lambda). (nil if not found && !err_raise)

    • phi (GPhys): latitude in radian (lambda). (nil if not found && !err_raise)

    • lond : Interger to show that longitude is the lon-th dim if found; (nil if not found && !err_raise)

    • latd : Interger to show that latitude is the lat-th dim iffound; (nil if not found && !err_raise)



264
265
266
267
268
269
270
271
272
273
 
# File 'lib/numru/ganalysis/planet.rb', line 264

def get_lambda_phi(gp, err_raise=true)
  lond, latd = find_lon_lat_dims(gp, err_raise)
  lam = lond && gp.axis(lond).to_gphys.convert_units("rad") # lon in rad
  lam.units = "" if lond                            # treat as non-dim
  phi = latd && gp.axis(latd).to_gphys.convert_units("rad") # lat in rad
  phi.units = "" if latd                            # treat as non-dim
  @@lam = lam
  @@phi = phi
  [lam, phi, lond, latd]
end

.get_phi(gp, err_raise = true) ⇒ Object 



283
284
285
286
287
288
289
 
# File 'lib/numru/ganalysis/planet.rb', line 283

def get_phi(gp, err_raise=true)
  latd = find_lat_dim(gp, err_raise)
  phi = latd && gp.axis(latd).to_gphys.convert_units("rad") # lat in rad
  phi.units = "" if latd                            # treat as non-dim
  @@phi = phi
  [phi, latd]
end

.grad_s(s) ⇒ Object 



188
189
190
191
192
193
194
195
196
197
198
 
# File 'lib/numru/ganalysis/planet.rb', line 188

def grad_s(s)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  xs = s.cderiv(lond,@@lonbc,lam) / (@@R*cos_phi)
  ys = s.cderiv(latd,latbc(phi),phi) / @@R
  xs.long_name = "xgrad(#{s.name})"
  xs.name = "xgrad"
  ys.long_name = "ygrad(#{s.name})"
  ys.name = "ygrad"
  [xs,ys]
end

.grad_sx(s) ⇒ Object 



200
201
202
203
204
205
206
207
 
# File 'lib/numru/ganalysis/planet.rb', line 200

def grad_sx(s)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  xs = s.cderiv(lond,@@lonbc,lam) / (@@R*cos_phi)
  xs.long_name = "xgrad(#{s.name})"
  xs.name = "xgrad"
  xs
end

.grad_sy(s) ⇒ Object 



209
210
211
212
213
214
215
216
 
# File 'lib/numru/ganalysis/planet.rb', line 209

def grad_sy(s)
  phi, latd  = get_phi(s)
  cos_phi = phi.cos
  ys = s.cderiv(latd,latbc(phi),phi) / @@R
  ys.long_name = "ygrad(#{s.name})"
  ys.name = "ygrad"
  ys
end

.grad_sy_cosphifact(s, cosphi_exponent) ⇒ Object 



218
219
220
221
222
223
224
225
226
 
# File 'lib/numru/ganalysis/planet.rb', line 218

def grad_sy_cosphifact(s,cosphi_exponent)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  cos_phi_factor = cos_phi**cosphi_exponent
  ys = (s*cos_phi_factor).cderiv(latd,latbc(phi),phi)/cos_phi_factor / @@R
  ys.long_name = "ygrad(#{s.name})"
  ys.name = "ygrad"
  ys
end

.init(planet) ⇒ Object 



19
20
21
22
23
24
25
26
27
 
# File 'lib/numru/ganalysis/planet.rb', line 19

def init(planet)
  case planet
  when Earth
    @@R = UNumeric[6.37e6, "m"]
    @@Ome = UNumeric[2*Math::PI/8.64e4,"s-1"]
  else
    raise "Unsupported predefined planet: #{planet}."
  end
end

.latbc(phi) ⇒ Object

< Differentian at the planets’s near surface. With suffix “_s” >



48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
 
# File 'lib/numru/ganalysis/planet.rb', line 48

def latbc(phi)
=begin
  # not so good
  pi2 = Math::PI/2
  eps = 0.1
  xs = phi[0..1].val 
  xe = phi[-2..-1].val
  if (pi2-xs[0].abs) / (xs[0]-xs[1]).abs < eps and
     (pi2-xe[-1].abs) / (xe[-1]-xe[-2]).abs < eps 
    GPhys::Derivative::MIRROR_B
  else
    GPhys::Derivative::LINEAR_EXT
  end
=end

  GPhys::Derivative::LINEAR_EXT
end

.omega ⇒ Object 



36
 
# File 'lib/numru/ganalysis/planet.rb', line 36

def omega; @@Ome; end

.omega=(o) ⇒ Object 



34
 
# File 'lib/numru/ganalysis/planet.rb', line 34

def omega=(o);  @@Ome = o; end

.radius ⇒ Object 



35
 
# File 'lib/numru/ganalysis/planet.rb', line 35

def radius; @@R; end

.radius=(r) ⇒ Object

< Adjustable planetary settings >



33
 
# File 'lib/numru/ganalysis/planet.rb', line 33

def radius=(r); @@R = r; end

.rot_s(u, v) ⇒ Object 



76
77
78
79
80
81
82
83
84
85
 
# File 'lib/numru/ganalysis/planet.rb', line 76

def rot_s(u,v)
  lam, phi, lond, latd  = get_lambda_phi(u)
  cos_phi = phi.cos
  dv_dlam = v.cderiv(lond,@@lonbc,lam)
  duc_dphi  = (u*cos_phi).cderiv(latd,latbc(phi),phi)
  rot = (dv_dlam - duc_dphi) / (@@R*cos_phi)
  rot.long_name = "rot(#{u.name},#{v.name})"
  rot.name = "rot"
  rot
end

.vor_s(u, v) ⇒ Object 



174
175
176
177
178
179
 
# File 'lib/numru/ganalysis/planet.rb', line 174

def vor_s(u,v)
  vor = rot_s(u,v)
  vor.long_name = "Relative Vorticity"
  vor.name = "vor"
  vor
end

.weight_cosphi(s, cos_exp, r_exp) ⇒ Object 



235
236
237
238
239
240
 
# File 'lib/numru/ganalysis/planet.rb', line 235

def weight_cosphi(s, cos_exp, r_exp)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  ys = s * (cos_phi**cos_exp * @@R**r_exp)
  ys
end

.weight_sinphi(s, sin_exp, r_exp) ⇒ Object 



242
243
244
245
246
247
 
# File 'lib/numru/ganalysis/planet.rb', line 242

def weight_sinphi(s, sin_exp, r_exp)
  lam, phi, lond, latd  = get_lambda_phi(s)
  sin_phi = phi.sin
  ys = s * (sin_phi**sin_exp * @@R**r_exp)
  ys
end

.weight_tanphi(s, tan_exp, r_exp) ⇒ Object 



228
229
230
231
232
233
 
# File 'lib/numru/ganalysis/planet.rb', line 228

def weight_tanphi(s, tan_exp, r_exp)
  lam, phi, lond, latd  = get_lambda_phi(s)
  tan_phi = phi.tan
  ys = s * (tan_phi**tan_exp * @@R**r_exp)
  ys
end