diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 3da021dbcc..ac2e668f8e 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 @@ -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 "//& @@ -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) & 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