Skip to content

Bugfix linear wave drag in MOM_barotropic #974

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
merged 3 commits into from
Sep 24, 2025
Merged

Bugfix linear wave drag in MOM_barotropic #974

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2941519
Bug fix for linear wave drag in MOM_barotropic
herrwang0 Sep 23, 2025
0c7b41a
Suppress warning message of negative eta in land
herrwang0 Sep 23, 2025
665d547
Flip the order of acceleration and velocity chksum
herrwang0 Sep 23, 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
26 changes: 16 additions & 10 deletions src/core/MOM_barotropic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1593,22 +1593,28 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
endif
if (CS%linear_wave_drag) then
!$OMP do
do j=js,je ; do I=is-1,ie ; if (CS%lin_drag_u(I,j) > 0.0) then
do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) * CS%lin_drag_u(I,j) > 0.0) then
Htot = 0.5 * (eta(i,j) + eta(i+1,j))
if (GV%Boussinesq) &
Htot = Htot + 0.5*GV%Z_to_H * (CS%bathyT(i,j) + CS%bathyT(i+1,j))
bt_rem_u(I,j) = bt_rem_u(I,j) * (Htot / (Htot + CS%lin_drag_u(I,j) * dtbt))

Rayleigh_u(I,j) = CS%lin_drag_u(I,j) / Htot
! If Htot==0., linear wave drag is not used and Rayleigh_u = 0.0 (from initialization)
! and bt_rem_u is unmodified.
if (Htot > 0.0) then
bt_rem_u(I,j) = bt_rem_u(I,j) * (Htot / (Htot + CS%lin_drag_u(I,j) * dtbt))
Rayleigh_u(I,j) = CS%lin_drag_u(I,j) / Htot
endif
endif ; enddo ; enddo
!$OMP do
do J=js-1,je ; do i=is,ie ; if (CS%lin_drag_v(i,J) > 0.0) then
do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) * CS%lin_drag_v(i,J) > 0.0) then
Htot = 0.5 * (eta(i,j) + eta(i,j+1))
if (GV%Boussinesq) &
Htot = Htot + 0.5*GV%Z_to_H * (CS%bathyT(i,j) + CS%bathyT(i,j+1))
bt_rem_v(i,J) = bt_rem_v(i,J) * (Htot / (Htot + CS%lin_drag_v(i,J) * dtbt))

Rayleigh_v(i,J) = CS%lin_drag_v(i,J) / Htot
! If Htot==0., linear wave drag is not used and Rayleigh_v = 0.0 (from initialization)
! and bt_rem_v is unmodified.
if (Htot > 0.0) then
bt_rem_v(i,J) = bt_rem_v(i,J) * (Htot / (Htot + CS%lin_drag_v(i,J) * dtbt))
Rayleigh_v(i,J) = CS%lin_drag_v(i,J) / Htot
endif
endif ; enddo ; enddo
endif

Expand Down Expand Up @@ -2614,7 +2620,7 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL
enddo ; enddo
else
do j=js,je ; do i=is,ie
if (eta_IC(i,j) < 0.0) then
if ((eta_IC(i,j) < 0.0) .and. (G%mask2dT(i,j) > 0.0)) then
write(mesg,'(" at ", ES12.4, ES12.4, i7, i7)') &
G%geoLonT(i,j), G%geoLatT(i,j), i + G%HI%idg_offset, j + G%HI%jdg_offset
call MOM_error(FATAL, "btstep: negative eta_IC at start of a non-Boussinesq barotropic solver "//&
Expand Down Expand Up @@ -2909,7 +2915,7 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL
enddo ; enddo
else
do j=js,je ; do i=is,ie
if (eta(i,j) < 0.0) then
if ((eta(i,j) < 0.0) .and. (G%mask2dT(i,j) > 0.0)) then
write(mesg,'(" at ", ES12.4, ES12.4, i7, i7)') &
G%geoLonT(i,j), G%geoLatT(i,j), i + G%HI%idg_offset, j + G%HI%jdg_offset
if (CS%bt_limit_integral_transport) &
Expand Down
15 changes: 8 additions & 7 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -707,13 +707,13 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
call cpu_clock_end(id_clock_mom_update)

if (CS%debug) then
call MOM_accel_chksum("Predictor accel", CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv, &
CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, symmetric=sym)
call uvchksum("Predictor 1 [uv]", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
call hchksum(h, "Predictor 1 h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
call uvchksum("Predictor 1 [uv]h", uh, vh, G%HI,haloshift=2, &
symmetric=sym, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
call MOM_accel_chksum("Predictor accel", CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv, &
CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, symmetric=sym)
! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
call MOM_state_chksum("Predictor 1 init", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1, &
symmetric=sym)
if (debug_redundant) then
Expand Down Expand Up @@ -976,14 +976,15 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
call cpu_clock_end(id_clock_mom_update)

if (CS%debug) then
call MOM_accel_chksum("Corrector accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, &
CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, &
symmetric=sym)
call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
call hchksum(h, "Corrector 1 h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
call uvchksum("Corrector 1 [uv]h", uh, vh, G%HI, haloshift=2, &
symmetric=sym, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)
call MOM_accel_chksum("Corrector accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, &
CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, &
symmetric=sym)
! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)

endif

! u <- u + dt d/dz visc d/dz u
Expand Down