From 73514f2dabfbf257b157eee182d52b27980942f4 Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 6 Sep 2024 19:35:15 -0300 Subject: [PATCH] Nodal modulation This commit fixes a few (potential) inconsistencies between open boundary tidal forcing and astronomical tidal forcing. 1. There was an inconsistency in the code that the nodal modulation can be calculated in OBC tidal forcing but not in the astronomical tidal forcing. This commit fixes this bug so that nodal modulation is applied to both forcings. 2. In the previous version of MOM_open_boundary.F90, a different set of tidal parameters can be set for open boundary tidal forcing, leading to potential inconsistency with astronomical tidal forcing. This commit obsoletes those for open boundary tidal forcing. 3. Another important bug fix is that the equilibrium phase is added to the SAL term, which was missing in the original code. --- src/core/MOM_open_boundary.F90 | 20 +++-- src/diagnostics/MOM_harmonic_analysis.F90 | 26 ++++--- .../lateral/MOM_tidal_forcing.F90 | 77 +++++++++++++++---- 3 files changed, 89 insertions(+), 34 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index fb833189ec..797f60bd9b 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1169,26 +1169,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.) + old_name="OBC_TIDE_ADD_EQ_PHASE", 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.) + old_name="OBC_TIDE_ADD_NODAL", 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.) + old_name="OBC_TIDE_REF_DATE", defaults=(/0, 0, 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., defaults=(/0, 0, 0/)) + old_name="OBC_TIDE_NODAL_REF_DATE", defaults=(/0, 0, 0/), do_not_log=tides) if (.not. OBC%add_eq_phase) then ! If equilibrium phase argument is not added, the input phases diff --git a/src/diagnostics/MOM_harmonic_analysis.F90 b/src/diagnostics/MOM_harmonic_analysis.F90 index 0013a8ab83..f2585d510a 100644 --- a/src/diagnostics/MOM_harmonic_analysis.F90 +++ b/src/diagnostics/MOM_harmonic_analysis.F90 @@ -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 @@ -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) 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, intent(in) :: freq(MAX_CONSTITUENTS) !< The frequency of a tidal constituent [T-1 ~> s-1] + real, intent(in) :: phase0(MAX_CONSTITUENTS) !< The phase of a tidal constituent at time 0 [rad] + real, intent(in) :: tide_fn(MAX_CONSTITUENTS) !< Amplitude modulation of tides by nodal cycle [nondim]. + real, intent(in) :: tide_un(MAX_CONSTITUENTS) !< 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 @@ -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 @@ -198,8 +204,8 @@ 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))) ! First column, corresponding to the zero frequency constituent (mean) CS%FtF(icos,1) = CS%FtF(icos,1) + cosomegat @@ -208,8 +214,8 @@ subroutine HA_accum_FtF(Time, CS) do cc=1,c 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))) ! Interior of the matrix, corresponding to the products of cosine and sine terms CS%FtF(icos,iccos) = CS%FtF(icos,iccos) + cosomegat * ccosomegat @@ -290,8 +296,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,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat ha1%FtSSH(i,j,isin) = ha1%FtSSH(i,j,isin) + (data(i,j) - ha1%ref(i,j)) * sinomegat diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 85c9b1ee81..5e624d4aea 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -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 @@ -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" @@ -385,11 +392,11 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) call get_param(param_file, mdl, "TIDE_REF_DATE", tide_ref_date, & "Year,month,day to use as reference date for tidal forcing. "//& "If not specified, defaults to 0.", & - defaults=(/0, 0, 0/)) + old_name="OBC_TIDE_REF_DATE", defaults=(/0, 0, 0/)) call get_param(param_file, mdl, "TIDE_USE_EQ_PHASE", CS%use_eq_phase, & "Correct phases by calculating equilibrium phase arguments for TIDE_REF_DATE. ", & - default=.false., fail_if_missing=.false.) + old_name="OBC_TIDE_ADD_EQ_PHASE", default=.false., fail_if_missing=.false.) if (sum(tide_ref_date) == 0) then ! tide_ref_date defaults to 0. CS%time_ref = set_date(1, 1, 1, 0, 0, 0) @@ -528,8 +535,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.", & + old_name="OBC_TIDE_ADD_NODAL", default=.false.) + call get_param(param_file, mdl, "TIDE_NODAL_REF_DATE", nodal_ref_date, & + "Fixed reference date to use for nodal modulation of astronomical tidal forcing.", & + old_name="OBC_TIDE_REF_DATE", fail_if_missing=.false., defaults=(/0, 0, 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) @@ -614,8 +659,8 @@ 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)) @@ -623,8 +668,8 @@ subroutine calc_tidal_forcing(Time, e_tide_eq, e_tide_sal, G, US, CS) 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)) @@ -632,8 +677,8 @@ subroutine calc_tidal_forcing(Time, e_tide_eq, e_tide_sal, G, US, CS) 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)) @@ -692,8 +737,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 @@ -702,8 +747,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)) @@ -713,8 +758,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))