Skip to content

Commit

Permalink
(*)Parenthesize thickness_diffuse for FMAs
Browse files Browse the repository at this point in the history 
  Added parentheses to 17 expressions in thickness_diffuse_full,
thickness_diffuse and thickness_diffuse_init to give rotationally consistent
solutions when fused-multiply-adds are enabled.  One comment was also added to
note that the calculation of PE_release_h is does not exhibit rotational
symmetry when MEKE_GM_SRC_ALT is set to true.  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 56d053a commit 49419f7
Showing 1 changed file with 31 additions and 30 deletions.
61 changes: 31 additions & 30 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,12 +218,12 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
!$OMP parallel do default(shared)
do j=js,je ; do I=is-1,ie
KH_u_CFL(I,j) = (0.25*CS%max_Khth_CFL) / &
(dt * (G%IdxCu(I,j)*G%IdxCu(I,j) + G%IdyCu(I,j)*G%IdyCu(I,j)))
(dt * ((G%IdxCu(I,j)*G%IdxCu(I,j)) + (G%IdyCu(I,j)*G%IdyCu(I,j))))
enddo ; enddo
!$OMP parallel do default(shared)
do j=js-1,je ; do I=is,ie
KH_v_CFL(i,J) = (0.25*CS%max_Khth_CFL) / &
(dt * (G%IdxCv(i,J)*G%IdxCv(i,J) + G%IdyCv(i,J)*G%IdyCv(i,J)))
(dt * ((G%IdxCv(i,J)*G%IdxCv(i,J)) + (G%IdyCv(i,J)*G%IdyCv(i,J))))
enddo ; enddo

! Calculates interface heights, e, in [Z ~> m].
Expand Down Expand Up @@ -535,8 +535,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
enddo ; enddo
! diagnose diffusivity at T-points
do j=js,je ; do i=is,ie
Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j) + hu(I,j)*KH_u_lay(I,j)) + &
(hv(i,J-1)*KH_v_lay(i,J-1) + hv(i,J)*KH_v_lay(i,J))) / &
Kh_t(i,j,k) = (((hu(I-1,j)*KH_u_lay(i-1,j)) + (hu(I,j)*KH_u_lay(I,j))) + &
((hv(i,J-1)*KH_v_lay(i,J-1)) + (hv(i,J)*KH_v_lay(i,J)))) / &
((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + 1.0e-20)
! Use this denominator instead if hu and hv are actual thicknesses rather than a 0/1 mask:
! ((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + h_neglect)
Expand Down Expand Up @@ -916,9 +916,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
drho_dS_u(I) * (S(i,j,k)-S(i,j,k-1)))
drdkR = (drho_dT_u(I) * (T(i+1,j,k)-T(i+1,j,k-1)) + &
drho_dS_u(I) * (S(i+1,j,k)-S(i+1,j,k-1)))
drdkDe_u(I,K) = drdkR * e(i+1,j,K) - drdkL * e(i,j,K)
drdkDe_u(I,K) = (drdkR * e(i+1,j,K)) - (drdkL * e(i,j,K))
elseif (find_work) then ! This is used in pure stacked SW mode
drdkDe_u(I,K) = drdkR * e(i+1,j,K) - drdkL * e(i,j,K)
drdkDe_u(I,K) = (drdkR * e(i+1,j,K)) - (drdkL * e(i,j,K))
endif
if (use_stanley) then
! Correction to the horizontal density gradient due to nonlinearity in
Expand Down Expand Up @@ -950,7 +950,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
! These unnormalized weights have been rearranged to minimize divisions.
wtL = hg2L*(haR*dzaR) ; wtR = hg2R*(haL*dzaL)

drdz = (wtL * drdkL + wtR * drdkR) / (dzaL*wtL + dzaR*wtR)
drdz = ((wtL * drdkL) + (wtR * drdkR)) / ((dzaL*wtL) + (dzaR*wtR))
! The expression for drdz above is mathematically equivalent to:
! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
! ((hg2L/haL) + (hg2R/haR))
Expand All @@ -963,7 +963,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
N2_unlim = drdz*G_rho0
else
N2_unlim = (GV%g_Earth*GV%RZ_to_H) * &
((wtL * drdkL + wtR * drdkR) / (haL*wtL + haR*wtR))
(((wtL * drdkL) + (wtR * drdkR)) / ((haL*wtL) + (haR*wtR)))
endif

dzg2A = dz(i,j,k-1)*dz(i+1,j,k-1) + dz_neglect2
Expand Down Expand Up @@ -1082,10 +1082,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

if (allocated(tv%SpV_avg) .and. (find_work .or. (k > nk_linear)) ) then
Rho_avg = ( ((h(i,j,k) + h(i,j,k-1)) + (h(i+1,j,k) + h(i+1,j,k-1))) + 4.0*hn_2 ) / &
( ((h(i,j,k)+hn_2) * tv%SpV_avg(i,j,k) + (h(i,j,k-1)+hn_2) * tv%SpV_avg(i,j,k-1)) + &
((h(i+1,j,k)+hn_2)*tv%SpV_avg(i+1,j,k) + (h(i+1,j,k-1)+hn_2)*tv%SpV_avg(i+1,j,k-1)) )
( (((h(i,j,k)+hn_2) * tv%SpV_avg(i,j,k)) + ((h(i,j,k-1)+hn_2) * tv%SpV_avg(i,j,k-1))) + &
(((h(i+1,j,k)+hn_2)*tv%SpV_avg(i+1,j,k)) + ((h(i+1,j,k-1)+hn_2)*tv%SpV_avg(i+1,j,k-1))) )
! Use an average density to convert the volume streamfunction estimate into a mass streamfunction.
Z_to_H = (GV%RZ_to_H*Rho_avg)
Z_to_H = GV%RZ_to_H*Rho_avg
else
Z_to_H = GV%Z_to_H
endif
Expand Down Expand Up @@ -1235,9 +1235,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
drho_dS_v(i) * (S(i,j,k)-S(i,j,k-1)))
drdkR = (drho_dT_v(i) * (T(i,j+1,k)-T(i,j+1,k-1)) + &
drho_dS_v(i) * (S(i,j+1,k)-S(i,j+1,k-1)))
drdkDe_v(i,K) = drdkR * e(i,j+1,K) - drdkL * e(i,j,K)
drdkDe_v(i,K) = (drdkR * e(i,j+1,K)) - (drdkL * e(i,j,K))
elseif (find_work) then ! This is used in pure stacked SW mode
drdkDe_v(i,K) = drdkR * e(i,j+1,K) - drdkL * e(i,j,K)
drdkDe_v(i,K) = (drdkR * e(i,j+1,K)) - (drdkL * e(i,j,K))
endif
if (use_stanley) then
! Correction to the horizontal density gradient due to nonlinearity in
Expand Down Expand Up @@ -1271,7 +1271,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
! These unnormalized weights have been rearranged to minimize divisions.
wtL = hg2L*(haR*dzaR) ; wtR = hg2R*(haL*dzaL)

