Skip to content

Baroclinic sea level calculation #1104

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
c2xu wants to merge 7 commits into
base: dev/gfdl
Choose a base branch
from
Open

Baroclinic sea level calculation #1104

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
41f2e9f
Baroclinic sea level calculation
c2xu May 1, 2026
cd016b8
Baroclinic sea level calculation
c2xu May 1, 2026
8ef7786
Baroclinic sea level calculation
c2xu May 2, 2026
650e20b
Baroclinic sea level calculation
chengz2 May 7, 2026
cec2d80
Baroclinic sea level calculation
chengz2 Jun 6, 2026
b31b374
Baroclinic sea level calculation
chengz2 Jul 20, 2026
1390983
Baroclinic sea level calculation
c2xu Jul 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 21 additions & 1 deletion src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ module MOM
use MOM_harmonic_analysis, only : HA_accum, harmonic_analysis_CS
use MOM_hor_index, only : hor_index_type, hor_index_init
use MOM_hor_index, only : rotate_hor_index
use MOM_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz
use MOM_interface_heights, only : find_eta, find_bsl, calc_derived_thermo, thickness_to_dz
use MOM_interface_filter, only : interface_filter, interface_filter_init, interface_filter_end
use MOM_interface_filter, only : interface_filter_CS
use MOM_internal_tides, only : int_tide_CS
Expand Down Expand Up @@ -342,6 +342,7 @@ module MOM
!! segment data with every call to update_OBC_tracer_data.
type(time_type) :: dt_obc_seg_interval !< A time_time representation of dt_obc_seg_period.
type(time_type) :: dt_obc_seg_time !< The next time OBC segment update is applied to OBGC tracers.
real :: rho_bsl !< Surface density for calculating the baroclinic sea level [R ~> kg m-3]

real, dimension(:,:), pointer :: frac_shelf_h => NULL() !< fraction of total area occupied
!! by ice shelf [nondim]
Expand Down Expand Up @@ -428,6 +429,8 @@ module MOM
!< Pointer to the control structure used for an alternate version of the mode-split RK2 dynamics
type(harmonic_analysis_CS), pointer :: HA_CSp => NULL()
!< Pointer to the control structure for harmonic analysis
logical :: HA_bsl
!< If true, perform harmonic analysis of baroclinic sea level
type(thickness_diffuse_CS) :: thickness_diffuse_CSp
!< Pointer to the control structure used for the isopycnal height diffusive transport.
!! This is also common referred to as Gent-McWilliams diffusion
Expand Down Expand Up @@ -616,6 +619,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1]
real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
ssh ! sea surface height, which may be based on eta_av [Z ~> m]
real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
bsl ! Baroclinic sea level [Z ~> m]
real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%GV)) :: &
dz ! Vertical distance across layers [Z ~> m]

Expand Down Expand Up @@ -1078,6 +1083,10 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
call post_tracer_integral_diagnostics(G, GV, US, CS%Tracer_reg, h, CS%tv, CS%diag)

call diag_copy_diag_to_storage(CS%diag_pre_sync, h, CS%diag)
if (associated(CS%HA_CSp) .and. CS%HA_bsl) then
call find_bsl(h, CS%tv, G, GV, US, CS%rho_bsl, bsl, dZref=G%Z_ref)
call HA_accum('bsl', bsl, Time_local, G, CS%HA_CSp)
endif
if (showCallTree) call callTree_waypoint("finished calculate_diagnostic_fields (step_MOM)")
call disable_averaging(CS%diag)
CS%t_dyn_rel_diag = 0.0
Expand Down Expand Up @@ -3716,6 +3725,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
endif
endif

if (associated(CS%HA_CSp)) then
call get_param(param_file, "MOM", "HA_BSL",CS%HA_bsl, &
"If true, perform harmonic analysis of baroclinic sea level.", &
default=.false., do_not_log=.true.)
call get_param(param_file, "MOM", "RHO_BSL", CS%rho_bsl, &
"Surface density for baroclinic sea level calculation.", &
units='kg m-3', default=1025.0, scale=US%kg_m3_to_R, do_not_log=.not.CS%HA_bsl)
else
CS%HA_bsl = .false.
endif

