+Fix dimensional rescaling with HARMONICS_SAL

  Corrected dimensional rescaling bugs in the spherical harmonics SAL code.  An
issue with horizontal length scaling was corrected by using G%Rad_Earth_L in
place of G%Rad_Earth in spherical_harmonics_init.  There are new optional
tmp_scale arguments to calc_SAL and spherical_harmonics_forward to allow the
rescaling to be undone before calling the reproducing sums.

  This commit also modifies the call to the reproducing sums in
spherical_harmonics_forward so that all real or imaginary components are
calculated with a single call, which reduces the cost of the SAL calculation
reproducing sums from about 6.7 times the cost with non-reproducing sums to just
5.5 times as much in testing with the tides_025 test case.

  There is also code added to avoid NaNs arising from a square root operating
on a negative argument from a 32-bit integer roll-over when a very large
number of harmonics components (more than 1024 x 1024) are unadvisedly being
used.

  While this commit corrects the dimensional scaling when HARMONICS_SAL is true,
all answers are bitwise identical when no rescaling is used or when the
spherical harmonics SAL is not used.  There are new optional arguments to two
publicly visible interfaces.

src/core/MOM_PressureForce_FV.F90

@@ -319,7 +319,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
       SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref &
                  - max(-G%bathyT(i,j)-G%Z_ref, 0.0)
     enddo ; enddo
-    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
 
     if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq) .or. (.not.CS%tides)) then
       !$OMP parallel do default(shared)
@@ -587,7 +587,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
           SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
         enddo ; enddo
       enddo
-      call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+      call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j)
@@ -618,7 +618,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
           SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
         enddo ; enddo
       enddo
-      call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+      call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
     else
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1

src/core/MOM_PressureForce_Montgomery.F90

@@ -216,7 +216,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
       enddo ; enddo ; enddo
     endif
 
-    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       geopot_bot(i,j) = geopot_bot(i,j) - GV%g_Earth*e_sal(i,j)
@@ -481,7 +481,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
         SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
       enddo ; enddo
     enddo
-    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j)

src/parameterizations/lateral/MOM_self_attr_load.F90

@@ -42,19 +42,21 @@ module MOM_self_attr_load
 !! be changed into bottom pressure anomaly in the future. Note that the SAL calculation applies to all motions
 !! across the spectrum. Tidal-specific methods that assume periodicity, i.e. iterative and read-in SAL, are
 !! stored in MOM_tidal_forcing module.
-subroutine calc_SAL(eta, eta_sal, G, CS)
+subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale)
   type(ocean_grid_type), intent(in)  :: G  !< The ocean's grid structure.
   real, dimension(SZI_(G),SZJ_(G)), intent(in)  :: eta     !< The sea surface height anomaly from
                                                            !! a time-mean geoid [Z ~> m].
   real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_sal !< The sea surface height anomaly from
                                                            !! self-attraction and loading [Z ~> m].
   type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init.
+  real,                   optional, intent(in)  :: tmp_scale !< A rescaling factor to temporarily convert eta
+                                                            !! to MKS units in reproducing sumes [m Z-1 ~> 1]
 
   ! Local variables
   integer :: n, m, l
   integer :: Isq, Ieq, Jsq, Jeq
   integer :: i, j
-  real :: eta_prop
+  real :: eta_prop ! The scalar constant of proportionality between eta and eta_sal [nondim]
 
   call cpu_clock_begin(id_clock_SAL)
 
@@ -69,7 +71,7 @@ subroutine calc_SAL(eta, eta_sal, G, CS)
 
   ! use the spherical harmonics method
   elseif (CS%use_sal_sht) then
-    call spherical_harmonics_forward(G, CS%sht, eta, CS%Snm_Re, CS%Snm_Im, CS%sal_sht_Nd)
+    call spherical_harmonics_forward(G, CS%sht, eta, CS%Snm_Re, CS%Snm_Im, CS%sal_sht_Nd, tmp_scale=tmp_scale)
 
     ! Multiply scaling factors to each mode
     do m = 0,CS%sal_sht_Nd
@@ -119,8 +121,8 @@ subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_Scaling)
   real, dimension(:), intent(out) :: Love_Scaling !< Scaling factors for inverse SHT [nondim]
 
   ! Local variables
