Skip to content

Separate self-attraction and loading as a new module #382

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

Merged
marshallward merged 9 commits into from
Aug 24, 2023
Merged

Separate self-attraction and loading as a new module #382

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
cc94945
Separate SAL from tidal_forcing
herrwang0 Mar 9, 2023
415395b
Decompose output from calc_tidal_forcing
herrwang0 Apr 26, 2023
47f6a0a
Fix SSH for calculating SAL with flooding points
herrwang0 Jun 6, 2023
c73236a
Renaming input parameter names for the SAL module
herrwang0 Jun 6, 2023
00bcbc6
Recover old answers with tides in Boussinesq mode
herrwang0 Jun 16, 2023
f1ba977
Refactor SAL and tides calls in Boussinesq mode
herrwang0 Jun 20, 2023
4ef6440
Change SAL related parameter names
herrwang0 Aug 17, 2023
12ee370
Modify reading SAL related parameters
herrwang0 Aug 17, 2023
121823d
*+Use TIDES_ANSWER_DATE with semi-Boussinesq tides
Hallberg-NOAA Aug 22, 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
10 changes: 6 additions & 4 deletions src/core/MOM_PressureForce.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ module MOM_PressureForce
use MOM_PressureForce_Mont, only : PressureForce_Mont_Bouss, PressureForce_Mont_nonBouss
use MOM_PressureForce_Mont, only : PressureForce_Mont_init
use MOM_PressureForce_Mont, only : PressureForce_Mont_CS
use MOM_self_attr_load, only : SAL_CS
use MOM_tidal_forcing, only : tidal_forcing_CS
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
Expand Down Expand Up @@ -80,15 +81,16 @@ subroutine PressureForce(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, e
end subroutine Pressureforce

!> Initialize the pressure force control structure
subroutine PressureForce_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
subroutine PressureForce_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, tides_CSp)
type(time_type), target, intent(in) :: Time !< Current model time
type(ocean_grid_type), intent(in) :: G !< Ocean 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 !< Parameter file handles
type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
type(PressureForce_CS), intent(inout) :: CS !< Pressure force control structure
type(tidal_forcing_CS), intent(inout), optional :: tides_CSp !< Tide control structure
type(SAL_CS), intent(in), optional :: SAL_CSp !< SAL control structure
type(tidal_forcing_CS), intent(in), optional :: tides_CSp !< Tide control structure
#include "version_variable.h"
character(len=40) :: mdl = "MOM_PressureForce" ! This module's name.

Expand All @@ -103,10 +105,10 @@ subroutine PressureForce_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)

if (CS%Analytic_FV_PGF) then
call PressureForce_FV_init(Time, G, GV, US, param_file, diag, &
CS%PressureForce_FV, tides_CSp)
CS%PressureForce_FV, SAL_CSp, tides_CSp)
else
call PressureForce_Mont_init(Time, G, GV, US, param_file, diag, &
CS%PressureForce_Mont, tides_CSp)
CS%PressureForce_Mont, SAL_CSp, tides_CSp)
endif
end subroutine PressureForce_init

Expand Down
241 changes: 191 additions & 50 deletions src/core/MOM_PressureForce_FV.F90
View file
Open in desktop

Large diffs are not rendered by default.

115 changes: 84 additions & 31 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ module MOM_PressureForce_Mont
use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_self_attr_load, only : calc_SAL, SAL_CS
use MOM_tidal_forcing, only : calc_tidal_forcing, tidal_forcing_CS
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
Expand All @@ -31,6 +32,7 @@ module MOM_PressureForce_Mont
!> Control structure for the Montgomery potential form of pressure gradient
type, public :: PressureForce_Mont_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.
real :: Rho0 !< The density used in the Boussinesq
!! approximation [R ~> kg m-3].
Expand All @@ -45,8 +47,10 @@ module MOM_PressureForce_Mont
real, allocatable :: PFv_bc(:,:,:) !< Meridional accelerations due to pressure gradients
!! deriving from density gradients within layers [L T-2 ~> m s-2].
!>@{ Diagnostic IDs
integer :: id_PFu_bc = -1, id_PFv_bc = -1, id_e_tidal = -1
integer :: id_PFu_bc = -1, id_PFv_bc = -1, id_e_sal = -1
integer :: id_e_tide = -1, id_e_tide_eq = -1, id_e_tide_sal = -1
!>@}
type(SAL_CS), pointer :: SAL_CSp => NULL() !< SAL control structure
type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< The tidal forcing control structure
end type PressureForce_Mont_CS

Expand Down Expand Up @@ -103,8 +107,10 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
! of a reduced gravity form of the equations [L2 T-2 ~> m2 s-2].
dp_star, & ! Layer thickness after compensation for compressibility [R L2 T-2 ~> Pa].
SSH, & ! The sea surface height anomaly, in depth units [Z ~> m].
e_tidal, & ! Bottom geopotential anomaly due to tidal forces from
! astronomical sources and self-attraction and loading [Z ~> m].
e_sal, & ! Bottom geopotential anomaly due to self-attraction and loading [Z ~> m].
e_tide_eq, & ! Bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m].
e_tide_sal, & ! Bottom geopotential anomaly due to harmonic self-attraction and loading
! specific to tides [Z ~> m].
geopot_bot ! Bottom geopotential relative to a temporally fixed reference value,
! including any tidal contributions [L2 T-2 ~> m2 s-2].
real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
Expand Down Expand Up @@ -180,12 +186,18 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
endif
endif

if (CS%tides) then
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
geopot_bot(i,j) = -GV%g_Earth * G%bathyT(i,j)
enddo ; enddo

! Calculate and add the self-attraction and loading geopotential anomaly.
if (CS%calculate_SAL) then
! Determine the sea surface height anomalies, to enable the calculation
! of self-attraction and loading.
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
SSH(i,j) = -G%bathyT(i,j) - G%Z_ref
SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0)
enddo ; enddo
if (use_EOS) then
!$OMP parallel do default(shared)
Expand All @@ -204,15 +216,19 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
enddo ; enddo ; enddo
endif

call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
geopot_bot(i,j) = -GV%g_Earth*(e_tidal(i,j) + G%bathyT(i,j))
geopot_bot(i,j) = geopot_bot(i,j) - GV%g_Earth*e_sal(i,j)
enddo ; enddo
else
endif

! 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)
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
geopot_bot(i,j) = -GV%g_Earth*G%bathyT(i,j)
geopot_bot(i,j) = geopot_bot(i,j) - GV%g_Earth*(e_tide_eq(i,j) + e_tide_sal(i,j))
enddo ; enddo
endif

Expand Down Expand Up @@ -348,7 +364,12 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb

if (CS%id_PFu_bc>0) call post_data(CS%id_PFu_bc, CS%PFu_bc, CS%diag)
if (CS%id_PFv_bc>0) call post_data(CS%id_PFv_bc, CS%PFv_bc, CS%diag)
if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag)
! 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_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)

end subroutine PressureForce_Mont_nonBouss

Expand Down Expand Up @@ -396,9 +417,11 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
real :: h_star(SZI_(G),SZJ_(G)) ! Layer thickness after compensation
! for compressibility [Z ~> m].
real :: SSH(SZI_(G),SZJ_(G)) ! The sea surface height anomaly, in depth units [Z ~> m].
real :: e_tidal(SZI_(G),SZJ_(G)) ! Bottom geopotential anomaly due to tidal
! forces from astronomical sources and self-
! attraction and loading, in depth units [Z ~> m].
real :: e_sal(SZI_(G),SZJ_(G)) ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m].
real :: e_tide_eq(SZI_(G),SZJ_(G)) ! Bottom geopotential anomaly due to tidal forces from astronomical sources
! [Z ~> m].
real :: e_tide_sal(SZI_(G),SZJ_(G)) ! Bottom geopotential anomaly due to harmonic self-attraction and loading
! specific to tides, in depth units [Z ~> m].
real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
! density [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
real :: I_Rho0 ! 1/Rho0 [R-1 ~> m3 kg-1].
Expand Down Expand Up @@ -440,33 +463,40 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
I_Rho0 = 1.0/CS%Rho0
G_Rho0 = GV%g_Earth / GV%Rho0

if (CS%tides) then
!$OMP parallel do default(shared)
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) = -G%bathyT(i,j) - G%Z_ref ; enddo
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_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
endif

! Here layer interface heights, e, are calculated.
if (CS%tides) then
call calc_SAL(SSH, e_sal, G, CS%SAL_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_tidal(i,j))
e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j)
enddo ; enddo
else
endif

! 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)
!$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(i,j,nz+1) = e(i,j,nz+1) - (e_tide_eq(i,j) + e_tide_sal(i,j))
enddo ; enddo
endif

!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do k=nz,1,-1 ; do i=Isq,Ieq+1
e(i,j,K) = e(i,j,K+1) + h(i,j,k)*GV%H_to_Z
Expand Down Expand Up @@ -582,25 +612,35 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
endif ! use_EOS

if (present(eta)) then
if (CS%tides) 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.
!$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) = e(i,j,1)*GV%Z_to_H + e_tidal(i,j)*GV%Z_to_H
eta(i,j) = eta(i,j) + (e_tide_eq(i,j)+e_tide_sal(i,j))*GV%Z_to_H
enddo ; enddo
else
endif
if (CS%calculate_SAL) 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
eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H
enddo ; enddo
endif
endif

if (CS%id_PFu_bc>0) call post_data(CS%id_PFu_bc, CS%PFu_bc, CS%diag)
if (CS%id_PFv_bc>0) call post_data(CS%id_PFv_bc, CS%PFv_bc, CS%diag)
if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag)
! 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_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)

end subroutine PressureForce_Mont_Bouss

Expand Down Expand Up @@ -821,14 +861,15 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
end subroutine Set_pbce_nonBouss

!> Initialize the Montgomery-potential form of PGF control structure
subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, tides_CSp)
type(time_type), target, intent(in) :: Time !< Current model time
type(ocean_grid_type), intent(in) :: G !< ocean 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 !< Parameter file handles
type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
type(PressureForce_Mont_CS), intent(inout) :: CS !< Montgomery PGF control structure
type(SAL_CS), intent(in), target, optional :: SAL_CSp !< SAL control structure
type(tidal_forcing_CS), intent(in), target, optional :: tides_CSp !< Tides control structure

! Local variables
Expand All @@ -841,6 +882,8 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_
CS%diag => diag ; CS%Time => Time
if (present(tides_CSp)) &
CS%tides_CSp => tides_CSp
if (present(SAL_CSp)) &
CS%SAL_CSp => SAL_CSp

mdl = "MOM_PressureForce_Mont"
call log_version(param_file, mdl, version, "")
Expand All @@ -852,6 +895,8 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_
units="kg m-3", default=1035.0, scale=US%R_to_kg_m3)
call get_param(param_file, mdl, "TIDES", CS%tides, &
"If true, apply tidal momentum forcing.", default=.false.)
call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, &
"If true, calculate self-attraction and loading.", default=CS%tides)
call get_param(param_file, mdl, "USE_EOS", use_EOS, default=.true., &
do_not_log=.true.) ! Input for diagnostic use only.

Expand All @@ -866,9 +911,17 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_
allocate(CS%PFv_bc(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.)
endif

if (CS%calculate_SAL) then
CS%id_e_sal = register_diag_field('ocean_model', 'e_sal', diag%axesT1, Time, &
'Self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m)
endif
if (CS%tides) then
CS%id_e_tidal = 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 = 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, &
'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, &
'Read-in tidal self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m)
endif

CS%GFS_scale = 1.0
Expand Down
Loading