Bugfix: Fix TC-RMW to correct the tangential and radial wind computations

Question

Bugfix: Fix TC-RMW to correct the tangential and radial wind computations

Closed this issue 6 months ago · 9 comments

Answer 1 · 2024-04-19T17:45:57.000Z

As discussed during the HAFS project meeting on April 19, 2024, a fix for this bug should be included in the MET-11.1.1 bugfix release, scheduled for May 1st.

Note that the existing tc_rmw use case computes tangential and radial winds output from HWRF GRIB2 Hurricane Gonzalo output.

@KathrynNewman will investigate running the HWRF GRIB2 data for this case through MetPy and/or DiaPost to generate "truth" data against which the TC-RMW output should be compared.

Answer 2 · 2024-04-25T19:26:17.000Z

@JohnHalleyGotway I used the script cross_section to plot tangential and radial winds from a lat/lon UKMO file. You may also look into MetPy's calculation.
Also, I have the Diapost script on Derecho (/glade/derecho/scratch/biswas/diapost.F90) if you want to take a look.

Answer 3 · 2024-05-01T23:11:26.000Z

I did find a seemingly unrelated bug.
When TC-RMW is configured to only write a single pressure level (e.g. "P500"), then no pressure dimension is added to the NetCDF output file. And the call to write_tc_pressure_level_data(...) writes data with the wrong dimensions. This behavior should be fixed but is not the underlying source of the bad output.
With multiple pressure levels:

	double pressure(pressure) ;
	double UGRD(track_point, pressure, range, azimuth) ;
	double VGRD(track_point, pressure, range, azimuth) ;
	double VR(track_point, pressure, range, azimuth) ;
	double VT(track_point, pressure, range, azimuth) ;

With a single pressure level:

	double pressure(pressure) ;
	double UGRD(track_point, range, azimuth) ;
	double VGRD(track_point, range, azimuth) ;
	double VR(track_point, range, azimuth) ;
	double VT(track_point, range, azimuth) ;

Recommend simplifying so that the pressure dimension is always used to define the data variables. This simplifies reading data from these files.

Answer 4 · 2024-05-03T17:15:48.000Z

Met with @KathrynNewman and @mrinalbiswas on 5/3/24.

John HG will test using a non-lat/lon grid to confirm that wind rotation for GRIB2 files really does work.

Here's a good reference from hafs-graphics:
https://github.com/hafs-community/hafs_graphics/blob/5cb6e5c6f9fa14e12c29b39e7443450f08b88491/ush/python/atmos/plot_azimuth_wind.py#L184

# Calculate tangential and radial wind
vt_p = np.ones((np.shape(XI)[0],np.shape(XI)[1],zsize))*np.nan
ur_p = np.ones((np.shape(XI)[0],np.shape(XI)[1],zsize))*np.nan
for j in range(np.shape(XI)[1]):
    for k in range(zsize):
        vt_p[:,j,k] = -u_p[:,j,k]*np.sin(theta)+v_p[:,j,k]*np.cos(theta)
        ur_p[:,j,k] = u_p[:,j,k]*np.cos(theta)+v_p[:,j,k]*np.sin(theta)

And this computation is very consistent with what's done by Diapost in derecho:/glade/derecho/scratch/biswas/diapost.F90:

 rcos = cos(phi) ; rsin = sin(phi) 
...
            wks(i,j,21) =   rcos*uur + rsin*vvr !*   radial   wind
            wks(i,j,22) = - rsin*uur + rcos*vvr !* tangential win

To see the output from the HAFS python plotting script, go to:
https://www.emc.ncep.noaa.gov/hurricane/HFSA/tcall.php

For Year 2023, Basin North Atlantic, select Storm IDALIA10L.
In the top right corner switch to Cycle IDALIA10L.2023082918, and under Storm-Scale-Atmos, select Storm: azimuthal mean wind.

Answer 5 · 2024-05-04T01:02:42.000Z

There seems to be a bug in how MET is indexing into the U and V values, as evidenced by the log messages below. For range = 0, (u, v) should remain constant for all azimuths. But, as shown below, they change:

DEBUG 4: wind_ne_to_rt() -> center (25.3, 84.8), range (km): 0, azimuth (deg): 0, point (25.3, 84.8), uv (-6.24, 21.6), radial wind: -6.24, tangential wind: 21.6
DEBUG 4: wind_ne_to_rt() -> center (25.3, 84.8), range (km): 0, azimuth (deg): 5, point (25.3, 84.8), uv (7.15, 6.89), radial wind: 7.72329, tangential wind: 6.24062
DEBUG 4: wind_ne_to_rt() -> center (25.3, 84.8), range (km): 0, azimuth (deg): 10, point (25.3, 84.8), uv (-9.06, 12.56), radial wind: -6.74134, tangential wind: 13.9424
DEBUG 4: wind_ne_to_rt() -> center (25.3, 84.8), range (km): 0, azimuth (deg): 15, point (25.3, 84.8), uv (7.47, 8.74), radial wind: 9.47754, tangential wind: 6.50881

