Skip to content

Commit

Permalink
(*)Parenthesize set_viscous_ML for FMAs
Browse files Browse the repository at this point in the history 
  Added parentheses to 19 expressions in set_viscous_ML, set_u_at_v and
set_v_at_u to treat the velocities at both edges of a tracer cell equivalently
when fused-multiply-adds are enabled, and thereby to exhibit exhibit
rotationally consistent solutions.  Also swapped the order of the u- and
v-components in the u-point calculation of Uh2 to mirror the order of the
corresponging v-point calculation for the same purpose.  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 ebf02a9 commit c0bef18
Showing 1 changed file with 28 additions and 28 deletions.
56 changes: 28 additions & 28 deletions src/parameterizations/vertical/MOM_set_viscosity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1840,8 +1840,8 @@ function set_v_at_u(v, h, G, GV, i, j, k, mask2dCv, OBC)
hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
set_v_at_u = 0.0
if (hwt_tot > 0.0) set_v_at_u = &
((hwt(0,0) * v(i,J,k) + hwt(1,-1) * v(i+1,J-1,k)) + &
(hwt(1,0) * v(i+1,J,k) + hwt(0,-1) * v(i,J-1,k))) / hwt_tot
(((hwt(0,0) * v(i,J,k)) + (hwt(1,-1) * v(i+1,J-1,k))) + &
((hwt(1,0) * v(i+1,J,k)) + (hwt(0,-1) * v(i,J-1,k)))) / hwt_tot

end function set_v_at_u

Expand Down Expand Up @@ -1885,8 +1885,8 @@ function set_u_at_v(u, h, G, GV, i, j, k, mask2dCu, OBC)
hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
set_u_at_v = 0.0
if (hwt_tot > 0.0) set_u_at_v = &
((hwt(0,0) * u(I,j,k) + hwt(-1,1) * u(I-1,j+1,k)) + &
(hwt(-1,0) * u(I-1,j,k) + hwt(0,1) * u(I,j+1,k))) / hwt_tot
(((hwt(0,0) * u(I,j,k)) + (hwt(-1,1) * u(I-1,j+1,k))) + &
((hwt(-1,0) * u(I-1,j,k)) + (hwt(0,1) * u(I,j+1,k)))) / hwt_tot

end function set_u_at_v

Expand Down Expand Up @@ -2154,8 +2154,8 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
if (associated(tv%p_surf)) press(I) = press(I) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i+1,j))
k2 = max(1,nkml)
I_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
T_EOS(I) = (h(i,j,k2)*tv%T(i,j,k2) + h(i+1,j,k2)*tv%T(i+1,j,k2)) * I_2hlay
S_EOS(I) = (h(i,j,k2)*tv%S(i,j,k2) + h(i+1,j,k2)*tv%S(i+1,j,k2)) * I_2hlay
T_EOS(I) = ((h(i,j,k2)*tv%T(i,j,k2)) + (h(i+1,j,k2)*tv%T(i+1,j,k2))) * I_2hlay
S_EOS(I) = ((h(i,j,k2)*tv%S(i,j,k2)) + (h(i+1,j,k2)*tv%S(i+1,j,k2))) * I_2hlay
enddo
call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, tv%eqn_of_state, &
(/Isq-G%IsdB+1,Ieq-G%IsdB+1/) )
Expand All @@ -2170,13 +2170,13 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
if (hlay > h_tiny) then ! Only consider non-vanished layers.
I_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
v_at_u = 0.5 * (h(i,j,k) * (v(i,J,k) + v(i,J-1,k)) + &
h(i+1,j,k) * (v(i+1,J,k) + v(i+1,J-1,k))) * I_2hlay
Uh2 = ((uhtot(I) - htot(I)*u(I,j,k))**2 + (vhtot(I) - htot(I)*v_at_u)**2)
v_at_u = 0.5 * ((h(i,j,k) * (v(i,J,k) + v(i,J-1,k))) + &
(h(i+1,j,k) * (v(i+1,J,k) + v(i+1,J-1,k)))) * I_2hlay
Uh2 = (uhtot(I) - htot(I)*u(I,j,k))**2 + (vhtot(I) - htot(I)*v_at_u)**2

if (use_EOS) then
T_lay = (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k)) * I_2hlay
S_lay = (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k)) * I_2hlay
T_lay = ((h(i,j,k)*tv%T(i,j,k)) + (h(i+1,j,k)*tv%T(i+1,j,k))) * I_2hlay
S_lay = ((h(i,j,k)*tv%S(i,j,k)) + (h(i+1,j,k)*tv%S(i+1,j,k))) * I_2hlay
if (nonBous_ML) then
gHprime = (GV%g_Earth * GV%H_to_RZ) * (dSpV_dT(I) * (Thtot(I) - T_lay*htot(I)) + &
dSpV_dS(I) * (Shtot(I) - S_lay*htot(I)))
Expand Down Expand Up @@ -2211,11 +2211,11 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
do I=Isq,Ieq ; if (do_i(I)) then
htot(I) = htot(I) + 0.5 * (h(i,j,k) + h(i+1,j,k))
uhtot(I) = uhtot(I) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(I,j,k)
vhtot(I) = vhtot(I) + 0.25 * (h(i,j,k) * (v(i,J,k) + v(i,J-1,k)) + &
h(i+1,j,k) * (v(i+1,J,k) + v(i+1,J-1,k)))
vhtot(I) = vhtot(I) + 0.25 * ((h(i,j,k) * (v(i,J,k) + v(i,J-1,k))) + &
(h(i+1,j,k) * (v(i+1,J,k) + v(i+1,J-1,k))))
if (use_EOS) then
Thtot(I) = Thtot(I) + 0.5 * (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k))
Shtot(I) = Shtot(I) + 0.5 * (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k))
Thtot(I) = Thtot(I) + 0.5 * ((h(i,j,k)*tv%T(i,j,k)) + (h(i+1,j,k)*tv%T(i+1,j,k)))
Shtot(I) = Shtot(I) + 0.5 * ((h(i,j,k)*tv%S(i,j,k)) + (h(i+1,j,k)*tv%S(i+1,j,k)))
else
Rhtot(i) = Rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * GV%Rlay(k)
endif
Expand Down Expand Up @@ -2379,8 +2379,8 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)

