Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

+Fix dimensional rescaling with HARMONICS_SAL #503

Merged
merged 2 commits into from
Oct 30, 2023
Merged

+Fix dimensional rescaling with HARMONICS_SAL #503

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b589181
+Fix dimensional rescaling with HARMONICS_SAL
Hallberg-NOAA Oct 13, 2023
ae5bc65
Merge branch 'dev/gfdl' into spherical_harm_scaling_fix
marshallward Oct 30, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions src/core/MOM_PressureForce_FV.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
14 changes: 8 additions & 6 deletions src/parameterizations/lateral/MOM_self_attr_load.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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, "")
Expand Down
42 changes: 28 additions & 14 deletions src/parameterizations/lateral/MOM_spherical_harmonics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)), &
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down