+Add support for tidal variations in open boundary condition data #1191
Conversation
# Conflicts: # src/core/MOM_open_boundary.F90
# Conflicts: # src/core/MOM_open_boundary.F90
Codecov Report
@@ Coverage Diff @@
## dev/gfdl #1191 +/- ##
============================================
- Coverage 46.08% 46.00% -0.08%
============================================
Files 214 224 +10
Lines 69399 70822 +1423
============================================
+ Hits 31984 32584 +600
- Misses 37415 38238 +823
Continue to review full report at Codecov.
|
…rence date was not provided
| !! For simplicity, the time associated with time_ref should | ||
| !! be at midnight. These formulas also only make sense if | ||
| !! the calendar is gregorian. | ||
| subroutine astro_longitudes_init(time_ref, longitudes_shpn) |
There was a problem hiding this comment.
We have been thinking that is would be preferable to do our trigonometric calculations in degrees, not radians. Given that this is being newly introduced here, perhaps it would make sense to write this routine to work in degrees. Regardless, the comment describing the longitudes_shpn structure needs to explicit describe their units, preferably using the same syntax (here it would be [rad] or [degrees]) used elsewhere in the MOM6 code.
There was a problem hiding this comment.
I can make the switch from radians to degrees, but there's two concerns I have:
- In addition to the new open boundary code, I would probably want to switch all of MOM_tidal_forcing to using degrees, including changing the tidal frequencies from rad/s to deg/s. Might change answers for any tests with tides=True?
- To switch to degree-based trig functions, I think we'd need to add
-fdec-mathwhen compiling with GNU. Not sure if flags are needed for other compilers.
I'll be happy to make these changes, but just wanted to double check since they may not be entirely trivial.
There was a problem hiding this comment.
Alternatively I could have eq_phase and nodal_fu do calculations with astronomical longitudes in degrees, but convert the result to radians before returning. This would alleviate 1. but not 2.
There was a problem hiding this comment.
Let's not change the frequencies from rad/sec to deg/sec, since this would very likely change answers. We usually put in answer-changing code changes with a flag to allow us to revert to the previous answers, or at least we do the changes in the smallest possible code change.
What would you think, @andrew-c-ross, about having eq_phase and nodal_fu return arguments in degrees that could then be converted (or not) as needed before they are used or stored, noting that these routines only seem to be called in one place?
src/core/MOM_open_boundary.F90
Outdated
| logical :: add_tide_constituents = .false. !< If true, add tidal constituents to the boundary elevation | ||
| !! and velocity. Will be set to true if n_tide_constituents > 0. | ||
| character(len=2), allocatable, dimension(:) :: tide_names !< Names of tidal constituents to add to the boundary data. | ||
| real, allocatable, dimension(:) :: tide_frequencies !< Angular frequencies of chosen tidal constituents. |
There was a problem hiding this comment.
The units of all real elements in public types are supposed to be documented in the comments describing them. In the cast of tide_frequencies, this might be [s-1] or [T-1 ~> s-1], depending on whether these units have been dimensionally rescaled. Similarly the units tide_eq_phases, tide_fn, tide_un, time_ref, and astro_shpn need to be documented here. If they are nondimensional, the convention is to describe their units as [nondim].
src/core/MOM_open_boundary.F90
Outdated
| character(len=200) :: config1 ! String for OBC_USER_CONFIG | ||
| real :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m] | ||
| integer, dimension(3) :: tide_ref_date, nodal_ref_date !< Reference date (t = 0) for tidal forcing |
There was a problem hiding this comment.
These dates should probably be stored using the standard time_type, rather than inventing a new 3-integer type.
There was a problem hiding this comment.
Is there a cleaner way you would recommend for reading in a date from MOM_input? My current strategy is to set the date in MOM_input with
OBC_TIDE_NODAL_REF_DATE = 1981,7,2
and then the subroutine get_param reads this in as a 3-integer array, which later gets does get converted to a time_type.
I could do the conversion to time_type immediately in this subroutine rather than passing the 3-int array over to initialize_obc_tides for conversion if that would make it cleaner.
src/core/MOM_open_boundary.F90
Outdated
| @@ -828,6 +987,66 @@ subroutine initialize_segment_data(G, OBC, PF) | |||
|
|
|||
| end subroutine initialize_segment_data | |||
|
|
|||
| subroutine initialize_obc_tides(OBC, tide_ref_date, nodal_ref_date, tide_constituent_str) | |||
| type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure | |||
| integer, dimension(3), intent(in) :: tide_ref_date !< Reference date (t = 0) for tidal forcing. | |||
There was a problem hiding this comment.
All of these date arguments should follow the convention throughout the rest of MOM6 and FMS models more generally and be cast as time_types.
src/core/MOM_open_boundary.F90
Outdated
| elseif (OBC%add_eq_phase) then | ||
| ! Astronomical longitudes were already calculated for use in equilibrium phases, | ||
| ! so use nodal longitude from that. | ||
| nodal_shpn = OBC%astro_shpn |
There was a problem hiding this comment.
This is an array syntax copy, which needs to use syntax that is consistent with what is happening here, namely this should be nodal_shpn(:) = OBC%astro_shpn(:). Also since this is an array of 4 different longitudes, perhaps this would be clearer if it were turned into a type with named elements and comments describing each one.
src/core/MOM_open_boundary.F90
Outdated
| allocate(OBC%tide_un(OBC%n_tide_constituents)) | ||
|
|
||
| do c=1,OBC%n_tide_constituents | ||
| OBC%tide_frequencies(c) = tidal_frequency(trim(OBC%tide_names(c))) |
There was a problem hiding this comment.
When tidal_frequency is called in tidal_forcing_init, this sets the default frequencies of the tidal constituents, but the user can override these at run-time. In the OBC implementation, there is not the ability to set these frequencies at run-time, and hence these are hard-coded. This would preclude using MOM6 for some paleoclimate applications, for example. I would recommend modifying tidal_frequency to also take a param_file argument and allow the user to override these frequencies at run-time.
There was a problem hiding this comment.
Would it be reasonable to have the frequency overrides for the body forcing also apply to the open boundary forcing? In other words the open boundary forcing would also use the frequency in TIDE_*_FREQ if it was specified. I can see pros and cons to doing it this way vs other ways, but this seems like a fairly clean and simple way of handling a niche case.
There was a problem hiding this comment.
What you propose does seem like an elegant solution. Certainly one would not want to be using different frequencies for the astronomical forcing and the boundary forcing, and this solution would avoid that.
There was a problem hiding this comment.
Overall these changes seem to properly address the needs of introducing tides as boundary forcing terms, which will be a valuable contribution to the MOM6 code. Thank you, @andrew-c-ross!
Given that changing interfaces or answers after the fact can be a challenge, there are a few things here that I would suggest should be considered with this pull request.
-
Please consider using the time_type used elsewhere in the MOM6 code for date and time variables.
-
Please document the units of all real variables where they are defined or passed in as arguments, preferably using the same syntax as in the rest of the MOM6 code.
-
Please do not use array-syntax arithmetic on arrays, and if you use array syntax for copies be sure to use syntax like
array1(:) = array2(:)to highlight what is being done. Please see https://github.com/NOAA-GFDL/MOM6/wiki/Code-style-guide for guidance. -
Please consider storing longitudes in degrees, not radians, as we have plans to convert most of the trigonometric functions to use degrees, and converting from degrees to radians and then back again introduces roundoff-level inaccuracies.
-
Some variables, such as the tidal frequencies used for the OBCs, have hard-coded values, and there should be the ability to use these values as defaults that can be overriden at run-time, as is the case in the MOM_tidal_forcing.F90 code.
There are specific comments describing each of these issues at various points in the code.
|
Thank you for the review, @Hallberg-NOAA . I will work on addressing these issues next week, and will let you know with a comment here once I think I have all of the issues addressed. My reasoning for storing the reference times as reals rather than time_types was that the reference times never change values or calendars once initialized, so I thought it would be more efficient to store and do calculations with reals. But I also realize that there are benefits to consistently using a time_type for them and so I will switch to using time_type. |
|
The latest version of this PR is being tested with https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/11310. |
Overview
This PR adds support for tidal variations in the open boundary conditions for velocity and elevation. In addition, some new subroutines needed for the boundary conditions are also used to improve the existing tidal forcing in MOM_tidal_forcing.F90. By default, all of these new options are disabled, and previous answers will not change.
Strategy for adding tidal variations
Tidal variations in the boundary data are added during each call to
update_OBC_segment_data. Boundary data updates in this subroutine generally proceed as they did before. However, when updating a field with tidal variation (U, V, or SSH), respective tidal fields are also read (e.g., Uamp and Uphase when updating U), the values of each harmonic constituent are calculated for the current model time, and the summed tidal forcing is added on top of the U, V, or SSH boundary data.Most parameters related to the boundary tidal forcing are attached to the open boundary control structure, rather than to each OBC segment, under the assumption that in most cases these parameters will not vary by segment.
New MOM_input parameters
OBC_TIDE_N_CONSTITUENTS: Sets the number of tidal constituents being added. By default, it is 0, in which case tidal variations are not included in the boundary conditions and none of the following
OBC_TIDE_*parameters will be used.OBC_TIDE_CONSTITUENTS: A list of the names of tidal constituents that will be added to the boundary. The names should be in the same order as they are in the values/files in
OBC_SEGMENT_*_DATA.OBC_TIDE_ADD_EQ_PHASE: If true, the equilibrium phase argument is included during the tidal calculations. The equilibrium phase argument gives the phase of the equilibrium tide on the Greenwich meridian at
time_ref=0. The tidal phases provided as input in the segment data should then be the local phase lag relative to the equilibrium tide at Greenwich. If false (default), the tidal phases provided as input should be the local phase lag relative totime_ref=0. Tide model products such as TPXO generally provide phases in the first form (Greenwich lag), soOBC_TIDE_ADD_EQ_PHASE = TRUEwill be a common usage.OBC_TIDE_REF_DATE: Year, month, and day of the tidal reference date (sets when
time_ref=0). Generally this will be a date close to the start of the model run.OBC_TIDE_ADD_NODAL: If true, the tidal boundary conditions are modulated by the 18.6 year nodal cycle. Defaults to false.
OBC_TIDE_NODAL_REF_DATE: Year, month, and day of the date to calculate the terms for the 18.6 year nodal modulation. In a typical usage, this option will be changed once per model year. Alternatively this could be done during every call to
update_OBC_segment_data, but the nodal modulation is slow and that would be excessive.TIDE_REF_DATE: Similar to
OBC_TIDE_REF_DATE, except it setstime_ref=0 for the existing forcing in MOM_tidal_forcing.F90. IfTIDE_REF_DATEis not specified, the reference date is set to 0 and previous answers are reproduced.TIDE_USE_EQ_PHASE: Similar to
OBC_TIDE_ADD_EQ_PHASE. If true, the equilibrium phase argument is included in the tidal forcing calculations. This should be set to true, and the calendar set to Gregorian, to get model tides that are in phase with real-world tides. If false (default), previous answers are reproduced.Modifications to MOM_input parameters
OBC_SEGMENT_*_DATA now includes six variables for tidal variation (Uamp, Uphase, Vamp, Vphase, SSHamp, SSHphase). If
OBC_TIDE_N_CONSTITUENTS> 0, all six variables are required to be provided. Like other variables, these can be specified by either value or file. If a file is provided, the file should have time as the first dimension. Time can be length-1 to keep the tidal amplitudes and phases constant in time (common), or time variation could be used as an alternative method of implementing nodal and other variations. The second dimension of the file is expected to be “constituent”, and the order of the constituents should matchOBC_TIDE_CONSTITUENTS. The remaining dimensions should be consistent with those for U, V, or SSH.New subroutines and functions
Added to MOM_tidal_forcing.F90
subroutine astro_longitudes_init: Calculates longitudes of the moon, sun, lunar perigee, and ascending node at a specified reference time. These longitudes are used by
eq_phaseandnodal_fu, below, to determine the equilibrium phase arguments and 18.6 year nodal modulation terms. There are a number of different equations for calculating these longitudes. The equation I use here is from Kowalik and Luick, "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019, who give equations based on Schureman, “Manual of Harmonic Analysis and Prediction of Tides”, 1958. Note that the equations use a measure of time since 1900-01-01 and are thus only accurate when following the Gregorian calendar.function eq_phase: Calculates the equilibrium phase argument for a given tidal constituent and reference time. This function is used to set the equilibrium phases for both the open boundary and body tidal forcing, if the respective options are enabled.
subroutine nodal_fu: Calculates the terms for amplitude and phase modulation by the 18.6 year nodal cycle for a given constituent.
function tidal_frequency: This function does a simple lookup of the angular frequency of a tidal constituent. The frequencies are the same as those previously used by MOM_tidal_forcing. Since both the body and boundary forcing will need access to these frequencies, I have placed them in a function.
Added to MOM_open_boundary.F90
subroutine initialize_obc_tides: Reads the MOM_input parameters related to OBC tides, allocates needed arrays, and calls the subroutines to calculate equilibrium phases and nodal modulation.
Other changes
Previously, if the boundary data for velocity was set by value rather than by file (e.g., U=value:2.0), the boundary velocity would erroneously be unset and remain at 0. When setting by value,
update_OBC_segment_datawould setsegment%field(m)%bt_vel(:,:) = segment%field(m)%value, which as far as I can tell is unused, and would not setnormal_trans,normal_vel,normal_vel_bt, etc. I have deleted all uses ofbt_veland instead havenormal_transand other variables updated the same way whether the data is from value or file.Because the equilibrium phase calculations are only meaningful under the Gregorian calendar, coupler_main.F90 in MOM6-examples/src/coupler will also need to be updated. My PR to FMScoupler to add the Gregorian calendar to this file has already been merged (NOAA-GFDL/FMScoupler@90098b1), so this should be easy.
Tidal variations are updated during each
update_OBC_segment_datacall, soupdate_OBC_datain MOM_boundary_update.F90 now checks if tides are being added and will callupdate_OBC_segment_dataeven if no other variables need updating.Several string lengths have been increased so that longer strings can be used for
OBC_SEGMENT_*_DATA.Possible future changes and improvements
OBC_TIDE_ADD_EQ_PHASEorTIDE_USE_EQ_PHASEare enabled but the calendar is not Gregorian; I have not added this yet since most of MOM is agnostic about the calendar.