Here a some partial solutions I got if ever anyone else is facing the same issue
Method 1 :
Compute interplevel at a given altitude,
Then compute 1st order derivative on the resulting 2D array
Repeat for every desired altitude
Method 2 :
Compute vertcross with z_levels being all the desired altitude, resulting in a (x, z) map.
Repeat for each Y value and concatenate in a 3D array
Then compute derivatives on each horizontal plane
These 2 methods are quite similar and should work but there are some drawbacks :
- They do not take into account the projection bias
- The resulting 3D array of gradient is computed on a new orthogonal mesh instead of on the existing hybrid vertical coordinate mesh
Here are some articles with potential methods :
Shchepetkin, Alexander, et J.C. McWilliams. « A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate ».
Journal of Geophysical Research C: Oceans 108 (15 mars 2003): 35‑1.
Lundquist, Katherine A. « Truly Horizontal Finite Difference Scheme in the WRF Model FY18 Q1 Update ». Lawrence Livermore National Lab. (LLNL), Livermore, CA (United States), 27 septembre 2018.
Truly horizontal finite difference scheme in the WRF model FY18 Q1 Update (Technical Report) | OSTI.GOV.