+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.
Hallberg-NOAA · Oct 13, 2023 · 0295d4e · 0295d4e
1 parent 89506fa
commit 0295d4e
Show file tree

Hide file tree

Showing 4 changed files with 41 additions and 21 deletions.
diff --git a/src/core/MOM_PressureForce_FV.F90 b/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

diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/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)

diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90
@@ -42,13 +42,15 @@ 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
@@ -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

diff --git a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 b/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,13 @@ 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
+    ! This expression 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)))
+    if (n <= 1024) then
+      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)))
+    else ! This branch for large n avoids NaNs due to 32-bit integer roll-over.
+      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)))
+    endif
   enddo ; enddo
 
   ! Calculate complex exponential factors
@@ -253,8 +271,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