From 43b41192a180392cc225840b2b23e413909effca Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Wed, 1 May 2024 11:19:41 -0600 Subject: [PATCH 01/15] Per #2841, unrelated, just removing a spurious character from a log message. --- src/libcode/vx_series_data/series_data.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libcode/vx_series_data/series_data.cc b/src/libcode/vx_series_data/series_data.cc index d7a0e587d9..d97015a967 100644 --- a/src/libcode/vx_series_data/series_data.cc +++ b/src/libcode/vx_series_data/series_data.cc @@ -49,7 +49,7 @@ void get_series_entry(int i_series, VarInfo* data_info, if(!found) { mlog << Error << "\nget_series_entry() -> " << "Could not find data for " << data_info->magic_time_str() - << " in file list:\n:" << write_css(search_files) << "\n\n"; + << " in file list:\n" << write_css(search_files) << "\n\n"; exit(1); } From 03e29fd66d6f677d41f048d8dba1ed51f89b4b19 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Wed, 1 May 2024 12:07:23 -0600 Subject: [PATCH 02/15] Per #2841, work in progress. Mostly just modifying whitespace and log messages so far --- src/libcode/vx_grid/tcrmw_grid.cc | 36 +++-- src/libcode/vx_grid/tcrmw_grid.h | 3 - src/tools/tc_utils/tc_rmw/tc_rmw.cc | 53 ++++--- .../tc_utils/tc_rmw/tc_rmw_wind_converter.cc | 135 ++++++++++-------- .../tc_utils/tc_rmw/tc_rmw_wind_converter.h | 29 ++-- 5 files changed, 128 insertions(+), 128 deletions(-) diff --git a/src/libcode/vx_grid/tcrmw_grid.cc b/src/libcode/vx_grid/tcrmw_grid.cc index 7ab34aa415..635f61ec8c 100644 --- a/src/libcode/vx_grid/tcrmw_grid.cc +++ b/src/libcode/vx_grid/tcrmw_grid.cc @@ -334,7 +334,6 @@ 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; @@ -344,18 +343,28 @@ 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); +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; @@ -387,28 +396,22 @@ void TcrmwGrid::range_azi_to_basis(const double range_deg, const double azi_deg, double u, v, w; - 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; - B_azi = u*Ir + v*Jr + w*Kr; - return; } @@ -427,7 +430,6 @@ x -= Nx*floor(x/Nx); x -= Nx*floor(x/Nx); - return; } @@ -500,7 +502,3 @@ return; //////////////////////////////////////////////////////////////////////// - - - - diff --git a/src/libcode/vx_grid/tcrmw_grid.h b/src/libcode/vx_grid/tcrmw_grid.h index c102332408..01d0a586f1 100644 --- a/src/libcode/vx_grid/tcrmw_grid.h +++ b/src/libcode/vx_grid/tcrmw_grid.h @@ -89,13 +89,10 @@ 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; - // // possibly toggles the signs of the radial and/or azimuthal components // diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw.cc b/src/tools/tc_utils/tc_rmw/tc_rmw.cc index 5479136ed0..daf4d82535 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw.cc @@ -67,22 +67,21 @@ static bool file_is_ok(const ConcatString &, const GrdFileType); static void process_rmw(); static void process_tracks(TrackInfoArray&); static void get_atcf_files(const StringArray&, - const StringArray&, StringArray&, StringArray&); + const StringArray&, StringArray&, StringArray&); static void process_track_files(const StringArray&, - const StringArray&, TrackInfoArray&); + const StringArray&, TrackInfoArray&); static bool is_keeper(const ATCFLineBase *); static void set_deck(const StringArray&); static void set_atcf_source(const StringArray&, - StringArray&, StringArray&); + StringArray&, StringArray&); static void set_data_files(const StringArray&); static void set_config(const StringArray&); static void set_out(const StringArray&); static void setup_grid(); static void setup_nc_file(); static void build_outfile_name(const ConcatString&, - const char*, ConcatString&); -static void compute_lat_lon(TcrmwGrid&, - double*, double*); + const char*, ConcatString&); +static void compute_lat_lon(TcrmwGrid&, double*, double*); static void process_fields(const TrackInfoArray&); //////////////////////////////////////////////////////////////////////// @@ -612,11 +611,11 @@ void setup_nc_file() { // Get VarInfo data_info = conf_info.data_info[i_var]; mlog << Debug(4) << "Processing field: " << data_info->magic_str() << "\n"; - string fname = data_info->name_attr(); + string fname = data_info->name_attr(); variable_levels[fname].push_back(data_info->level_attr()); variable_long_names[fname] = data_info->long_name_attr(); variable_units[fname] = data_info->units_attr(); - wind_converter.update_input(fname, data_info->units_attr()); + wind_converter.update_input(fname, data_info->units_attr()); } // Define pressure levels @@ -650,7 +649,7 @@ void compute_lat_lon(TcrmwGrid& tcrmw_grid, ia * tcrmw_grid.azimuth_delta_deg(), lat, lon); lat_arr[i] = lat; - lon_arr[i] = - lon; + lon_arr[i] = -lon; // switch degrees west to east } } } @@ -736,31 +735,29 @@ void process_fields(const TrackInfoArray& tracks) { // Find data for this track point get_series_entry(i_point, data_info, data_files, ftype, data_dp, latlon_arr); - // Check data range - double data_min, data_max; - data_dp.data_range(data_min, data_max); - mlog << Debug(4) << "data_min:" << data_min << "\n"; - mlog << Debug(4) << "data_max:" << data_max << "\n"; - - // Regrid data + // Regrid data and log the range of values before and after + double dmin, dmax, dmin_rgd, dmax_rgd; + data_dp.data_range(dmin, dmax); data_dp = met_regrid(data_dp, latlon_arr, grid, data_info->regrid()); - data_dp.data_range(data_min, data_max); - mlog << Debug(4) << "data_min:" << data_min << "\n"; - mlog << Debug(4) << "data_max:" << data_max << "\n"; + data_dp.data_range(dmin_rgd, dmax_rgd); + + mlog << Debug(4) << data_info->magic_str() + << " input range (" << dmin << ", " << dmax + << "), regrid range (" << dmin_rgd << ", " << dmax_rgd << ")\n"; - // if this is "U", setup everything for matching "V" and compute the radial/tangential - if(wind_converter.compute_winds_if_input_is_u(i_point, sname, slevel, valid_time, data_files, ftype, + // if this is "U", setup everything for matching "V" and compute the radial/tangential + if(wind_converter.compute_winds_if_input_is_u(i_point, sname, slevel, valid_time, data_files, ftype, latlon_arr, lat_arr, lon_arr, grid, data_dp, tcrmw_grid)) { write_tc_pressure_level_data(nc_out, tcrmw_grid, pressure_level_indices, data_info->level_attr(), i_point, - data_3d_vars[conf_info.radial_velocity_field_name.string()], - wind_converter.get_wind_r_arr()); - write_tc_pressure_level_data(nc_out, tcrmw_grid, - pressure_level_indices, data_info->level_attr(), i_point, - data_3d_vars[conf_info.tangential_velocity_field_name.string()], - wind_converter.get_wind_t_arr()); + data_3d_vars[conf_info.radial_velocity_field_name.string()], + wind_converter.get_wind_r_arr()); + write_tc_pressure_level_data(nc_out, tcrmw_grid, + pressure_level_indices, data_info->level_attr(), i_point, + data_3d_vars[conf_info.tangential_velocity_field_name.string()], + wind_converter.get_wind_t_arr()); } - + // Write data if(variable_levels[data_info->name_attr()].size() > 1) { write_tc_pressure_level_data(nc_out, tcrmw_grid, diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc index f94994b314..b48f572021 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc @@ -29,13 +29,13 @@ using namespace std; //////////////////////////////////////////////////////////////////////// static void wind_ne_to_ra(const TcrmwGrid&, - const DataPlane&, const DataPlane&, - const double*, const double*, double*, double*); + const DataPlane&, const DataPlane&, + const double*, const double*, + double*, double*); //////////////////////////////////////////////////////////////////////// -void TCRMW_WindConverter::_free_winds_arrays(void) -{ +void TCRMW_WindConverter::_free_winds_arrays(void) { if (_windR != nullptr) { delete [] _windR; _windR = nullptr; @@ -94,6 +94,7 @@ void TCRMW_WindConverter::init(const TCRMWConfInfo *conf) { _vIndexMap[varlevel] = i_var; } } + // test for consistency if (_uIndexMap.size() != _vIndexMap.size()) { mlog << Warning << "Uneven number of u/v wind inputs, no wind conversion will be done:\n" @@ -104,7 +105,8 @@ void TCRMW_WindConverter::init(const TCRMWConfInfo *conf) { map::const_iterator iu, iv; for (iu=_uIndexMap.begin(), iv=_vIndexMap.begin(); iu!=_uIndexMap.end(); ++iu, ++iv) { if (iu->first != iv->first) { - mlog << Warning << "Ordering of u/v wind input levels not the same, not implemented, no wind conversions will be done:\n" + mlog << Warning << "Ordering of u/v wind input levels not the same, " + << "not implemented, no wind conversions will be done:\n" << " " << iu->first << " " << iv->first << "\n"; _computeWinds = false; } @@ -128,11 +130,9 @@ void TCRMW_WindConverter::update_input(const string &variableName, const string //////////////////////////////////////////////////////////////////////// void TCRMW_WindConverter::append_nc_output_vars(map > &variable_levels, - map &variable_long_names, - map &variable_units) { - if (!_computeWinds) { - return; - } + map &variable_long_names, + map &variable_units) { + if (!_computeWinds) return; if (_foundUInInput && _foundVInInput) { variable_levels[_conf->tangential_velocity_field_name] = variable_levels[_conf->u_wind_field_name.string()]; @@ -145,13 +145,14 @@ void TCRMW_WindConverter::append_nc_output_vars(map > &va else { if (!_foundUInInput) { mlog << Warning << "\nTCWRMW_WindConverter::checkInputs() -> " - << "field not found in input \"" << _conf->u_wind_field_name << "\"\n\n"; + << "field not found in input \"" << _conf->u_wind_field_name << "\"\n\n"; } if (!_foundVInInput) { mlog << Warning << "\nTCWRMW_WindConverter::checkInputs() -> " - << "field not found in input \"" << _conf->v_wind_field_name << "\"\n\n"; + << "field not found in input \"" << _conf->v_wind_field_name << "\"\n\n"; } - mlog << Warning << "\nNot computing radial and tangential winds\n\n"; + mlog << Warning << "\nTCWRMW_WindConverter::checkInputs() -> " + << "Not computing radial and tangential winds\n\n"; _computeWinds = false; } } @@ -159,17 +160,17 @@ void TCRMW_WindConverter::append_nc_output_vars(map > &va //////////////////////////////////////////////////////////////////////// bool TCRMW_WindConverter::compute_winds_if_input_is_u(int i_point, - const string &varName, - const string &varLevel, - unixtime valid_time, - const StringArray &data_files, - const GrdFileType &ftype, - const Grid &latlon_arr, - const double *lat_arr, - const double *lon_arr, - const Grid &grid, - const DataPlane &data_dp, - const TcrmwGrid &tcrmw_grid) { + const string &varName, + const string &varLevel, + unixtime valid_time, + const StringArray &data_files, + const GrdFileType &ftype, + const Grid &latlon_arr, + const double *lat_arr, + const double *lon_arr, + const Grid &grid, + const DataPlane &data_dp, + const TcrmwGrid &tcrmw_grid) { if (!_computeWinds) { return false; } @@ -190,53 +191,61 @@ bool TCRMW_WindConverter::compute_winds_if_input_is_u(int i_point, DataPlane data_dpV; Grid latlon_arrV; - get_series_entry(i_point, data_infoV, data_files, ftype, data_dpV, - latlon_arrV); - double data_min, data_max; - data_dpV.data_range(data_min, data_max); - mlog << Debug(4) << "V data_min:" << data_min << "\n"; - mlog << Debug(4) << "V data_max:" << data_max << "\n"; + get_series_entry(i_point, data_infoV, data_files, ftype, + data_dpV, latlon_arrV); + double dmin, dmax, dmin_rgd, dmax_rgd; + data_dpV.data_range(dmin, dmax); data_dpV = met_regrid(data_dpV, latlon_arr, grid, data_infoV->regrid()); - data_dpV.data_range(data_min, data_max); - mlog << Debug(4) << "V data_min:" << data_min << "\n"; - mlog << Debug(4) << "V data_max:" << data_max << "\n"; + data_dpV.data_range(dmin_rgd, dmax_rgd); + + mlog << Debug(4) << data_infoV->magic_str() + << " input range (" << dmin << ", " << dmax + << "), regrid range (" << dmin_rgd << ", " << dmax_rgd << ")\n"; - // here's the conversion, at last + // 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, - _windR, _windT); - // _windR and _windT now set + _windR, _windT); + return true; } - //////////////////////////////////////////////////////////////////////// void wind_ne_to_ra(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) - 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]; - 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 << ", rt:" << wind_r << "," << wind_t <<"\n"; - } - wind_r_arr[i] = wind_r; - wind_t_arr[i] = wind_t; - } - } + 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) + 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; + } // end for ia + } // end for ir + + return; } +//////////////////////////////////////////////////////////////////////// diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h index 8ad0642a17..f87d86193a 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h @@ -15,7 +15,7 @@ // // Mod# Date Name Description // ---- ---- ---- ----------- -// 000 05/11/22 DaveAlbo New +// 000 05/11/22 Albo New // //////////////////////////////////////////////////////////////////////// @@ -32,7 +32,6 @@ using std::map; using std::string; - //////////////////////////////////////////////////////////////////////// // // Constants @@ -98,8 +97,8 @@ class TCRMW_WindConverter { // if configured to compute winds, but didn't find U or V, turn off // the wind computations and report an error void append_nc_output_vars(std::map > &variable_levels, - std::map &variable_long_names, - std::map &variable_units); + std::map &variable_long_names, + std::map &variable_units); // Check input varName against U, and if it's a match, lookup V using the // map members, and then compute tangential and radial winds if it is so @@ -109,17 +108,17 @@ class TCRMW_WindConverter { // If true if returned, the winds can be accessed by calls to // get_wind_t_arr() and get_wind_r_arr() bool compute_winds_if_input_is_u(int i_point, - const string &varName, - const string &varLevel, - unixtime valid_time, - const StringArray &data_files, - const GrdFileType &ftype, - const Grid &latlon_arr, - const double *lat_arr, - const double *lon_arr, - const Grid &grid, - const DataPlane &data_dp, - const TcrmwGrid &tcrmw_grid); + const string &varName, + const string &varLevel, + unixtime valid_time, + const StringArray &data_files, + const GrdFileType &ftype, + const Grid &latlon_arr, + const double *lat_arr, + const double *lon_arr, + const Grid &grid, + const DataPlane &data_dp, + const TcrmwGrid &tcrmw_grid); }; From 6c3b9428a32074fa5da9fa15f685a8613808858a Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 2 May 2024 17:38:21 +0000 Subject: [PATCH 03/15] Per #2841, running tc_rmw with only UGRD input causes the tool to hang. Adding break statements to prevent the hang. --- .../tc_utils/tc_rmw/tc_rmw_wind_converter.cc | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc index b48f572021..f14665ea41 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc @@ -102,13 +102,16 @@ void TCRMW_WindConverter::init(const TCRMWConfInfo *conf) { << _conf->v_wind_field_name.string() << " has " << _vIndexMap.size() << " inputs\n"; _computeWinds = false; } - map::const_iterator iu, iv; - for (iu=_uIndexMap.begin(), iv=_vIndexMap.begin(); iu!=_uIndexMap.end(); ++iu, ++iv) { - if (iu->first != iv->first) { - mlog << Warning << "Ordering of u/v wind input levels not the same, " - << "not implemented, no wind conversions will be done:\n" - << " " << iu->first << " " << iv->first << "\n"; - _computeWinds = false; + if (_computeWinds) { + map::const_iterator iu, iv; + for (iu=_uIndexMap.begin(), iv=_vIndexMap.begin(); iu!=_uIndexMap.end(); ++iu, ++iv) { + if (iu->first != iv->first) { + mlog << Warning << "Ordering of u/v wind input levels not the same, " + << "not implemented, no wind conversions will be done:\n" + << " " << iu->first << " " << iv->first << "\n"; + _computeWinds = false; + break; + } } } } From 2d7d8c1661f6d2a631aaefd95d7b3d03c7e004a5 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 2 May 2024 17:53:42 +0000 Subject: [PATCH 04/15] Per #2841, add new has_pressure_level() utility function and update tc_rmw to use it to deteremine whether or not the pressure dimension should be written the output rather than basing that off the number of levels being > 1. --- src/libcode/vx_tc_util/vx_tc_nc_util.cc | 17 ++++++++++++++++- src/libcode/vx_tc_util/vx_tc_nc_util.h | 2 ++ src/tools/tc_utils/tc_rmw/tc_rmw.cc | 2 +- 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/libcode/vx_tc_util/vx_tc_nc_util.cc b/src/libcode/vx_tc_util/vx_tc_nc_util.cc index 4c5297d8a1..f2f545cdef 100644 --- a/src/libcode/vx_tc_util/vx_tc_nc_util.cc +++ b/src/libcode/vx_tc_util/vx_tc_nc_util.cc @@ -158,6 +158,21 @@ void write_tc_rmw(NcFile* nc_out, //////////////////////////////////////////////////////////////////////// +bool has_pressure_level(vector levels) { + + bool status = false; + + for (int j = 0; j < levels.size(); j++) { + if (levels[j].substr(0, 1) == "P") { + status = true; + break; + } + } + + return status; +} + +//////////////////////////////////////////////////////////////////////// set get_pressure_level_strings( map > variable_levels) { @@ -538,7 +553,7 @@ void def_tc_variables(NcFile* nc_out, string long_name = variable_long_names[i->first]; string units = variable_units[i->first]; - if (levels.size() > 1) { + if (has_pressure_level(levels)) { data_var = nc_out->addVar( var_name, ncDouble, dims_3d); add_att(&data_var, "long_name", long_name); diff --git a/src/libcode/vx_tc_util/vx_tc_nc_util.h b/src/libcode/vx_tc_util/vx_tc_nc_util.h index 981d55e08e..e0239b531f 100644 --- a/src/libcode/vx_tc_util/vx_tc_nc_util.h +++ b/src/libcode/vx_tc_util/vx_tc_nc_util.h @@ -35,6 +35,8 @@ extern void write_tc_track_point(netCDF::NcFile*, extern void write_tc_rmw(netCDF::NcFile*, const netCDF::NcDim&, const TrackInfo&); +extern bool has_pressure_level(std::vector); + extern std::set get_pressure_level_strings( std::map >); diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw.cc b/src/tools/tc_utils/tc_rmw/tc_rmw.cc index daf4d82535..728e11a28f 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw.cc @@ -759,7 +759,7 @@ void process_fields(const TrackInfoArray& tracks) { } // Write data - if(variable_levels[data_info->name_attr()].size() > 1) { + if(has_pressure_level(variable_levels[data_info->name_attr()])) { write_tc_pressure_level_data(nc_out, tcrmw_grid, pressure_level_indices, data_info->level_attr(), i_point, data_3d_vars[data_info->name_attr()], data_dp.data()); From 776dd72e1232122183eb8b38876397ab0290b720 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Fri, 3 May 2024 20:51:16 +0000 Subject: [PATCH 05/15] Per #1849, fix to radial and tangential winds. --- src/libcode/vx_grid/tcrmw_grid.cc | 86 ++++--------------- src/libcode/vx_grid/tcrmw_grid.h | 22 ++--- .../tc_utils/tc_rmw/tc_rmw_wind_converter.cc | 49 +++++------ 3 files changed, 48 insertions(+), 109 deletions(-) diff --git a/src/libcode/vx_grid/tcrmw_grid.cc b/src/libcode/vx_grid/tcrmw_grid.cc index 635f61ec8c..6e9b105729 100644 --- a/src/libcode/vx_grid/tcrmw_grid.cc +++ b/src/libcode/vx_grid/tcrmw_grid.cc @@ -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; @@ -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; diff --git a/src/libcode/vx_grid/tcrmw_grid.h b/src/libcode/vx_grid/tcrmw_grid.h index 01d0a586f1..3757111e9e 100644 --- a/src/libcode/vx_grid/tcrmw_grid.h +++ b/src/libcode/vx_grid/tcrmw_grid.h @@ -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 @@ -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; - }; diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc index f14665ea41..7aec340e56 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc @@ -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 // //////////////////////////////////////////////////////////////////////// @@ -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*); @@ -206,7 +207,7 @@ 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; @@ -214,37 +215,33 @@ bool TCRMW_WindConverter::compute_winds_if_input_is_u(int i_point, //////////////////////////////////////////////////////////////////////// -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 From fb0fab0963148395ad4164c6a92675c5f5f74601 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Mon, 6 May 2024 19:20:51 +0000 Subject: [PATCH 06/15] Per #2841, just changing a code comment about the longitude swap --- src/tools/tc_utils/tc_diag/tc_diag.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tools/tc_utils/tc_diag/tc_diag.cc b/src/tools/tc_utils/tc_diag/tc_diag.cc index c808faa29a..11bb4c04e8 100644 --- a/src/tools/tc_utils/tc_diag/tc_diag.cc +++ b/src/tools/tc_utils/tc_diag/tc_diag.cc @@ -847,7 +847,7 @@ void compute_lat_lon(TcrmwGrid& grid, ia * grid.azimuth_delta_deg(), lat, lon); lat_arr[i] = lat; - lon_arr[i] = -lon; // degrees east to west + lon_arr[i] = -lon; // degrees west to east } } From a086c2455b0cf59bcdc8cff4fc4126d3ab1f6ecc Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Mon, 6 May 2024 19:36:58 +0000 Subject: [PATCH 07/15] Per #2841, update variable names for consistency and understanding. --- src/tools/tc_utils/tc_rmw/tc_rmw.cc | 34 ++++++++++++++--------------- src/tools/tc_utils/tc_rmw/tc_rmw.h | 12 +++------- 2 files changed, 20 insertions(+), 26 deletions(-) diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw.cc b/src/tools/tc_utils/tc_rmw/tc_rmw.cc index 728e11a28f..5c69775609 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw.cc @@ -549,19 +549,19 @@ void set_out(const StringArray& a) { void setup_grid() { - grid_data.name = "TCRMW"; - grid_data.range_n = conf_info.n_range; - grid_data.azimuth_n = conf_info.n_azimuth; + tcrmw_data.name = "TCRMW"; + tcrmw_data.range_n = conf_info.n_range; + tcrmw_data.azimuth_n = conf_info.n_azimuth; // Define the maximum range in km based on the fixed increment if(is_bad_data(conf_info.rmw_scale)) { - grid_data.range_max_km = + tcrmw_data.range_max_km = conf_info.delta_range_km * (conf_info.n_range - 1); } - tcrmw_grid.set_from_data(grid_data); - grid.set(grid_data); + tcrmw_grid.set_from_data(tcrmw_data); + grid_out.set(tcrmw_data); } //////////////////////////////////////////////////////////////////////// @@ -648,8 +648,8 @@ void compute_lat_lon(TcrmwGrid& tcrmw_grid, ir * tcrmw_grid.range_delta_km(), ia * tcrmw_grid.azimuth_delta_deg(), lat, lon); - lat_arr[i] = lat; - lon_arr[i] = -lon; // switch degrees west to east + lat_arr[i] = lat; + lon_arr[i] = -lon; // degrees west to east } } } @@ -690,12 +690,12 @@ void process_fields(const TrackInfoArray& tracks) { << point.lon() << ").\n"; // Set grid center - grid_data.lat_center = point.lat(); - grid_data.lon_center = -1.0*point.lon(); // internal sign change + tcrmw_data.lat_center = point.lat(); + tcrmw_data.lon_center = -1.0*point.lon(); // internal sign change // Define the maximum range in km relative to the radius of maximum winds if(!is_bad_data(conf_info.rmw_scale)) { - grid_data.range_max_km = + tcrmw_data.range_max_km = conf_info.rmw_scale * point.mrd() * tc_km_per_nautical_miles * (conf_info.n_range - 1); @@ -703,9 +703,9 @@ void process_fields(const TrackInfoArray& tracks) { // Re-define the range/azimuth grid tcrmw_grid.clear(); - tcrmw_grid.set_from_data(grid_data); - grid.clear(); - grid.set(grid_data); + tcrmw_grid.set_from_data(tcrmw_data); + grid_out.clear(); + grid_out.set(tcrmw_data); // Compute lat and lon coordinate arrays compute_lat_lon(tcrmw_grid, lat_arr, lon_arr); @@ -733,12 +733,12 @@ void process_fields(const TrackInfoArray& tracks) { data_info->set_valid(valid_time); // Find data for this track point - get_series_entry(i_point, data_info, data_files, ftype, data_dp, latlon_arr); + get_series_entry(i_point, data_info, data_files, ftype, data_dp, grid_in); // Regrid data and log the range of values before and after double dmin, dmax, dmin_rgd, dmax_rgd; data_dp.data_range(dmin, dmax); - data_dp = met_regrid(data_dp, latlon_arr, grid, data_info->regrid()); + data_dp = met_regrid(data_dp, grid_in, grid_out, data_info->regrid()); data_dp.data_range(dmin_rgd, dmax_rgd); mlog << Debug(4) << data_info->magic_str() @@ -747,7 +747,7 @@ void process_fields(const TrackInfoArray& tracks) { // if this is "U", setup everything for matching "V" and compute the radial/tangential if(wind_converter.compute_winds_if_input_is_u(i_point, sname, slevel, valid_time, data_files, ftype, - latlon_arr, lat_arr, lon_arr, grid, data_dp, tcrmw_grid)) { + grid_in, grid_out, data_dp, tcrmw_grid)) { write_tc_pressure_level_data(nc_out, tcrmw_grid, pressure_level_indices, data_info->level_attr(), i_point, data_3d_vars[conf_info.radial_velocity_field_name.string()], diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw.h b/src/tools/tc_utils/tc_rmw/tc_rmw.h index 87e37b74be..85ffef1c64 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw.h +++ b/src/tools/tc_utils/tc_rmw/tc_rmw.h @@ -84,7 +84,6 @@ static TCRMWConfInfo conf_info; static GrdFileType ftype; static TCRMW_WindConverter wind_converter; - // Optional arguments static ConcatString out_dir; static ConcatString out_prefix; @@ -136,20 +135,15 @@ static std::map pressure_level_indices; // //////////////////////////////////////////////////////////////////////// -static DataPlane dp; -static Grid latlon_arr; -static TcrmwData grid_data; +static Grid grid_in; +static TcrmwData tcrmw_data; static TcrmwGrid tcrmw_grid; -static Grid grid; +static Grid grid_out; // Grid coordinate arrays static double* lat_arr; static double* lon_arr; -// Wind arrays -/* static double* wind_r_arr; */ -/* static double* wind_t_arr; */ - //////////////////////////////////////////////////////////////////////// #endif // __TC_RMW_H__ From 34c9beac9846b7b8e8b98ff2765f58ec10d10d19 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Mon, 6 May 2024 19:45:11 +0000 Subject: [PATCH 08/15] Per #2841, the substantive bug fixed here is how we index into the U and V data in wind_ne_to_rt() using i_rev instead of i. Also, update the variable names for consistency and clarity.. --- .../tc_utils/tc_rmw/tc_rmw_wind_converter.cc | 63 ++++++++++--------- .../tc_utils/tc_rmw/tc_rmw_wind_converter.h | 8 +-- 2 files changed, 35 insertions(+), 36 deletions(-) diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc index 7aec340e56..f6afa06079 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc @@ -31,7 +31,6 @@ using namespace std; static void wind_ne_to_rt(const TcrmwGrid&, const DataPlane&, const DataPlane&, - const double*, const double*, double*, double*); //////////////////////////////////////////////////////////////////////// @@ -169,11 +168,9 @@ bool TCRMW_WindConverter::compute_winds_if_input_is_u(int i_point, unixtime valid_time, const StringArray &data_files, const GrdFileType &ftype, - const Grid &latlon_arr, - const double *lat_arr, - const double *lon_arr, - const Grid &grid, - const DataPlane &data_dp, + const Grid &grid_in, + const Grid &grid_out, + const DataPlane &u_wind_dp, const TcrmwGrid &tcrmw_grid) { if (!_computeWinds) { return false; @@ -181,34 +178,33 @@ bool TCRMW_WindConverter::compute_winds_if_input_is_u(int i_point, int uIndex = -1; int vIndex = -1; - VarInfo *data_infoV = (VarInfo *) 0; + VarInfo *v_wind_info = (VarInfo *) 0; if (varName == _conf->u_wind_field_name.string()) { uIndex = _uIndexMap[varLevel]; vIndex = _vIndexMap[varLevel]; - data_infoV = _conf->data_info[vIndex]; - data_infoV->set_valid(valid_time); + v_wind_info = _conf->data_info[vIndex]; + v_wind_info->set_valid(valid_time); } else { // not the U input return false; } - DataPlane data_dpV; - Grid latlon_arrV; - get_series_entry(i_point, data_infoV, data_files, ftype, - data_dpV, latlon_arrV); + DataPlane v_wind_dp; + Grid v_wind_grid; + get_series_entry(i_point, v_wind_info, data_files, ftype, + v_wind_dp, v_wind_grid); double dmin, dmax, dmin_rgd, dmax_rgd; - data_dpV.data_range(dmin, dmax); - data_dpV = met_regrid(data_dpV, latlon_arr, grid, data_infoV->regrid()); - data_dpV.data_range(dmin_rgd, dmax_rgd); + v_wind_dp.data_range(dmin, dmax); + v_wind_dp = met_regrid(v_wind_dp, v_wind_grid, grid_out, v_wind_info->regrid()); + v_wind_dp.data_range(dmin_rgd, dmax_rgd); - mlog << Debug(4) << data_infoV->magic_str() + mlog << Debug(4) << v_wind_info->magic_str() << " input range (" << dmin << ", " << dmax << "), regrid range (" << dmin_rgd << ", " << dmax_rgd << ")\n"; // Compute the radial and tangential winds and store in _windR and _windT - wind_ne_to_rt(tcrmw_grid, data_dp, data_dpV, lat_arr, lon_arr, - _windR, _windT); + wind_ne_to_rt(tcrmw_grid, u_wind_dp, v_wind_dp, _windR, _windT); return true; } @@ -217,31 +213,36 @@ bool TCRMW_WindConverter::compute_winds_if_input_is_u(int i_point, 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) { + int n_rng = tcrmw_grid.range_n(); + int n_azi = tcrmw_grid.azimuth_n(); + // 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++) { + for(int ir = 0; ir < n_rng; ir++) { + for(int ia = 0; ia < n_azi; ia++) { - int i = ir * tcrmw_grid.azimuth_n() + ia; + // Store data in reverse order + int i_rev = (n_rng - ir - 1) * n_azi + ia; 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]); + double lat, lon; + tcrmw_grid.range_azi_to_latlon(range_km, azi_deg, lat, lon); + + tcrmw_grid.wind_ne_to_rt(azi_deg, u_dp.data()[i_rev], v_dp.data()[i_rev], + wind_r_arr[i_rev], wind_t_arr[i_rev]); mlog << Debug(4) << "wind_ne_to_rt() -> " - << "center (" << tcrmw_grid.lat_center_deg() + << "center lat/lon (" << 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"; - + << ", point lat/lon (" << lat << ", " << lon + << "), uv (" << u_dp.data()[i_rev] << ", " << v_dp.data()[i_rev] + << "), radial wind: " << wind_r_arr[i_rev] + << ", tangential wind: " << wind_t_arr[i_rev] << "\n"; } // end for ia } // end for ir diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h index f87d86193a..d80b8ed15d 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.h @@ -113,11 +113,9 @@ class TCRMW_WindConverter { unixtime valid_time, const StringArray &data_files, const GrdFileType &ftype, - const Grid &latlon_arr, - const double *lat_arr, - const double *lon_arr, - const Grid &grid, - const DataPlane &data_dp, + const Grid &grid_in, + const Grid &grid_out, + const DataPlane &u_wind_dp, const TcrmwGrid &tcrmw_grid); }; From 99d7d3d2c813fa18bb5a938397261a5641bf2595 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Mon, 13 May 2024 16:19:12 -0600 Subject: [PATCH 09/15] Per #2841, patch apparent copy/paste bug in tcrmw_grid.cc code. --- src/libcode/vx_grid/tcrmw_grid.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/libcode/vx_grid/tcrmw_grid.cc b/src/libcode/vx_grid/tcrmw_grid.cc index 6e9b105729..14be0516ea 100644 --- a/src/libcode/vx_grid/tcrmw_grid.cc +++ b/src/libcode/vx_grid/tcrmw_grid.cc @@ -378,7 +378,7 @@ RotatedLatLonGrid::latlon_to_xy(true_lat, true_lon, x, y); x -= Nx*floor(x/Nx); -x -= Nx*floor(x/Nx); +y -= Ny*floor(y/Ny); return; @@ -394,7 +394,7 @@ void TcrmwGrid::xy_to_latlon(double x, double y, double & true_lat, double & tru x -= Nx*floor(x/Nx); -x -= Nx*floor(x/Nx); +y -= Ny*floor(y/Ny); RotatedLatLonGrid::xy_to_latlon(x, y, true_lat, true_lon); From afc4e3053ef2226fd614eccbf83bfdec03f3adf8 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 20 Jun 2024 17:27:09 +0000 Subject: [PATCH 10/15] Per #2841, only changing whitespace. --- src/libcode/vx_grid/tcrmw_grid.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libcode/vx_grid/tcrmw_grid.cc b/src/libcode/vx_grid/tcrmw_grid.cc index 14be0516ea..9c47fef682 100644 --- a/src/libcode/vx_grid/tcrmw_grid.cc +++ b/src/libcode/vx_grid/tcrmw_grid.cc @@ -136,7 +136,7 @@ Ir.set_xyz(1.0, 0.0, 0.0); Jr.set_xyz(0.0, 1.0, 0.0); Kr.set_xyz(0.0, 0.0, 1.0); -Range_n = 0; +Range_n = 0; Azimuth_n = 0; Range_max_km = 0.0; From 33a00f9c0b4c68fc2ddb0460f089555022a8f9c0 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 20 Jun 2024 18:26:35 +0000 Subject: [PATCH 11/15] Per #2841, update logic in EarthRotation::set_tcrmw() to change the rotation of TCRMW grids from pointing north to pointing east. --- src/libcode/vx_grid/earth_rotation.cc | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/libcode/vx_grid/earth_rotation.cc b/src/libcode/vx_grid/earth_rotation.cc index 35049e0076..ec8b855379 100644 --- a/src/libcode/vx_grid/earth_rotation.cc +++ b/src/libcode/vx_grid/earth_rotation.cc @@ -179,7 +179,13 @@ M23 = clat*clon; M33 = slat; */ -set_np(lat_center, lon_center, lon_center - 180.0); + // + // MET #2841 define the rotation by subtracting 90 degrees + // instead of 180 to define TCRMW grids as pointing east + // instead of north. + // + +set_np(lat_center, lon_center, lon_center - 90.0); // // From 090028133ac5446162109a9a495c5e0d3390ca45 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Thu, 20 Jun 2024 23:24:40 +0000 Subject: [PATCH 12/15] Per #2841, fix warning message --- src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc index f6afa06079..76b54037c7 100644 --- a/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc +++ b/src/tools/tc_utils/tc_rmw/tc_rmw_wind_converter.cc @@ -147,14 +147,14 @@ void TCRMW_WindConverter::append_nc_output_vars(map > &va } else { if (!_foundUInInput) { - mlog << Warning << "\nTCWRMW_WindConverter::checkInputs() -> " + mlog << Warning << "\nTCWRMW_WindConverter::append_nc_output_vars() -> " << "field not found in input \"" << _conf->u_wind_field_name << "\"\n\n"; } if (!_foundVInInput) { - mlog << Warning << "\nTCWRMW_WindConverter::checkInputs() -> " + mlog << Warning << "\nTCWRMW_WindConverter::append_nc_output_vars() -> " << "field not found in input \"" << _conf->v_wind_field_name << "\"\n\n"; } - mlog << Warning << "\nTCWRMW_WindConverter::checkInputs() -> " + mlog << Warning << "\nTCWRMW_WindConverter::append_nc_output_vars() -> " << "Not computing radial and tangential winds\n\n"; _computeWinds = false; } From e5bb35af867cb482d59c8e726ba8d68fd32260d9 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Fri, 21 Jun 2024 21:52:45 +0000 Subject: [PATCH 13/15] Replace tab with spaces --- src/libcode/vx_tc_util/vx_tc_nc_util.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libcode/vx_tc_util/vx_tc_nc_util.cc b/src/libcode/vx_tc_util/vx_tc_nc_util.cc index f2f545cdef..5a13616c6f 100644 --- a/src/libcode/vx_tc_util/vx_tc_nc_util.cc +++ b/src/libcode/vx_tc_util/vx_tc_nc_util.cc @@ -165,7 +165,7 @@ bool has_pressure_level(vector levels) { for (int j = 0; j < levels.size(); j++) { if (levels[j].substr(0, 1) == "P") { status = true; - break; + break; } } From 9e251a5a637719604084f9d4f7eb22661c06f82e Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Sat, 22 Jun 2024 02:38:57 +0000 Subject: [PATCH 14/15] Per #2841, correct the units for the azimuth netcdf output variable --- src/libcode/vx_tc_util/vx_tc_nc_util.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libcode/vx_tc_util/vx_tc_nc_util.cc b/src/libcode/vx_tc_util/vx_tc_nc_util.cc index 5a13616c6f..9278fd5118 100644 --- a/src/libcode/vx_tc_util/vx_tc_nc_util.cc +++ b/src/libcode/vx_tc_util/vx_tc_nc_util.cc @@ -338,7 +338,7 @@ void def_tc_range_azimuth(NcFile* nc_out, add_att(&range_var, "_FillValue", bad_data_double); add_att(&azimuth_var, "long_name", "azimuth"); - add_att(&azimuth_var, "units", "degrees_clockwise_from_north"); + add_att(&azimuth_var, "units", "degrees_clockwise_from_east"); add_att(&azimuth_var, "standard_name", "azimuth"); add_att(&azimuth_var, "_FillValue", bad_data_double); From 2bde096fefc05344bda864108b87b942ba7343ee Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Sat, 22 Jun 2024 11:16:27 +0000 Subject: [PATCH 15/15] Per #2841, reverse the x dimension of the rotated latlon grid to effectively switch from counterclockwise rotation to clockwise. --- docs/Users_Guide/tc-diag.rst | 2 +- docs/Users_Guide/tc-rmw.rst | 2 +- src/libcode/vx_grid/tcrmw_grid.cc | 7 +++++++ 3 files changed, 9 insertions(+), 2 deletions(-) diff --git a/docs/Users_Guide/tc-diag.rst b/docs/Users_Guide/tc-diag.rst index 53c1ee1b0d..1ef12703d2 100644 --- a/docs/Users_Guide/tc-diag.rst +++ b/docs/Users_Guide/tc-diag.rst @@ -24,7 +24,7 @@ Originally developed for the Statistical Hurricane Intensity Prediction Scheme ( TC-Diag is run once for each initialization time to produce diagnostics for each user-specified combination of TC tracks and model fields. The user provides track data (such as one or more ATCF a-deck track files), along with track filtering criteria as needed, to select one or more tracks to be processed. The user also provides gridded model data from which diagnostics should be computed. Gridded data can be provided for multiple concurrent storms, multiple models, and/or multiple domains (i.e. parent and nest) in a single run. -TC-Diag first determines the list of valid times that appear in any one of the tracks. For each valid time, it processes all track points for that time. For each track point, it reads the gridded model fields requested in the configuration file and transforms the gridded data to a range-azimuth cylindrical coordinates grid. For each domain, it writes the range-azimuth data to a temporary NetCDF file. +TC-Diag first determines the list of valid times that appear in any one of the tracks. For each valid time, it processes all track points for that time. For each track point, it reads the gridded model fields requested in the configuration file and transforms the gridded data to a range-azimuth cylindrical coordinates grid, as described for the TC-RMW tool in :numref:`tc-rmw`. For each domain, it writes the range-azimuth data to a temporary NetCDF file. .. note:: The current version of the tool does not yet include the capabilities described in the next three paragraphs. These additional capabilities are planned to be added in the MET v12.0.0 release later in 2023. diff --git a/docs/Users_Guide/tc-rmw.rst b/docs/Users_Guide/tc-rmw.rst index d37ea66e42..9c67d56312 100644 --- a/docs/Users_Guide/tc-rmw.rst +++ b/docs/Users_Guide/tc-rmw.rst @@ -7,7 +7,7 @@ TC-RMW Tool Introduction ============ -The TC-RMW tool regrids tropical cyclone model data onto a moving range-azimuth grid centered on points along the storm track provided in ATCF format, most likely the adeck generated from the file. The radial grid spacing may be set as a factor of the radius of maximum winds (RMW). If wind fields are specified in the configuration file the radial and tangential wind components will be computed. Any regridding method available in MET can be used to interpolate data on the model output grid to the specified range-azimuth grid. The regridding will be done separately on each vertical level. The model data files must coincide with track points in a user provided ATCF formatted track file. +The TC-RMW tool regrids tropical cyclone model data onto a moving range-azimuth grid centered on points along the storm track provided in ATCF format, most likely the adeck generated from the file. The radial grid spacing can be defined in kilometers or as a factor of the radius of maximum winds (RMW). The azimuthal grid spacing is defined in degrees clockwise from due east. If wind vector fields are specified in the configuration file the radial and tangential wind components will be computed. Any regridding method available in MET can be used to interpolate data on the model output grid to the specified range-azimuth grid. The regridding will be done separately on each vertical level. The model data files must coincide with track points in a user provided ATCF formatted track file. Practical information ===================== diff --git a/src/libcode/vx_grid/tcrmw_grid.cc b/src/libcode/vx_grid/tcrmw_grid.cc index 9c47fef682..d0c32d4e51 100644 --- a/src/libcode/vx_grid/tcrmw_grid.cc +++ b/src/libcode/vx_grid/tcrmw_grid.cc @@ -288,6 +288,7 @@ y = (lat_rot - RData.rot_lat_ll)/(RData.delta_rot_lat); x = lon_rot/(RData.delta_rot_lon); +x = Nx - x; // MET #2841 switch from counterclockwise to clockwise RotatedLatLonGrid::xy_to_latlon(x, y, lat, lon); @@ -310,6 +311,8 @@ const double range_max_deg = deg_per_km*Range_max_km; RotatedLatLonGrid::latlon_to_xy(lat, lon, x, y); +x = Nx - x; // MET #2841 switch from counterclockwise to clockwise + azi_deg = x*(RData.delta_rot_lon); range_deg = range_max_deg - y*(RData.delta_rot_lat); @@ -378,6 +381,8 @@ RotatedLatLonGrid::latlon_to_xy(true_lat, true_lon, x, y); x -= Nx*floor(x/Nx); +x = Nx - x; // MET #2841 switch from counterclockwise to clockwise + y -= Ny*floor(y/Ny); return; @@ -394,6 +399,8 @@ void TcrmwGrid::xy_to_latlon(double x, double y, double & true_lat, double & tru x -= Nx*floor(x/Nx); +x = Nx - x; // MET #2841 switch from counterclockwise to clockwise + y -= Ny*floor(y/Ny); RotatedLatLonGrid::xy_to_latlon(x, y, true_lat, true_lon);