Hey all. Not a reviewer, but I did just refactor these impvmix routines for GPUs (see the ocean/gpuimpvmix branch in my philipwjones fork) and have a few comments. For both performance reasons and eliminating the need for all these separate routines, it's actually better to compute the bottom drag term as an (nEdges) array in a separate loop before computing the ABC coefficients.

Second, as a general rule, don't use real numbers as exponents if they are in fact integers. The algorithm for computing real exponentiation is more expensive than the integer form so you want the compiler to choose the right form. This is especially true if the exponent is 2 or 3 where the compiler will often replace with a simple multiply (or even better, for a simple square, do the replacement yourself and avoid the exponentiation).

For fixed constant real numbers, be sure to add the _RKIND suffix (eg 0.5_RKIND) to ensure correct precision - you may only get single precision otherwise. Finally, try to combine/precompute some of these factors (especially divisions). Although flops are cheap, they are not free and precomputing common factors almost always wins.

Rebased to head of ocean/develop. Passes nightly test suite.

The PR has been verified against previous studies (Pringle et al 2021) and also with TPXO8 data.
Attached are plots to show such comparisons.

The following plot can be compared with Pringle et al 2021 (Figs 5 and 6)

Below is the plot showing code agreement before and after rebasing

The changes here were absorbed into #607, which is also on drag, so they could be tested together.