Last active
January 9, 2023 01:20
-
-
Save tenomoto/f072310cc49cf83c15cb7a2a9fe15c2b to your computer and use it in GitHub Desktop.
Calculate tangential and radial winds by locating the typhoon centre at the North Pole
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment