Bugfix #2841 develop tang_rad_winds (#2921)

* Per #2841, port over fixes from bugfix_2841_main_v11.1_tang_rad_winds for the develop branch

* Per #2841, clarify in the docs that azimuths are defined in degrees counter-clockwise from due east.

* Per #2841, just updating with output from enum_to_string.

* Per #2841, tweak the documentation.

* Per #2841, correct the location of using namespace lines.

* Per #2841, update compute_tc_diag.py to no longer skip writing the radial and tangential wind diagnostics.

* Per #2841, update compute_tc_diag.py to no longer skip writing radial and tangential wind diagnostics.

* Revert "Per #2841, update compute_tc_diag.py to no longer skip writing radial and tangential wind diagnostics."

This reverts commit f097345.

* Revert "Per #2841, update compute_tc_diag.py to no longer skip writing the radial and tangential wind diagnostics."

This reverts commit c040215.

* Per #2841, update comp_dir.sh logic to include .dat in the files that are diffed

* Replace tab with spaces

* Per #2841, correct the units for the azimuth netcdf output variable

* Per #2841, reverse the x dimension of the rotated latlon grid to effectively switch from counterclockwise rotation to clockwise.

---------

Co-authored-by: MET Tools Test Account <[email protected]>

docs/Users_Guide/tc-diag.rst

@@ -15,7 +15,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, as described in :numref:`Contributor's Guide Section %s <tmp_files_tc_diag>`.
+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, as described in :numref:`Contributor's Guide Section %s <tmp_files_tc_diag>`.
 
 Once the input data have been processed into the temporary NetCDF files, TC-Diag then calls one or more Python diagnostics scripts, as specified in the configuration file, to compute tropical cyclone diagnostic values. The computed diagnostics values are retrieved from the Python script and stored in memory.
 

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
 =====================

internal/test_unit/R_test/comp_dir.R

@@ -51,9 +51,9 @@ strDir1 = gsub("/$", "", listArgs[1]);
 strDir2 = gsub("/$", "", listArgs[2]);
 
 # build a list of files in each stat folder
-listTest1 = system(paste("find", strDir1, "| egrep '\\.stat$|\\.txt$|\\.tcst|\\.nc$|\\.out$|\\.ps$|\\.png$' | sort"), intern=T);
+listTest1 = system(paste("find", strDir1, "| egrep '\\.stat$|\\.txt$|\\.tcst|\\.nc$|\\.out$|\\.ps$|\\.png$|\\.dat$' | sort"), intern=T);
 listTest1Files = gsub(paste(strDir1, "/", sep=""), "", listTest1);
-listTest2 = system(paste("find", strDir2, "| egrep '\\.stat$|\\.txt$|\\.tcst|\\.nc$|\\.out$|\\.ps$|\\.png$' | sort"), intern=T);
+listTest2 = system(paste("find", strDir2, "| egrep '\\.stat$|\\.txt$|\\.tcst|\\.nc$|\\.out$|\\.ps$|\\.png$|\\.dat$' | sort"), intern=T);
 listTest2Files = gsub(paste(strDir2, "/", sep=""), "", listTest2);
 
 if( 1 <= verb ){ cat("dir1:", strDir1, "contains", length(listTest1Files), "files\n");
@@ -103,10 +103,11 @@ for(strFile in listTest1Files[ listTest1Files %in% listTest2Files ]){
 		compareNc(strFile1, strFile2, verb, strict, file_size_delta, compare_nc_var);
 	}
 
-	# if the files are PostScript, PNG, or end in .out, compare accordingly
+	# if the files are PostScript, PNG, or end in .out or .dat, compare accordingly
 	else if( TRUE == grepl("\\.out$", strFile1, perl=T) ||
 				TRUE == grepl("\\.ps$",  strFile1, perl=T) ||
-				TRUE == grepl("\\.png$", strFile1, perl=T) ){
+				TRUE == grepl("\\.png$", strFile1, perl=T) || 
+				TRUE == grepl("\\.dat$", strFile1, perl=T) ){
 		if( 1 <= verb ){ cat("file1: ", strFile1, "\nfile2: ", strFile2, "\n", sep=""); }
 		compareDiff(strFile1, strFile2, verb);
 	}

src/basic/vx_util/main.cc

@@ -1,5 +1,5 @@
 // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
-// ** Copyright UCAR (c) 2022 - 2023
+// ** Copyright UCAR (c) 2022 - 2024
 // ** University Corporation for Atmospheric Research (UCAR)
 // ** National Center for Atmospheric Research (NCAR)
 // ** Research Applications Lab (RAL)

src/libcode/vx_gis/shapetype_to_string.cc

@@ -28,7 +28,6 @@
 
 #include "shapetype_to_string.h"
 
-
 using namespace std;
 
 

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);
 
    //
    //

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;
@@ -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);
@@ -324,54 +327,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);
-
-
-
-
-
-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;
 
