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

Jump to bottom

Nodal modulation #725

Open
wants to merge 2 commits into
base: dev/gfdl
Choose a base branch
from
+92 −32
Open

Nodal modulation #725

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
211cab5
Nodal modulation
c2xu Sep 6, 2024
405213a
Merge branch 'NOAA-GFDL:dev/gfdl' into c2xu/nodal_modulation
c2xu Sep 14, 2024
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
20 changes: 12 additions & 8 deletions src/core/MOM_open_boundary.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1127,26 +1127,30 @@ subroutine initialize_obc_tides(OBC, US, param_file)
type(astro_longitudes) :: nodal_longitudes !< Solar and lunar longitudes for tidal forcing
type(time_type) :: nodal_time !< Model time to calculate nodal modulation for.
integer :: c !< Index to tidal constituent.
logical :: tides !< True if astronomical tides are also used.

call get_param(param_file, mdl, "OBC_TIDE_CONSTITUENTS", tide_constituent_str, &
"Names of tidal constituents being added to the open boundaries.", &
fail_if_missing=.true.)

call get_param(param_file, mdl, "OBC_TIDE_ADD_EQ_PHASE", OBC%add_eq_phase, &
call get_param(param_file, mdl, "TIDES", tides, &
"If true, apply tidal momentum forcing.", default=.false., do_not_log=.true.)

call get_param(param_file, mdl, "TIDE_USE_EQ_PHASE", OBC%add_eq_phase, &
"If true, add the equilibrium phase argument to the specified tidal phases.", &
default=.false., fail_if_missing=.false.)
default=.false., do_not_log=tides)

call get_param(param_file, mdl, "OBC_TIDE_ADD_NODAL", OBC%add_nodal_terms, &
call get_param(param_file, mdl, "TIDE_ADD_NODAL", OBC%add_nodal_terms, &
"If true, include 18.6 year nodal modulation in the boundary tidal forcing.", &
default=.false.)
default=.false., do_not_log=tides)

call get_param(param_file, mdl, "OBC_TIDE_REF_DATE", tide_ref_date, &
call get_param(param_file, mdl, "TIDE_REF_DATE", tide_ref_date, &
"Reference date to use for tidal calculations and equilibrium phase.", &
fail_if_missing=.true.)
default=0, do_not_log=tides)

call get_param(param_file, mdl, "OBC_TIDE_NODAL_REF_DATE", nodal_ref_date, &
call get_param(param_file, mdl, "TIDE_NODAL_REF_DATE", nodal_ref_date, &
"Fixed reference date to use for nodal modulation of boundary tides.", &
fail_if_missing=.false., default=0)
default=0, do_not_log=tides)

