From b58918192da3ea42627509b0561e105718fb5b6a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 13 Oct 2023 08:34:37 -0400 Subject: [PATCH] +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 | 6 +-- src/core/MOM_PressureForce_Montgomery.F90 | 4 +- .../lateral/MOM_self_attr_load.F90 | 14 ++++--- .../lateral/MOM_spherical_harmonics.F90 | 42 ++++++++++++------- 4 files changed, 41 insertions(+), 25 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 64df200f31..5fb3ade634 100644 --- 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 index 3de713c801..6d982bc7e3 100644 --- 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 index 20d239eb53..7f7215c9d8 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/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, "") diff --git a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 b/src/parameterizations/lateral/MOM_spherical_harmonics.F90 index 2a72d26a20..26258e6b8e 100644 --- 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,12 +88,13 @@ 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 @@ -94,8 +102,8 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd) 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