Skip to content

Commit

Permalink
Alternative method for SAL and tides in Boussinesq
Browse files Browse the repository at this point in the history 
Add an alternative method for SAL and tides in Boussinesq mode. The
current method adjusts interface heights with geopotential height
anomaly for SAL and tides. For non-Boussinesq mode, the current method
is algebraically the same as taking the gradient of SAL and tide
geopotential (body forcing). For Boussinesq mode, there is a baroclinic
component of tidal forcing and SAL. The alternative method is added to
calculate the gradient of tidal forcing and SAL directly at the cost of
additional multiplications. The new method is controlled by runtime
parameter BOUSSINESQ_SAL_TIDES.
  • Loading branch information
@herrwang0
herrwang0 committed Sep 29, 2024
1 parent 55ad253 commit aca2291
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 40 deletions.
120 changes: 82 additions & 38 deletions src/core/MOM_PressureForce_FV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,8 @@ 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
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_tidal_eq = -1 !< Diagnostic identifier
Expand Down Expand Up @@ -821,7 +823,12 @@ 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_tidal_eq+e_tidal_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_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)
Expand Down Expand Up @@ -1189,43 +1196,39 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
enddo ; enddo
enddo

! Self-attraction and loading (SAL) and tides (new answers)
if (CS%tides_answer_date>20230630) then
! SAL
if (CS%calculate_SAL) then
if (CS%sal_use_bpa) then
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
pbot_anom(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1)
enddo ; enddo
else
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
pbot_anom(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topogrpahy above sea level
enddo ; enddo
endif
call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
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
! Self-attraction and loading (SAL) (new answers)
if (CS%calculate_SAL .and. CS%tides_answer_date>20230630) then
if (CS%sal_use_bpa) then
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
pbot_anom(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1)
enddo ; enddo
else
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
pbot_anom(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topography above sea level
enddo ; enddo
endif
call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
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

! Tides
if (CS%tides) then
call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
! Tides (new answers)
if (CS%tides .and. CS%tides_answer_date>20230630) 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 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
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
Expand Down Expand Up @@ -1604,6 +1607,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
Expand Down Expand Up @@ -1649,7 +1680,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
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
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
Expand All @@ -1662,7 +1693,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
enddo ; enddo
endif
endif
if (CS%calculate_SAL .and. (CS%tides_answer_date>20230630)) then
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
Expand Down Expand Up @@ -1709,7 +1740,12 @@ 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+e_tidal_eq+e_tidal_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_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)
Expand Down Expand Up @@ -1780,6 +1816,14 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
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.)
Expand Down
4 changes: 2 additions & 2 deletions src/parameterizations/lateral/MOM_self_attr_load.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ module MOM_self_attr_load
real, allocatable :: Love_scaling(:)
!< Dimensional coefficients for harmonic SAL, which are functions of Love numbers
!! [nondim or L2 Z-1 T-2 ~> m s-2 or R-1 ~> m3 kg-1 or L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
real, allocatable :: Snm_Re(:), Snm_Im(:)
!< Real and imaginary coefficients for harmonic SAL [Z ~> m or L2 T-2 ~> m2 s-2]
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

integer :: id_clock_SAL !< CPU clock for self-attraction and loading
Expand Down

0 comments on commit aca2291

Please sign in to comment.