Skip to content

Commit

Permalink
Per #1849, fix to radial and tangential winds.
Browse files Browse the repository at this point in the history
  • Loading branch information
@JohnHalleyGotway
JohnHalleyGotway committed May 3, 2024
1 parent 2d7d8c1 commit 776dd72
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 109 deletions.
86 changes: 18 additions & 68 deletions src/libcode/vx_grid/tcrmw_grid.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -324,63 +324,23 @@ return;
////////////////////////////////////////////////////////////////////////


void TcrmwGrid::wind_ne_to_ra (const double lat, const double lon,
const double east_component, const double north_component,
double & radial_component, double & azimuthal_component) const
void TcrmwGrid::wind_ne_to_rt (const double azi_deg,
const double u_wind, const double v_wind,
double & radial_wind, double & tangential_wind) const

{

Vector E, N, V;
Vector B_range, B_azi;
double azi_deg, range_deg, range_km;

latlon_to_range_azi(lat, lon, range_km, azi_deg);

range_deg = deg_per_km*range_km;

E = latlon_to_east (lat, lon);
N = latlon_to_north (lat, lon);

V = east_component*E + north_component*N;

range_azi_to_basis(range_deg, azi_deg, B_range, B_azi);

radial_component = dot(V, B_range);

azimuthal_component = dot(V, B_azi);

/* JHG From diapost.F90
xxx = wcos(icen(k),jcen(k))
yyy = wsin(icen(k),jcen(k))
IF ( xxx > sqrt2_ ) THEN
alfa = Asin(yyy)
ELSEIF ( xxx < -sqrt2_ ) THEN
alfa = pi_cnst*Sign(1.,yyy) - Asin(yyy)
ELSE ! ( |xxx| <= sqrt2_ ) !
alfa = Acos(xxx)*Sign(1.,yyy)
ENDIF !* SOUTH-POLAR PROJECTION
phi = pi3_2 - alfa - i*dphi
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 wind
*/

return;
double rcos = cosd(azi_deg);
double rsin = sind(azi_deg);

if (is_bad_data(u_wind) || is_bad_data(v_wind)) {
radial_wind = bad_data_double;
tangential_wind = bad_data_double;
}
else {
radial_wind = rcos*u_wind + rsin*v_wind;
tangential_wind = -1.0*rsin*u_wind + rcos*v_wind;
}


////////////////////////////////////////////////////////////////////////


void TcrmwGrid::wind_ne_to_ra_conventional (const double lat, const double lon,
const double east_component, const double north_component,
double & radial_component, double & azimuthal_component) const

{

wind_ne_to_ra(lat, lon, east_component, north_component, radial_component, azimuthal_component);

return;

Expand All @@ -390,27 +350,17 @@ return;
////////////////////////////////////////////////////////////////////////


void TcrmwGrid::range_azi_to_basis(const double range_deg, const double azi_deg, Vector & B_range, Vector & B_azi) const
void TcrmwGrid::wind_ne_to_rt (const double lat, const double lon,
const double u_wind, const double v_wind,
double & radial_wind, double & tangential_wind) const

{

double u, v, w;
double range_km, azi_deg;

u = cosd(range_deg)*sind(azi_deg);

v = cosd(range_deg)*cosd(azi_deg);

w = -sind(range_deg);

B_range = u*Ir + v*Jr + w*Kr;

u = cosd(azi_deg);

v = -sind(azi_deg);

w = 0.0;
latlon_to_range_azi(lat, lon, range_km, azi_deg);

B_azi = u*Ir + v*Jr + w*Kr;
wind_ne_to_rt(azi_deg, u_wind, v_wind, radial_wind, tangential_wind);

return;

Expand Down
22 changes: 7 additions & 15 deletions src/libcode/vx_grid/tcrmw_grid.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,11 +35,8 @@ class TcrmwGrid : public RotatedLatLonGrid {

void calc_ijk(); // calculate rotated basis vectors

void range_azi_to_basis(const double range_deg, const double azi_deg, Vector & B_range, Vector & B_azi) const;

TcrmwData TData;


Vector Ir, Jr, Kr;

int Range_n, Azimuth_n; // # of points in the radial and azimuthal directions
Expand Down Expand Up @@ -89,20 +86,15 @@ class TcrmwGrid : public RotatedLatLonGrid {

void xy_to_latlon(double x, double y, double & true_lat, double & true_lon) const;

void wind_ne_to_ra(const double lat, const double lon,
const double east_component, const double north_component,
double & radial_component, double & azimuthal_component) const;
void wind_ne_to_rt(const double azi_deg,
const double u_wind, const double v_wind,
double & radial_wind, double & tangential_wind) const;

void wind_ne_to_rt(const double lat, const double lon,
const double u_wind, const double v_wind,
double & radial_wind, double & tangential_wind) const;

//
// possibly toggles the signs of the radial and/or azimuthal components
//
// to align with the conventions used in the TC community
//

void wind_ne_to_ra_conventional (const double lat, const double lon,
const double east_component, const double north_component,
double & radial_component, double & azimuthal_component) const;

};


Expand Down
49 changes: 23 additions & 26 deletions src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
// ---- ---- ---- -----------
// 000 05/11/22 Albo Pulled the wind conversion into a class
// 001 09/28/22 Prestopnik MET #2227 Remove namespace std from header files
// 002 05/03/24 Halley Gotway MET #2841 Fix radial and tangential winds
//
////////////////////////////////////////////////////////////////////////

Expand All @@ -28,7 +29,7 @@ using namespace std;

////////////////////////////////////////////////////////////////////////

static void wind_ne_to_ra(const TcrmwGrid&,
static void wind_ne_to_rt(const TcrmwGrid&,
const DataPlane&, const DataPlane&,
const double*, const double*,
double*, double*);
Expand Down Expand Up @@ -206,45 +207,41 @@ bool TCRMW_WindConverter::compute_winds_if_input_is_u(int i_point,
<< "), regrid range (" << dmin_rgd << ", " << dmax_rgd << ")\n";

// Compute the radial and tangential winds and store in _windR and _windT
wind_ne_to_ra(tcrmw_grid, data_dp, data_dpV, lat_arr, lon_arr,
wind_ne_to_rt(tcrmw_grid, data_dp, data_dpV, lat_arr, lon_arr,
_windR, _windT);

return true;
}

////////////////////////////////////////////////////////////////////////

void wind_ne_to_ra(const TcrmwGrid& tcrmw_grid,
void wind_ne_to_rt(const TcrmwGrid& tcrmw_grid,
const DataPlane& u_dp, const DataPlane& v_dp,
const double* lat_arr, const double* lon_arr,
double* wind_r_arr, double* wind_t_arr) {

// Transform (u, v) to (radial, azimuthal)
// Transform (u, v) to (radial, tangential) winds
for(int ir = 0; ir < tcrmw_grid.range_n(); ir++) {
for(int ia = 0; ia < tcrmw_grid.azimuth_n(); ia++) {

int i = ir * tcrmw_grid.azimuth_n() + ia;
double lat = lat_arr[i];
double lon = -lon_arr[i]; // convert degrees east to west
double u = u_dp.data()[i];
double v = v_dp.data()[i];
double wind_r;
double wind_t;
if(is_bad_data(u) || is_bad_data(v)) {
mlog << Debug(4) << "wind_ne_to_ra() -> "
<< "latlon (" << lat << "," << lon
<< ") winds are missing\n";
wind_r = bad_data_double;
wind_t = bad_data_double;
} else {
tcrmw_grid.wind_ne_to_ra(lat, lon, u, v, wind_r, wind_t);
mlog << Debug(4) << "wind_ne_to_ra() -> "
<< "latlon (" << lat << ", " << lon
<< "), uv (" << u << ", " << v
<< "), radial wind: " << wind_r
<< ", tangential wind: " << wind_t << "\n";
}
wind_r_arr[i] = wind_r;
wind_t_arr[i] = wind_t;

double azi_deg = ia * tcrmw_grid.azimuth_delta_deg();
double range_km = ir * tcrmw_grid.range_delta_km();

tcrmw_grid.wind_ne_to_rt(azi_deg, u_dp.data()[i], v_dp.data()[i],
wind_r_arr[i], wind_t_arr[i]);

mlog << Debug(4) << "wind_ne_to_rt() -> "
<< "center (" << tcrmw_grid.lat_center_deg()
<< ", " << tcrmw_grid.lon_center_deg()
<< "), range (km): " << range_km
<< ", azimuth (deg): " << azi_deg
<< ", point (" << lat_arr[i] << ", " << -lon_arr[i]
<< "), uv (" << u_dp.data()[i] << ", " << v_dp.data()[i]
<< "), radial wind: " << wind_r_arr[i]
<< ", tangential wind: " << wind_t_arr[i] << "\n";

} // end for ia
} // end for ir

Expand Down

0 comments on commit 776dd72

Please sign in to comment.