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

[JohnHalleyGotway]

Describe the Problem

The computation of radial and tangential winds in the TC-RMW tool are incorrect and should be fixed.

Expected Behavior

Provide a clear and concise description of what you expected to happen here.

Environment

Describe your runtime environment:
1. Machine: (e.g. HPC name, Linux Workstation, Mac Laptop)
2. OS: (e.g. RedHat Linux, MacOS)
3. Software version number(s)

To Reproduce

Describe the steps to reproduce the behavior:
1. Go to '...'
2. Click on '....'
3. Scroll down to '....'
4. See error
Post relevant sample data following these instructions:
https://dtcenter.org/community-code/model-evaluation-tools-met/met-help-desk#ftp

Relevant Deadlines

List relevant project deadlines here or state NONE.

Funding Source

@JohnHalleyGotway should charge 2784603 for this bugfix.

Define the Metadata

Assignee

• Select engineer(s) or no engineer required
• Select scientist(s) or no scientist required

Labels

• Review default alert labels
• Select component(s)
• Select priority
• Select requestor(s)

Milestone and Projects

• Select Milestone as the next bugfix version
• Select Coordinated METplus-X.Y Support project for support of the current coordinated release
• Select MET-X.Y.Z Development project for development toward the next official release

Define Related Issue(s)

Consider the impact to the other METplus components.

• METplus, MET, METdataio, METviewer, METexpress, METcalcpy, METplotpy

Bugfix Checklist

See the METplus Workflow for details.

• Complete the issue definition above, including the Time Estimate and Funding Source.
• Fork this repository or create a branch of main_<Version>.
  Branch name: bugfix_<Issue Number>_main_<Version>_<Description>
• Fix the bug and test your changes.
• Add/update log messages for easier debugging.
• Add/update unit tests.
• Add/update documentation.
• Push local changes to GitHub.
• Submit a pull request to merge into main_<Version>.
  Pull request: bugfix <Issue Number> main_<Version> <Description>
• Define the pull request metadata, as permissions allow.
  Select: Reviewer(s) and Development issue
  Select: Milestone as the next bugfix version
  Select: Coordinated METplus-X.Y Support project for support of the current coordinated release
• Iterate until the reviewer(s) accept and merge your changes.
• Delete your fork or branch.
• Complete the steps above to fix the bug on the develop branch.
  Branch name: bugfix_<Issue Number>_develop_<Description>
  Pull request: bugfix <Issue Number> develop <Description>
  Select: Reviewer(s) and Development issue
  Select: Milestone as the next official version
  Select: MET-X.Y.Z Development project for development toward the next official release
• Close this issue.

[JohnHalleyGotway]

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.

[mrinalbiswas]

@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.

[JohnHalleyGotway]

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.

[JohnHalleyGotway]

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.
  

[JohnHalleyGotway]

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

[JohnHalleyGotway]

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

[JohnHalleyGotway]

I tested using Hurricane Idalia on seneca with following commands:

1. 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

3. Plot the result from MET:

bash
source setup_env.bash
./test_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.

[JohnHalleyGotway]

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

1. 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.

2. 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).