diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90
index 42e6514ab9..cccfa1c2e4 100644
--- a/src/core/MOM_PressureForce_FV.F90
+++ b/src/core/MOM_PressureForce_FV.F90
@@ -39,8 +39,10 @@ module MOM_PressureForce_FV
!> Finite volume pressure gradient control structure
type, public :: PressureForce_FV_CS ; private
logical :: initialized = .false. !< True if this control structure has been initialized.
- logical :: calculate_SAL !< If true, calculate self-attraction and loading.
- logical :: tides !< If true, apply tidal momentum forcing.
+ logical :: calculate_SAL = .false. !< If true, calculate self-attraction and loading.
+ logical :: sal_use_bpa = .false. !< If true, use bottom pressure anomaly instead of SSH
+ !! to calculate SAL.
+ logical :: tides = .false. !< If true, apply tidal momentum forcing.
real :: Rho0 !< The density used in the Boussinesq
!! approximation [R ~> kg m-3].
real :: GFS_scale !< A scaling of the surface pressure gradients to
@@ -80,10 +82,12 @@ module MOM_PressureForce_FV
!! equation of state is 0 to account for the displacement of the sea
!! surface including adjustments for atmospheric or sea-ice pressure.
logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF
- integer :: tides_answer_date !< Recover old answers with tides in Boussinesq mode
+ logical :: bq_sal_tides = .false. !< If true, use an alternative method for SAL and tides
+ !! in Boussinesq mode
+ integer :: tides_answer_date = 99991231 !< Recover old answers with tides
integer :: id_e_tide = -1 !< Diagnostic identifier
- integer :: id_e_tide_eq = -1 !< Diagnostic identifier
- integer :: id_e_tide_sal = -1 !< Diagnostic identifier
+ integer :: id_e_tidal_eq = -1 !< Diagnostic identifier
+ integer :: id_e_tidal_sal = -1 !< Diagnostic identifier
integer :: id_e_sal = -1 !< Diagnostic identifier
integer :: id_rho_pgf = -1 !< Diagnostic identifier
integer :: id_rho_stanley_pgf = -1 !< Diagnostic identifier
@@ -141,12 +145,16 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
! the pressure anomaly at the top of the layer [R L4 T-4 ~> Pa m2 s-2].
real, dimension(SZI_(G),SZJ_(G)) :: &
dp, & ! The (positive) change in pressure across a layer [R L2 T-2 ~> Pa].
- SSH, & ! The sea surface height anomaly, in depth units [Z ~> m].
+ SSH, & ! Sea surfae height anomaly for self-attraction and loading. Used if
+ ! CALCULATE_SAL is True and SAL_USE_BPA is False [Z ~> m].
+ pbot, & ! Total bottom pressure for self-attraction and loading. Used if
+ ! CALCULATE_SAL is True and SAL_USE_BPA is True [R L2 T-2 ~> Pa].
e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m].
- e_tide_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m].
- e_tide_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
+ e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m].
+ e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
! specific to tides [Z ~> m].
- e_sal_tide, & ! The summation of self-attraction and loading and tidal forcing [Z ~> m].
+ e_sal_and_tide, & ! The summation of self-attraction and loading and tidal forcing, used for recovering
+ ! old answers only [Z ~> m].
dM ! The barotropic adjustment to the Montgomery potential to
! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
@@ -399,15 +407,24 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
enddo ; enddo
enddo
- ! Calculate and add the self-attraction and loading geopotential anomaly.
+ ! Calculate and add self-attraction and loading (SAL) geopotential height anomaly to interface height.
if (CS%calculate_SAL) then
- !$OMP parallel do default(shared)
- do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- SSH(i,j) = (za(i,j,1) - 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, tmp_scale=US%Z_to_m)
+ if (CS%sal_use_bpa) then
+ !$OMP parallel do default(shared)
+ do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ pbot(i,j) = p(i,j,nz+1)
+ enddo ; enddo
+ call calc_SAL(pbot, 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
+ SSH(i,j) = (za(i,j,1) - 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, tmp_scale=US%Z_to_m)
+ endif
+ ! This gives new answers after the change of separating SAL from tidal forcing module.
if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq) .or. (.not.CS%tides)) then
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
@@ -416,21 +433,21 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
endif
endif
- ! Calculate and add the tidal geopotential anomaly.
+ ! Calculate and add tidal geopotential height anomaly to interface height.
if (CS%tides) then
if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq)) then
- call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp)
+ call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- za(i,j,1) = za(i,j,1) - GV%g_Earth * (e_tide_eq(i,j) + e_tide_sal(i,j))
+ za(i,j,1) = za(i,j,1) - GV%g_Earth * (e_tidal_eq(i,j) + e_tidal_sal(i,j))
enddo ; enddo
else ! This block recreates older answers with tides.
if (.not.CS%calculate_SAL) e_sal(:,:) = 0.0
- call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_tide, e_tide_eq, e_tide_sal, &
+ call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, &
G, US, CS%tides_CSp)
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- za(i,j,1) = za(i,j,1) - GV%g_Earth * e_sal_tide(i,j)
+ za(i,j,1) = za(i,j,1) - GV%g_Earth * e_sal_and_tide(i,j)
enddo ; enddo
endif
endif
@@ -808,10 +825,15 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
! To be consistent with old runs, tidal forcing diagnostic also includes total SAL.
! New diagnostics are given for each individual field.
- if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal+e_tide_eq+e_tide_sal, CS%diag)
+ if (CS%id_e_tide>0) then
+ if (CS%tides_answer_date>20230630) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ e_sal_and_tide(i,j) = e_sal(i,j) + e_tidal_eq(i,j) + e_tidal_sal(i,j)
+ enddo ; enddo ; endif
+ call post_data(CS%id_e_tide, e_sal_and_tide, CS%diag)
+ endif
if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag)
- if (CS%id_e_tide_eq>0) call post_data(CS%id_e_tide_eq, e_tide_eq, CS%diag)
- if (CS%id_e_tide_sal>0) call post_data(CS%id_e_tide_sal, e_tide_sal, CS%diag)
+ if (CS%id_e_tidal_eq>0) call post_data(CS%id_e_tidal_eq, e_tidal_eq, CS%diag)
+ if (CS%id_e_tidal_sal>0) call post_data(CS%id_e_tidal_sal, e_tidal_sal, CS%diag)
if (CS%id_MassWt_u>0) call post_data(CS%id_MassWt_u, MassWt_u, CS%diag)
if (CS%id_MassWt_v>0) call post_data(CS%id_MassWt_v, MassWt_v, CS%diag)
@@ -846,14 +868,17 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
! Local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! Interface height in depth units [Z ~> m].
real, dimension(SZI_(G),SZJ_(G)) :: &
- e_sal_tide, & ! The summation of self-attraction and loading and tidal forcing [Z ~> m].
+ e_sal_and_tide, & ! The summation of self-attraction and loading and tidal forcing [Z ~> m].
e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m].
- e_tide_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources
+ e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources
! [Z ~> m].
- e_tide_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
+ e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
! specific to tides [Z ~> m].
Z_0p, & ! The height at which the pressure used in the equation of state is 0 [Z ~> m]
- SSH, & ! The sea surface height anomaly, in depth units [Z ~> m].
+ SSH, & ! Sea surfae height anomaly for self-attraction and loading. Used if
+ ! CALCULATE_SAL is True and SAL_USE_BPA is False [Z ~> m].
+ pbot, & ! Total bottom pressure for self-attraction and loading. Used if
+ ! CALCULATE_SAL is True and SAL_USE_BPA is True [R L2 T-2 ~> Pa].
dM ! The barotropic adjustment to the Montgomery potential to
! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
real, dimension(SZI_(G)) :: &
@@ -950,7 +975,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m].
real :: I_Rho0 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1].
- real :: G_Rho0 ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
+ real :: G_Rho0 ! G_Earth / Rho_0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1]
real :: rho_ref ! The reference density [R ~> kg m-3].
real :: dz_neglect ! A minimal thickness [Z ~> m], like e.
@@ -1003,77 +1028,51 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0
endif
- if (CS%tides_answer_date>20230630) then
- do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- e(i,j,nz+1) = -G%bathyT(i,j)
- enddo ; enddo
+ do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ e(i,j,nz+1) = -G%bathyT(i,j)
+ enddo ; enddo
- ! Calculate and add the self-attraction and loading geopotential anomaly.
- if (CS%calculate_SAL) then
- ! Determine the surface height anomaly for calculating self attraction
- ! and loading. This should really be based on bottom pressure anomalies,
- ! but that is not yet implemented, and the current form is correct for
- ! barotropic tides.
- !$OMP parallel do default(shared)
- do j=Jsq,Jeq+1
- do i=Isq,Ieq+1
- SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0)
- enddo
- do k=1,nz ; do i=Isq,Ieq+1
- SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
- enddo ; enddo
+ ! The following two if-blocks are used to recover old answers for self-attraction and loading
+ ! (SAL) and tides only. The old algorithm moves interface heights before density calculations,
+ ! and therefore is incorrect without SSH_IN_EOS_PRESSURE_FOR_PGF=True (added in August 2024).
+ ! See the code right after Pa calculation loop for the new algorithm.
+
+ ! Calculate and add SAL geopotential anomaly to interface height (old answers)
+ if (CS%calculate_SAL .and. CS%tides_answer_date<=20250131) then
+ !$OMP parallel do default(shared)
+ do j=Jsq,Jeq+1
+ do i=Isq,Ieq+1
+ SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0)
enddo
- 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)
+ do k=1,nz ; do i=Isq,Ieq+1
+ SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
enddo ; enddo
- endif
+ enddo
+ call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
- ! Calculate and add the tidal geopotential anomaly.
- if (CS%tides) then
- call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp)
+ if (CS%tides_answer_date>20230630) then ! answers_date between [20230701, 20250131]
!$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_tide_eq(i,j) + e_tide_sal(i,j))
- enddo ; enddo
- endif
- else ! Old answers
- ! Calculate and add the self-attraction and loading geopotential anomaly.
- if (CS%calculate_SAL) then
- ! Determine the surface height anomaly for calculating self attraction
- ! and loading. This should really be based on bottom pressure anomalies,
- ! but that is not yet implemented, and the current form is correct for
- ! barotropic tides.
- !$OMP parallel do default(shared)
- do j=Jsq,Jeq+1
- do i=Isq,Ieq+1
- SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0)
- enddo
- do k=1,nz ; do i=Isq,Ieq+1
- 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, tmp_scale=US%Z_to_m)
- else
- !$OMP parallel do default(shared)
- do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- e_sal(i,j) = 0.0
+ e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j)
enddo ; enddo
endif
+ endif
- ! Calculate and add the tidal geopotential anomaly.
- if (CS%tides) then
- call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_tide, e_tide_eq, e_tide_sal, &
- G, US, CS%tides_CSp)
- !$OMP parallel do default(shared)
+ ! Calculate and add tidal geopotential anomaly to interface height (old answers)
+ if (CS%tides .and. CS%tides_answer_date<=20250131) then
+ if (CS%tides_answer_date>20230630) then ! answers_date between [20230701, 20250131]
+ call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
+ !$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- e(i,j,nz+1) = -(G%bathyT(i,j) + e_sal_tide(i,j))
+ e(i,j,nz+1) = e(i,j,nz+1) - (e_tidal_eq(i,j) + e_tidal_sal(i,j))
enddo ; enddo
- else
+ else ! answers_date before 20230701
+ if (.not.CS%calculate_SAL) e_sal(:,:) = 0.0
+ call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, &
+ G, US, CS%tides_CSp)
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- e(i,j,nz+1) = -(G%bathyT(i,j) + e_sal(i,j))
+ e(i,j,nz+1) = e(i,j,nz+1) - e_sal_and_tide(i,j)
enddo ; enddo
endif
endif
@@ -1223,6 +1222,42 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
enddo ; enddo
enddo
+ ! Calculate and add SAL geopotential anomaly to interface height (new answers)
+ if (CS%calculate_SAL .and. CS%tides_answer_date>20250131) then
+ if (CS%sal_use_bpa) then
+ !$OMP parallel do default(shared)
+ do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ pbot(i,j) = pa(i,j,nz+1) - GxRho * e(i,j,nz+1)
+ enddo ; enddo
+ call calc_SAL(pbot, 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
+ SSH(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topography above sea level
+ enddo ; enddo
+ call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
+ endif
+ if (.not.CS%bq_sal_tides) then ; do K=1,nz+1
+ !$OMP parallel do default(shared)
+ do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ e(i,j,K) = e(i,j,K) - e_sal(i,j)
+ pa(i,j,K) = pa(i,j,K) - GxRho * e_sal(i,j)
+ enddo ; enddo
+ enddo ; endif
+ endif
+
+ ! Calculate and add tidal geopotential anomaly to interface height (new answers)
+ if (CS%tides .and. CS%tides_answer_date>20250131) then
+ call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
+ if (.not.CS%bq_sal_tides) then ; do K=1,nz+1
+ !$OMP parallel do default(shared)
+ do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ e(i,j,K) = e(i,j,K) - (e_tidal_eq(i,j) + e_tidal_sal(i,j))
+ pa(i,j,K) = pa(i,j,K) - GxRho * (e_tidal_eq(i,j) + e_tidal_sal(i,j))
+ enddo ; enddo
+ enddo ; endif
+ endif
+
if (CS%correction_intxpa .or. CS%reset_intxpa_integral) then
! Determine surface temperature and salinity for use in the pressure gradient corrections
if (use_ALE .and. (CS%Recon_Scheme > 0)) then
@@ -1599,6 +1634,34 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
((h(i,j,k) + h(i,j+1,k)) + h_neglect))
enddo ; enddo ; enddo
+ ! Calculate SAL geopotential anomaly and add its gradient to pressure gradient force
+ if (CS%calculate_SAL .and. CS%tides_answer_date>20230630 .and. CS%bq_sal_tides) then
+ !$OMP parallel do default(shared)
+ do k=1,nz
+ do j=js,je ; do I=Isq,Ieq
+ PFu(I,j,k) = PFu(I,j,k) + (e_sal(i+1,j) - e_sal(i,j)) * GV%g_Earth * G%IdxCu(I,j)
+ enddo ; enddo
+ do J=Jsq,Jeq ; do i=is,ie
+ PFv(i,J,k) = PFv(i,J,k) + (e_sal(i,j+1) - e_sal(i,j)) * GV%g_Earth * G%IdyCv(i,J)
+ enddo ; enddo
+ enddo
+ endif
+
+ ! Calculate tidal geopotential anomaly and add its gradient to pressure gradient force
+ if (CS%tides .and. CS%tides_answer_date>20230630 .and. CS%bq_sal_tides) then
+ !$OMP parallel do default(shared)
+ do k=1,nz
+ do j=js,je ; do I=Isq,Ieq
+ PFu(I,j,k) = PFu(I,j,k) + ((e_tidal_eq(i+1,j) + e_tidal_sal(i+1,j)) &
+ - (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * GV%g_Earth * G%IdxCu(I,j)
+ enddo ; enddo
+ do J=Jsq,Jeq ; do i=is,ie
+ PFv(i,J,k) = PFv(i,J,k) + ((e_tidal_eq(i,j+1) + e_tidal_sal(i,j+1)) &
+ - (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * GV%g_Earth * G%IdyCv(i,J)
+ enddo ; enddo
+ enddo
+ endif
+
if (CS%GFS_scale < 1.0) then
! Adjust the Montgomery potential to make this a reduced gravity model.
if (use_EOS) then
@@ -1640,36 +1703,29 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
if (present(eta)) then
! eta is the sea surface height relative to a time-invariant geoid, for comparison with
! what is used for eta in btstep. See how e was calculated about 200 lines above.
- if (CS%tides_answer_date>20230630) then
- !$OMP parallel do default(shared)
- do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- eta(i,j) = e(i,j,1)*GV%Z_to_H
- enddo ; enddo
- if (CS%tides) then
- !$OMP parallel do default(shared)
- do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- eta(i,j) = eta(i,j) + (e_tide_eq(i,j)+e_tide_sal(i,j))*GV%Z_to_H
- enddo ; enddo
- endif
- if (CS%calculate_SAL) then
- !$OMP parallel do default(shared)
- do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H
- enddo ; enddo
- endif
- else ! Old answers
- if (CS%tides) then
+ !$OMP parallel do default(shared)
+ do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ eta(i,j) = e(i,j,1)*GV%Z_to_H
+ enddo ; enddo
+ if (CS%tides .and. (.not.CS%bq_sal_tides)) then
+ if (CS%tides_answer_date>20230630) then
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- eta(i,j) = e(i,j,1)*GV%Z_to_H + (e_sal_tide(i,j))*GV%Z_to_H
+ eta(i,j) = eta(i,j) + (e_tidal_eq(i,j)+e_tidal_sal(i,j))*GV%Z_to_H
enddo ; enddo
else
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- eta(i,j) = (e(i,j,1) + e_sal(i,j))*GV%Z_to_H
+ eta(i,j) = eta(i,j) + e_sal_and_tide(i,j)*GV%Z_to_H
enddo ; enddo
endif
endif
+ if (CS%calculate_SAL .and. (CS%tides_answer_date>20230630) .and. (.not.CS%bq_sal_tides)) then
+ !$OMP parallel do default(shared)
+ do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H
+ enddo ; enddo
+ endif
endif
if (CS%use_stanley_pgf) then
@@ -1711,10 +1767,15 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
! To be consistent with old runs, tidal forcing diagnostic also includes total SAL.
! New diagnostics are given for each individual field.
- if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal_tide, CS%diag)
+ if (CS%id_e_tide>0) then
+ if (CS%tides_answer_date>20230630) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ e_sal_and_tide(i,j) = e_sal(i,j) + e_tidal_eq(i,j) + e_tidal_sal(i,j)
+ enddo ; enddo ; endif
+ call post_data(CS%id_e_tide, e_sal_and_tide, CS%diag)
+ endif
if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag)
- if (CS%id_e_tide_eq>0) call post_data(CS%id_e_tide_eq, e_tide_eq, CS%diag)
- if (CS%id_e_tide_sal>0) call post_data(CS%id_e_tide_sal, e_tide_sal, CS%diag)
+ if (CS%id_e_tidal_eq>0) call post_data(CS%id_e_tidal_eq, e_tidal_eq, CS%diag)
+ if (CS%id_e_tidal_sal>0) call post_data(CS%id_e_tidal_sal, e_tidal_sal, CS%diag)
if (CS%id_MassWt_u>0) call post_data(CS%id_MassWt_u, MassWt_u, CS%diag)
if (CS%id_MassWt_v>0) call post_data(CS%id_MassWt_v, MassWt_v, CS%diag)
@@ -1769,19 +1830,29 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "TIDES", CS%tides, &
"If true, apply tidal momentum forcing.", default=.false.)
- if (CS%tides) then
- call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
- "This sets the default value for the various _ANSWER_DATE parameters.", &
- default=99991231)
- call get_param(param_file, mdl, "TIDES_ANSWER_DATE", CS%tides_answer_date, &
- "The vintage of self-attraction and loading (SAL) and tidal forcing calculations in "//&
- "Boussinesq mode. Values below 20230701 recover the old answers in which the SAL is "//&
- "part of the tidal forcing calculation. The change is due to a reordered summation "//&
- "and the difference is only at bit level.", default=20230630)
- endif
+ call get_param(param_file, '', "DEFAULT_ANSWER_DATE", default_answer_date, default=99991231)
+ if (CS%tides) &
+ call get_param(param_file, mdl, "TIDES_ANSWER_DATE", CS%tides_answer_date, "The vintage of "//&
+ "self-attraction and loading (SAL) and tidal forcing calculations. Setting "//&
+ "dates before 20230701 recovers old answers (Boussinesq and non-Boussinesq "//&
+ "modes) when SAL is part of the tidal forcing calculation. The answer "//&
+ "difference is only at bit level and due to a reordered summation. Setting "//&
+ "dates before 20250201 recovers answers (Boussinesq mode) that interface "//&
+ "heights are modified before pressure force integrals are calculated.", &
+ default=20230630, do_not_log=(.not.CS%tides))
call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, &
"If true, calculate self-attraction and loading.", default=CS%tides)
-
+ if (CS%calculate_SAL) &
+ call get_param(param_file, '', "SAL_USE_BPA", CS%sal_use_bpa, default=.false., &
+ do_not_log=.true.)
+ if ((CS%tides .or. CS%calculate_SAL) .and. GV%Boussinesq) &
+ call get_param(param_file, mdl, "BOUSSINESQ_SAL_TIDES", CS%bq_sal_tides, "If true, "//&
+ "in Boussinesq mode, use an alternative method to include self-attraction "//&
+ "and loading (SAL) and tidal forcings in pressure gradient, in which their "//&
+ "gradients are calculated separately, instead of adding geopotential "//&
+ "anomalies as corrections to the interface height. This alternative method "//&
+ "elimates a baroclinic component of the SAL and tidal forcings.", &
+ default=.false.)
call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &
"If true, Temperature and salinity are used as state variables.", &
default=.true., do_not_log=.true.)
@@ -1796,6 +1867,9 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
"pressure used in the equation of state calculations for the Boussinesq pressure "//&
"gradient forces, including adjustments for atmospheric or sea-ice pressure.", &
default=.false., do_not_log=.not.GV%Boussinesq)
+ if (CS%tides .and. CS%tides_answer_date<=20250131 .and. CS%use_SSH_in_Z0p) &
+ call MOM_error(FATAL, trim(mdl) // ", PressureForce_FV_init: SSH_IN_EOS_PRESSURE_FOR_PGF "//&
+ "needs to be FALSE to recover tide answers before 20250131.")
call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, &
"If True, use the ALE algorithm (regridding/remapping). "//&
@@ -1880,9 +1954,9 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
if (CS%tides) then
CS%id_e_tide = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, Time, &
'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=US%Z_to_m)
- CS%id_e_tide_eq = register_diag_field('ocean_model', 'e_tide_eq', diag%axesT1, Time, &
+ CS%id_e_tidal_eq = register_diag_field('ocean_model', 'e_tide_eq', diag%axesT1, Time, &
'Equilibrium tides height anomaly', 'meter', conversion=US%Z_to_m)
- CS%id_e_tide_sal = register_diag_field('ocean_model', 'e_tide_sal', diag%axesT1, Time, &
+ CS%id_e_tidal_sal = register_diag_field('ocean_model', 'e_tide_sal', diag%axesT1, Time, &
'Read-in tidal self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m)
endif
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index f602b65240..2978bd46b9 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -1506,7 +1506,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p
call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
cont_stencil = continuity_stencil(CS%continuity_CSp)
call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
- if (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp)
+ if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp)
if (CS%use_tides) then
call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp)
HA_CSp => CS%HA_CSp
diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90
index 87e46795b5..31ee923aae 100644
--- a/src/core/MOM_dynamics_split_RK2b.F90
+++ b/src/core/MOM_dynamics_split_RK2b.F90
@@ -1419,7 +1419,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US,
call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
cont_stencil = continuity_stencil(CS%continuity_CSp)
call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
- if (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp)
+ if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp)
if (CS%use_tides) then
call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp)
HA_CSp => CS%HA_CSp
diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90
index 72c7dbe6cd..c4b7c89edc 100644
--- a/src/core/MOM_dynamics_unsplit.F90
+++ b/src/core/MOM_dynamics_unsplit.F90
@@ -707,7 +707,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS
call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
cont_stencil = continuity_stencil(CS%continuity_CSp)
call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
- if (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp)
+ if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp)
if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
CS%SAL_CSp, CS%tides_CSp)
diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90
index cc37f1c2bc..7ca2c700f7 100644
--- a/src/core/MOM_dynamics_unsplit_RK2.F90
+++ b/src/core/MOM_dynamics_unsplit_RK2.F90
@@ -671,7 +671,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag
call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
cont_stencil = continuity_stencil(CS%continuity_CSp)
call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
- if (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp)
+ if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp)
if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
CS%SAL_CSp, CS%tides_CSp)
diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90
index c15b6bd54b..fb833189ec 100644
--- a/src/core/MOM_open_boundary.F90
+++ b/src/core/MOM_open_boundary.F90
@@ -5470,7 +5470,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
endif
I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale
if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
- ! Calculate weights. Both a and u_L are nodim. Adding them together has no meaning.
+ ! Calculate weights. Both a and u_L are nondim. Adding them together has no meaning.
! However, since they cannot be both non-zero, adding them works like a switch.
! When InvLscale_out is 0 and outflow, only interior data is applied to reservoirs
! When InvLscale_in is 0 and inflow, only nudged data is applied to reservoirs
diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90
index f72b6f513e..9aa325cfb8 100644
--- a/src/parameterizations/lateral/MOM_self_attr_load.F90
+++ b/src/parameterizations/lateral/MOM_self_attr_load.F90
@@ -1,16 +1,18 @@
module MOM_self_attr_load
-use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_MODULE
-use MOM_domains, only : pass_var
-use MOM_error_handler, only : MOM_error, FATAL, WARNING
-use MOM_file_parser, only : read_param, get_param, log_version, param_file_type
+use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_MODULE
+use MOM_domains, only : pass_var
+use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
+use MOM_grid, only : ocean_grid_type
+use MOM_io, only : slasher, MOM_read_data
+use MOM_load_love_numbers, only : Love_Data
use MOM_obsolete_params, only : obsolete_logical, obsolete_int
-use MOM_grid, only : ocean_grid_type
-use MOM_unit_scaling, only : unit_scale_type
use MOM_spherical_harmonics, only : spherical_harmonics_init, spherical_harmonics_end
use MOM_spherical_harmonics, only : spherical_harmonics_forward, spherical_harmonics_inverse
use MOM_spherical_harmonics, only : sht_CS, order2index, calc_lmax
-use MOM_load_love_numbers, only : Love_Data
+use MOM_unit_scaling, only : unit_scale_type
+use MOM_verticalGrid, only : verticalGrid_type
implicit none ; private
@@ -27,12 +29,21 @@ module MOM_self_attr_load
logical :: use_tidal_sal_prev = .false.
!< If true, read the tidal SAL from the previous iteration of the tides to
!! facilitate convergence.
- real :: sal_scalar_value !< The constant of proportionality between sea surface height
- !! (really it should be bottom pressure) anomalies and bottom
- !! geopotential anomalies [nondim].
- type(sht_CS), allocatable :: sht !< Spherical harmonic transforms (SHT) control structure
- integer :: sal_sht_Nd !< Maximum degree for SHT [nodim]
- real, allocatable :: Love_Scaling(:) !< Love number for each SHT mode [nodim]
+ logical :: use_bpa = .false.
+ !< If true, use bottom pressure anomaly instead of SSH to calculate SAL.
+ real :: eta_prop
+ !< The partial derivative of eta_sal with the local value of eta [nondim].
+ real :: linear_scaling
+ !< Dimensional coefficients for scalar SAL [nondim or Z T2 L-2 R-1 ~> m Pa-1]
+ type(sht_CS), allocatable :: sht
+ !< Spherical harmonic transforms (SHT) control structure
+ integer :: sal_sht_Nd
+ !< Maximum degree for spherical harmonic transforms [nondim]
+ real, allocatable :: ebot_ref(:,:)
+ !< Reference bottom pressure scaled by Rho_0 and G_Earth[Z ~> m]
+ real, allocatable :: Love_scaling(:)
+ !< Dimensional coefficients for harmonic SAL, which are functions of Love numbers
+ !! [nondim] or [Z T2 L-2 R-1 ~> m Pa-1], depending on the value of use_ppa.
real, allocatable :: Snm_Re(:), & !< Real SHT coefficient for SHT SAL [Z ~> m]
Snm_Im(:) !< Imaginary SHT coefficient for SHT SAL [Z ~> m]
end type SAL_CS
@@ -41,47 +52,54 @@ module MOM_self_attr_load
contains
-!> This subroutine calculates seawater self-attraction and loading based on sea surface height. This should
-!! 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.
+!> This subroutine calculates seawater self-attraction and loading based on either sea surface height (SSH) or bottom
+!! pressure anomaly. 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.
+!! The input field can be either SSH [Z ~> m] or total bottom pressure [R L2 T-2 ~> Pa]. If total bottom pressure is
+!! used, bottom pressure anomaly is first calculated by subtracting a reference bottom pressure from an input file.
+!! The output field is expressed as geopotential height anomaly, and therefore has the unit of [Z ~> m].
subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale)
- type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
+ 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].
+ !! a time-mean geoid or total bottom pressure [Z ~> m] or [R L2 T-2 ~> Pa].
+ real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_sal !< The geopotential 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]
+ 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
+ real, dimension(SZI_(G),SZJ_(G)) :: bpa ! SSH or bottom pressure anomaly [Z ~> m] or [R L2 T-2 ~> Pa]
integer :: n, m, l
integer :: Isq, Ieq, Jsq, Jeq
integer :: i, j
- real :: eta_prop ! The scalar constant of proportionality between eta and eta_sal [nondim]
call cpu_clock_begin(id_clock_SAL)
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
+ if (CS%use_bpa) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ bpa(i,j) = eta(i,j) - CS%ebot_ref(i,j)
+ enddo ; enddo ; else ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+ bpa(i,j) = eta(i,j)
+ enddo ; enddo ; endif
+
! use the scalar approximation and/or iterative tidal SAL
if (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then
- call scalar_SAL_sensitivity(CS, eta_prop)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
- eta_sal(i,j) = eta_prop*eta(i,j)
+ eta_sal(i,j) = CS%linear_scaling * bpa(i,j)
enddo ; enddo
! 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, tmp_scale=tmp_scale)
+ call spherical_harmonics_forward(G, CS%sht, bpa, 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
l = order2index(m, CS%sal_sht_Nd)
do n = m,CS%sal_sht_Nd
- CS%Snm_Re(l+n-m) = CS%Snm_Re(l+n-m) * CS%Love_Scaling(l+n-m)
- CS%Snm_Im(l+n-m) = CS%Snm_Im(l+n-m) * CS%Love_Scaling(l+n-m)
+ CS%Snm_Re(l+n-m) = CS%Snm_Re(l+n-m) * CS%Love_scaling(l+n-m)
+ CS%Snm_Im(l+n-m) = CS%Snm_Im(l+n-m) * CS%Love_scaling(l+n-m)
enddo
enddo
@@ -98,38 +116,36 @@ subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale)
call cpu_clock_end(id_clock_SAL)
end subroutine calc_SAL
-!> This subroutine calculates the partial derivative of the local geopotential height with the input
-!! sea surface height due to the scalar approximation of self-attraction and loading.
+!> This subroutine returns eta_prop member of SAL_CS type, which is the non-dimensional partial
+!! derivative of the local geopotential height with the input sea surface height due to the scalar
+!! approximation of self-attraction and loading.
subroutine scalar_SAL_sensitivity(CS, deta_sal_deta)
type(SAL_CS), intent(in) :: CS !< The control structure returned by a previous call to SAL_init.
real, intent(out) :: deta_sal_deta !< The partial derivative of eta_sal with
!! the local value of eta [nondim].
-
- if (CS%use_sal_scalar .and. CS%use_tidal_sal_prev) then
- deta_sal_deta = 2.0*CS%sal_scalar_value
- elseif (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then
- deta_sal_deta = CS%sal_scalar_value
- else
- deta_sal_deta = 0.0
- endif
+ deta_sal_deta = CS%eta_prop
end subroutine scalar_SAL_sensitivity
!> This subroutine calculates coefficients of the spherical harmonic modes for self-attraction and loading.
!! The algorithm is based on the SAL implementation in MPAS-ocean, which was modified by Kristin Barton from
!! routine written by K. Quinn (March 2010) and modified by M. Schindelegger (May 2017).
-subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_Scaling)
- integer, intent(in) :: nlm !< Maximum spherical harmonics degree [nondim]
- real, intent(in) :: rhoW !< The average density of sea water [R ~> kg m-3]
- real, intent(in) :: rhoE !< The average density of Earth [R ~> kg m-3]
- real, dimension(:), intent(out) :: Love_Scaling !< Scaling factors for inverse SHT [nondim]
+subroutine calc_love_scaling(rhoW, rhoE, grav, CS)
+ real, intent(in) :: rhoW !< The average density of sea water [R ~> kg m-3]
+ real, intent(in) :: rhoE !< The average density of Earth [R ~> kg m-3]
+ real, intent(in) :: grav !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
+ type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init.
! Local variables
+ real :: coef_rhoE ! A scaling coefficient of solid Earth density. coef_rhoE = rhoW / rhoE with USE_BPA=False
+ ! and coef_rhoE = 1.0 / (rhoE * grav) with USE_BPA=True. [nondim] or [Z T2 L-2 R-1 ~> m Pa-1]
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_tot ! Size of the stored Love numbers [nondim]
+ integer :: nlm ! Maximum spherical harmonics degree [nondim]
integer :: n, m, l
n_tot = size(Love_Data, dim=2)
+ nlm = CS%sal_sht_Nd
if (nlm+1 > n_tot) call MOM_error(FATAL, "MOM_tidal_forcing " // &
"calc_love_scaling: maximum spherical harmonics degree is larger than " // &
@@ -146,36 +162,47 @@ subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_Scaling)
KDat(2) = (-1.0 / 3.0) * H1 - (2.0 / 3.0) * L1 - 1.0
endif
+ if (CS%use_bpa) then
+ coef_rhoE = 1.0 / (rhoE * grav) ! [Z T2 L-2 R-1 ~> m Pa-1]
+ else
+ coef_rhoE = rhoW / rhoE ! [nondim]
+ endif
+
do m=0,nlm ; do n=m,nlm
- l = order2index(m,nlm)
- Love_Scaling(l+n-m) = (3.0 / real(2*n+1)) * (rhoW / rhoE) * (1.0 + KDat(n+1) - HDat(n+1))
+ l = order2index(m, nlm)
+ ! Love_scaling has the same as coef_rhoE.
+ CS%Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * coef_rhoE * (1.0 + KDat(n+1) - HDat(n+1))
enddo ; enddo
end subroutine calc_love_scaling
!> This subroutine initializes the self-attraction and loading control structure.
-subroutine SAL_init(G, US, param_file, CS)
- type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
- type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
- type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
+subroutine SAL_init(G, GV, US, param_file, CS)
+ type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
+ type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
+ type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
+ type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
type(SAL_CS), intent(inout) :: CS !< Self-attraction and loading control structure
! Local variables
# include "version_variable.h"
character(len=40) :: mdl = "MOM_self_attr_load" ! This module's name.
integer :: lmax ! Total modes of the real spherical harmonics [nondim]
- real :: rhoW ! The average density of sea water [R ~> kg m-3].
real :: rhoE ! The average density of Earth [R ~> kg m-3].
+ character(len=200) :: filename, ebot_ref_file, inputdir ! Strings for file/path
+ character(len=200) :: ebot_ref_varname ! Variable name in file
+ logical :: calculate_sal=.false.
+ logical :: tides=.false., use_tidal_sal_file=.false., bq_sal_tides_bug=.false.
+ integer :: tides_answer_date=99991231 ! Recover old answers with tides
+ real :: sal_scalar_value, tide_sal_scalar_value ! Scaling SAL factors [nondim]
+ integer :: isd, ied, jsd, jed
- logical :: calculate_sal
- logical :: tides, use_tidal_sal_file
- real :: tide_sal_scalar_value ! Scaling SAL factor [nondim]
+ isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
! Read all relevant parameters and write them to the model log.
call log_version(param_file, mdl, version, "")
call get_param(param_file, '', "TIDES", tides, default=.false., do_not_log=.True.)
- call get_param(param_file, mdl, "CALCULATE_SAL", calculate_sal, "If true, calculate "//&
- " self-attraction and loading.", default=tides, do_not_log=.True.)
+ call get_param(param_file, '', "CALCULATE_SAL", calculate_sal, default=tides, do_not_log=.True.)
if (.not. calculate_sal) return
if (tides) then
@@ -183,17 +210,47 @@ subroutine SAL_init(G, US, param_file, CS)
default=.false., do_not_log=.True.)
call get_param(param_file, '', "TIDAL_SAL_FROM_FILE", use_tidal_sal_file, &
default=.false., do_not_log=.True.)
+ call get_param(param_file, '', "TIDES_ANSWER_DATE", tides_answer_date, &
+ default=20230630, do_not_log=.True.)
endif
+ call get_param(param_file, mdl, "SAL_USE_BPA", CS%use_bpa, &
+ "If true, use bottom pressure anomaly to calculate self-attraction and "// &
+ "loading (SAL). Otherwise sea surface height anomaly is used, which is "// &
+ "only correct for homogenous flow.", default=.False.)
+ if (CS%use_bpa) then
+ call get_param(param_file, '', "INPUTDIR", inputdir, default=".", do_not_log=.True.)
+ inputdir = slasher(inputdir)
+ call get_param(param_file, mdl, "REF_BOT_PRES_FILE", ebot_ref_file, &
+ "Reference bottom pressure file used by self-attraction and loading (SAL).", &
+ default="pbot.nc")
+ call get_param(param_file, mdl, "REF_BOT_PRES_VARNAME", ebot_ref_varname, &
+ "The name of the variable in REF_BOT_PRES_FILE with reference bottom "//&
+ "pressure. The variable should have the unit of Pa.", &
+ default="pbot")
+ filename = trim(inputdir)//trim(ebot_ref_file)
+ call log_param(param_file, mdl, "INPUTDIR/REF_BOT_PRES_FILE", filename)
+
+ allocate(CS%ebot_ref(isd:ied, jsd:jed), source=0.0)
+ call MOM_read_data(filename, trim(ebot_ref_varname), CS%ebot_ref, G%Domain,&
+ scale=US%Pa_to_RL2_T2)
+ call pass_var(CS%ebot_ref, G%Domain)
+ endif
+ if (tides_answer_date<=20250131 .and. CS%use_bpa) &
+ call MOM_error(FATAL, trim(mdl) // ", SAL_init: SAL_USE_BPA needs to be false to recover "//&
+ "tide answers before 20250131.")
call get_param(param_file, mdl, "SAL_SCALAR_APPROX", CS%use_sal_scalar, &
"If true, use the scalar approximation to calculate self-attraction and "//&
"loading.", default=tides .and. (.not. use_tidal_sal_file))
+ if (CS%use_sal_scalar .and. CS%use_bpa) &
+ call MOM_error(WARNING, trim(mdl) // ", SAL_init: Using bottom pressure anomaly for scalar "//&
+ "approximation SAL is unsubstantiated.")
call get_param(param_file, '', "TIDE_SAL_SCALAR_VALUE", tide_sal_scalar_value, &
units="m m-1", default=0.0, do_not_log=.True.)
if (tide_sal_scalar_value/=0.0) &
call MOM_error(WARNING, "TIDE_SAL_SCALAR_VALUE is a deprecated parameter. "//&
"Use SAL_SCALAR_VALUE instead." )
- call get_param(param_file, mdl, "SAL_SCALAR_VALUE", CS%sal_scalar_value, &
+ call get_param(param_file, mdl, "SAL_SCALAR_VALUE", sal_scalar_value, &
"The constant of proportionality between sea surface "//&
"height (really it should be bottom pressure) anomalies "//&
"and bottom geopotential anomalies. This is only used if "//&
@@ -207,20 +264,35 @@ subroutine SAL_init(G, US, param_file, CS)
"The maximum degree of the spherical harmonics transformation used for "// &
"calculating the self-attraction and loading term.", &
default=0, do_not_log=(.not. CS%use_sal_sht))
- call get_param(param_file, '', "RHO_0", rhoW, default=1035.0, scale=US%kg_m3_to_R, &
- units="kg m-3", do_not_log=.True.)
call get_param(param_file, mdl, "RHO_SOLID_EARTH", rhoE, &
"The mean solid earth density. This is used for calculating the "// &
"self-attraction and loading term.", units="kg m-3", &
default=5517.0, scale=US%kg_m3_to_R, do_not_log=(.not. CS%use_sal_sht))
+ ! Set scaling coefficients for scalar approximation
+ if (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then
+ if (CS%use_sal_scalar .and. CS%use_tidal_sal_prev) then
+ CS%eta_prop = 2.0 * sal_scalar_value
+ else
+ CS%eta_prop = sal_scalar_value
+ endif
+ if (CS%use_bpa) then
+ CS%linear_scaling = CS%eta_prop / (GV%Rho0 * GV%g_Earth)
+ else
+ CS%linear_scaling = CS%eta_prop
+ endif
+ else
+ CS%eta_prop = 0.0 ; CS%linear_scaling = 0.0
+ endif
+
+ ! Set scaling coefficients for spherical harmonics
if (CS%use_sal_sht) then
lmax = calc_lmax(CS%sal_sht_Nd)
- allocate(CS%Snm_Re(lmax)); CS%Snm_Re(:) = 0.0
- allocate(CS%Snm_Im(lmax)); CS%Snm_Im(:) = 0.0
+ allocate(CS%Snm_Re(lmax), source=0.0)
+ allocate(CS%Snm_Im(lmax), source=0.0)
- allocate(CS%Love_Scaling(lmax)); CS%Love_Scaling(:) = 0.0
- call calc_love_scaling(CS%sal_sht_Nd, rhoW, rhoE, CS%Love_Scaling)
+ allocate(CS%Love_scaling(lmax), source=0.0)
+ call calc_love_scaling(GV%Rho0, rhoE, GV%g_Earth, CS)
allocate(CS%sht)
call spherical_harmonics_init(G, param_file, CS%sht)
@@ -234,8 +306,11 @@ end subroutine SAL_init
subroutine SAL_end(CS)
type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call
!! to SAL_init; it is deallocated here.
+
+ if (allocated(CS%ebot_ref)) deallocate(CS%ebot_ref)
+
if (CS%use_sal_sht) then
- if (allocated(CS%Love_Scaling)) deallocate(CS%Love_Scaling)
+ if (allocated(CS%Love_scaling)) deallocate(CS%Love_scaling)
if (allocated(CS%Snm_Re)) deallocate(CS%Snm_Re)
if (allocated(CS%Snm_Im)) deallocate(CS%Snm_Im)
call spherical_harmonics_end(CS%sht)
@@ -247,20 +322,20 @@ end subroutine SAL_end
!!
!! \section section_SAL Self attraction and loading
!!
-!! This module contains methods to calculate self-attraction and loading (SAL) as a function of sea surface height (SSH)
-!! (rather, it should be bottom pressure anomaly). SAL is primarily used for fast evolving processes like tides or
-!! storm surges, but the effect applies to all motions.
+!! This module contains methods to calculate self-attraction and loading (SAL) as a function of sea surface height or
+!! bottom pressure anomaly. SAL is primarily used for fast evolving processes like tides or storm surges, but the
+!! effect applies to all motions.
!!
!! If SAL_SCALAR_APPROX is true, a scalar approximation is applied (\cite Accad1978) and the SAL is simply
-!! a fraction (set by SAL_SCALAR_VALUE, usually around 10% for global tides) of local SSH.
-!! For tides, the scalar approximation can also be used to iterate the SAL to convergence [see
-!! USE_PREVIOUS_TIDES in MOM_tidal_forcing, \cite Arbic2004].
+!! a fraction (set by SAL_SCALAR_VALUE, usually around 10% for global tides) of local SSH. For tides, the
+!! scalar approximation can also be used to iterate the SAL to convergence [see USE_PREVIOUS_TIDES in
+!! MOM_tidal_forcing, \cite Arbic2004].
!!
!! If SAL_HARMONICS is true, a more accurate online spherical harmonic transforms are used to calculate
-!! SAL. Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is
-!! set by SAL_HARMONICS_DEGREE. The algorithm is based on SAL calculation in Model for Prediction Across
-!! Scales (MPAS)-Ocean
-!! developed by Los Alamos National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023].
+!! SAL. Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is set
+!! by SAL_HARMONICS_DEGREE. The algorithm is based on SAL calculation in Model for Prediction Across
+!! Scales (MPAS)-Ocean developed by Los Alamos National Laboratory and University of Michigan
+!! [\cite Barton2022 and \cite Brus2023].
!!
!! References:
!!
diff --git a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 b/src/parameterizations/lateral/MOM_spherical_harmonics.F90
index 4f9cee03aa..067e2e3e94 100644
--- a/src/parameterizations/lateral/MOM_spherical_harmonics.F90
+++ b/src/parameterizations/lateral/MOM_spherical_harmonics.F90
@@ -1,12 +1,12 @@
!> Laplace's spherical harmonic transforms (SHT)
module MOM_spherical_harmonics
+use MOM_coms_infra, only : sum_across_PEs
+use MOM_coms, only : reproducing_sum
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, &
CLOCK_MODULE, CLOCK_ROUTINE, CLOCK_LOOP
use MOM_error_handler, only : MOM_error, FATAL
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
-use MOM_coms_infra, only : sum_across_PEs
-use MOM_coms, only : reproducing_sum
implicit none ; private
diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
index 7ca432fea4..8ed75a9f44 100644
--- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
+++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
@@ -73,7 +73,7 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int)
Kd, & ! A column of diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
h_top, h_bot ! Distances from the top or bottom [H ~> m or kg m-2].
real :: dz_h_int ! The ratio of the vertical distances across the layers surrounding an interface
- ! over the layer thicknesses [H Z-1 ~> nonodim or kg m-3]
+ ! over the layer thicknesses [H Z-1 ~> nondim or kg m-3]
real :: ustar ! The local friction velocity [Z T-1 ~> m s-1]
real :: absf ! The absolute value of the Coriolis parameter [T-1 ~> s-1]
real :: htot ! The sum of the thicknesses [H ~> m or kg m-2].