From 41f2e9ff4050a1211d0cec2472f2b3317acfaaa6 Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 1 May 2026 13:27:57 -0400 Subject: [PATCH 1/7] Baroclinic sea level calculation Algorithms implemented in MOM_interface_heights.F90 --- src/core/MOM_interface_heights.F90 | 99 +++++++++++++++++++++++++++++- 1 file changed, 98 insertions(+), 1 deletion(-) diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index 8821c0cb16..5097760eb1 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -19,7 +19,7 @@ module MOM_interface_heights #include -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 @@ -290,6 +290,103 @@ 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, 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] + 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] + dp, & ! Pressure change across a layer in Boussinesq mode [R L2 T-2 ~> Pa] + ! or 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 L T-2 ~> Pa m] + ! or layer-integrated geopotential change in non-Boussinesq mode [L2 T-2 Z ~> m3 s-2] + p_int ! Vertical integral of pressure at the bottom of a layer [R L2 T-2 Z ~> Pa m] + ! or that scaled by GV%g_Earth in 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 :: Rho0 ! Reference density, which must be the surface density [R ~> kg m-3] + 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 + + Rho0 = GV%Rlay(1) + + 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 + pt(i,j) = 0.0 ; pb(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 (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), 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) + (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 (G%mask2dT(i,j) > 0.0) then + bsl(i,j) = - (p_int(i,j) * I_gEarth) / (Rho0 * bathyT(i,j)) + endif ; enddo ; enddo + else ! (.not. GV%Boussinesq) + do k=1,nz + !$OMP do + do j=js,je ; do i=is,ie + pb(i,j) = pt(i,j) + (GV%g_Earth * GV%H_to_RZ) * h(i,j,k) + enddo ; enddo + call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), pt, pb, 0.0, G%HI, & + tv%eqn_of_state, US, dp, dp_int) + !$OMP do + do j=js,je ; do i=is,ie + p_int(i,j) = p_int(i,j) + (pt(i,j) * dp(i,j) + dp_int(i,j)) + pt(i,j) = pb(i,j) + enddo ; enddo + enddo + !$OMP do + do j=js,je ; do i=is,ie + bsl(i,j) = (eta(i,j,1) - ((p_int(i,j) * (I_gEarth * I_gEarth)) / & + (Rho0 * bathyT(i,j)) - 0.5 * bathyT(i,j))) & + + 0.5 * (eta(i,j,1) * eta(i,j,1)) / BathyT(i,j) + enddo ; enddo + endif ! (GV%Boussinesq) + else ! (.not. associated(tv%eqn_of_state)) + !$OMP do + do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) 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) / (Rho0 * 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) From cd016b834d3dcafde73f884b8809f3339e4619d4 Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 1 May 2026 19:26:45 -0400 Subject: [PATCH 2/7] Baroclinic sea level calculation Implemented in MOM_diagnostics.F90 --- src/diagnostics/MOM_diagnostics.F90 | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index de44b07c49..5c9ad968fb 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -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 @@ -126,6 +126,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 @@ -920,6 +921,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 @@ -983,6 +985,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, 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 @@ -2294,6 +2302,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. From 8ef778615e44c6169a8993a02b4af4c51d9ada27 Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 1 May 2026 21:12:34 -0400 Subject: [PATCH 3/7] Baroclinic sea level calculation Implemented for harmonic analysis --- src/core/MOM.F90 | 18 +++++++++++++++++- src/diagnostics/MOM_harmonic_analysis.F90 | 10 ++++++++-- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index b51627b08e..dee8196169 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -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 @@ -428,6 +428,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 @@ -616,6 +618,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] @@ -1078,6 +1082,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, 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 @@ -3716,6 +3724,14 @@ 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.) + 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, & diff --git a/src/diagnostics/MOM_harmonic_analysis.F90 b/src/diagnostics/MOM_harmonic_analysis.F90 index 45d2c76146..db1ebff146 100644 --- a/src/diagnostics/MOM_harmonic_analysis.F90 +++ b/src/diagnostics/MOM_harmonic_analysis.F90 @@ -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 @@ -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) From 650e20b734df34c5f8236ba3b878fb0bd682ec73 Mon Sep 17 00:00:00 2001 From: William Xu Date: Thu, 7 May 2026 11:38:17 -0600 Subject: [PATCH 4/7] Baroclinic sea level calculation * Rewrite algorithm for the EOS and non-Boussinesq case. * Skip calculations at grid points with potentially zero depth. --- src/core/MOM_interface_heights.F90 | 44 ++++++++++++++++-------------- 1 file changed, 24 insertions(+), 20 deletions(-) diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index 5097760eb1..50e80a9d27 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -308,10 +308,11 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) 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] - ! or geopotential change across a layer in non-Boussinesq mode [L2 T-2 ~> m2 s-2] + 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 L T-2 ~> Pa m] - ! or layer-integrated geopotential change in non-Boussinesq mode [L2 T-2 Z ~> m3 s-2] + dg_int, & ! Layer-integrated geopotential change in non-Boussinesq mode [L2 T-2 Z ~> m3 s-2] p_int ! Vertical integral of pressure at the bottom of a layer [R L2 T-2 Z ~> Pa m] ! or that scaled by GV%g_Earth in 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] @@ -319,6 +320,8 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) ! 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 :: Rho0 ! Reference density, which must be the surface density [R ~> kg m-3] + real :: SpV0 ! Reference specific volume, which must be the surface specific volume [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 @@ -327,15 +330,16 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) I_gEarth = 1.0 / GV%g_Earth - Rho0 = GV%Rlay(1) + Rho0 = GV%Rlay(1) ; SpV0 = 1.0 / GV%Rlay(1) 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 - pt(i,j) = 0.0 ; pb(i,j) = 0.0 ; p_int(i,j) = 0.0 ; bsl(i,j) = 0.0 + 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 @@ -344,39 +348,39 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), eta(:,:,k), eta(:,:,k+1), 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 + 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 (G%mask2dT(i,j) > 0.0) then + do j=js,je ; do i=is,ie ; if (maskT(i,j)) then bsl(i,j) = - (p_int(i,j) * I_gEarth) / (Rho0 * bathyT(i,j)) endif ; enddo ; enddo else ! (.not. GV%Boussinesq) do k=1,nz !$OMP do - do j=js,je ; do i=is,ie - pb(i,j) = pt(i,j) + (GV%g_Earth * GV%H_to_RZ) * h(i,j,k) - enddo ; enddo - call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), pt, pb, 0.0, G%HI, & - tv%eqn_of_state, US, dp, dp_int) + 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, SpV0, G%HI, & + tv%eqn_of_state, US, dg, dg_int) !$OMP do - do j=js,je ; do i=is,ie - p_int(i,j) = p_int(i,j) + (pt(i,j) * dp(i,j) + dp_int(i,j)) + 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) - enddo ; enddo + endif ; enddo ; enddo enddo !$OMP do - do j=js,je ; do i=is,ie - bsl(i,j) = (eta(i,j,1) - ((p_int(i,j) * (I_gEarth * I_gEarth)) / & - (Rho0 * bathyT(i,j)) - 0.5 * bathyT(i,j))) & - + 0.5 * (eta(i,j,1) * eta(i,j,1)) / BathyT(i,j) - enddo ; enddo + do j=js,je ; do i=is,ie ; if (maskT(i,j)) then + bsl(i,j) = (p_int(i,j) * (I_gEarth * I_gEarth)) / (Rho0 * bathyT(i,j)) + endif ;enddo ; enddo endif ! (GV%Boussinesq) else ! (.not. associated(tv%eqn_of_state)) !$OMP do - do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then + 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))) From cec2d801df9d724e68065ae41feaeda655b3474b Mon Sep 17 00:00:00 2001 From: William Xu Date: Sat, 6 Jun 2026 16:05:51 -0600 Subject: [PATCH 5/7] Baroclinic sea level calculation * Fix dimensions of variables * Set the surface density as a run-time parameter, RHO_BSL, with the default value being RHO_BSL = 1025 kg m-3. --- src/core/MOM.F90 | 6 +++++- src/core/MOM_interface_heights.F90 | 28 +++++++++++++--------------- src/diagnostics/MOM_diagnostics.F90 | 6 +++++- 3 files changed, 23 insertions(+), 17 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index dee8196169..adcdd814ac 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -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] @@ -1083,7 +1084,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS 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 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)") @@ -3728,6 +3729,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & 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 diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index 50e80a9d27..8157cad7a3 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -291,10 +291,11 @@ 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, bsl, dZref) +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 @@ -311,16 +312,15 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) 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] 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 L T-2 ~> Pa m] - dg_int, & ! Layer-integrated geopotential change in non-Boussinesq mode [L2 T-2 Z ~> m3 s-2] - p_int ! Vertical integral of pressure at the bottom of a layer [R L2 T-2 Z ~> Pa m] - ! or that scaled by GV%g_Earth in Boussinesq and EOS mode [R L4 T-4 ~> Pa 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 :: Rho0 ! Reference density, which must be the surface density [R ~> kg m-3] - real :: SpV0 ! Reference specific volume, which must be the surface specific volume [R-1 ~> m3 kg-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 @@ -328,9 +328,7 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref - I_gEarth = 1.0 / GV%g_Earth - - Rho0 = GV%Rlay(1) ; SpV0 = 1.0 / GV%Rlay(1) + 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) @@ -345,7 +343,7 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) 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), Rho0, & + 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 @@ -355,7 +353,7 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) 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) / (Rho0 * bathyT(i,j)) + bsl(i,j) = - (p_int(i,j) * I_gEarth) / (rho_s * bathyT(i,j)) endif ; enddo ; enddo else ! (.not. GV%Boussinesq) do k=1,nz @@ -364,7 +362,7 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) 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, SpV0, G%HI, & + 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 @@ -375,7 +373,7 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) 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)) / (Rho0 * bathyT(i,j)) + bsl(i,j) = (p_int(i,j) * (I_gEarth * I_gEarth)) / (rho_s * bathyT(i,j)) endif ;enddo ; enddo endif ! (GV%Boussinesq) else ! (.not. associated(tv%eqn_of_state)) @@ -385,7 +383,7 @@ subroutine find_bsl(h, tv, G, GV, US, bsl, dZref) 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) / (Rho0 * bathyT(i,j))) + 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 diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 5c9ad968fb..f8aafa5ad9 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -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. @@ -987,7 +988,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) endif if (CS%id_bsl > 0) then - call find_bsl(h, tv, G, GV, US, bsl, dZref=G%Z_ref) + 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 @@ -2305,6 +2306,9 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 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. From b31b374bee58c4b85777f6fb5a051ec983796648 Mon Sep 17 00:00:00 2001 From: William Xu Date: Mon, 20 Jul 2026 12:34:51 -0600 Subject: [PATCH 6/7] Baroclinic sea level calculation * Clarify the description of local variables in find_bsl(), and the description of RHO_BSL in MOM_parameter_doc * Move find_eta so that it is not called in non-Boussinesq mode --- src/core/MOM.F90 | 4 +++- src/core/MOM_interface_heights.F90 | 27 +++++++++++++++------------ src/diagnostics/MOM_diagnostics.F90 | 6 ++++-- 3 files changed, 22 insertions(+), 15 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index adcdd814ac..a4bbc0b7b9 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3730,7 +3730,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & "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.", & + "Reference density for calculating the baroclinic sea level based on "//& + "integrals of the differences between the actual layer density and this "//& + "value. By defulat, it should be the surface density.", & 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. diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index 8157cad7a3..feef02d30a 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -290,12 +290,13 @@ 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 +!> Calculates the baroclinic sea level, following Xu et al., 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] + !! which is the reference density for all anomalies defined below 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 @@ -307,14 +308,16 @@ subroutine find_bsl(h, tv, G, GV, US, rho_s, bsl, dZref) 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] - 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] + pt, & ! Pressure anomaly at the top of a layer in Boussinesq mode [R L2 T-2 ~> Pa] + ! Or pressure at the top of a layer in non-Boussinesq mode [R L2 T-2 ~> Pa] + pb, & ! Pressure at the bottom of a layer in non-Boussinesq mode[R L2 T-2 ~> Pa] + gz, & ! Geopotential anomaly at the bottom of a layer [L2 T-2 ~> m2 s-2] + dp, & ! Pressure anomaly change across a layer in Boussinesq mode [R L2 T-2 ~> Pa] + ! Or pressure change across a layer in non-Boussinesq mode [R L2 T-2 ~> Pa] + dg, & ! Geopotential anomaly change across a layer in non-Boussinesq mode [L2 T-2 ~> m2 s-2] + dp_int, & ! Layer-integrated pressure anomaly change in Boussinesq mode [R L2 Z T-2 ~> Pa m] + dg_int, & ! Layer-integrated geopotential anomaly change in non-Boussinesq mode [R L4 T-4 ~> Pa m2 s-2] + p_int ! Vertical integral of pressure anomaly 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] @@ -330,8 +333,6 @@ subroutine find_bsl(h, tv, G, GV, US, rho_s, bsl, 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) - !$OMP parallel default(shared) !$OMP do do j=js,je ; do i=is,ie @@ -342,6 +343,7 @@ subroutine find_bsl(h, tv, G, GV, US, rho_s, bsl, dZref) if (associated(tv%eqn_of_state)) then if (GV%Boussinesq) then + call find_eta(h, tv, G, GV, US, eta, halo_size=1, dZref=dZ_ref) 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) @@ -374,9 +376,10 @@ subroutine find_bsl(h, tv, G, GV, US, rho_s, bsl, dZref) !$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)) - endif ;enddo ; enddo + endif ; enddo ; enddo endif ! (GV%Boussinesq) else ! (.not. associated(tv%eqn_of_state)) + call find_eta(h, tv, G, GV, US, eta, halo_size=1, dZref=dZ_ref) !$OMP do do j=js,je ; do i=is,ie ; if (maskT(i,j)) then do k=2,nz diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index f8aafa5ad9..8f075c271a 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -2306,8 +2306,10 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 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.", & + call get_param(param_file, "MOM", "RHO_BSL", CS%rho_bsl, & + "Reference density for calculating the baroclinic sea level based on "//& + "integrals of the differences between the actual layer density and this "//& + "value. By defulat, it should be the surface density.", & 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 From 13909838673298cbaa693b696fe55581c5f9d913 Mon Sep 17 00:00:00 2001 From: Chengzhu Xu Date: Wed, 29 Jul 2026 17:09:28 -0600 Subject: [PATCH 7/7] Baroclinic sea level calculation Split find_bsl() into three subroutines, depending on whether an equation of state is used and whether the Boussinesq approximation is used. --- src/core/MOM_interface_heights.F90 | 209 ++++++++++++++++++++--------- 1 file changed, 147 insertions(+), 62 deletions(-) diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index feef02d30a..95b8c0c2ec 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -290,7 +290,17 @@ 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., submitted to JPO +!> Baroclinic sea level calculation, following Xu et al., submitted to JPO +!! This subroutine calculates the baroclinic sea level if an equation of state is not used, +!! and calls bsl_boussinesq or bsl_non_boussinesq otherwise. +!! Note: 1) The barotropic-baroclinic sea level decomposition is useful only in the deep ocean, +!! where baroclinic waves are "large", yielding sea level fluctuations of a few centimeters. +!! 2) The baroclinic sea level represents the temporal anomaly to a mean state around which the +!! dynamics are linearzed. Thus, it must always be interpreted with its temporal average removed, +!! and the "raw" baroclinic sea level calculated from this subroutine is not a useful quantity. +!! 3) The baroclinic sea level calculated in this subroutine is based on the linear decomposition. +!! It should thus not be used to interprete sea level fluctuations where nonlinearity dominates, +!! such as in shallow waters or over steep bottom topography. 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 @@ -308,22 +318,10 @@ subroutine find_bsl(h, tv, G, GV, US, rho_s, bsl, dZref) 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 anomaly at the top of a layer in Boussinesq mode [R L2 T-2 ~> Pa] - ! Or pressure at the top of a layer in non-Boussinesq mode [R L2 T-2 ~> Pa] - pb, & ! Pressure at the bottom of a layer in non-Boussinesq mode[R L2 T-2 ~> Pa] - gz, & ! Geopotential anomaly at the bottom of a layer [L2 T-2 ~> m2 s-2] - dp, & ! Pressure anomaly change across a layer in Boussinesq mode [R L2 T-2 ~> Pa] - ! Or pressure change across a layer in non-Boussinesq mode [R L2 T-2 ~> Pa] - dg, & ! Geopotential anomaly change across a layer in non-Boussinesq mode [L2 T-2 ~> m2 s-2] - dp_int, & ! Layer-integrated pressure anomaly change in Boussinesq mode [R L2 Z T-2 ~> Pa m] - dg_int, & ! Layer-integrated geopotential anomaly change in non-Boussinesq mode [R L4 T-4 ~> Pa m2 s-2] - p_int ! Vertical integral of pressure anomaly 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] + p_int ! Vertical integral of pressure anomaly at the bottom of a layer, + ! normalized by GV%g_Earth [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 @@ -331,67 +329,154 @@ subroutine find_bsl(h, tv, G, GV, US, rho_s, bsl, dZref) dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref + if (associated(tv%eqn_of_state)) then + if (GV%Boussinesq) then + call bsl_boussinesq(h, tv, G, GV, US, rho_s, bsl, dZ_ref) + else ! (.not. GV%Boussinesq) + call bsl_non_boussinesq(h, tv, G, GV, US, rho_s, bsl, dZ_ref) + endif ! (GV%Boussinesq) + else ! (.not. associated(tv%eqn_of_state)) + 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_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 + 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 + !$OMP end parallel + endif ! (associated(tv%eqn_of_state)) + +end subroutine find_bsl + +!> BSL calculation when an EOS is used and the Boussinesq approximation is used +subroutine bsl_boussinesq(h, tv, G, GV, US, rho_s, bsl, dZ_ref) + 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] + !! which is the reference density for all anomalies defined below + 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, intent(in) :: dZ_ref !< 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 anomaly at the top of a layer [R L2 T-2 ~> Pa] + dp, & ! Pressure anomaly change across a layer [R L2 T-2 ~> Pa] + dp_int, & ! Layer-integrated pressure anomaly change [R L2 Z T-2 ~> Pa m] + p_int ! Vertical integral of pressure anomaly at the bottom of a layer [R L2 Z T-2 ~> Pa m] + real :: I_gEarth ! The inverse of the gravitational acceleration [T2 Z L-2 ~> s2 m-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 + + I_gEarth = 1.0 / GV%g_Earth; + + 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 + pt(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 + + 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)) + endif ; enddo ; enddo + !$OMP end parallel + +end subroutine bsl_boussinesq + +!> BSL calculation when an EOS is used and the Boussinesq approximation is not used +subroutine bsl_non_boussinesq(h, tv, G, GV, US, rho_s, bsl, dZ_ref) + 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] + !! which is the reference density for all anomalies defined below + 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, intent(in) :: dZ_ref !< 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 anomaly at the bottom of a layer [L2 T-2 ~> m2 s-2] + dp, & ! Pressure change across a layer [R L2 T-2 ~> Pa] + dg, & ! Geopotential anomaly change across a layer [L2 T-2 ~> m2 s-2] + dg_int, & ! Layer-integrated geopotential anomaly change [R L4 T-4 ~> Pa m2 s-2] + g_int ! Vertically-integrated geopotential anomaly at the bottom of a layer [R L4 T-4 ~> Pa m2 s-2] + 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 + I_gEarth = 1.0 / GV%g_Earth; SpV_s = 1.0 / rho_s !$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 + pt(i,j) = 0.0 ; pb(i,j) = 0.0 ; gz(i,j) = 0.0 ; g_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 - call find_eta(h, tv, G, GV, US, eta, halo_size=1, dZref=dZ_ref) - 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)) - 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)) - endif ; enddo ; enddo - endif ! (GV%Boussinesq) - else ! (.not. associated(tv%eqn_of_state)) - call find_eta(h, tv, G, GV, US, eta, halo_size=1, dZref=dZ_ref) + do k=1,nz !$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))) + 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 - endif ! (associated(tv%eqn_of_state)) + 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) + g_int(i,j) = g_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) = (g_int(i,j) * (I_gEarth * I_gEarth)) / (rho_s * bathyT(i,j)) + endif ; enddo ; enddo !$OMP end parallel -end subroutine find_bsl +end subroutine bsl_non_boussinesq !> Calculate derived thermodynamic quantities for re-use later. subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)