drdz = (wtL * drdkL + wtR * drdkR) / (dzaL*wtL + dzaR*wtR)
drdz = ((wtL * drdkL) + (wtR * drdkR)) / ((dzaL*wtL) + (dzaR*wtR))
! The expression for drdz above is mathematically equivalent to:
! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
! ((hg2L/haL) + (hg2R/haR))
Expand All @@ -1284,7 +1284,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
N2_unlim = drdz*G_rho0
else
N2_unlim = (GV%g_Earth*GV%RZ_to_H) * &
((wtL * drdkL + wtR * drdkR) / (haL*wtL + haR*wtR))
(((wtL * drdkL) + (wtR * drdkR)) / ((haL*wtL) + (haR*wtR)))
endif

dzg2A = dz(i,j,k-1)*dz(i,j+1,k-1) + dz_neglect2
Expand Down Expand Up @@ -1401,10 +1401,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
do i=is,ie
if (allocated(tv%SpV_avg) .and. (find_work .or. (k > nk_linear)) ) then
Rho_avg = ( ((h(i,j,k) + h(i,j,k-1)) + (h(i,j+1,k) + h(i,j+1,k-1))) + 4.0*hn_2 ) / &
( ((h(i,j,k)+hn_2) * tv%SpV_avg(i,j,k) + (h(i,j,k-1)+hn_2) * tv%SpV_avg(i,j,k-1)) + &
((h(i,j+1,k)+hn_2)*tv%SpV_avg(i,j+1,k) + (h(i,j+1,k-1)+hn_2)*tv%SpV_avg(i,j+1,k-1)) )
( (((h(i,j,k)+hn_2) * tv%SpV_avg(i,j,k)) + ((h(i,j,k-1)+hn_2) * tv%SpV_avg(i,j,k-1))) + &
(((h(i,j+1,k)+hn_2)*tv%SpV_avg(i,j+1,k)) + ((h(i,j+1,k-1)+hn_2)*tv%SpV_avg(i,j+1,k-1))) )
! Use an average density to convert the volume streamfunction estimate into a mass streamfunction.
Z_to_H = (GV%RZ_to_H*Rho_avg)
Z_to_H = GV%RZ_to_H*Rho_avg
else
Z_to_H = GV%Z_to_H
endif
Expand Down Expand Up @@ -1510,7 +1510,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
drho_dS_u(I) * (S(i+1,j,1)-S(i,j,1))
if (allocated(tv%SpV_avg)) then
G_scale = GV%H_to_RZ * GV%g_Earth * &
( ((h(i,j,1)+hn_2) * tv%SpV_avg(i,j,1) + (h(i+1,j,1)+hn_2) * tv%SpV_avg(i+1,j,1)) / &
( ( ((h(i,j,1)+hn_2) * tv%SpV_avg(i,j,1)) + ((h(i+1,j,1)+hn_2) * tv%SpV_avg(i+1,j,1)) ) / &
( (h(i,j,1) + h(i+1,j,1)) + 2.0*hn_2 ) )
endif
endif
Expand Down Expand Up @@ -1547,7 +1547,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
drho_dS_v(i) * (S(i,j+1,1)-S(i,j,1))
if (allocated(tv%SpV_avg)) then
G_scale = GV%H_to_RZ * GV%g_Earth * &
( ((h(i,j,1)+hn_2) * tv%SpV_avg(i,j,1) + (h(i,j+1,1)+hn_2) * tv%SpV_avg(i,j+1,1)) / &
( ( ((h(i,j,1)+hn_2) * tv%SpV_avg(i,j,1)) + ((h(i,j+1,1)+hn_2) * tv%SpV_avg(i,j+1,1)) ) / &
( (h(i,j,1) + h(i,j+1,1)) + 2.0*hn_2 ) )
endif
endif
Expand All @@ -1572,22 +1572,23 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
if (CS%MEKE_src_answer_date >= 20240601) then
do j=js,je ; do i=is,ie ; do k=nz,1,-1
PE_release_h = -0.25 * GV%H_to_RZ * &
( (KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + &
Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k)) + &
(Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + &
Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k)) )
( ((KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k)) + &
(Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k))) + &
((Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k)) + &
(Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k))) )
MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h
enddo ; enddo ; enddo
else
do j=js,je ; do i=is,ie ; do k=nz,1,-1
PE_release_h = -0.25 * GV%H_to_RZ * &
(KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k) + &
Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k) + &
Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k) + &
Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k))
((KH_u(I,j,k)*(Slope_x_PE(I,j,k)**2) * hN2_x_PE(I,j,k)) + &
(Kh_u(I-1,j,k)*(Slope_x_PE(I-1,j,k)**2) * hN2_x_PE(I-1,j,k)) + &
(Kh_v(i,J,k)*(Slope_y_PE(i,J,k)**2) * hN2_y_PE(i,J,k)) + &
(Kh_v(i,J-1,k)*(Slope_y_PE(i,J-1,k)**2) * hN2_y_PE(i,J-1,k)))
MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h
enddo ; enddo ; enddo
endif