if (.not. OBC%add_eq_phase) then
! If equilibrium phase argument is not added, the input phases
Expand Down
26 changes: 16 additions & 10 deletions src/diagnostics/MOM_harmonic_analysis.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,9 @@ module MOM_harmonic_analysis
time_ref !< Reference time (t = 0) used to calculate tidal forcing
real, dimension(MAX_CONSTITUENTS) :: &
freq, & !< The frequency of a tidal constituent [T-1 ~> s-1]
phase0 !< The phase of a tidal constituent at time 0 [rad]
phase0, & !< The phase of a tidal constituent at time 0 [rad]
tide_fn, & !< Amplitude modulation of tides by nodal cycle [nondim].
tide_un !< Phase modulation of tides by nodal cycle [rad].
real, allocatable :: FtF(:,:) !< Accumulator of (F' * F) for all fields [nondim]
integer :: nc !< The number of tidal constituents in use
integer :: length !< Number of fields of which harmonic analysis is to be performed
Expand All @@ -60,13 +62,15 @@ module MOM_harmonic_analysis

!> This subroutine sets static variables used by this module and initializes CS%list.
!! THIS MUST BE CALLED AT THE END OF tidal_forcing_init.
subroutine HA_init(Time, US, param_file, time_ref, nc, freq, phase0, const_name, CS)
subroutine HA_init(Time, US, param_file, time_ref, nc, freq, phase0, const_name, tide_fn, tide_un, CS)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be more concise if tidal_forcing_CS is used as a single input, rather than its members?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is because tidal_forcing_CS is not transparent and MOM_harmonic_analysis cannot see what's inside. @Hallberg-NOAA commented on this in his review of the initial PR for inline harmonic analysis.

type(time_type), intent(in) :: Time !< The current model time
type(time_type), intent(in) :: time_ref !< Reference time (t = 0) used to calculate tidal forcing
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
real, dimension(MAX_CONSTITUENTS), intent(in) :: freq !< The frequency of a tidal constituent [T-1 ~> s-1]
real, dimension(MAX_CONSTITUENTS), intent(in) :: phase0 !< The phase of a tidal constituent at time 0 [rad]
real, dimension(MAX_CONSTITUENTS), intent(in) :: freq !< The frequency of a tidal constituent [T-1 ~> s-1]
real, dimension(MAX_CONSTITUENTS), intent(in) :: phase0 !< The phase of a tidal constituent at time 0 [rad]
real, dimension(MAX_CONSTITUENTS), intent(in) :: tide_fn !< Amplitude modulation of tides by nodal cycle [nondim].
real, dimension(MAX_CONSTITUENTS), intent(in) :: tide_un !< Phase modulation of tides by nodal cycle [rad].
integer, intent(in) :: nc !< The number of tidal constituents in use
character(len=16), intent(in) :: const_name(MAX_CONSTITUENTS) !< The name of each constituent
type(harmonic_analysis_CS), intent(out) :: CS !< Control structure of the MOM_harmonic_analysis module
Expand Down Expand Up @@ -135,6 +139,8 @@ subroutine HA_init(Time, US, param_file, time_ref, nc, freq, phase0, const_name,
CS%time_ref = time_ref
CS%freq = freq
CS%phase0 = phase0
CS%tide_fn = tide_fn
CS%tide_un = tide_un
CS%nc = nc
CS%const_name = const_name
CS%length = 0
Expand Down Expand Up @@ -196,17 +202,17 @@ subroutine HA_accum_FtF(Time, CS)
do c=1,nc
icos = 2*c
isin = 2*c+1
cosomegat = cos(CS%freq(c) * now + CS%phase0(c))
sinomegat = sin(CS%freq(c) * now + CS%phase0(c))
cosomegat = CS%tide_fn(c) * cos(CS%freq(c) * now + CS%phase0(c) + CS%tide_un(c))
sinomegat = CS%tide_fn(c) * sin(CS%freq(c) * now + CS%phase0(c) + CS%tide_un(c))
CS%FtF(icos,1) = CS%FtF(icos,1) + cosomegat
CS%FtF(isin,1) = CS%FtF(isin,1) + sinomegat
CS%FtF(1,icos) = CS%FtF(icos,1)
CS%FtF(1,isin) = CS%FtF(isin,1)
do cc=c,nc
iccos = 2*cc
issin = 2*cc+1
ccosomegat = cos(CS%freq(cc) * now + CS%phase0(cc))
ssinomegat = sin(CS%freq(cc) * now + CS%phase0(cc))
ccosomegat = CS%tide_fn(cc) * cos(CS%freq(cc) * now + CS%phase0(cc) + CS%tide_un(cc))
ssinomegat = CS%tide_fn(cc) * sin(CS%freq(cc) * now + CS%phase0(cc) + CS%tide_un(cc))
CS%FtF(icos,iccos) = CS%FtF(icos,iccos) + cosomegat * ccosomegat
CS%FtF(icos,issin) = CS%FtF(icos,issin) + cosomegat * ssinomegat
CS%FtF(isin,iccos) = CS%FtF(isin,iccos) + sinomegat * ccosomegat
Expand Down Expand Up @@ -280,8 +286,8 @@ subroutine HA_accum_FtSSH(key, data, Time, G, CS)
do c=1,nc
icos = 2*c
isin = 2*c+1
cosomegat = cos(CS%freq(c) * now + CS%phase0(c))
sinomegat = sin(CS%freq(c) * now + CS%phase0(c))
cosomegat = CS%tide_fn(c) * cos(CS%freq(c) * now + CS%phase0(c) + CS%tide_un(c))
sinomegat = CS%tide_fn(c) * sin(CS%freq(c) * now + CS%phase0(c) + CS%tide_un(c))
do j=js,je ; do i=is,ie
ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j))
ha1%FtSSH(i,j,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat
Expand Down
5 changes: 5 additions & 0 deletions src/diagnostics/MOM_obsolete_params.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,6 +155,11 @@ subroutine find_obsolete_params(param_file)
call obsolete_logical(param_file, "USE_GRID_SPACE_DIAGNOSTIC_AXES", &
hint="Instead use USE_INDEX_DIAGNOSTIC_AXIS.")

call obsolete_logical(param_file, "OBC_TIDE_ADD_EQ_PHASE", hint="Instead use TIDE_USE_EQ_PHASE.")
call obsolete_logical(param_file, "OBC_TIDE_ADD_NODAL", hint="Instead use TIDE_ADD_NODAL.")
call obsolete_int(param_file, "OBC_TIDE_REF_DATE", hint="Instead use TIDE_REF_DATE.")
call obsolete_int(param_file, "OBC_TIDE_NODAL_REF_DATE", hint="Instead use TIDE_NODAL_REF_DATE.")

! Write the file version number to the model log.
call log_version(param_file, mdl, version)

Expand Down
73 changes: 59 additions & 14 deletions src/parameterizations/lateral/MOM_tidal_forcing.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,9 @@ module MOM_tidal_forcing
ampsal(:,:,:), & !< The amplitude of the SAL [Z ~> m].
cosphase_prev(:,:,:), & !< The cosine of the phase of the amphidromes in the previous tidal solutions [nondim].
sinphase_prev(:,:,:), & !< The sine of the phase of the amphidromes in the previous tidal solutions [nondim].
amp_prev(:,:,:) !< The amplitude of the previous tidal solution [Z ~> m].
amp_prev(:,:,:), & !< The amplitude of the previous tidal solution [Z ~> m].
tide_fn(:), & !< Amplitude modulation of tides by nodal cycle [nondim].
tide_un(:) !< Phase modulation of tides by nodal cycle [rad].
end type tidal_forcing_CS

integer :: id_clock_tides !< CPU clock for tides
Expand Down Expand Up @@ -251,9 +253,14 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS)
real, dimension(MAX_CONSTITUENTS) :: amp_def ! Default amplitude for each tidal constituent [m]
real, dimension(MAX_CONSTITUENTS) :: love_def ! Default love number for each constituent [nondim]
integer, dimension(3) :: tide_ref_date !< Reference date (t = 0) for tidal forcing.
integer, dimension(3) :: nodal_ref_date !< Reference date for calculating nodal modulation for tidal forcing.
logical :: use_M2, use_S2, use_N2, use_K2, use_K1, use_O1, use_P1, use_Q1
logical :: use_MF, use_MM
logical :: tides ! True if a tidal forcing is to be used.
logical :: add_nodal_terms = .false. !< If true, insert terms for the 18.6 year modulation when
!! calculating tidal forcing.
type(time_type) :: nodal_time !< Model time to calculate nodal modulation for.
type(astro_longitudes) :: nodal_longitudes !< Solar and lunar longitudes for tidal forcing
logical :: HA_ssh, HA_ubt, HA_vbt
! This include declares and sets the variable "version".
# include "version_variable.h"
Expand Down Expand Up @@ -527,8 +534,46 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS)
enddo
endif