-  real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames
-  real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers
+  real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames [nondim]
+  real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers [nondim]
   integer :: n_tot ! Size of the stored Love numbers
   integer :: n, m, l
 
@@ -163,7 +165,7 @@ subroutine SAL_init(G, US, param_file, CS)
 
   logical :: calculate_sal
   logical :: tides, use_tidal_sal_file
-  real :: tide_sal_scalar_value
+  real :: tide_sal_scalar_value ! Scaling SAL factor [nondim]
 
   ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")

src/parameterizations/lateral/MOM_spherical_harmonics.F90

@@ -42,7 +42,7 @@ module MOM_spherical_harmonics
 contains
 
 !> Calculates forward spherical harmonics transforms
-subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
+subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd, tmp_scale)
   type(ocean_grid_type), intent(in)    :: G            !< The ocean's grid structure.
   type(sht_CS),          intent(inout) :: CS           !< Control structure for SHT
   real, dimension(SZI_(G),SZJ_(G)), &
@@ -51,13 +51,20 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
   real,                  intent(out)   :: Snm_Im(:)    !< SHT coefficients for the imaginary modes (sine) [A]
   integer,     optional, intent(in)    :: Nd           !< Maximum degree of the spherical harmonics
                                                        !! overriding ndegree in the CS [nondim]
+  real,        optional, intent(in)    :: tmp_scale    !< A temporary rescaling factor to convert
+                                                       !! var to MKS units during the reproducing
+                                                       !! sums [a A-1 ~> 1]
   ! local variables
-  integer :: Nmax ! Local copy of the maximum degree of the spherical harmonics [nondim]
-  integer :: Ltot ! Local copy of the number of spherical harmonics [nondim]
+  integer :: Nmax ! Local copy of the maximum degree of the spherical harmonics
+  integer :: Ltot ! Local copy of the number of spherical harmonics
   real, dimension(SZI_(G),SZJ_(G)) :: &
     pmn,   & ! Current associated Legendre polynomials of degree n and order m [nondim]
     pmnm1, & ! Associated Legendre polynomials of degree n-1 and order m [nondim]
     pmnm2    ! Associated Legendre polynomials of degree n-2 and order m [nondim]
+  real :: scale ! A rescaling factor to temporarily convert var to MKS units during the
+                ! reproducing sums [a A-1 ~> 1]
+  real :: I_scale ! The inverse of scale [A a-1 ~> 1]
+  real :: sum_tot ! The total of all components output by the reproducing sum in arbitrary units [a]
   integer :: i, j, k
   integer :: is, ie, js, je, isd, ied, jsd, jed
   integer :: m, n, l
@@ -81,21 +88,22 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
   do l=1,Ltot ; Snm_Re(l) = 0.0; Snm_Im(l) = 0.0 ; enddo
 
   if (CS%reprod_sum) then
+    scale = 1.0 ; if (present(tmp_scale)) scale = tmp_scale
     do m=0,Nmax
       l = order2index(m, Nmax)
 
       do j=js,je ; do i=is,ie
-        CS%Snm_Re_raw(i,j,l) = var(i,j) * CS%Pmm(i,j,m+1) * CS%cos_lonT_wtd(i,j,m+1)
-        CS%Snm_Im_raw(i,j,l) = var(i,j) * CS%Pmm(i,j,m+1) * CS%sin_lonT_wtd(i,j,m+1)
+        CS%Snm_Re_raw(i,j,l) = (scale*var(i,j)) * CS%Pmm(i,j,m+1) * CS%cos_lonT_wtd(i,j,m+1)
+        CS%Snm_Im_raw(i,j,l) = (scale*var(i,j)) * CS%Pmm(i,j,m+1) * CS%sin_lonT_wtd(i,j,m+1)
         pmnm2(i,j) = 0.0
         pmnm1(i,j) = CS%Pmm(i,j,m+1)
       enddo ; enddo
 
       do n = m+1, Nmax ; do j=js,je ; do i=is,ie
         pmn(i,j) = &
           CS%a_recur(n+1,m+1) * CS%cos_clatT(i,j) * pmnm1(i,j) - CS%b_recur(n+1,m+1) * pmnm2(i,j)
