Baroclinic sea level calculation

src/core/MOM.F90

@@ -99,7 +99,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
@@ -405,6 +405,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
@@ -591,6 +593,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]
 
@@ -1047,6 +1051,10 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
                           G, GV, US, CS%diagnostics_CSp)
       call post_tracer_diagnostics_at_sync(CS%Tracer_reg, h, CS%diag_pre_sync, CS%diag, G, GV, CS%t_dyn_rel_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, 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
@@ -3553,6 +3561,14 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
     CS%dt_obc_seg_time = Time + CS%dt_obc_seg_interval
   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.)
+  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, &

src/core/MOM_interface_heights.F90

@@ -17,7 +17,7 @@ module MOM_interface_heights
 
 #include <MOM_memory.h>
 
-public find_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
+public find_eta, find_bsl, 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
@@ -263,6 +263,106 @@ 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 (BSL) following Equation (4) of Zaron & Ray (2023), JPO,
+!! 53, 2591-2596. DOI: 10.1175/JPO-D-23-0073.1. The baroclinic pressure at z = 0 is defined by
+!!   p'(0) = p(0) - (1/H) int_{-H}^0 p(z) dz,
+!! and the baroclinic sea level is defined by
+!!   bsl = p'(0) / (g rho_0).
+!! Here, p'(0) can be computed using the full density/pressure fields, or the density/pressure
+!! anomalies. The results differ by an offset of 0.5 H, but the physical interpretation of BSL
+!! is unaffected. Because the BSL is a perturbation from the mean state, the interpretation is
+!! meaningful only if the mean value is removed. In this subroutine, the anomalies are used in
+!! the calculation, so that the BSL represents deviations from zero, instead of 0.5 H.
+subroutine find_bsl(h, tv, G, GV, US, 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, 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]
+    p, &            ! Pressure at the top of a layer [R L2 T-2 ~> Pa]
+    dp, &           ! Pressure change across a layer [R L2 T-2 ~> Pa]
+    p_int, &        ! Vertical integral of pressure at the bottom of a layer [R L T-2 ~> Pa m]
+                    ! or that normalized by (g rho_0) in Boussinesq and non-EOS mode [Z ~> m2]
+    dp_int          ! Vertical integral of pressure difference across a layer [R L T-2 ~> Pa m]
+  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
+  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
+
+  call find_eta(h, tv, G, GV, US, eta, halo_size=1, dZref=dZ_ref)
+
+  !$OMP parallel default(shared)
+  !$OMP do
+  do j=js,je ; do i=is,ie
+    p(i,j) = 0.0 ; p_int(i,j) = 0.0 ; bsl(i,j) = 0.0
+    bathyT(i,j) = G%bathyT(i,j) + dZ_ref
+  enddo ; enddo
+
+  if (GV%Boussinesq) then
+    if (associated(tv%eqn_of_state)) then
+      ! In this case, the algorithm follows
+      !   bsl = (rho_s/rho_0) eta - (int_{-H}^{eta} p(z) dz) / (g rho_0 H),
+      ! where rho_0 is the reference density, rho_s is the surface density anomaly, and p(z) is
+      ! the pressure anomaly. The depth integral is evaluated layer-by-layer. Within each layer,
+      !   int_{z_{k+1}}^{z_k} p(z) dz
+      !       = p(z_k) h_k + int_{z_{k+1}}^{z_k} [p(z) - p(z_k)] dz,
+      ! where the integrated pressure difference is given by dp_int.
+      do k=1,nz
+        call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), eta(:,:,k), eta(:,:,k+1), GV%Rho0, &
+                            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 (G%mask2dT(i,j) > 0.0) then
+          p_int(i,j) = p_int(i,j) + dp_int(i,j) + p(i,j) * (GV%H_to_Z * h(i,j,k))
+          p(i,j) = p(i,j) + dp(i,j)
+        endif ; enddo ; enddo
+      enddo
+      !$OMP do
+      do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
+        bsl(i,j) = (GV%g_prime(1) / GV%g_Earth - 1.0) * eta(i,j,1) - &
+                   p_int(i,j) / ((GV%g_Earth * GV%Rho0) * bathyT(i,j))
+      endif ; enddo ; enddo
+    else ! (.not.associated(tv%eqn_of_state))
+      ! The integration by parts is used to represent the depth integral by the density anomaly,
+      !   int_{-H}^{eta} p(z) dz = int_{-H}^{eta} g rho(z) (H + z) dz.
+      ! Within the k-th layer, the density is constant, i.e., rho(z) = rho_k, so that
+      !   int_{z_{k+1}}^{z_k} rho(z) (H + z) dz
+      !       = rho_k [H (z_k - z_{k+1}) + 0.5 (z_k^2 - z_{k+1}^2)].
+      ! With some algebra, the full depth integral can be expressed as
+      !   int_{-H}^eta rho(z) (H + z) dz
+      !       = 0.5 sum_{k=1}^N (rho_k - rho_{k-1}) (z_k + H)^2
+      ! where rho_0 = 0. Note that this expression remains unchanged if computed using the full
+      ! density field, except that rho_0 must be replaced by the reference density.
+      !$OMP do
+      do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
+        do k = 1,nz
+          p_int(i,j) = p_int(i,j) + (GV%g_prime(k) / GV%g_Earth) * &
+                                    ((eta(i,j,k) + bathyT(i,j)) * (eta(i,j,k) + bathyT(i,j)))
+        enddo
+        bsl(i,j) = (GV%g_prime(1) / GV%g_Earth - 1.0) * eta(i,j,1) - &
+                   0.5 * (p_int(i,j) - (eta(i,j,1) + bathyT(i,j)) * &
+                                       (eta(i,j,1) + bathyT(i,j))) / bathyT(i,j)
+      endif ; enddo ; enddo
+    endif ! (associated(tv%eqn_of_state))
+  else ! (.not.GV%Boussiensq)
+    ! I'm not sure if the same definition of baroclinic pressure is applicable to non-Boussinesq
+    ! case, since the pressure and velocity modes are not the same. I'll leave it to the future.
+    call MOM_error(FATAL, "BSL calculation is not yet available in non-Boussinesq mode.")
+  endif
+  !$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)