! visc%tbl_thick_shelf_u(I,j) = max(CS%Htbl_shelf_min, &
! dztot(I) / (0.5 + sqrt(0.25 + &
! (htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))**2 / &
! (ustar(i))**2 )) )
! ((htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))**2) / &
! (ustar(i)**2) )) )
ustar1 = ustar(i)
h2f2 = (htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)) + h_neglect*CS%omega)**2
tbl_thick = max(CS%Htbl_shelf_min, &
Expand Down Expand Up @@ -2433,8 +2433,8 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
if (associated(tv%p_surf)) press(i) = press(i) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i,j+1))
k2 = max(1,nkml)
I_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
T_EOS(i) = (h(i,j,k2)*tv%T(i,j,k2) + h(i,j+1,k2)*tv%T(i,j+1,k2)) * I_2hlay
S_EOS(i) = (h(i,j,k2)*tv%S(i,j,k2) + h(i,j+1,k2)*tv%S(i,j+1,k2)) * I_2hlay
T_EOS(i) = ((h(i,j,k2)*tv%T(i,j,k2)) + (h(i,j+1,k2)*tv%T(i,j+1,k2))) * I_2hlay
S_EOS(i) = ((h(i,j,k2)*tv%S(i,j,k2)) + (h(i,j+1,k2)*tv%S(i,j+1,k2))) * I_2hlay
enddo
call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, &
tv%eqn_of_state, (/is-G%IsdB+1,ie-G%IsdB+1/) )
Expand All @@ -2449,13 +2449,13 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
if (hlay > h_tiny) then ! Only consider non-vanished layers.
I_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
u_at_v = 0.5 * (h(i,j,k) * (u(I-1,j,k) + u(I,j,k)) + &
h(i,j+1,k) * (u(I-1,j+1,k) + u(I,j+1,k))) * I_2hlay
Uh2 = ((uhtot(I) - htot(I)*u_at_v)**2 + (vhtot(I) - htot(I)*v(i,J,k))**2)
u_at_v = 0.5 * ((h(i,j,k) * (u(I-1,j,k) + u(I,j,k))) + &
(h(i,j+1,k) * (u(I-1,j+1,k) + u(I,j+1,k)))) * I_2hlay
Uh2 = (vhtot(i) - htot(i)*v(i,J,k))**2 + (uhtot(i) - htot(i)*u_at_v)**2

if (use_EOS) then
T_lay = (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k)) * I_2hlay
S_lay = (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k)) * I_2hlay
T_lay = ((h(i,j,k)*tv%T(i,j,k)) + (h(i,j+1,k)*tv%T(i,j+1,k))) * I_2hlay
S_lay = ((h(i,j,k)*tv%S(i,j,k)) + (h(i,j+1,k)*tv%S(i,j+1,k))) * I_2hlay
if (nonBous_ML) then
gHprime = (GV%g_Earth * GV%H_to_RZ) * (dSpV_dT(i) * (Thtot(i) - T_lay*htot(i)) + &
dSpV_dS(i) * (Shtot(i) - S_lay*htot(i)))
Expand Down Expand Up @@ -2490,11 +2490,11 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
do i=is,ie ; if (do_i(i)) then
htot(i) = htot(i) + 0.5 * (h(i,J,k) + h(i,j+1,k))
vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,J,k)
uhtot(i) = uhtot(i) + 0.25 * (h(i,j,k) * (u(I-1,j,k) + u(I,j,k)) + &
h(i,j+1,k) * (u(I-1,j+1,k) + u(I,j+1,k)))
uhtot(i) = uhtot(i) + 0.25 * ((h(i,j,k) * (u(I-1,j,k) + u(I,j,k))) + &
(h(i,j+1,k) * (u(I-1,j+1,k) + u(I,j+1,k))))
if (use_EOS) then
Thtot(i) = Thtot(i) + 0.5 * (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k))
Shtot(i) = Shtot(i) + 0.5 * (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k))
Thtot(i) = Thtot(i) + 0.5 * ((h(i,j,k)*tv%T(i,j,k)) + (h(i,j+1,k)*tv%T(i,j+1,k)))
Shtot(i) = Shtot(i) + 0.5 * ((h(i,j,k)*tv%S(i,j,k)) + (h(i,j+1,k)*tv%S(i,j+1,k)))
else
Rhtot(i) = Rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * GV%Rlay(k)
endif
Expand Down

0 comments on commit c0bef18

Please sign in to comment.