call get_param(param_file, mdl, "TIDE_ADD_NODAL", add_nodal_terms, &
"If true, include 18.6 year nodal modulation in the astronomical tidal forcing.", &
default=.false.)
call get_param(param_file, mdl, "TIDE_NODAL_REF_DATE", nodal_ref_date, &
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reasonable scenario that nodal correction should use a date other than TIDE_REF_DATE?
I appreciate the flexibility for introducing a different reference date for nodal correction, but it would also seem rather confusing.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is what has been done in MOM_open_boundary.F90. I think the reason is that these two dates could be different when tides are forced at the open boundaries by some external data. @andrew-c-ross probably has a better answer.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When I wrote the original OBC tide code I wasn't aware of an easy way for MOM6 to know the current date/time when the OBC and tidal forcing were being initialized. If Time_init is available though I think it would make sense to use that for both the tide and nodal correction ref dates. The only exception would be if, when restarting the model, the init time was something other than the time of the restart, since we would want to update the date used for the nodal correction ~yearly.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the answers!

"Fixed reference date to use for nodal modulation of astronomical tidal forcing.", &
fail_if_missing=.false., default=0)

! If the nodal correction is based on a different time, initialize that.
! Otherwise, it can use N from the time reference.
if (add_nodal_terms) then
if (sum(nodal_ref_date) /= 0) then
! A reference date was provided for the nodal correction
nodal_time = set_date(nodal_ref_date(1), nodal_ref_date(2), nodal_ref_date(3))
call astro_longitudes_init(nodal_time, nodal_longitudes)
elseif (CS%use_eq_phase) then
! Astronomical longitudes were already calculated for use in equilibrium phases,
! so use nodal longitude from that.
nodal_longitudes = CS%tidal_longitudes
else
! Tidal reference time is a required parameter, so calculate the longitudes from that.
call astro_longitudes_init(CS%time_ref, nodal_longitudes)
endif
endif

allocate(CS%tide_fn(nc))
allocate(CS%tide_un(nc))

do c=1,nc
! Find nodal corrections if needed
if (add_nodal_terms) then
call nodal_fu(trim(CS%const_name(c)), nodal_longitudes%N, CS%tide_fn(c), CS%tide_un(c))
else
CS%tide_fn(c) = 1.0
CS%tide_un(c) = 0.0
endif
enddo