src/diagnostics/MOM_diagnostics.F90

@@ -24,7 +24,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_col_mass
+use MOM_interface_heights, only : find_eta, find_bsl, 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
@@ -114,6 +114,7 @@ module MOM_diagnostics
   integer :: id_drho_dT        = -1, id_drho_dS        = -1
   integer :: id_h_pre_sync     = -1
   integer :: id_tosq           = -1, id_sosq           = -1
+  integer :: id_bsl            = -1
 
   !>@}
   type(wave_speed_CS) :: wave_speed  !< Wave speed control struct
@@ -902,8 +903,9 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
     btm_pres,&! The pressure at the ocean bottom, or CMIP variable 'pbo'.
               ! This is the column mass multiplied by gravity plus the pressure
               ! at the ocean surface [R L2 T-2 ~> Pa].
-    tr_int    ! vertical integral of a tracer times density,
+    tr_int, & ! vertical integral of a tracer times density,
               ! (Rho_0 in a Boussinesq model) [Conc R Z ~> Conc kg m-2].
+    bsl       ! Baroclinic sea level [Z ~> m].
   real    :: IG_Earth  ! Inverse of gravitational acceleration [T2 Z L-2 ~> s2 m-1].
 
   integer :: i, j, k, is, ie, js, je, nz
@@ -951,6 +953,11 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
     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, bsl, dZref=G%Z_ref)
+    call post_data(CS%id_bsl, bsl, CS%diag)
+  endif
+
 end subroutine calculate_vertical_integrals
 
 !> This subroutine calculates terms in the mechanical energy budget.
@@ -2101,6 +2108,9 @@ 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)
 
   ! Register time derivatives and allocate memory for diagnostics that need
   ! access from across several modules.

src/diagnostics/MOM_harmonic_analysis.F90

@@ -89,7 +89,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
@@ -266,6 +266,9 @@ subroutine HA_init(Time, US, param_file, nc, CS)
   call get_param(param_file, mdl, "HA_SSH", HA_ssh, &
                  "If true, perform harmonic analysis of sea serface 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) call HA_register('bsl', 'h', CS)
   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)