(*)Parenthesize initialization squares for FMAs

  Added parentheses to 20 sums of squares of x- and y- distances or velocity
components used for initialization in 8 modules to give 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.

src/initialization/MOM_state_initialization.F90

@@ -1595,7 +1595,7 @@ real function my_psi(ig,jg)
 
     x = 2.0*(G%geoLonBu(ig,jg)-G%west_lon) / G%len_lon - 1.0  ! -1<x<1
     y = 2.0*(G%geoLatBu(ig,jg)-G%south_lat) / G%len_lat - 1.0 ! -1<y<1
-    r = sqrt( x**2 + y**2 ) ! Circular stream function is a function of radius only
+    r = sqrt( (x**2) + (y**2) ) ! Circular stream function is a function of radius only
     r = min(1.0, r) ! Flatten stream function in corners of box
     my_psi = 0.5*(1.0 - cos(dpi*r))
     my_psi = my_psi * (circular_max_u * G%US%m_to_L*G%len_lon*1e3 / dpi) ! len_lon is in km

src/tracer/advection_test_tracer.F90

@@ -226,13 +226,13 @@ subroutine initialize_advection_test_tracer(restart, day, G, GV, h,diag, OBC, CS
       do j=js,je ; do i=is,ie
         locx = abs(G%geoLonT(i,j)-CS%x_origin)/CS%x_width
         locy = abs(G%geoLatT(i,j)-CS%y_origin)/CS%y_width
-        if (locx**2+locy**2<=1.0) CS%tr(i,j,k,m) = 1.0
+        if ((locx**2) + (locy**2) <= 1.0) CS%tr(i,j,k,m) = 1.0
       enddo ; enddo
       k=5 ! Cut cylinder
       do j=js,je ; do i=is,ie
         locx = (G%geoLonT(i,j)-CS%x_origin)/CS%x_width
         locy = (G%geoLatT(i,j)-CS%y_origin)/CS%y_width
-        if (locx**2+locy**2<=1.0) CS%tr(i,j,k,m) = 1.0
+        if ((locx**2) + (locy**2) <= 1.0) CS%tr(i,j,k,m) = 1.0
         if (locx>0.0 .and. abs(locy)<0.2) CS%tr(i,j,k,m) = 0.0
       enddo ; enddo
 

src/user/Idealized_Hurricane.F90

@@ -368,7 +368,7 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx
 
   ! Implementing Holland (1980) parameteric wind profile
 
-  radius = SQRT(XX**2 + YY**2)
+  radius = SQRT((XX**2) + (YY**2))
 
   !/ BGR
   ! rkm - r converted to km for Holland prof.
@@ -451,7 +451,7 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx
   dV = U10*cos(Adir-Alph) - Vocn + V_TS
 
   !  Use a simple drag coefficient as a function of U10 (from Sullivan et al., 2010)
-  du10 = sqrt(du**2+dv**2)
+  du10 = sqrt((du**2) + (dv**2))
   if (dU10 < 11.0*US%m_s_to_L_T) then
     Cd = 1.2e-3
   elseif (dU10 < 20.0*US%m_s_to_L_T) then
@@ -465,8 +465,8 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx
   endif
 
   ! Compute stress vector
-  TX = US%L_to_Z * CS%rho_a * Cd * sqrt(dU**2 + dV**2) * dU
-  TY = US%L_to_Z * CS%rho_a * Cd * sqrt(dU**2 + dV**2) * dV
+  TX = US%L_to_Z * CS%rho_a * Cd * sqrt((dU**2) + (dV**2)) * dU
+  TY = US%L_to_Z * CS%rho_a * Cd * sqrt((dU**2) + (dV**2)) * dV
 
 end subroutine idealized_hurricane_wind_profile
 
@@ -541,7 +541,7 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C
   !/ BR
   ! Calculate x position as a function of time.
   xx = US%s_to_T*( t0 - time_type_to_real(day)) * CS%hurr_translation_spd * cos(transdir)
-  rad = sqrt(xx**2 + CS%dy_from_center**2)
+  rad = sqrt((xx**2) + (CS%dy_from_center**2))
   !/ BR
   ! rkm - rad converted to km for Holland prof.
   !       used in km due to error, correct implementation should
@@ -619,7 +619,7 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C
     !BR
     !  Add a simple drag coefficient as a function of U10 |
     !/----------------------------------------------------|
-    du10 = sqrt(du**2+dv**2)
+    du10 = sqrt((du**2) + (dv**2))
     if (dU10 < 11.0*US%m_s_to_L_T) then
       Cd = 1.2e-3
     elseif (dU10 < 20.0*US%m_s_to_L_T) then
@@ -641,7 +641,7 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C
     Vocn = 0. ! sfc_state%v(i,J)
     dU = U10*sin(Adir-pie-Alph) - Uocn + U_TS
     dV = U10*cos(Adir-Alph) - Vocn + V_TS
-    du10=sqrt(du**2+dv**2)
+    du10 = sqrt((du**2) + (dv**2))
     if (dU10 < 11.0*US%m_s_to_L_T) then
       Cd = 1.2e-3
     elseif (dU10 < 20.0*US%m_s_to_L_T) then

src/user/Neverworld_initialization.F90

@@ -157,7 +157,7 @@ real function dist_line_fixed_x(x, y, x0, y0, y1)
   dx = x - x0
   yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1
   dy = y - yr ! =0 within y0<y<y1, =y0-y for y<y0, =y-y1 for y>y1
-  dist_line_fixed_x = sqrt( dx*dx + dy*dy )
+  dist_line_fixed_x = sqrt( (dx*dx) + (dy*dy) )
 end function dist_line_fixed_x
 
 !> Distance between points x,y and a line segment (x0,y0) and (x1,y0).
@@ -229,7 +229,7 @@ real function circ_ridge(lon, lat, lon0, lat0, ring_radius, ring_thickness, ridg
   real :: r ! A relative position [degrees]
   real :: frac_ht ! The fractional height of the topography [nondim]
 
-  r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point
+  r = sqrt( ((lon - lon0)**2) + ((lat - lat0)**2) ) ! Pseudo-distance from a point
   r = abs( r - ring_radius) ! Pseudo-distance from a circle
   frac_ht = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height
   circ_ridge = 1. - frac_ht ! Fractional depths (1-frac_ridge_height) .. 1
@@ -292,8 +292,8 @@ subroutine Neverworld_initialize_thickness(h, depth_tot, G, GV, US, param_file,
       h(i,j,k) = e0(k) - e_interface ! Nominal thickness
       x = (G%geoLonT(i,j)-G%west_lon)/G%len_lon
       y = (G%geoLatT(i,j)-G%south_lat)/G%len_lat
-      r1 = sqrt((x-0.7)**2+(y-0.2)**2)
-      r2 = sqrt((x-0.3)**2+(y-0.25)**2)
+      r1 = sqrt(((x-0.7)**2) + ((y-0.2)**2))
+      r2 = sqrt(((x-0.3)**2) + ((y-0.25)**2))
       h(i,j,k) = h(i,j,k) + pert_amp * (e0(k) - e0(nz+1)) * &
                             (spike(r1,0.15)-spike(r2,0.15)) ! Prescribed perturbation
       if (h_noise /= 0.) then

src/user/SCM_CVMix_tests.F90

@@ -217,7 +217,7 @@ subroutine SCM_CVMix_tests_wind_forcing(sfc_state, forces, day, G, US, CS)
   enddo ; enddo
   call pass_vector(forces%taux, forces%tauy, G%Domain, To_All)
 
-  mag_tau = sqrt(CS%tau_x*CS%tau_x + CS%tau_y*CS%tau_y)
+  mag_tau = sqrt((CS%tau_x*CS%tau_x) + (CS%tau_y*CS%tau_y))
   if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
     forces%ustar(i,j) = sqrt( US%L_to_Z * mag_tau / CS%Rho0 )
   enddo ; enddo ; endif

src/user/basin_builder.F90

@@ -208,7 +208,7 @@ real function dist_line_fixed_x(x, y, x0, y0, y1)
   dx = x - x0
   yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1
   dy = y - yr ! =0 within y0<y<y1, =y0-y for y<y0, =y-y1 for y>y1
-  dist_line_fixed_x = sqrt( dx*dx + dy*dy )
+  dist_line_fixed_x = sqrt( (dx*dx) + (dy*dy) )
 end function dist_line_fixed_x
 
 !> Distance between points x,y and a line segment (x0,y0) and (x1,y0).
@@ -310,7 +310,7 @@ real function circ_conic_ridge(lon, lat, lon0, lat0, ring_radius, ring_thickness
   real :: r  ! A relative position [degrees]
   real :: frac_ht ! The fractional height of the topography [nondim]
 
-  r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point
+  r = sqrt( ((lon - lon0)**2) + ((lat - lat0)**2) ) ! Pseudo-distance from a point
   r = abs( r - ring_radius) ! Pseudo-distance from a circle
   frac_ht = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height
   circ_conic_ridge = 1. - frac_ht ! nondim depths (1-frac_ridge_height) .. 1
@@ -329,7 +329,7 @@ real function circ_scurve_ridge(lon, lat, lon0, lat0, ring_radius, ring_thicknes
   real :: s  ! A function of the normalized position [nondim]
   real :: frac_ht ! The fractional height of the topography [nondim]
 
-  r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point
+  r = sqrt( ((lon - lon0)**2) + ((lat - lat0)**2) ) ! Pseudo-distance from a point
   r = abs( r - ring_radius) ! Pseudo-distance from a circle
   s = 1. - scurve(r, 0., ring_thickness) ! 0 .. 1
   frac_ht = s * ridge_height ! 0 .. frac_ridge_height

src/user/circle_obcs_initialization.F90

@@ -102,7 +102,7 @@ subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, US, param_file,
   latC = G%south_lat + 0.5*G%len_lat
   lonC = G%west_lon + 0.5*G%len_lon + xOffset
   do j=js,je ; do i=is,ie
-    rad = sqrt((G%geoLonT(i,j)-lonC)**2+(G%geoLatT(i,j)-latC)**2)/(diskrad)
+    rad = sqrt(((G%geoLonT(i,j)-lonC)**2) + ((G%geoLatT(i,j)-latC)**2)) / diskrad
     ! if (rad <= 6.*diskrad) h(i,j,k) = h(i,j,k)+10.0*exp( -0.5*( rad**2 ) )
     rad = min( rad, 1. ) ! Flatten outside radius of diskrad
     rad = rad*(2.*asin(1.)) ! Map 0-1 to 0-pi

src/user/seamount_initialization.F90

@@ -72,7 +72,7 @@ subroutine seamount_initialize_topography( D, G, param_file, max_depth )
     ! Compute normalized zonal coordinates (x,y=0 at center of domain)
     x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5
     y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5
-    D(i,j) = G%max_depth * ( 1.0 - delta * exp(-(rLx*x)**2 -(rLy*y)**2) )
+    D(i,j) = G%max_depth * ( 1.0 - delta * exp(-((rLx*x)**2) - ((rLy*y)**2)) )
   enddo ; enddo
 
 end subroutine seamount_initialize_topography