Skip to content

+Fix OBC reproducibility across restarts #907

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

Merged
Hallberg-NOAA merged 6 commits into from
Jun 6, 2025
Merged

+Fix OBC reproducibility across restarts #907

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
55017b7
Refactoring for OBC restart fix
Hallberg-NOAA May 23, 2025
580f3a2
+Add the runtime parameter VERTEX_SHEAR_OBC_BUG
Hallberg-NOAA May 27, 2025
89a0e5d
+Add calc_isoneutral_slopes optional argument OBC_N2
Hallberg-NOAA May 29, 2025
f3b6cf3
+Add the runtime parameter MIXING_COEFS_OBC_BUG
Hallberg-NOAA May 29, 2025
8689af6
+Add OBC_RESERVOIR_INIT_BUG and fix OBC restarts
Hallberg-NOAA May 29, 2025
97934a9
(*)Correct rotation with OBC tracer reservoirs
Hallberg-NOAA May 30, 2025
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
33 changes: 24 additions & 9 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,9 +113,11 @@ module MOM
use MOM_obsolete_diagnostics, only : register_obsolete_diagnostics
use MOM_open_boundary, only : ocean_OBC_type, OBC_registry_type
use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs
use MOM_open_boundary, only : setup_OBC_tracer_reservoirs, fill_temp_salt_segments
use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields
use MOM_open_boundary, only : open_boundary_setup_vert, update_OBC_segment_data
use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init, write_OBC_info, chksum_OBC_segments
use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init, open_boundary_halo_update
use MOM_open_boundary, only : write_OBC_info, chksum_OBC_segments
use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init
use MOM_porous_barriers, only : porous_barrier_CS
use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS
Expand Down Expand Up @@ -777,7 +779,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS

if (CS%VarMix%use_variable_mixing) then
call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag)
call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt)
call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, CS%OBC, dt)
call calc_depth_function(G, CS%VarMix)
call disable_averaging(CS%diag)
endif
Expand Down Expand Up @@ -1284,7 +1286,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, &
endif
endif
endif
if (CS%debug_OBCs .and. associated(CS%OBC)) call chksum_OBC_segments(CS%OBC, G, GV, US, 3)
! if (CS%debug_OBCs .and. associated(CS%OBC)) call chksum_OBC_segments(CS%OBC, G, GV, US, 3)

if (CS%do_dynamics .and. CS%split) then !--------------------------- start SPLIT
! This section uses a split time stepping scheme for the dynamic equations,
Expand Down Expand Up @@ -2067,7 +2069,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
if (.not. skip_diffusion) then
if (CS%VarMix%use_variable_mixing) then
call pass_var(CS%h, G%Domain)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, CS%OBC, dt_offline)
call calc_depth_function(G, CS%VarMix)
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
endif
Expand All @@ -2094,7 +2096,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
if (.not. skip_diffusion) then
if (CS%VarMix%use_variable_mixing) then
call pass_var(CS%h, G%Domain)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, CS%OBC, dt_offline)
call calc_depth_function(G, CS%VarMix)
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
endif
Expand Down Expand Up @@ -3260,7 +3262,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
endif

if (associated(OBC_in)) then
call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, CS%OBC)
if (use_temperature) then
call pass_var(CS%tv%T, G%Domain)
call pass_var(CS%tv%S, G%Domain)
call fill_temp_salt_segments(G, GV, US, CS%OBC, CS%tv)
endif

call rotate_OBC_init(OBC_in, G, CS%OBC)
call calc_derived_thermo(CS%tv, CS%h, G, GV, US)
call update_OBC_segment_data(G, GV, US, CS%OBC, CS%tv, CS%h, Time)
endif
Expand All @@ -3274,7 +3282,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
endif
if (use_ice_shelf) &
deallocate(frac_shelf_in,mass_shelf_in)
else
else ! The model is being run without grid rotation. This is true of all production runs.
if (use_ice_shelf) then
call initialize_ice_shelf(param_file, G, Time, ice_shelf_CSp, diag_ptr, Time_init, &
dirs%output_directory, calve_ice_shelf_bergs=point_calving)
Expand Down Expand Up @@ -3343,8 +3351,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &

