Skip to content

Doing bicubic interpolation on data with inverted latitudes leads to incorrect interpolation #1889

Description

@uwagura

Describe the bug
Performing bicubic interpolation on data where the latitude values are listed in descending order rather than ascending order leads to an incorrect interpolation. I believe this is because of the following lines in horiz_interp/include/horiz_interp_bicubic.inc:

<----- NOTE - following lines largely written by AI ---->

if( yz .le. lat_in(1) ) then
   jcl = 1
   jcu = 1
else if( yz .ge. lat_in(nlat_in) ) then
   ...

This logic assumes lat_in(1) is the minimum (southernmost) latitude. When the array is descending, lat_in(1) is instead the maximum (northernmost) latitude, so the condition yz .le. lat_in(1) is TRUE for every destination point. As a result:

  • jcl = jcu = 1 for all points.
  • rat_y = 0.0 for all points.
  • All interpolated values equal the field value at row 1 of the source grid (the northernmost latitude).

The helper functions INDL_ and INDU_ used in the else branch also assume an ascending array, but they are never reached because the first if branch catches everything.
<---- END AI assisted section ---->

CEFI configurations use bicubic interpolation for the winds. When I was preparing new forcing files using cdo, the latitude values in my winds got flipped , and as a result of this bug, I pretty much lost all of my wind stress in my experiment:

Image

Left shows u* in my experiment with the ascending latitude values ( i.e the expected value) and right shows the same with descending latitude values.

Note that billinearly interpolated variables don't seem to have this issue, as there seems to be logic there to handle descending coordinates in horiz_interp/include/horiz_interp_bilinear.inc:

  44 
  45     decreasing_lat = .false.
  46     if (lat_in(1) > lat_in(2)) decreasing_lat = .true.
  47 

To Reproduce
Run any experiment with bicubic interpolation for a given forcing file that has descending latitude values

Expected behavior
Consistent results regardless of the order of the latitude coordinate in a forcing file

System Environment
This was on Gaea c6, so

  • OS: SUSE Linux Enterprise Server 15.6
  • Compiler(s): Type and version ifort 2021.10.0
  • MPI type, and version mpich 8.1.21
  • netCDF Version: netCDF 4.9.0 , netCDF-Fortran 4.5.3

Additional context
As noted above, I was heavily assisted by Claude Sonnet 4.6 via Github Copilot to discover and describe this issue. I also asked for Claude's assistance in developing a fix, which I am currently testing and can open as PR if developers agree that this is and issue and the fix works.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions