Skip to content

Commit 9762fdf

Browse files
herrwang0herrwang0
committed
Alternative method for SAL and tides in Boussinesq
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.
1 parent 55ad253 commit 9762fdf

File tree

1 file changed

+70
-36
lines changed

1 file changed

+70
-36
lines changed

‎src/core/MOM_PressureForce_FV.F90‎

Lines changed: 70 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,8 @@ module MOM_PressureForce_FV
8282
!! equation of state is 0 to account for the displacement of the sea
8383
!! surface including adjustments for atmospheric or sea-ice pressure.
8484
logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF
85+
logical :: bq_sal_tides = .false. !< If true, use an alternative method for SAL and tides
86+
!! in Boussinesq mode
8587
integer :: tides_answer_date = 99991231 !< Recover old answers with tides
8688
integer :: id_e_tide = -1 !< Diagnostic identifier
8789
integer :: id_e_tidal_eq = -1 !< Diagnostic identifier
@@ -1189,43 +1191,39 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
11891191
enddo ; enddo
11901192
enddo
11911193

1192-
! Self-attraction and loading (SAL) and tides (new answers)
1193-
if (CS%tides_answer_date>20230630) then
1194-
! SAL
1195-
if (CS%calculate_SAL) then
1196-
if (CS%sal_use_bpa) then
1197-
!$OMP parallel do default(shared)
1198-
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1199-
pbot_anom(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1)
1200-
enddo ; enddo
1201-
else
1202-
!$OMP parallel do default(shared)
1203-
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1204-
pbot_anom(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topogrpahy above sea level
1205-
enddo ; enddo
1206-
endif
1207-
call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
1208-
do K=1,nz+1
1209-
!$OMP parallel do default(shared)
1210-
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1211-
e(i,j,K) = e(i,j,K) - e_sal(i,j)
1212-
pa(i,j,K) = pa(i,j,K) - GxRho * e_sal(i,j)
1213-
enddo ; enddo
1214-
enddo
1194+
! Self-attraction and loading (SAL) (new answers)
1195+
if (CS%calculate_SAL .and. CS%tides_answer_date>20230630) then
1196+
if (CS%sal_use_bpa) then
1197+
!$OMP parallel do default(shared)
1198+
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1199+
pbot_anom(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1)
1200+
enddo ; enddo
1201+
else
1202+
!$OMP parallel do default(shared)
1203+
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1204+
pbot_anom(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topogrpahy above sea level
1205+
enddo ; enddo
12151206
endif
1207+
call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
1208+
if (.not.CS%bq_sal_tides) then ; do K=1,nz+1
1209+
!$OMP parallel do default(shared)
1210+
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1211+
e(i,j,K) = e(i,j,K) - e_sal(i,j)
1212+
pa(i,j,K) = pa(i,j,K) - GxRho * e_sal(i,j)
1213+
enddo ; enddo
1214+
enddo ; endif
1215+
endif
12161216

1217-
! Tides
1218-
if (CS%tides) then
1219-
call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
1217+
! Tides (new answers)
1218+
if (CS%tides .and. CS%tides_answer_date>20230630) then
1219+
call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
1220+
if (.not.CS%bq_sal_tides) then ; do K=1,nz+1
12201221
!$OMP parallel do default(shared)
1221-
do K=1,nz+1
1222-
!$OMP parallel do default(shared)
1223-
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1224-
e(i,j,K) = e(i,j,K) - (e_tidal_eq(i,j) + e_tidal_sal(i,j))
1225-
pa(i,j,K) = pa(i,j,K) - GxRho * (e_tidal_eq(i,j) + e_tidal_sal(i,j))
1226-
enddo ; enddo
1227-
enddo
1228-
endif
1222+
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1223+
e(i,j,K) = e(i,j,K) - (e_tidal_eq(i,j) + e_tidal_sal(i,j))
1224+
pa(i,j,K) = pa(i,j,K) - GxRho * (e_tidal_eq(i,j) + e_tidal_sal(i,j))
1225+
enddo ; enddo
1226+
enddo ; endif
12291227
endif
12301228

12311229
if (CS%correction_intxpa .or. CS%reset_intxpa_integral) then
@@ -1604,6 +1602,34 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
16041602
((h(i,j,k) + h(i,j+1,k)) + h_neglect))
16051603
enddo ; enddo ; enddo
16061604

1605+
! Calculate SAL geopotential anomaly and add its gradient to pressure gradient force
1606+
if (CS%calculate_SAL .and. CS%tides_answer_date>20230630 .and. CS%bq_sal_tides) then
1607+
!$OMP parallel do default(shared)
1608+
do k=1,nz
1609+
do j=js,je ; do I=Isq,Ieq
1610+
PFu(I,j,k) = PFu(I,j,k) + (e_sal(i+1,j) - e_sal(i,j)) * GV%g_Earth * G%IdxCu(I,j)
1611+
enddo ; enddo
1612+
do J=Jsq,Jeq ; do i=is,ie
1613+
PFv(i,J,k) = PFv(i,J,k) + (e_sal(i,j+1) - e_sal(i,j)) * GV%g_Earth * G%IdyCv(i,J)
1614+
enddo ; enddo
1615+
enddo
1616+
endif
1617+
1618+
! Calculate tidal geopotential anomaly and add its gradient to pressure gradient force
1619+
if (CS%tides .and. CS%tides_answer_date>20230630 .and. CS%bq_sal_tides) then
1620+
!$OMP parallel do default(shared)
1621+
do k=1,nz
1622+
do j=js,je ; do I=Isq,Ieq
1623+
PFu(I,j,k) = PFu(I,j,k) + ((e_tidal_eq(i+1,j) + e_tidal_sal(i+1,j)) &
1624+
- (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * GV%g_Earth * G%IdxCu(I,j)
1625+
enddo ; enddo
1626+
do J=Jsq,Jeq ; do i=is,ie
1627+
PFv(i,J,k) = PFv(i,J,k) + ((e_tidal_eq(i,j+1) + e_tidal_sal(i,j+1)) &
1628+
- (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * GV%g_Earth * G%IdyCv(i,J)
1629+
enddo ; enddo
1630+
enddo
1631+
endif
1632+
16071633
if (CS%GFS_scale < 1.0) then
16081634
! Adjust the Montgomery potential to make this a reduced gravity model.
16091635
if (use_EOS) then
@@ -1649,7 +1675,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
16491675
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
16501676
eta(i,j) = e(i,j,1)*GV%Z_to_H
16511677
enddo ; enddo
1652-
if (CS%tides) then
1678+
if (CS%tides .and. (.not.CS%bq_sal_tides)) then
16531679
if (CS%tides_answer_date>20230630) then
16541680
!$OMP parallel do default(shared)
16551681
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
@@ -1662,7 +1688,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
16621688
enddo ; enddo
16631689
endif
16641690
endif
1665-
if (CS%calculate_SAL .and. (CS%tides_answer_date>20230630)) then
1691+
if (CS%calculate_SAL .and. (CS%tides_answer_date>20230630) .and. (.not.CS%bq_sal_tides)) then
16661692
!$OMP parallel do default(shared)
16671693
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
16681694
eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H
@@ -1780,6 +1806,14 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
17801806
if (CS%calculate_SAL) &
17811807
call get_param(param_file, '', "SAL_USE_BPA", CS%sal_use_bpa, default=.false., &
17821808
do_not_log=.true.)
1809+
if ((CS%tides .or. CS%calculate_SAL) .and. GV%Boussinesq) &
1810+
call get_param(param_file, mdl, "BOUSSINESQ_SAL_TIDES", CS%bq_sal_tides, "If true, "//&
1811+
"in Boussinesq mode, use an alternative method to include self-attraction "//&
1812+
"and loading (SAL) and tidal forcings in pressure gradient, in which their "//&
1813+
"gradients are calculated separately, instead of adding geopotential "//&
1814+
"anomalies as corrections to the interface height. This alternative method "//&
1815+
"elimates a baroclinic component of the SAL and tidal forcings.", &
1816+
default=.false.)
17831817
call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &
17841818
"If true, Temperature and salinity are used as state variables.", &
17851819
default=.true., do_not_log=.true.)

0 commit comments

Comments
 (0)