call callTree_waypoint("dynamics initialized (initialize_MOM)")

CS%mixedlayer_restrat = mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, &
Expand Down
101 changes: 100 additions & 1 deletion src/core/MOM_interface_heights.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ module MOM_interface_heights

#include <MOM_memory.h>

public find_eta, find_dz_for_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
public find_eta, find_bsl, find_dz_for_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
public calc_derived_thermo
public convert_MLD_to_ML_thickness
public find_rho_bottom, find_col_avg_SpV, find_col_mass
Expand Down Expand Up @@ -290,6 +290,105 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref)

end subroutine find_eta_2d

!> Calculates the baroclinic sea level, following Xu et al., to be submitted to JPO
subroutine find_bsl(h, tv, G, GV, US, rho_s, bsl, dZref)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, intent(in) :: rho_s !< Surface density [R ~> kg m-3]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
!! thermodynamic variables
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: bsl !< Baroclinic sea level [Z ~> m]
real, optional, intent(in) :: dZref !< The difference in the
!! reference height between G%bathyT and eta [Z ~> m]. The default is 0

! Local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! layer interface heights [Z ~> m]
real, dimension(SZI_(G),SZJ_(G)) :: &
bathyT, & ! Bathymetry at T points plus dZ_ref [Z ~> m]
pt, & ! Pressure at the top of a layer [R L2 T-2 ~> Pa]
pb, & ! Pressure at the bottom of a layer [R L2 T-2 ~> Pa]
gz, & ! Geopotential at the bottom of a layer [L2 T-2 ~> m2 s-2]
dp, & ! Pressure change across a layer in Boussinesq mode [R L2 T-2 ~> Pa]

@Hallberg-NOAA Hallberg-NOAA Jul 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The comment here should make it clear that this is the anomaly in the pressure change across a layer relative to what would be found if the density of the layer were rho_s. Without this clarification, it is very hard to see how the code in this routine would generate the baroclinic sea level anomaly. This comment also pertains to the descriptions of dg, dp_int, dg_int and p_int in the next few lines, all of which are also anomalies relative to the reference state.

dg, & ! Geopotential change across a layer in non-Boussinesq mode [L2 T-2 ~> m2 s-2]
dp_int, & ! Layer-integrated pressure change in Boussinesq mode [R L2 Z T-2 ~> Pa m]
dg_int, & ! Layer-integrated geopotential change in non-Boussinesq mode [R L4 T-4 ~> Pa m2 s-2]
p_int ! Vertical integral of pressure at the bottom of a layer [R L2 Z T-2 ~> Pa m]
! or that scaled by GV%g_Earth in non-Boussinesq and EOS mode [R L4 T-4 ~> Pa m2 s-2]
! or that normalized by GV%g_Earth in non-EOS mode [R Z2 ~> Pa s2]
real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]
! dZ_ref is 0 unless the optional argument dZref is present
real :: I_gEarth ! The inverse of the gravitational acceleration [T2 Z L-2 ~> s2 m-1]
real :: SpV_s ! Specific volume of the surface layer [R-1 ~> m3 kg-1]
logical, dimension(SZI_(G),SZJ_(G)) :: maskT ! Mask at T points for skipping land points in calculations
integer :: i, j, k, is, ie, js, je, nz

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref

I_gEarth = 1.0 / GV%g_Earth; SpV_s = 1.0 / rho_s

call find_eta(h, tv, G, GV, US, eta, halo_size=1, dZref=dZ_ref)

@Hallberg-NOAA Hallberg-NOAA Jul 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This call to find_eta is only needed in Boussinesq mode. Because this can be an expensive call in non-Boussinesq mode, it should be avoided when it is not needed. In fact, given the substantial differences in the algorithms and the meanings of the variables between Boussinesq and non-Boussinesq modes, I would be tempted to go so far as to split the body of this routine into two separate routines that are called depending on whether the Boussinsesq approximation is being made.


