search

dtcenter / MET

Model Evaluation Tools
https://dtcenter.org/community-code/model-evaluation-tools-met
Apache License 2.0
74 stars 22 forks source link

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

Closed JohnHalleyGotway closed 3 weeks ago

JohnHalleyGotway commented 4 months ago

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

Labels

Milestone and Projects

Define Related Issue(s)

Consider the impact to the other METplus components.

Bugfix Checklist

See the METplus Workflow for details.

JohnHalleyGotway commented 2 months ago

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 commented 2 months ago

@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 commented 2 months ago

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 commented 2 months ago

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

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

JohnHalleyGotway commented 2 months ago

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 commented 2 months ago

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:

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 commented 2 months ago

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

JohnHalleyGotway commented 2 months ago

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

JohnHalleyGotway commented 3 weeks ago

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