-        CS%Snm_Re_raw(i,j,l+n-m) = var(i,j) * pmn(i,j) * CS%cos_lonT_wtd(i,j,m+1)
-        CS%Snm_Im_raw(i,j,l+n-m) = var(i,j) * pmn(i,j) * CS%sin_lonT_wtd(i,j,m+1)
+        CS%Snm_Re_raw(i,j,l+n-m) = (scale*var(i,j)) * pmn(i,j) * CS%cos_lonT_wtd(i,j,m+1)
+        CS%Snm_Im_raw(i,j,l+n-m) = (scale*var(i,j)) * pmn(i,j) * CS%sin_lonT_wtd(i,j,m+1)
         pmnm2(i,j) = pmnm1(i,j)
         pmnm1(i,j) = pmn(i,j)
       enddo ; enddo ; enddo
@@ -125,10 +133,15 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
   if (id_clock_sht_global_sum>0) call cpu_clock_begin(id_clock_sht_global_sum)
 
   if (CS%reprod_sum) then
-    do l=1,Ltot
-      Snm_Re(l) = reproducing_sum(CS%Snm_Re_raw(:,:,l))
-      Snm_Im(l) = reproducing_sum(CS%Snm_Im_raw(:,:,l))
-    enddo
+    sum_tot = reproducing_sum(CS%Snm_Re_raw(:,:,1:Ltot), sums=Snm_Re(1:Ltot))
+    sum_tot = reproducing_sum(CS%Snm_Im_raw(:,:,1:Ltot), sums=Snm_Im(1:Ltot))
+    if (scale /= 1.0) then
+      I_scale = 1.0 / scale
+      do l=1,Ltot
+        Snm_Re(l) = I_scale * Snm_Re(l)
+        Snm_Im(l) = I_scale * Snm_Im(l)
+      enddo
+    endif
   else
     call sum_across_PEs(Snm_Re, Ltot)
     call sum_across_PEs(Snm_Im, Ltot)
@@ -240,8 +253,9 @@ subroutine spherical_harmonics_init(G, param_file, CS)
   allocate(CS%a_recur(CS%ndegree+1, CS%ndegree+1)); CS%a_recur(:,:) = 0.0
   allocate(CS%b_recur(CS%ndegree+1, CS%ndegree+1)); CS%b_recur(:,:) = 0.0
   do m=0,CS%ndegree ; do n=m+1,CS%ndegree
+    ! These expressione will give NaNs with 32-bit integers for n > 23170, but this is trapped elsewhere.
     CS%a_recur(n+1,m+1) = sqrt(real((2*n-1) * (2*n+1)) / real((n-m) * (n+m)))
-    CS%b_recur(n+1,m+1) = sqrt(real((2*n+1) * (n+m-1) * (n-m-1)) / real((n-m) * (n+m) * (2*n-3)))
+    CS%b_recur(n+1,m+1) = sqrt((real(2*n+1) * real((n+m-1) * (n-m-1))) / (real((n-m) * (n+m)) * real(2*n-3)))
   enddo ; enddo
 
   ! Calculate complex exponential factors
@@ -253,8 +267,8 @@ subroutine spherical_harmonics_init(G, param_file, CS)
     do j=js,je ; do i=is,ie
       CS%cos_lonT(i,j,m+1)     = cos(real(m) * (G%geolonT(i,j)*RADIAN))
       CS%sin_lonT(i,j,m+1)     = sin(real(m) * (G%geolonT(i,j)*RADIAN))
-      CS%cos_lonT_wtd(i,j,m+1) = CS%cos_lonT(i,j,m+1) * G%areaT(i,j) / G%Rad_Earth**2
-      CS%sin_lonT_wtd(i,j,m+1) = CS%sin_lonT(i,j,m+1) * G%areaT(i,j) / G%Rad_Earth**2
+      CS%cos_lonT_wtd(i,j,m+1) = CS%cos_lonT(i,j,m+1) * G%areaT(i,j) / G%Rad_Earth_L**2
+      CS%sin_lonT_wtd(i,j,m+1) = CS%sin_lonT(i,j,m+1) * G%areaT(i,j) / G%Rad_Earth_L**2
     enddo ; enddo
   enddo