| 22 Jun 2026
A High-Fidelity CUDA Implementation of the Shchepetkin Density-Jacobian Pressure Gradient Scheme for ROMS
Abstract. Baroclinic pressure gradients in terrain-following ocean models require careful numerical treatment to avoid generating spurious currents over steep topography. The Shchepetkin density Jacobian algorithm achieves high accuracy through harmonic mean density slope reconstruction and fourth-order monotonized cubic polynomial corrections, but has largely remained on CPU architectures despite growing GPU adoption in Earth system models. This work presents the first publicly available CUDA implementation of the complete Shchepetkin density-Jacobian pressure gradient scheme as used in the Regional Ocean Modeling System (ROMS). No algorithmic simplifications were introduced during GPU porting. Validation against CPU reference solutions on the "Tall Isolated Seamount" benchmark (Beckmann and Haidvogel, 1993) with steep topography (r ≈ 0.4) demonstrates maximum relative error of πͺ (10−6) across the 54 × 51 × 13 computational domain. A "lake at rest" test using a horizontally uniform density field over the same Seamount bathymetry confirms the well-balanced property: the kernel produces exactly zero pressure gradient force to machine precision. The small discrepancies in the Seamount comparison reflect floating-point precision and architectural differences between GPU and CPU rather than algorithmic deficiencies. The standalone kernel architecture enables integration into existing ocean models. Source code and validation data are made publicly available under an open-source license.
This is a strong piece of engineering work showing a faithful CUDA implementation of a numerically delicate ROMS pressure-gradient scheme. The thread-per-column strategy is well suited to the vertical recurrence, and the standalone binary-extraction framework provides a practical way to verify the port through point-by-point comparison. I also appreciated the discussion in Sect. 3.3 of the indexing error identified during development. That example clearly shows why this type of quantitative verification is important.
My main concern is not with the implementation itself, but with how some of the results are framed. At several points, the manuscript moves from showing that the GPU reproduces the CPU calculation to broader claims about validation, well-balanced behavior, and practical performance. Those claims either need additional evidence or should be stated more cautiously. The following issues should be addressed before publication.
1. Verification and validation should be distinguished more clearly
The study shows that, for identical inputs, the CUDA implementation closely reproduces the existing Fortran routine. It shows that the translation from the CPU implementation to CUDA has been performed correctly and that no obvious numerical or indexing errors remain in the tested configuration. However, this comparison does not independently validate the pressure-gradient scheme itself. Any error already present in the CPU implementation would also be reproduced by the GPU version and would not be revealed by the comparison. I would suggest revising the wording (validation -> verification) in the abstract, Sect. 2.3, Sect. 3, and the conclusions so that the main result is framed as numerical fidelity to the existing ROMS implementation rather than independent confirmation of the schemeβs correctness.
2. The seamount experiment is a single-call comparison rather than a dynamic benchmark
The introduction emphasizes that pressure-gradient errors can generate spurious currents that persist or grow over time. This is also the main reason the seamount problem is widely used as a benchmark. In the present study, however, the comparison appears to be based on one time step. It does not show how the GPU implementation behaves during a time integration. A stronger experiment would integrate the CUDA kernel within ROMS over many baroclinic time steps and compare diagnostics such as maximum spurious velocity and domain-integrated kinetic energy. If full integration is beyond the scope of the current study, the manuscript should state more explicitly that the present experiment is a single-call verification test and should avoid broader claims about long-term stability or circulation behavior.
3. The lake-at-rest test should be described more narrowly
The current interpretation of the lake-at-rest test appears stronger than the test supports. Because the density anomaly is set to zero everywhere, the density-dependent baroclinic terms vanish directly. The result therefore confirms that the kernel returns zero for a zero-density-anomaly input, but it does not test cancellation among nonzero density, pressure, and geometric terms. I would suggest narrowing the associated well-balanced claim. A more demanding alternative would be to add a case in which the expected tendency is independently known and nonzero intermediate terms must cancel numerically.
4. The performance comparison should be placed in a more representative context
The reported 15Γ speedup is based on comparison with one Intel Xeon Gold 6154 core on a 54x51x13 domain and is not representative of how ROMS is typically used.Β For practical applications, the relevant CPU runs would involve multiple cores and MPI domain decomposition rather than a single core. The current result should be presented as a preliminary kernel-level comparison for the specific grid and hardware tested. The manuscript should clearly distinguish single-core speedup from node-level or resource-equivalent performance. It would also be useful to report the CPU compiler, optimization flags, and vectorization settings, since the baseline walltime may depend strongly on the build configuration.
The timing methodology also needs to state whether the reported walltime includes memory allocation, host-to-device and device-to-host transfers, kernel launches, synchronization, warm-up iterations, and CUDA context initialization, since the manuscript itself notes that launch and transfer latency may be significant for the small seamount grid.
5. The MAX_N = 64 limitation needs clarification
The statement that βMAX_N = 64 accommodates the vertical resolution of any current ROMS configurationβ appears too broad. Please clarify whether MAX_N is a compile-time parameter that users can modify, what happens when N>64, and whether the code fails with a clear error.Β
6. Figure 2a is described as showing absolute error, but the color bar contains both positive and negative values. Please revise either the plotted quantity or the caption and color-bar label.
7. The βfirst publicly availableβ claim is stated definitively in the abstract and conclusions but qualified as βto our knowledgeβ in Sect. 1.2. Please use consistent wording throughout and briefly indicate the basis for the claim.