The description in the Fortran code says:
This subroutine performs interpolation of a bivariate function,
z(x,y), on a rectangular grid in the x-y plane. It is based on
the revised Akima method.
In this subroutine, the interpolating function is a piecewise
function composed of a set of bicubic (bivariate third-degree)
polynomials, each applicable to a rectangle of the input grid
in the x-y plane. Each polynomial is determined locally.
This subroutine has the accuracy of a bicubic polynomial, i.e.,
it interpolates accurately when all data points lie on a
surface of a bicubic polynomial.
The grid lines can be unevenly spaced.