if (CS%debug) then
call hchksum(MEKE%GM_src, 'MEKE%GM_src', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
call uvchksum("KH_[uv]", Kh_u, Kh_v, G%HI, scale=US%L_to_m**2*US%s_to_T, &
Expand Down Expand Up @@ -2198,11 +2199,11 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
allocate(CS%Kh_eta_u(G%IsdB:G%IedB, G%jsd:G%jed), source=0.)
allocate(CS%Kh_eta_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.)
do j=G%jsc,G%jec ; do I=G%isc-1,G%iec
grid_sp = sqrt((2.0*G%dxCu(I,j)**2 * G%dyCu(I,j)**2) / (G%dxCu(I,j)**2 + G%dyCu(I,j)**2))
grid_sp = sqrt((2.0*G%dxCu(I,j)**2 * G%dyCu(I,j)**2) / ((G%dxCu(I,j)**2) + (G%dyCu(I,j)**2)))
CS%Kh_eta_u(I,j) = G%OBCmaskCu(I,j) * MAX(0.0, CS%Kh_eta_bg + CS%Kh_eta_vel * grid_sp)
enddo ; enddo
do J=G%jsc-1,G%jec ; do i=G%isc,G%iec
grid_sp = sqrt((2.0*G%dxCv(i,J)**2 * G%dyCv(i,J)**2) / (G%dxCv(i,J)**2 + G%dyCv(i,J)**2))
grid_sp = sqrt((2.0*G%dxCv(i,J)**2 * G%dyCv(i,J)**2) / ((G%dxCv(i,J)**2) + (G%dyCv(i,J)**2)))
CS%Kh_eta_v(i,J) = G%OBCmaskCv(i,J) * MAX(0.0, CS%Kh_eta_bg + CS%Kh_eta_vel * grid_sp)
enddo ; enddo
endif
Expand Down

0 comments on commit 49419f7

Please sign in to comment.