I am trying to compute a decomposition of the pressure field in MOM6. My goal is to compute the baroclinic pressure at the z=0 level, the "mean" ocean surface, and then convert this to an equivalent sea surface height anomaly by dividing by g*rho0. By "baroclinic pressure", I am referring to the hydrostatic pressure anomaly minus its vertical average.

My current setup is using isopycnal layers (based on MOM6-examples/ocean_only/global_ALE/layer/ example) with the Boussinesq approximation, using 50 layers and 1/12-deg horizontal resolution. I'm studying the pressure gradient code in MOM_PressureForce_Montgomery.F90, and it looks like the pressure gradient force is being computed from the Montgomery function with, grad M + z grad \rho, (PFu_bc and PFu around lines 557-561, MOM6/src/core/MOM_PressureForce_Montgomery.F90 at main · mom-ocean/MOM6).

Is the z grad \rho term non-zero when using layer coordinates in the vertical? I ask, because it looks like the non-EOS code option neglects this term.

If so, then I wonder if it is possible to initialize the model using realistic (WOA) temperature and salinity, but using a non-EOS-based setting? I'm not sure of the right question to ask. I'd like to run the model in non-EOS isopycnal mode (so that grad rho is zero in each layer), but I'd like the layers to be configured to represent the actual ocean. I only want to model the quasi-linear evolution of waves for 30-days or so; the whole ocean does not need to be in realistic thermodynamic equilibrium.

-Ed