@@ -381,33 +353,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;
-
-
-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;
-
+double range_km, azi_deg;
 
-u =  cosd(azi_deg);
-
-v = -sind(azi_deg);
-
-w =  0.0;
-
-
-B_azi = u*Ir + v*Jr + w*Kr;
+latlon_to_range_azi(lat, lon, range_km, azi_deg);
 
+wind_ne_to_rt(azi_deg, u_wind, v_wind, radial_wind, tangential_wind);
 
 return;
 
@@ -425,8 +381,9 @@ RotatedLatLonGrid::latlon_to_xy(true_lat, true_lon, x, y);
 
 x -= Nx*floor(x/Nx);
 
-x -= Nx*floor(x/Nx);
+x = Nx - x; // MET #2841 switch from counterclockwise to clockwise
 
+y -= Ny*floor(y/Ny);
 
 return;
 
@@ -442,7 +399,9 @@ void TcrmwGrid::xy_to_latlon(double x, double y, double & true_lat, double & tru
 
 x -= Nx*floor(x/Nx);
 
-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);
 
@@ -500,7 +459,3 @@ return;
 
 ////////////////////////////////////////////////////////////////////////
 
-
-
-
-

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,23 +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_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;
 
-      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
-         //
-         //      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;
-
 };
 
 

src/libcode/vx_tc_util/vx_tc_nc_util.cc

@@ -8,7 +8,6 @@
 
 ////////////////////////////////////////////////////////////////////////
 
-
 #include<netcdf>
 
 #include "vx_tc_nc_util.h"
@@ -159,6 +158,21 @@ void write_tc_rmw(NcFile* nc_out,
 
 ////////////////////////////////////////////////////////////////////////
 
+bool has_pressure_level(vector<string> 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<string> get_pressure_level_strings(
     map<string, vector<string> > variable_levels) {
@@ -289,7 +303,7 @@ void def_tc_pressure(NcFile* nc_out,
     put_nc_data(&pressure_var, &pressure_data[0]);
 
     // Cleanup
-    if(pressure_data) { delete [] pressure_data; pressure_data = (double *) nullptr; }
+    if(pressure_data) { delete [] pressure_data; pressure_data = (double *) 0; }
 
     return;
 }
@@ -324,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);
 
@@ -342,8 +356,8 @@ void def_tc_range_azimuth(NcFile* nc_out,
     put_nc_data(&azimuth_var, &azimuth_data[0]);
 
     // Cleanup
-    if(range_data)   { delete [] range_data;   range_data   = (double *) nullptr; }
-    if(azimuth_data) { delete [] azimuth_data; azimuth_data = (double *) nullptr; }
+    if(range_data)   { delete [] range_data;   range_data   = (double *) 0; }
+    if(azimuth_data) { delete [] azimuth_data; azimuth_data = (double *) 0; }
 
     return;
 }
@@ -539,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);

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<std::string>);
+
 extern std::set<std::string> get_pressure_level_strings(
     std::map<std::string, std::vector<std::string> >);
 

src/tools/tc_utils/tc_diag/tc_diag.cc

@@ -867,7 +867,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
       }
    }