NOAA-GFDL / FRE-NCtools

Tools for manipulating and creating netCDF inputs for FMS managed models
GNU Lesser General Public License v3.0
19 stars 28 forks source link

bronx-19 fregrid is broken for stretched grid sources #180

Closed nikizadehgfdl closed 2 years ago

nikizadehgfdl commented 2 years ago

Describe the bug fregrid in bronx-19 produces huge numbers all over the map when the source data is on stretched grid CS.

To Reproduce

cd /work/Niki.Zadeh/mosaic_generation/exchange_grid_toolset/test_stretched_grid_lwh/c256r25L49_tlat32N_am4p0_AM4p1chem_camsEm_dustddep_nsp5f1800_dp0.5
fregrid --standard_dimension --input_mosaic C256_mosaic.nc --input_file 20080101.atmos_month --interp_method conserve_order2 --nlon 576 --nlat 360 --scalar_field olr,t_ref --output_file out.bronx-19.nc

Expected behavior No huge numbers in output file. bronx-18 was doing it right.

System Environment PAN with module load fre/bronx-19

Additional context The huge numbers stem from an update to poly_area here: https://github.com/NOAA-GFDL/FRE-NCtools/blob/master/tools/libfrencutils/mosaic_util.c#L412

Although this same update fixed an issue with generating the source stretched grid at the South Pole (issue #45 ). I think more work needs to be done on poly_area to understand the algorithm for a general fix. At the least I suggest bypassing the fix in PR #60 for stretched grid generation by introducing a switch:

-    if(dx <= -M_PI) dx = dx + 2.0*M_PI;
+    if(dx < -M_PI) dx = dx + 2.0*M_PI;
+    if(sign_bug && dx == -M_PI) dx = dx + 2.0*M_PI;

I found out that the rest of the updates to poly_area in PR #60 were only cosmetic and not actually needed for stretched grid.

nikizadehgfdl commented 2 years ago

@ceblanton or @ngs333 can you please assign this issue to me since I am working on it.

nikizadehgfdl commented 2 years ago

I found the root cause. The fix I made for poly_area in issue #45 produces the corrected value for area in grid cells that straddle the South Pole that would otherwise be 25% of the area of the Earth. But, the function poly_ctrlat() in create_xgrid.c which is called by frergid --conserve_order2 has a similar bug for those grid cells and comes up with a vastly wrong latitude for the centroid of those grid cells. However since ctrlat is devided by cell_area in fregrid calculations the two nonsense huge values (for incorrect area and ctrlat) used to cancel each other out and fregrid would end up with a finite value. After my fix to poly_area the cancellation would not happen and we end up with very wrong (huge negative) values for otherwise positive quantities (e.g., olr) for the few grid cells around the South Pole. The fix is simply to do the same fix I did for poly_area() in poly_ctrlat().

In my opinion the fregrid estimates are not accurate at the South Pole cells for stretched grids, with or without any of my fixes. It's only that the values make sense after I put the fix in both poly_area and ctrlat, or without either of the fixes (as in bronx-18).