!$OMP parallel default(shared)
!$OMP do
do j=js,je ; do i=is,ie
pt(i,j) = 0.0 ; pb(i,j) = 0.0 ; gz(i,j) = 0.0 ; p_int(i,j) = 0.0 ; bsl(i,j) = 0.0
bathyT(i,j) = G%bathyT(i,j) + dZ_ref
maskT(i,j) = G%mask2dT(i,j) > 0.0 .and. bathyT(i,j) > 0.0
enddo ; enddo

if (associated(tv%eqn_of_state)) then
if (GV%Boussinesq) then
do k=1,nz
call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), eta(:,:,k), eta(:,:,k+1), rho_s, &
GV%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dp, dp_int)
!$OMP do
do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
p_int(i,j) = p_int(i,j) + (pt(i,j) * (GV%H_to_Z * h(i,j,k)) + dp_int(i,j))
pt(i,j) = pt(i,j) + dp(i,j)
endif ; enddo ; enddo
enddo
!$OMP do
do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
bsl(i,j) = - (p_int(i,j) * I_gEarth) / (rho_s * bathyT(i,j))

@Hallberg-NOAA Hallberg-NOAA Jul 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suspect that there are cases where dividing by bathyT here will not work because in some cases (such as the Great Lakes) bathyT can be equal to or above the reference geoid and still have water. Would it make more sense to use the vertical sum of the layer thicknesses (perhaps in a max with GV%H_subroundoff) in the denominator instead?

endif ; enddo ; enddo
else ! (.not. GV%Boussinesq)
do k=1,nz
!$OMP do
do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
dp(i,j) = GV%g_Earth * (GV%H_to_RZ * h(i,j,k))
pb(i,j) = pt(i,j) + dp(i,j)
endif ; enddo ; enddo
call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), pt, pb, SpV_s, G%HI, &
tv%eqn_of_state, US, dg, dg_int)
!$OMP do
do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
gz(i,j) = gz(i,j) + dg(i,j)
p_int(i,j) = p_int(i,j) + (gz(i,j) * dp(i,j) + dg_int(i,j))
pt(i,j) = pb(i,j)
endif ; enddo ; enddo
enddo
!$OMP do
do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
bsl(i,j) = (p_int(i,j) * (I_gEarth * I_gEarth)) / (rho_s * bathyT(i,j))

@Hallberg-NOAA Hallberg-NOAA Jul 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that the denominator here should be the sum of the layer thicknesses.

Also, there should be spaces around the semicolon on the next line to follow the pattern used elsewhere.

endif ;enddo ; enddo
endif ! (GV%Boussinesq)
else ! (.not. associated(tv%eqn_of_state))
!$OMP do
do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
do k=2,nz
p_int(i,j) = p_int(i,j) + (GV%Rlay(k) - GV%Rlay(k-1)) * &
((eta(i,j,k) + bathyT(i,j)) * (eta(i,j,k) + bathyT(i,j)))
enddo
bsl(i,j) = - 0.5 * (p_int(i,j) / (rho_s * bathyT(i,j)))
endif ; enddo ; enddo
endif ! (associated(tv%eqn_of_state))
!$OMP end parallel

end subroutine find_bsl

!> Calculate derived thermodynamic quantities for re-use later.
subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)
Expand Down
17 changes: 16 additions & 1 deletion src/diagnostics/MOM_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ module MOM_diagnostics
use MOM_error_handler, only : MOM_error, FATAL, WARNING
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_interface_heights, only : find_eta, find_dz_for_eta, find_col_mass
use MOM_interface_heights, only : find_eta, find_bsl, find_dz_for_eta, find_col_mass
use MOM_spatial_means, only : global_area_mean, global_layer_mean
use MOM_spatial_means, only : global_volume_mean, global_area_integral
use MOM_tracer_registry, only : tracer_registry_type, post_tracer_transport_diagnostics
Expand Down Expand Up @@ -64,6 +64,7 @@ module MOM_diagnostics
!! non-Boussinesq layer thicknesses as are used to find the free
!! surface height, instead of using an approximate thickness
!! based on division by the mid-layer density.
real :: rho_bsl !< Surface density for calculating the baroclinic sea level [R ~> kg m-3]