if (ALE_remap_init_conds(CS%ALE_CSp) .and. .not. query_initialized(CS%h,"h",restart_CSp)) then
! This block is controlled by the ALE parameter REMAP_AFTER_INITIALIZATION.
! \todo This block exists for legacy reasons and we should phase it out of
! all examples. !###
! \todo This block exists for legacy reasons and we should phase it out of all examples. !###
if (CS%debug) then
call uvchksum("Pre ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1, unscale=US%L_T_to_m_s)
call hchksum(CS%h,"Pre ALE adjust init cond h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
Expand Down Expand Up @@ -3657,6 +3664,14 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
call register_diags_offline_transport(Time, CS%diag, CS%offline_CSp, GV, US)
endif

if (associated(CS%OBC)) then
! At this point any information related to the tracer reservoirs has either been read from
! the restart file or has been specified in the segments. Initialize the tracer reservoir
! values from the segments if they have not been set via the restart file.
call setup_OBC_tracer_reservoirs(G, GV, CS%OBC, restart_CSp)
call open_boundary_halo_update(G, CS%OBC)
endif

call register_obsolete_diagnostics(param_file, CS%diag)

if (use_frazil) then
Expand Down
84 changes: 76 additions & 8 deletions src/core/MOM_isopycnal_slopes.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,8 @@ module MOM_isopycnal_slopes

!> Calculate isopycnal slopes, and optionally return other stratification dependent functions such as N^2
!! and dz*S^2*g-prime used, or calculable from factors used, during the calculation.
subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, &
slope_x, slope_y, N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC)
subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, slope_x, slope_y, &
N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC, OBC_N2)
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
Expand Down Expand Up @@ -61,6 +61,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
!! Eady growth rate at v-points. [Z T-1 ~> m s-1]
integer, optional, intent(in) :: halo !< Halo width over which to compute
type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
logical, optional, intent(in) :: OBC_N2 !< If present and true, use interior data
!! to calculate stratification at open boundary
!! condition faces.

! Local variables
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
Expand Down Expand Up @@ -127,6 +130,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan

logical :: present_N2_u, present_N2_v
logical :: local_open_u_BC, local_open_v_BC ! True if u- or v-face OBCs exist anywhere in the global domain.
logical :: OBC_friendly ! If true, open boundary conditions are in use and only interior data should
! be used to calculate N2 at OBC faces.
integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of
! state calculations at u-points.
integer, dimension(2) :: EOSdom_v ! The shifted i-computational domain to use for equation of
Expand Down Expand Up @@ -155,9 +160,11 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan

local_open_u_BC = .false.
local_open_v_BC = .false.
OBC_friendly = .false.
if (present(OBC)) then ; if (associated(OBC)) then
local_open_u_BC = OBC%open_u_BCs_exist_globally
local_open_v_BC = OBC%open_v_BCs_exist_globally
if (present(OBC_N2)) OBC_friendly = OBC_N2
endif ; endif

use_EOS = associated(tv%eqn_of_state)
Expand Down Expand Up @@ -241,12 +248,12 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
enddo