Answer 6 · 2024-05-06T17:49:56.000Z

There appears to be an issue in the ordering of the data.

Using Gonzalo data for 2014101312 initialization, 6-hour forecast for testing:

AL, 08, 2014101312, 03, HCLT, 006, 176N,  626W,  59,  994, XX,  34, NEQ, 0068, 0044, 0033, 0065,  -99,  -99,  15,   0,   0,    ,   0,    ,   0,   0,           ,  ,   ,    ,   0,   0,   0,   0,       THERMO PARAMS,     -50,    1730,     366, N, 10, DT, -999

Storm center lat/lon is (17.6, -62.6).
Checking UGRD and VGRD at 850mb with ncview:

(lat, lon) (17.602, -62.6) = (i, j) (220, 254)
UGRD (220, 253) = -4.95
VGRD (220, 253) = 0.11

And that's very close to the regridded (u, v) value of (-4.707, 0.071). However, the problem is that this value is reported for range = 200 rather than range = 0, where it should exist:

DEBUG 4: wind_ne_to_rt() -> center lat/lon (17.6, 62.6), range (km): 200, azimuth (deg): 315, point lat/lon (18.8658, 61.2576), uv (-4.707, 0.071), radial wind: -3.37856, tangential wind: -3.27815

Answer 7 · 2024-05-07T19:31:52.000Z

I tested using Hurricane Idalia on seneca with following commands:

Run plot_azimuth_wind.py using the plot_atmos.yml config file in that directory:

cd /d1/projects/MET/MET_issues/bugfix_2841 # on seneca
conda activate plot_tcrmw
python plot_azimuth_wind.py

That creates the plot shown in the comment above. But I've modified it to dump out debug info.
2. Run TC-RMW using these same inputs.

./run_tc_rmw_idalia.sh

Plot the result from MET:

bash
source setup_env.bash
./run_plot_cross_section.sh

I think that MET and the HAFS plotting script differ by 90 degrees in their rotation angles.
For MET, 0 degrees is north (same longitude, latitude to the north):

DEBUG 4: wind_ne_to_rt() -> center lat/lon (25.3, 84.8), range (km): 200, azimuth (deg): 0, point lat/lon (27.0966, 84.8), uv (-9.84, 1.95), radial wind: -9.84, tangential wind: 1.95

In MET, the 200 km lat - center lat = 1.7966 deg * 111 km/deg = 200 km.

For HAFS 0, degrees is east (same latitude, longitude to the east):

rng: 200.0 km, lev: P 500
theta (deg): [  0.  ... ]
lat: [25.27999878 ... ]
lon: [83.90307905 ... ]

For HAFS, the 200 km lon - center lon = 0.896921 deg * 111 km/deg = 100km.

So MET is specifying the diameter of the circle while HAFS is specifying the radius. Need to modify the configuration for better testing.

Switching HAFS from 200km every 50km to 400km every 100km, I get:

rng: 400.0 km, lev: P 500
theta (deg): [  0.  ... ]
lat: [25.27999878 ... ]
lon: [83.00615346 ... ]

The "400km" lon - center lon = 1.7938466 deg * 111 km/deg = 200km.

Answer 8 · 2024-05-10T17:29:50.000Z

@JohnHalleyGotway and @KathrynNewman met on 5/10/24 to debug and found the following:

The Python diag_lib used for the MET TC-Diag tool, assumes that a rotation angle of 0 means due east based on the following lines of code:

    # The x,y coordinates of each cylindrical grid point in km relative to the
    # center of the grid.
    cyl_x_km = rad_2d_km * np.cos(theta_2d_radians)
    cyl_y_km = rad_2d_km * np.sin(theta_2d_radians)

theta_2d_radians = 0 yields a cylindrical point due east, not due north. If the sines and cosines were switched, we'd get the opposite.
So MET should change it definition of the cylindrical coordinates grid to match.
@JohnHalleyGotway will fix this which will change all existing TC-RMW and TC-Diag output.

Strongly suspect that the HAFS plot actually covers 200km rather than 400. Recommend starting a GitHub discussion here. @KathrynNewman will post a new discussion and also clarify that tangential winds are generally positive in the northern hemisphere and negative in the southern hemisphere (rather than using different signs for them).

Answer 9 · 2024-06-21T21:50:32.000Z

@JohnHalleyGotway and @KathrynNewman met to review this PRs #2920 and #2921 and found the following:

The lat/lon location for azimuth = 0 degrees is correct, meaning that it now points due east which matches the HAFS python plotting script logic.
However, while MET rotates through azimuths in a counter-clockwise direction, the HAFS plotting script does it clockwise.
Recommend updating MET to also rotate through the azimuths in a clockwise direction.

Bugfix: Fix TC-RMW to correct the tangential and radial wind computations

Describe the Problem

Expected Behavior

Environment

To Reproduce

Relevant Deadlines

Funding Source

Define the Metadata

Assignee

Labels

Milestone and Projects

Define Related Issue(s)

Bugfix Checklist