if (present(HA_CS)) then
call HA_init(Time, US, param_file, CS%time_ref, CS%nc, CS%freq, CS%phase0, CS%const_name, HA_CS)
call HA_init(Time, US, param_file, CS%time_ref, CS%nc, CS%freq, CS%phase0, CS%const_name, &
CS%tide_fn, CS%tide_un, HA_CS)
call get_param(param_file, mdl, "HA_SSH", HA_ssh, &
"If true, perform harmonic analysis of sea serface height.", default=.false.)
if (HA_ssh) call HA_register('ssh', 'h', HA_CS)
Expand Down Expand Up @@ -613,26 +658,26 @@ subroutine calc_tidal_forcing(Time, e_tide_eq, e_tide_sal, G, US, CS)

do c=1,CS%nc
m = CS%struct(c)
amp_cosomegat = CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c))
amp_sinomegat = CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c))
amp_cosomegat = CS%amp(c)*CS%love_no(c)*CS%tide_fn(c) * cos(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
amp_sinomegat = CS%amp(c)*CS%love_no(c)*CS%tide_fn(c) * sin(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
e_tide_eq(i,j) = e_tide_eq(i,j) + (amp_cosomegat*CS%cos_struct(i,j,m) + &
amp_sinomegat*CS%sin_struct(i,j,m))
enddo ; enddo
enddo

if (CS%use_tidal_sal_file) then ; do c=1,CS%nc
cosomegat = cos(CS%freq(c)*now)
sinomegat = sin(CS%freq(c)*now)
cosomegat = CS%tide_fn(c) * cos(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
sinomegat = CS%tide_fn(c) * sin(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
e_tide_sal(i,j) = e_tide_sal(i,j) + CS%ampsal(i,j,c) * &
(cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c))
enddo ; enddo
enddo ; endif

if (CS%use_tidal_sal_prev) then ; do c=1,CS%nc
cosomegat = cos(CS%freq(c)*now)
sinomegat = sin(CS%freq(c)*now)
cosomegat = CS%tide_fn(c) * cos(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
sinomegat = CS%tide_fn(c) * sin(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
e_tide_sal(i,j) = e_tide_sal(i,j) - CS%sal_scalar * CS%amp_prev(i,j,c) * &
(cosomegat*CS%cosphase_prev(i,j,c) + sinomegat*CS%sinphase_prev(i,j,c))
Expand Down Expand Up @@ -691,8 +736,8 @@ subroutine calc_tidal_forcing_legacy(Time, e_sal, e_sal_tide, e_tide_eq, e_tide_

do c=1,CS%nc
m = CS%struct(c)
amp_cosomegat = CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c))
amp_sinomegat = CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c))
amp_cosomegat = CS%amp(c)*CS%love_no(c)*CS%tide_fn(c) * cos(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
amp_sinomegat = CS%amp(c)*CS%love_no(c)*CS%tide_fn(c) * sin(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
amp_cossin = (amp_cosomegat*CS%cos_struct(i,j,m) + amp_sinomegat*CS%sin_struct(i,j,m))
e_sal_tide(i,j) = e_sal_tide(i,j) + amp_cossin
Expand All @@ -701,8 +746,8 @@ subroutine calc_tidal_forcing_legacy(Time, e_sal, e_sal_tide, e_tide_eq, e_tide_
enddo

if (CS%use_tidal_sal_file) then ; do c=1,CS%nc
cosomegat = cos(CS%freq(c)*now)
sinomegat = sin(CS%freq(c)*now)
cosomegat = CS%tide_fn(c) * cos(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
sinomegat = CS%tide_fn(c) * sin(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
amp_cossin = CS%ampsal(i,j,c) &
* (cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c))
Expand All @@ -712,8 +757,8 @@ subroutine calc_tidal_forcing_legacy(Time, e_sal, e_sal_tide, e_tide_eq, e_tide_
enddo ; endif

if (CS%use_tidal_sal_prev) then ; do c=1,CS%nc
cosomegat = cos(CS%freq(c)*now)
sinomegat = sin(CS%freq(c)*now)
cosomegat = CS%tide_fn(c) * cos(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
sinomegat = CS%tide_fn(c) * sin(CS%freq(c)*now + CS%phase0(c) + CS%tide_un(c))
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
amp_cossin = -CS%sal_scalar * CS%amp_prev(i,j,c) &
* (cosomegat*CS%cosphase_prev(i,j,c) + sinomegat*CS%sinphase_prev(i,j,c))
Expand Down
Loading