The field of computational fluid dynamics (CFD) is dedicated to solving the Navier-Stokes equations (or some simplification of them). A subset of CFD, ocean and atmospheric models numerically solve the same equations for realistic applications. What are the differences and trade-offs between the general CFD approaches and the applied realistic cases?
Atmosphere and ocean have highly-stratified flows in which the Coriolis force is a major source of dynamics. Maintaining geostrophic balance is extremely important and many numerical schemes are intended to be exactly compatible (at least in the absence of topography) to avoid radiating energy in gravity waves. Due to the stratification, limiting vertical numerical diffusion is extremely important and special grids are often used (especially in the ocean) for that purpose. Many methods are effectively 2.5-dimensional formulations.
For climate simulation over long time periods, conservation of energy and other fluxes (like salt) are often considered to be critical for statistically meaningful results. Methods that are less accurate and have certain numerical artifacts may be chosen in order to avoid damping out the dynamics. Note that long-term dynamics may not homogenize at continental scales averaged over multiple decades.
Industrial CFD solvers tend to be used for flows that are more isotropic (genuinely 3D) and often neglect Coriolis. They often have stronger forcing and thus less critical energy conservation requirements. It is common to deal with strong shocks, in which case nonlinear spatial discretizations must be used, despite being more dissipative.
Because lab experiments can actually be performed for most industrial applications, the software experiences more validation. Weather models also have constant validation, but climate models are almost impossible to validate due to the time scales involved and unavoidable over-fitting.
Jed Brown described the traditional approach used in the mesoscale and larger scale models. Actually, in microscale many atmospheric models are very close to to traditional CFD codes, use similar finite volume discretizations, similar 3D grids where vertical is treated similarly as horizontal, and so on. Depending on the resolutions even features like buildings are resolved with the same approaches known from engineering CFD, like the immersed boundary methods or body fitted grids.
You can encounter all the discretization techniques you know from the engineering CFD, like finite differences, finite volumes, pseudo-spectral and even finite elements. The same pressure correction (fractional-step) methods are often used to solve the incompressible Navier-Stokes equations (with the Boussinesq or anelastic terms for buoyancy).
Of course, different parametrization for the heat and momentum fluxes near the surface are commonly used, taking into account the specifics of the land-surface interactions like the Monin-Obukhov similarity or other semi-empirical relations.
The whole method of large-eddy simulation (LES), now very popular in engineering, actually originates in boundary layer meteorology. I would even say that many atmospheric modelers at this scale wouldn't hesitate at all to call their work CFD.
In many (but not all) applications you also have to add the Coriolis force. The schemes do not have to be well-balanced however, it is just one additional volume force. If you also compute the processes like cloud formation, precipitation and radiation, things get more complicated, but the same holds for engineering models which solve reaction kinetics, combustion and similar.
This class of models also includes those accounting for the ocean-atmosphere interactions you asked for, see for example https://ams.confex.com/ams/pdfpapers/172658.pdf
Difference between weather prediction software and "casual CFD solver" is how weather prediction works with transition of water. Water is being treated as second component, so model becomes 3 dimensional with 2 components.
Second (not the least important) is that weather (and other convection) models work directly with $\omega$ from $d \omega/dt = (\omega \nabla) u + \nu \nabla^2 \omega$.