From 29415196b450973c32de26d6aa96a71e6e59c89a Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 22 Sep 2025 22:42:48 -0400 Subject: [PATCH 1/3] Bug fix for linear wave drag in MOM_barotropic * Linear wave drag is limited to be only applied to land points, using velocity point masks mask2dC[uv]. * Rayleigh_[uv] calculation and bt_rem_[uv] update from linear wave drag is limited for Htot>0 only. This patch eliminates potential NaN in Rayleigh_[uv] in an unusual scenario that Htot==0.0 and lin_drag_[uv]/=0. The changes do not change answers: bt_rem_[uv] is zero at land points regardless. Rayleigh_[uv] is added to [uv]_accel_bt which is masked before updating velocity. --- src/core/MOM_barotropic.F90 | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 3da021dbcc..eccd6e4a4e 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -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 From 0c7b41a60993b344b78df87de294f6cf4dc53e5b Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 22 Sep 2025 22:52:10 -0400 Subject: [PATCH 2/3] Suppress warning message of negative eta in land In MOM_barotropic and non-Boussinesq mode, warning message on negative eta is now only issued at wet points, consistently with Boussinesq. --- src/core/MOM_barotropic.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index eccd6e4a4e..ac2e668f8e 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2620,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 "//& @@ -2915,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) & From 665d54749eec317450c945bbcac1985a29169cc1 Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 22 Sep 2025 23:05:26 -0400 Subject: [PATCH 3/3] Flip the order of acceleration and velocity chksum In MOM_dynamics_split_RK2, now accleration chksum is printed before velocity with debug on, so that we could know which accleration term is responsible for a NaN in velocity. --- src/core/MOM_dynamics_split_RK2.F90 | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 0df9dedc74..2f345ff9e9 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -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 @@ -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