do I=is-1,ie
GxSpV_u(I) = G_Rho0 !This will be changed if both use_EOS and allocated(tv%SpV_avg) are true
GxSpV_u(I) = G_Rho0 ! This will be changed if both use_EOS and allocated(tv%SpV_avg) are true
enddo
!$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, &
!$OMP h_neglect,dz_neglect,h_neglect2, &
!$OMP present_N2_u,G_Rho0,N2_u,slope_x,dzSxN,EOSdom_u,EOSdom_h1, &
!$OMP local_open_u_BC,dzu,OBC,use_stanley) &
!$OMP local_open_u_BC,dzu,OBC,use_stanley,OBC_friendly) &
!$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, &
!$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
!$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, &
Expand All @@ -266,6 +273,22 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
T_u(I) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1)))
S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1)))
enddo
if (OBC_friendly) then
do I=is-1,ie
l_seg = OBC%segnum_u(I,j)
if (l_seg /= OBC_NONE) then
if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then
pres_u(I) = pres(i,j,K)
T_u(I) = 0.5*(T(i,j,k) + T(i,j,k-1))
S_u(I) = 0.5*(S(i,j,k) + S(i,j,k-1))
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then
pres_u(I) = pres(i+1,j,K)
T_u(I) = 0.5*(T(i+1,j,k) + T(i+1,j,k-1))
S_u(I) = 0.5*(S(i+1,j,k) + S(i+1,j,k-1))
endif
endif
enddo
endif
call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, &
tv%eqn_of_state, EOSdom_u)
if (present_N2_u .or. (present(dzSxN))) then
Expand Down Expand Up @@ -338,8 +361,21 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
! The expression for drdz above is mathematically equivalent to:
! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
! ((hg2L/haL) + (hg2R/haR))
! This is the gradient of density along geopotentials.
! which is an estimate of the gradient of density across geopotentials.
if (present_N2_u) then
if (OBC_friendly) then ; if (OBC%segnum_u(I,j) /= OBC_NONE) then
l_seg = OBC%segnum_u(I,j)
if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then
drdz = drdkL / dzaL ! Note that drdz is not used for slopes at OBC faces.
if (use_EOS .and. allocated(tv%SpV_avg)) &
GxSpV_u(I) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then
drdz = drdkR / dzaR
if (use_EOS .and. allocated(tv%SpV_avg)) &
GxSpV_u(I) = GV%g_Earth * 0.5 * (tv%SpV_avg(i+1,j,k) + tv%SpV_avg(i+1,j,k-1))
endif
endif ; endif

N2_u(I,j,K) = GxSpV_u(I) * drdz * G%mask2dCu(I,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
endif

Expand Down Expand Up @@ -375,6 +411,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
endif
slope = slope * max(G%mask2dT(i,j), G%mask2dT(i+1,j))
endif

slope_x(I,j,K) = slope
if (present(dzSxN)) &
dzSxN(I,j,K) = sqrt( GxSpV_u(I) * max(0., (wtL * ( dzaL * drdkL )) &
Expand All @@ -391,7 +428,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
!$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, &
!$OMP h,h_neglect,e,dz_neglect, &
!$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,dzSyN,EOSdom_v, &
!$OMP dzv,local_open_v_BC,OBC,use_stanley) &
!$OMP dzv,local_open_v_BC,OBC,use_stanley,OBC_friendly) &
!$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, &
!$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
!$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, &
Expand All @@ -411,6 +448,22 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1)))
S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1)))
enddo
if (OBC_friendly) then
do i=is,ie
l_seg = OBC%segnum_v(i,J)
if (l_seg /= OBC_NONE) then
if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then
pres_v(i) = pres(i,j,K)
T_v(i) = 0.5*(T(i,j,k) + T(i,j,k-1))
S_v(i) = 0.5*(S(i,j,k) + S(i,j,k-1))
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then
pres_v(i) = pres(i,j+1,K)
T_v(i) = 0.5*(T(i,j+1,k) + T(i,j+1,k-1))
S_v(i) = 0.5*(S(i,j+1,k) + S(i,j+1,k-1))
endif
endif
enddo
endif
call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, &
tv%eqn_of_state, EOSdom_v)

Expand Down Expand Up @@ -490,8 +543,23 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
! The expression for drdz above is mathematically equivalent to:
! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
! ((hg2L/haL) + (hg2R/haR))
! This is the gradient of density along geopotentials.
if (present_N2_v) N2_v(i,J,K) = GxSpV_v(i) * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
! which is an estimate of the gradient of density across geopotentials.
if (present_N2_v) then
if (OBC_friendly) then ; if (OBC%segnum_v(i,J) /= OBC_NONE) then
l_seg = OBC%segnum_v(i,J)
if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then
drdz = drdkL / dzaL ! Note that drdz is not used for slopes at OBC faces.
if (use_EOS .and. allocated(tv%SpV_avg)) &
GxSpV_v(i) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then
drdz = drdkL / dzaL
if (use_EOS .and. allocated(tv%SpV_avg)) &
GxSpV_v(i) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j+1,k) + tv%SpV_avg(i,j+1,k-1))
endif
endif ; endif

N2_v(i,J,K) = GxSpV_v(i) * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
endif

if (use_EOS) then
drdy = ((wtA * drdjA + wtB * drdjB) / (wtA + wtB) - &
Expand Down
Loading