Calculate tangential and radial winds by locating the typhoon centre at the North Pole

rotate.ncl

pi = acos(-1.0)
tau = 2 * pi
deg2rad = pi / 180.0
rad2deg = 180.0 / pi

procedure lonlat2xyz(lon, lat, x, y, z)
begin
  x = cos(lon) * cos(lat)
  y = sin(lon) * cos(lat)
  z = sin(lat)
end

procedure xyz2lonlat(x, y, z, lon, lat)
begin
  lon = where(x .ne. 0.0, mod(atan2(y, x) + tau, tau), 0.0)
  lat= asin(z)
end

procedure uv2xyzd(u, v, lon, lat, xd, yd, zd)
local coslon, sinlon, sinlat
begin
  coslon = cos(lon)
  sinlon = sin(lon)
  sinlat = sin(lat)
  xd = -u * sinlon - v * coslon * sinlat
  yd =  u * coslon - v * sinlon * sinlat
  zd =               v * cos(lat)
end

procedure xyzd2uv(xd, yd, zd, lon, lat, u, v)
local coslon, sinlon
begin
  coslon = cos(lon)
  sinlon = sin(lon)
  coslat = cos(lat)
  sinlat = sin(lat)
  u = (/yd * coslon - xd * sinlon/)
  v = (/-(xd * coslon + yd * sinlon) * sinlat + zd * coslat/)
end

procedure generate_points(nlon, nlat, dlat, lon, lat)
local dlon
begin
  dlon = 2 * pi / nlon
  lon = dlon * ispan(0, nlon - 1, 1)
  lat = 0.5 * pi - dlat * ispan(0, nlat - 1, 1)
end

procedure np2tc(lonc, latc, x, y, z, xx, yy, zz)
local coslonc, sinlonc, coslatc, sinlatc
begin
  coslonc = cos(lonc)
  sinlonc = sin(lonc)
  coslatc = cos(latc)
  sinlatc = sin(latc)
  xx = coslonc * sinlatc * x - sinlonc * y + coslonc * coslatc * z
  yy = sinlonc * sinlatc * x + coslonc * y + sinlonc * coslatc * z
  zz =          -coslatc * x               +           sinlatc * z
end

procedure tc2np(lonc, latc, x, y, z, xx, yy, zz)
local coslonc, sinlonc, coslatc, sinlatc
begin
  coslonc = cos(lonc)
  sinlonc = sin(lonc)
  coslatc = cos(latc)
  sinlatc = sin(latc)
  xx = coslonc * sinlatc * x + sinlonc * sinlatc * y - coslatc * z
  yy =          -sinlonc * x +           coslonc * y
  zz = coslonc * coslatc * x + sinlonc * coslatc * y + sinlatc * z
end

procedure rotate_lonlat(lonc, latc, lon, lat, lonout, latout)
local x, y, z, xx, yy, zz
begin
  nlon = dimsizes(lon)
  nlat = dimsizes(lat)
  ij = 0
  do j = 0, nlat(0) - 1
    do i = 0, nlon(0) - 1
      x = 0.0
      y = 0.0
      z = 0.0
      lonlat2xyz(lon(i), lat(j), x, y, z)
      xx = 0.0
      yy = 0.0
      zz = 0.0
      np2tc(lonc, latc, x, y, z, xx, yy, zz)
      xyz2lonlat(xx, yy, zz, lonout(ij), latout(ij))
      ij = ij + 1
    end do
  end do
end

; # of polar grid points
nlon = 100
nlat =  11 
; input
infile="http://www.jamstec.go.jp/esc/fes/dods/alera2/stream2010/plev/fcst/mean/yyyymmddhh"
yyyymmddhh=2012091400
lonc = 128.5 * deg2rad
latc = 17.5 * deg2rad
uvarname = "u10"
vvarname = "v10"
; output
outfile = "outuv.grd"
begin
; generate polar grid
  dlat =   1.0 * deg2rad
  lonin = new(nlon, "float")
  latin = new(nlat, "float")
  generate_points(nlon, nlat, dlat, lonin, latin)
; (lonin, lonlat) => TC(lonout, latout)
  lonout = new(nlon * nlat , "float")
  latout = new(nlon * nlat , "float")
  rotate_lonlat(lonc, latc, lonin, latin, lonout, latout)
; read input
  f = addfile(infile, "r")
  time = cd_calendar(f->time, -3)
  t = closest_val(yyyymmddhh, time)
  lon = doubletofloat(f->lon) * deg2rad
  lat = doubletofloat(f->lat) * deg2rad
  u = f->$uvarname$(t, :, :)
  v = f->$vvarname$(t, :, :)
; (u, v) => (xdot, ydot, zdot)
  xd = u
  yd = u
  zd = u
  lon2d = conform(u, lon, 1)
  lat2d = conform(u, lat, 0)
  uv2xyzd(u, v, lon2d, lat2d, xd, yd, zd)
; bilinear interpolation of winds
  xdtc = linint2_points(lon, lat, xd, True, lonout, latout, 0)
  ydtc = linint2_points(lon, lat, yd, True, lonout, latout, 0)
  zdtc = linint2_points(lon, lat, zd, True, lonout, latout, 0)
; rotate coordinates to polar grid
  xdnp = xdtc
  ydnp = xdtc
  zdnp = xdtc
  lonc2d = conform(xdtc, lonc, -1)
  latc2d = conform(xdtc, latc, -1)
  tc2np(lonc2d, latc2d, xdtc, ydtc, zdtc, xdnp, ydnp, zdnp)
; (xdot, ydot, zdot) => (u, v)
  unp = xdnp
  vnp = xdnp
  lonnp = ndtooned(conform_dims((/nlat, nlon/), lonin, 1))
  xyzd2uv(xdnp, ydnp, zdnp, lonnp, latnp, unp, vnp)
; save rotated winds
  system("rm -f "+outfile) 
  setfileoption("bin", "WriteByteOrder", "BigEndian")
  fbindirwrite(outfile, unp)
  fbindirwrite(outfile, vnp)
end