Skip to content

Commit

Permalink
(*)Parenthesize MOM_barotropic for FMAs
Browse files Browse the repository at this point in the history 
  Added parentheses to 18 expressions in btstep, and one more each in set_dtbt
and barotropic_init to exhibit rotationally consistent solutions when
fused-multiply-adds are enabled.  All answers are bitwise identical in cases
without FMAs, but answers could change when FMAs are enabled.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Jul 29, 2024
1 parent 64b851c commit 46e8b66
Showing 1 changed file with 47 additions and 47 deletions.
94 changes: 47 additions & 47 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -915,10 +915,10 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
do J=js-1,je ; do I=is-1,ie
q(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * &
((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / &
(max((G%areaT(i,j) * max(GV%Z_to_H*G%bathyT(i,j) + eta_in(i,j), 0.0) + &
G%areaT(i+1,j+1) * max(GV%Z_to_H*G%bathyT(i+1,j+1) + eta_in(i+1,j+1), 0.0)) + &
(G%areaT(i+1,j) * max(GV%Z_to_H*G%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + &
G%areaT(i,j+1) * max(GV%Z_to_H*G%bathyT(i,j+1) + eta_in(i,j+1), 0.0)), h_neglect) )
(max(((G%areaT(i,j) * max(GV%Z_to_H*G%bathyT(i,j) + eta_in(i,j), 0.0)) + &
(G%areaT(i+1,j+1) * max(GV%Z_to_H*G%bathyT(i+1,j+1) + eta_in(i+1,j+1), 0.0))) + &
((G%areaT(i+1,j) * max(GV%Z_to_H*G%bathyT(i+1,j) + eta_in(i+1,j), 0.0)) + &
(G%areaT(i,j+1) * max(GV%Z_to_H*G%bathyT(i,j+1) + eta_in(i,j+1), 0.0))), h_neglect) )
enddo ; enddo
else
!$OMP parallel do default(shared)
Expand All @@ -933,8 +933,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
do J=js-1,je ; do I=is-1,ie
q(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * &
((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / &
(max((G%areaT(i,j) * eta_in(i,j) + G%areaT(i+1,j+1) * eta_in(i+1,j+1)) + &
(G%areaT(i+1,j) * eta_in(i+1,j) + G%areaT(i,j+1) * eta_in(i,j+1)), h_neglect) )
(max(((G%areaT(i,j) * eta_in(i,j)) + (G%areaT(i+1,j+1) * eta_in(i+1,j+1))) + &
((G%areaT(i+1,j) * eta_in(i+1,j)) + (G%areaT(i,j+1) * eta_in(i,j+1))), h_neglect) )
enddo ; enddo
endif

Expand Down Expand Up @@ -1477,14 +1477,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
!$OMP parallel do default(shared)
do j=js,je ; do I=is-1,ie
Cor_ref_u(I,j) = &
((azon(I,j) * vbt_Cor(i+1,j) + czon(I,j) * vbt_Cor(i ,j-1)) + &
(bzon(I,j) * vbt_Cor(i ,j) + dzon(I,j) * vbt_Cor(i+1,j-1)))
(((azon(I,j) * vbt_Cor(i+1,j)) + (czon(I,j) * vbt_Cor(i ,j-1))) + &
((bzon(I,j) * vbt_Cor(i ,j)) + (dzon(I,j) * vbt_Cor(i+1,j-1))))
enddo ; enddo
!$OMP parallel do default(shared)
do J=js-1,je ; do i=is,ie
Cor_ref_v(i,J) = -1.0 * &
((amer(I-1,j) * ubt_Cor(I-1,j) + cmer(I ,j+1) * ubt_Cor(I ,j+1)) + &
(bmer(I ,j) * ubt_Cor(I ,j) + dmer(I-1,j+1) * ubt_Cor(I-1,j+1)))
(((amer(I-1,j) * ubt_Cor(I-1,j)) + (cmer(I ,j+1) * ubt_Cor(I ,j+1))) + &
((bmer(I ,j) * ubt_Cor(I ,j)) + (dmer(I-1,j+1) * ubt_Cor(I-1,j+1))))
enddo ; enddo

! Now start new halo updates.
Expand Down Expand Up @@ -1645,16 +1645,16 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
! gravity waves, but it is a conservative estimate since it ignores the
! stabilizing effect of the bottom drag.
Idt_max2 = 0.5 * (dgeo_de * (1.0 + 2.0*bebt)) * (G%IareaT(i,j) * &
((gtot_E(i,j) * (Datu(I,j)*G%IdxCu(I,j)) + &
gtot_W(i,j) * (Datu(I-1,j)*G%IdxCu(I-1,j))) + &
(gtot_N(i,j) * (Datv(i,J)*G%IdyCv(i,J)) + &
gtot_S(i,j) * (Datv(i,J-1)*G%IdyCv(i,J-1)))) + &
(((gtot_E(i,j) * (Datu(I,j)*G%IdxCu(I,j))) + &
(gtot_W(i,j) * (Datu(I-1,j)*G%IdxCu(I-1,j)))) + &
((gtot_N(i,j) * (Datv(i,J)*G%IdyCv(i,J))) + &
(gtot_S(i,j) * (Datv(i,J-1)*G%IdyCv(i,J-1))))) + &
((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
(G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1))) * CS%BT_Coriolis_scale**2 )
H_eff_dx2 = max(H_min_dyn * ((G%IdxT(i,j))**2 + (G%IdyT(i,j))**2), &
H_eff_dx2 = max(H_min_dyn * ((G%IdxT(i,j)**2) + (G%IdyT(i,j)**2)), &
G%IareaT(i,j) * &
((Datu(I,j)*G%IdxCu(I,j) + Datu(I-1,j)*G%IdxCu(I-1,j)) + &
(Datv(i,J)*G%IdyCv(i,J) + Datv(i,J-1)*G%IdyCv(i,J-1)) ) )
(((Datu(I,j)*G%IdxCu(I,j)) + (Datu(I-1,j)*G%IdxCu(I-1,j))) + &
((Datv(i,J)*G%IdyCv(i,J)) + (Datv(i,J-1)*G%IdyCv(i,J-1))) ) )
dyn_coef_max = CS%const_dyn_psurf * max(0.0, 1.0 - dtbt**2 * Idt_max2) / &
(dtbt**2 * H_eff_dx2)

Expand Down Expand Up @@ -1974,10 +1974,10 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
! On odd-steps, update v first.
!$OMP do schedule(static)
do J=jsv-1,jev ; do i=isv-1,iev+1
Cor_v(i,J) = -1.0*((amer(I-1,j) * ubt(I-1,j) + cmer(I,j+1) * ubt(I,j+1)) + &
(bmer(I,j) * ubt(I,j) + dmer(I-1,j+1) * ubt(I-1,j+1))) - Cor_ref_v(i,J)
PFv(i,J) = ((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_N(i,j) - &
(eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1)) * &
Cor_v(i,J) = -1.0*(((amer(I-1,j) * ubt(I-1,j)) + (cmer(I,j+1) * ubt(I,j+1))) + &
((bmer(I,j) * ubt(I,j)) + (dmer(I-1,j+1) * ubt(I-1,j+1)))) - Cor_ref_v(i,J)
PFv(i,J) = (((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_N(i,j)) - &
((eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1))) * &
dgeo_de * CS%IdyCv(i,J)
enddo ; enddo
!$OMP end do nowait
Expand Down Expand Up @@ -2049,11 +2049,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
! Now update the zonal velocity.
!$OMP do schedule(static)
do j=jsv,jev ; do I=isv-1,iev
Cor_u(I,j) = ((azon(I,j) * vbt(i+1,J) + czon(I,j) * vbt(i,J-1)) + &
(bzon(I,j) * vbt(i,J) + dzon(I,j) * vbt(i+1,J-1))) - &
Cor_u(I,j) = (((azon(I,j) * vbt(i+1,J)) + (czon(I,j) * vbt(i,J-1))) + &
((bzon(I,j) * vbt(i,J)) + (dzon(I,j) * vbt(i+1,J-1)))) - &
Cor_ref_u(I,j)
PFu(I,j) = ((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_E(i,j) - &
(eta_PF_BT(i+1,j)-eta_PF(i+1,j))*gtot_W(i+1,j)) * &
PFu(I,j) = (((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_E(i,j)) - &
((eta_PF_BT(i+1,j)-eta_PF(i+1,j))*gtot_W(i+1,j))) * &
dgeo_de * CS%IdxCu(I,j)
enddo ; enddo
!$OMP end do nowait
Expand Down Expand Up @@ -2128,11 +2128,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
! On even steps, update u first.
!$OMP do schedule(static)
do j=jsv-1,jev+1 ; do I=isv-1,iev
Cor_u(I,j) = ((azon(I,j) * vbt(i+1,J) + czon(I,j) * vbt(i,J-1)) + &
(bzon(I,j) * vbt(i,J) + dzon(I,j) * vbt(i+1,J-1))) - &
Cor_u(I,j) = (((azon(I,j) * vbt(i+1,J)) + (czon(I,j) * vbt(i,J-1))) + &
((bzon(I,j) * vbt(i,J)) + (dzon(I,j) * vbt(i+1,J-1)))) - &
Cor_ref_u(I,j)
PFu(I,j) = ((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_E(i,j) - &
(eta_PF_BT(i+1,j)-eta_PF(i+1,j))*gtot_W(i+1,j)) * &
PFu(I,j) = (((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_E(i,j)) - &
((eta_PF_BT(i+1,j)-eta_PF(i+1,j))*gtot_W(i+1,j))) * &
dgeo_de * CS%IdxCu(I,j)
enddo ; enddo
!$OMP end do nowait
Expand Down Expand Up @@ -2206,20 +2206,20 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
if (CS%use_old_coriolis_bracket_bug) then
!$OMP do schedule(static)
do J=jsv-1,jev ; do i=isv,iev
Cor_v(i,J) = -1.0*((amer(I-1,j) * ubt(I-1,j) + bmer(I,j) * ubt(I,j)) + &
(cmer(I,j+1) * ubt(I,j+1) + dmer(I-1,j+1) * ubt(I-1,j+1))) - Cor_ref_v(i,J)
PFv(i,J) = ((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_N(i,j) - &
(eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1)) * &
Cor_v(i,J) = -1.0*(((amer(I-1,j) * ubt(I-1,j)) + (bmer(I,j) * ubt(I,j))) + &
((cmer(I,j+1) * ubt(I,j+1)) + (dmer(I-1,j+1) * ubt(I-1,j+1)))) - Cor_ref_v(i,J)
PFv(i,J) = (((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_N(i,j)) - &
((eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1))) * &
dgeo_de * CS%IdyCv(i,J)
enddo ; enddo
!$OMP end do nowait
else
!$OMP do schedule(static)
do J=jsv-1,jev ; do i=isv,iev
Cor_v(i,J) = -1.0*((amer(I-1,j) * ubt(I-1,j) + cmer(I,j+1) * ubt(I,j+1)) + &
(bmer(I,j) * ubt(I,j) + dmer(I-1,j+1) * ubt(I-1,j+1))) - Cor_ref_v(i,J)
PFv(i,J) = ((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_N(i,j) - &
(eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1)) * &
Cor_v(i,J) = -1.0*(((amer(I-1,j) * ubt(I-1,j)) + (cmer(I,j+1) * ubt(I,j+1))) + &
((bmer(I,j) * ubt(I,j)) + (dmer(I-1,j+1) * ubt(I-1,j+1)))) - Cor_ref_v(i,J)
PFv(i,J) = (((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_N(i,j)) - &
((eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1))) * &
dgeo_de * CS%IdyCv(i,J)
enddo ; enddo
!$OMP end do nowait
Expand Down Expand Up @@ -2576,14 +2576,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
do k=1,nz
do j=js,je ; do I=is-1,ie
accel_layer_u(I,j,k) = (u_accel_bt(I,j) - &
((pbce(i+1,j,k) - gtot_W(i+1,j)) * e_anom(i+1,j) - &
(pbce(i,j,k) - gtot_E(i,j)) * e_anom(i,j)) * CS%IdxCu(I,j) )
(((pbce(i+1,j,k) - gtot_W(i+1,j)) * e_anom(i+1,j)) - &
((pbce(i,j,k) - gtot_E(i,j)) * e_anom(i,j))) * CS%IdxCu(I,j) )
if (abs(accel_layer_u(I,j,k)) < accel_underflow) accel_layer_u(I,j,k) = 0.0
enddo ; enddo
do J=js-1,je ; do i=is,ie
accel_layer_v(i,J,k) = (v_accel_bt(i,J) - &
((pbce(i,j+1,k) - gtot_S(i,j+1)) * e_anom(i,j+1) - &
(pbce(i,j,k) - gtot_N(i,j)) * e_anom(i,j)) * CS%IdyCv(i,J) )
(((pbce(i,j+1,k) - gtot_S(i,j+1)) * e_anom(i,j+1)) - &
((pbce(i,j,k) - gtot_N(i,j)) * e_anom(i,j))) * CS%IdyCv(i,J) )
if (abs(accel_layer_v(i,J,k)) < accel_underflow) accel_layer_v(i,J,k) = 0.0
enddo ; enddo
enddo
Expand Down Expand Up @@ -2904,8 +2904,8 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
! This is pretty accurate for gravity waves, but it is a conservative
! estimate since it ignores the stabilizing effect of the bottom drag.
Idt_max2 = 0.5 * (1.0 + 2.0*CS%bebt) * (G%IareaT(i,j) * &
((gtot_E(i,j)*Datu(I,j)*G%IdxCu(I,j) + gtot_W(i,j)*Datu(I-1,j)*G%IdxCu(I-1,j)) + &
(gtot_N(i,j)*Datv(i,J)*G%IdyCv(i,J) + gtot_S(i,j)*Datv(i,J-1)*G%IdyCv(i,J-1))) + &
(((gtot_E(i,j)*Datu(I,j)*G%IdxCu(I,j)) + (gtot_W(i,j)*Datu(I-1,j)*G%IdxCu(I-1,j))) + &
((gtot_N(i,j)*Datv(i,J)*G%IdyCv(i,J)) + (gtot_S(i,j)*Datv(i,J-1)*G%IdyCv(i,J-1)))) + &
((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
(G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1))) * CS%BT_Coriolis_scale**2 )
if (Idt_max2 * min_max_dt2 > 1.0) min_max_dt2 = 1.0 / Idt_max2
Expand Down Expand Up @@ -4850,10 +4850,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
if (G%mask2dT(i,j)+G%mask2dT(i,j+1)+G%mask2dT(i+1,j)+G%mask2dT(i+1,j+1)>0.) then
CS%q_D(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * &
((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / &
(Z_to_H * max(((G%areaT(i,j) * max(Mean_SL+G%bathyT(i,j),0.0) + &
G%areaT(i+1,j+1) * max(Mean_SL+G%bathyT(i+1,j+1),0.0)) + &
(G%areaT(i+1,j) * max(Mean_SL+G%bathyT(i+1,j),0.0) + &
G%areaT(i,j+1) * max(Mean_SL+G%bathyT(i,j+1),0.0))), GV%H_subroundoff) )
(Z_to_H * max((((G%areaT(i,j) * max(Mean_SL+G%bathyT(i,j),0.0)) + &
(G%areaT(i+1,j+1) * max(Mean_SL+G%bathyT(i+1,j+1),0.0))) + &
((G%areaT(i+1,j) * max(Mean_SL+G%bathyT(i+1,j),0.0)) + &
(G%areaT(i,j+1) * max(Mean_SL+G%bathyT(i,j+1),0.0)))), GV%H_subroundoff) )
else ! All four h points are masked out so q_D(I,J) will is meaningless
CS%q_D(I,J) = 0.
endif
Expand Down

0 comments on commit 46e8b66

Please sign in to comment.