type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
!! regulate the timing of diagnostic output.
Expand Down Expand Up @@ -126,6 +127,7 @@ module MOM_diagnostics
integer :: id_h_pre_sync = -1
integer :: id_tosq = -1, id_sosq = -1
integer :: id_t20d = -1, id_t17d = -1
integer :: id_bsl = -1

!>@}
type(wave_speed_CS) :: wave_speed !< Wave speed control struct
Expand Down Expand Up @@ -920,6 +922,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
z_top, & ! Height of the top of a layer or the ocean [Z ~> m].
z_bot, & ! Height of the bottom of a layer (for id_mass) or the
! (positive) depth of the ocean (for id_col_ht) [Z ~> m].
bsl, & ! Baroclinic sea level [Z ~> m]
mass, & ! integrated mass of the water column [R Z ~> kg m-2]. For
! non-Boussinesq models this is rho*dz. For Boussinesq
! models, this is either the integral of in-situ density
Expand Down Expand Up @@ -983,6 +986,12 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
endif
if (CS%id_col_mass > 0) call post_data(CS%id_col_mass, mass, CS%diag)
endif

if (CS%id_bsl > 0) then
call find_bsl(h, tv, G, GV, US, CS%rho_bsl, bsl, dZref=G%Z_ref)
call post_data(CS%id_bsl, bsl, CS%diag)
endif

if (CS%id_t20d > 0 .or. CS%id_t17d > 0) then
call PPM%init(GV%ke, h_neglect=0.)
do j=js,je ; do i=is,ie
Expand Down Expand Up @@ -2294,6 +2303,12 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
CS%id_pbo = register_diag_field('ocean_model', 'pbo', diag%axesT1, Time, &
long_name='Sea Water Pressure at Sea Floor', standard_name='sea_water_pressure_at_sea_floor', &
units='Pa', conversion=US%RL2_T2_to_Pa)
CS%id_bsl = register_diag_field('ocean_model', 'bsl', diag%axesT1, Time, &
long_name='Baroclinic Sea Level', standard_name='baroclinic_sea_level', &
units='m', conversion=US%Z_to_m)
call get_param(param_file, mdl, "RHO_BSL", CS%rho_bsl, &
"Surface density for baroclinic sea level calculation.", &
units='kg m-3', default=1025.0, scale=US%kg_m3_to_R, do_not_log=CS%id_bsl<=0)

! Register time derivatives and allocate memory for diagnostics that need
! access from across several modules.
Expand Down
10 changes: 8 additions & 2 deletions src/diagnostics/MOM_harmonic_analysis.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ subroutine HA_init(Time, US, param_file, nc, CS)
type(HA_type) :: ha1 !< A temporary, null field used for initializing CS%list
real :: HA_start_time !< Start time of harmonic analysis [T ~> s]
real :: HA_end_time !< End time of harmonic analysis [T ~> s]
logical :: HA_ssh, HA_ubt, HA_vbt
logical :: HA_ssh, HA_bsl, HA_ubt, HA_vbt
character(len=40) :: mdl="MOM_harmonic_analysis" !< This module's name
character(len=255) :: mesg
integer :: year, month, day, hour, minute, second
Expand Down Expand Up @@ -268,8 +268,14 @@ subroutine HA_init(Time, US, param_file, nc, CS)

! Register variables/fields to be analyzed
call get_param(param_file, mdl, "HA_SSH", HA_ssh, &
"If true, perform harmonic analysis of sea serface height.", default=.false.)
"If true, perform harmonic analysis of sea surface height.", default=.false.)
if (HA_ssh) call HA_register('ssh', 'h', CS)
call get_param(param_file, mdl, "HA_BSL", HA_bsl, &
"If true, perform harmonic analysis of baroclinic sea level.", default=.false.)
if (HA_bsl) then
call HA_register('bsl', 'h', CS)
endif

call get_param(param_file, mdl, "HA_UBT", HA_ubt, &
"If true, perform harmonic analysis of zonal barotropic velocity.", default=.false.)
if (HA_ubt) call HA_register('ubt', 'u', CS)
Expand Down
Loading