Baroclinic sea level calculation

[c2xu]

Reference: Zaron, E. D., and Ray, R. D. (2023). Clarifying the distinction between steric and baroclinic sea surface height. Journal of Physical Oceanography, 53, 2591-2596.
DOI: 10.1175/JPO-D-23-0073.1

[c2xu]

1. In line 1048 of MOM.F90, I'm not 100% sure if Time_local is the correct time to use.
2. I think eta calculated in line 973 of MOM_diagnostics.F90 should be the same as ssh calculated in line 1025 of MOM.F90. But then in line 1075 of MOM.F90, ssh is calculated differently. I don't fully understand the algorithm for the latter, but it seems more accurate based on the harmonic analysis and the comparison with TPXO data. I'm not sure if this is the SSH that should be used for calculating the baroclinic sea level.

[awallcraft]

On line 270 of src/diagnostics/MOM_harmonic_analysis.F90 "sea serface height." should be "sea surface height.".

[awallcraft]

Baroclinic sea level is a useful new diagnostic, but the current code has unacceptable limitations:

a) This field is valid, and potentially useful, with or without tides, but it is currently only calculated if harmonic analysis is active which in turn requires tides and RK2 splitting. It needs to be reconfigured so that registering baroclinic_sea_level is all that is required to get this diagnostic.

b) bsl is only calculated if (present(Time) .and. GV%Boussinesq) - and Time implies HA_CSp. The requirement for an active harmonic analysis needs to be removed, but in addition bsl needs to be calculated in the nonBoussinesq case.

[herrwang0]

I agree with Alan that BSL is a field we would like both as a diagnostic for output and an internal field for runtime harmonic analysis.

Therefore, I am not sure subroutine calculate_vertical_integrals is the right place for this particular calculation. I think MOM_diagnostics is best reserved for "end products". If there is other usage of BSL other than as a diagnostic, e.g. harmonic analysis, it should be somewhere else (in this case, MOM.F90 would be my choice). It would also save a bit of overhead for calling a function from a different module (for the scenario that you only want to calculate the harmonics while don't need the diagnostic BSL itself).

[herrwang0]

Could you add the equation for BSL in the comment? I trust you must have tested both EoS and layer-mode versions. But it would help me understand the code better with a little bit of context (is the equation $\eta_{BSL} = -\frac{1}{2\rho_0 h}\int^{\rho_b}_{\rho_s} (z+D)^2 d\rho$ ?).

[herrwang0]

This is a nice contribution. And I agree with Alan's comments that it BSL shouldn't depend on harmonic analysis and a non-Boussinesq version would be nice.

In addition, I think the calculation of BSL should be moved out of MOM_diagnostics, as it is used a part of harmonic analysis.

Also, I am having trouble understanding the calculation, especially how the implementation here is connected to Zaron & Ray (2023). Related, I did a little bit of derivation on this particular calculation and I think it could be more economic to use the PBCE term (in split mode). It could also help with non-Boussinesq case. I understand that this is probably beyond the scope of this PR, so I am happy to have a sidebar conversation or share some notes if there is interest.

[herrwang0]

1. In line 1048 of MOM.F90, I'm not 100% sure if Time_local is the correct time to use.
2. I think eta calculated in line 973 of MOM_diagnostics.F90 should be the same as ssh calculated in line 1025 of MOM.F90. But then in line 1075 of MOM.F90, ssh is calculated differently. I don't fully understand the algorithm for the latter, but it seems more accurate based on the harmonic analysis and the comparison with TPXO data. I'm not sure if this is the SSH that should be used for calculating the baroclinic sea level.

1. Time_local is the timestamp at the end of each time step. This is the timestamp of most diagnostics outside of barotropic solver.
2. In split mode, SSH and eta are two different variables. SSH is the average of SSH_inst over a cycle (could be a coupling time step or a forcing time step). And SSH_inst itself is the averaged sea surface height over a baroclinic timestep from the BT solver, it is consistent with the sea surface height used by pressure gradient anomaly in BT solver. On the other hand, eta calculated from find_eta is the final sea surface height at the end of each time step.

[c2xu]

Thanks @awallcraft and @herrwang0. I revised this PR so that the BSL calculation is done in its own subroutine in MOM_interface_heights.F90, and the activation of BSL calculation no longer relies on the harmonic analysis. I also made a minor change to the harmonic analysis of SSH in PR #999, so that the time used for the harmonic analyses of SSH and BSL is consistent. I also left a few comments on the derivation of the algorithm, but the full derivation can be found in this notes. I'm not sure about the algorithm in non-Boussinesq case, but I'll implement it in a separate PR in the future once I figure it out.