MOM_tidal_forcing.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Tidal contributions to geopotential
7
8use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, &
9 clock_module, clock_routine
10use mom_domains, only : pass_var
11use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
12use mom_file_parser, only : get_param, log_version, param_file_type
13use mom_grid, only : ocean_grid_type
14use mom_io, only : field_exists, file_exists, mom_read_data
15use mom_time_manager, only : set_date, time_type, time_minus_signed
17
18implicit none ; private
19
22! MOM_open_boundary uses the following to set tides on the boundary.
24
25#include <MOM_memory.h>
26
27integer, parameter :: max_constituents = 10 !< The maximum number of tidal
28 !! constituents that could be used.
29!> Simple type to store astronomical longitudes used to calculate tidal phases.
30type, public :: astro_longitudes
31 real :: s !< Mean longitude of moon [rad]
32 real :: h !< Mean longitude of sun [rad]
33 real :: p !< Mean longitude of lunar perigee [rad]
34 real :: n !< Longitude of ascending node [rad]
35end type astro_longitudes
36
37!> The control structure for the MOM_tidal_forcing module
38type, public :: tidal_forcing_cs ; private
39 logical :: use_tidal_sal_file !< If true, Read the tidal self-attraction
40 !! and loading from input files, specified
41 !! by TIDAL_INPUT_FILE.
42 logical :: use_tidal_sal_prev !< If true, use the SAL from the previous
43 !! iteration of the tides to facilitate convergence.
44 logical :: use_eq_phase !< If true, tidal forcing is phase-shifted to match
45 !! equilibrium tide. Set to false if providing tidal phases
46 !! that have already been shifted by the
47 !! astronomical/equilibrium argument.
48 real :: sal_scalar = 0.0 !< The constant of proportionality between self-attraction and
49 !! loading (SAL) geopotential anomaly and total geopotential geopotential
50 !! anomalies. This is only used if USE_PREVIOUS_TIDES is true. [nondim].
51 integer :: nc !< The number of tidal constituents in use.
52 real, dimension(MAX_CONSTITUENTS) :: &
53 freq, & !< The frequency of a tidal constituent [rad T-1 ~> rad s-1].
54 phase0, & !< The phase of a tidal constituent at time 0 [rad].
55 amp, & !< The amplitude of a tidal constituent at time 0 [Z ~> m].
56 love_no !< The Love number of a tidal constituent at time 0 [nondim].
57 integer :: struct(max_constituents) !< An encoded spatial structure for each constituent
58 character (len=16) :: const_name(max_constituents) !< The name of each constituent
59
60 type(time_type) :: time_ref !< Reference time (t = 0) used to calculate tidal forcing.
61 type(astro_longitudes) :: tidal_longitudes !< Astronomical longitudes used to calculate
62 !! tidal phases at t = 0.
63 real, allocatable :: &
64 sin_struct(:,:,:), & !< The sine based structures that can be associated with
65 !! the astronomical forcing [nondim].
66 cos_struct(:,:,:), & !< The cosine based structures that can be associated with
67 !! the astronomical forcing [nondim].
68 cosphasesal(:,:,:), & !< The cosine of the phase of the self-attraction and loading amphidromes [nondim].
69 sinphasesal(:,:,:), & !< The sine of the phase of the self-attraction and loading amphidromes [nondim].
70 ampsal(:,:,:), & !< The amplitude of the SAL [Z ~> m].
71 cosphase_prev(:,:,:), & !< The cosine of the phase of the amphidromes in the previous tidal solutions [nondim].
72 sinphase_prev(:,:,:), & !< The sine of the phase of the amphidromes in the previous tidal solutions [nondim].
73 amp_prev(:,:,:), & !< The amplitude of the previous tidal solution [Z ~> m].
74 tide_fn(:), & !< Amplitude modulation of tides by nodal cycle [nondim].
75 tide_un(:) !< Phase modulation of tides by nodal cycle [rad].
76end type tidal_forcing_cs
77
78integer :: id_clock_tides !< CPU clock for tides
79
80contains
81
82!> Finds astronomical longitudes s, h, p, and N,
83!! the mean longitude of the moon, sun, lunar perigee, and ascending node, respectively,
84!! at the specified reference time time_ref.
85!! These formulas were obtained from
86!! Kowalik and Luick, "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019
87!! (their Equation I.71), which are based on Schureman, 1958.
88!! For simplicity, the time associated with time_ref should
89!! be at midnight. These formulas also only make sense if
90!! the calendar is Gregorian.
91subroutine astro_longitudes_init(time_ref, longitudes)
92 type(time_type), intent(in) :: time_ref !> Time to calculate longitudes for.
93 type(astro_longitudes), intent(out) :: longitudes !> Lunar and solar longitudes at time_ref.
94
95 ! Local variables
96 real :: d !> Time since the reference date [days]
97 real :: t !> Time in Julian centuries [centuries]
98 real, parameter :: pi = 4.0 * atan(1.0) !> 3.14159... [nondim]
99
100 ! Find date at time_ref in days since midnight at the start of 1900-01-01
101 d = time_minus_signed(time_ref, set_date(1900, 1, 1, 0, 0, 0)) / (24.0 * 3600.0)
102 ! Time since 1900-01-01 in Julian centuries
103 ! Kowalik and Luick use 36526, but Schureman uses 36525 which I think is correct.
104 t = d / 36525.0
105 ! Calculate longitudes, including converting to radians on [0, 2pi)
106 ! s: Mean longitude of moon
107 longitudes%s = mod((277.0248 + 481267.8906 * t) + 0.0011 * (t**2), 360.0) * pi / 180.0
108 ! h: Mean longitude of sun
109 longitudes%h = mod((280.1895 + 36000.7689 * t) + 3.0310e-4 * (t**2), 360.0) * pi / 180.0
110 ! p: Mean longitude of lunar perigee
111 longitudes%p = mod((334.3853 + 4069.0340 * t) - 0.0103 * (t**2), 360.0) * pi / 180.0
112 ! n: Longitude of ascending node
113 longitudes%N = mod((259.1568 - 1934.142 * t) + 0.0021 * (t**2), 360.0) * pi / 180.0
114end subroutine astro_longitudes_init
115
116!> Calculates the equilibrium phase argument for the given tidal
117!! constituent constit and the astronomical longitudes and the reference time.
118!! These formulas follow Table I.4 of Kowalik and Luick,
119!! "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019.
120function eq_phase(constit, longitudes)
121 character (len=2), intent(in) :: constit !> Name of constituent (e.g., M2).
122 type(astro_longitudes), intent(in) :: longitudes !> Mean longitudes calculated using astro_longitudes_init
123 real, parameter :: pi = 4.0 * atan(1.0) !> 3.14159... [nondim]
124 real :: eq_phase !> The equilibrium phase argument for the constituent [rad].
125
126 select case (constit)
127 case ("M2")
128 eq_phase = 2 * (longitudes%h - longitudes%s)
129 case ("S2")
130 eq_phase = 0.0
131 case ("N2")
132 eq_phase = (- 3 * longitudes%s + 2 * longitudes%h) + longitudes%p
133 case ("K2")
134 eq_phase = 2 * longitudes%h
135 case ("K1")
136 eq_phase = longitudes%h + pi / 2.0
137 case ("O1")
138 eq_phase = (- 2 * longitudes%s + longitudes%h) - pi / 2.0
139 case ("P1")
140 eq_phase = - longitudes%h - pi / 2.0
141 case ("Q1")
142 eq_phase = ((- 3 * longitudes%s + longitudes%h) + longitudes%p) - pi / 2.0
143 case ("MF")
144 eq_phase = 2 * longitudes%s
145 case ("MM")
146 eq_phase = longitudes%s - longitudes%p
147 case default
148 call mom_error(fatal, "eq_phase: unrecognized constituent")
149 end select
150end function eq_phase
151
152!> Looks up angular frequencies for the main tidal constituents.
153!! Values used here are from previous versions of MOM.
154function tidal_frequency(constit)
155 character (len=2), intent(in) :: constit !> Constituent to look up
156 real :: tidal_frequency !> Angular frequency [rad s-1]
157
158 select case (constit)
159 case ("M2")
160 tidal_frequency = 1.4051890e-4
161 case ("S2")
162 tidal_frequency = 1.4544410e-4
163 case ("N2")
164 tidal_frequency = 1.3787970e-4
165 case ("K2")
166 tidal_frequency = 1.4584234e-4
167 case ("K1")
168 tidal_frequency = 0.7292117e-4
169 case ("O1")
170 tidal_frequency = 0.6759774e-4
171 case ("P1")
172 tidal_frequency = 0.7252295e-4
173 case ("Q1")
174 tidal_frequency = 0.6495854e-4
175 case ("MF")
176 tidal_frequency = 0.053234e-4
177 case ("MM")
178 tidal_frequency = 0.026392e-4
179 case default
180 call mom_error(fatal, "tidal_frequency: unrecognized constituent")
181 end select
182end function tidal_frequency
183
184!> Find amplitude (f) and phase (u) modulation of tidal constituents by the 18.6
185!! year nodal cycle. Values here follow Table I.6 in Kowalik and Luick,
186!! "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019.
187subroutine nodal_fu(constit, nodelon, fn, un)
188 character (len=2), intent(in) :: constit !> Tidal constituent to find modulation for.
189 real, intent(in) :: nodelon !> Longitude of ascending node [rad], which
190 !! can be calculated using astro_longitudes_init.
191 real, intent(out) :: fn !> Amplitude modulation [nondim]
192 real, intent(out) :: un !> Phase modulation [rad]
193
194 real, parameter :: radians = 4.0 * atan(1.0) / 180.0 !> Converts degrees to radians [nondim]
195
196 select case (constit)
197 case ("M2")
198 fn = 1.0 - 0.037 * cos(nodelon)
199 un = -2.1 * radians * sin(nodelon)
200 case ("S2")
201 fn = 1.0 ! Solar S2 has no amplitude modulation.
202 un = 0.0 ! S2 has no phase modulation.
203 case ("N2")
204 fn = 1.0 - 0.037 * cos(nodelon)
205 un = -2.1 * radians * sin(nodelon)
206 case ("K2")
207 fn = 1.024 + 0.286 * cos(nodelon)
208 un = -17.7 * radians * sin(nodelon)
209 case ("K1")
210 fn = 1.006 + 0.115 * cos(nodelon)
211 un = -8.9 * radians * sin(nodelon)
212 case ("O1")
213 fn = 1.009 + 0.187 * cos(nodelon)
214 un = 10.8 * radians * sin(nodelon)
215 case ("P1")
216 fn = 1.0 ! P1 has no amplitude modulation.
217 un = 0.0 ! P1 has no phase modulation.
218 case ("Q1")
219 fn = 1.009 + 0.187 * cos(nodelon)
220 un = 10.8 * radians * sin(nodelon)
221 case ("MF")
222 fn = 1.043 + 0.414 * cos(nodelon)
223 un = -23.7 * radians * sin(nodelon)
224 case ("MM")
225 fn = 1.0 - 0.130 * cos(nodelon)
226 un = 0.0 ! MM has no phase modulation.
227 case default
228 call mom_error(fatal, "nodal_fu: unrecognized constituent")
229 end select
230
231end subroutine nodal_fu
232
233!> This subroutine allocates space for the static variables used
234!! by this module. The metrics may be effectively 0, 1, or 2-D arrays,
235!! while fields like the background viscosities are 2-D arrays.
236!! ALLOC is a macro defined in MOM_memory.h for allocate or nothing with
237!! static memory.
238subroutine tidal_forcing_init(Time, G, US, param_file, CS)
239 type(time_type), intent(in) :: time !< The current model time.
240 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
241 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
242 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
243 type(tidal_forcing_cs), intent(inout) :: cs !< Tidal forcing control structure
244
245 ! Local variables
246 real, dimension(SZI_(G), SZJ_(G)) :: &
247 phase, & ! The phase of some tidal constituent [radians].
248 lat_rad, lon_rad ! Latitudes and longitudes of h-points [radians].
249 real :: deg_to_rad ! A conversion factor from degrees to radians [radian degree-1]
250 real, dimension(MAX_CONSTITUENTS) :: freq_def ! Default frequency for each tidal constituent [rad s-1]
251 real, dimension(MAX_CONSTITUENTS) :: phase0_def ! Default reference phase for each tidal constituent [rad]
252 real, dimension(MAX_CONSTITUENTS) :: amp_def ! Default amplitude for each tidal constituent [m]
253 real, dimension(MAX_CONSTITUENTS) :: love_def ! Default love number for each constituent [nondim]
254 integer, dimension(3) :: tide_ref_date !< Reference date (t = 0) for tidal forcing.
255 integer, dimension(3) :: nodal_ref_date !< Reference date for calculating nodal modulation for tidal forcing.
256 logical :: use_m2, use_s2, use_n2, use_k2, use_k1, use_o1, use_p1, use_q1
257 logical :: use_mf, use_mm
258 logical :: tides ! True if a tidal forcing is to be used.
259 logical :: add_nodal_terms = .false. !< If true, insert terms for the 18.6 year modulation when
260 !! calculating tidal forcing.
261 type(time_type) :: nodal_time !< Model time to calculate nodal modulation for.
262 type(astro_longitudes) :: nodal_longitudes !< Solar and lunar longitudes for tidal forcing
263 ! This include declares and sets the variable "version".
264# include "version_variable.h"
265 character(len=40) :: mdl = "MOM_tidal_forcing" ! This module's name.
266 character(len=128) :: mesg
267 character(len=200) :: tidal_input_files(4*max_constituents)
268 integer :: i, j, c, is, ie, js, je, isd, ied, jsd, jed, nc
269
270 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
271 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
272
273 ! Read all relevant parameters and write them to the model log.
274 call log_version(param_file, mdl, version, "")
275 call get_param(param_file, mdl, "TIDES", tides, &
276 "If true, apply tidal momentum forcing.", default=.false.)
277
278 if (.not.tides) return
279
280 ! Set up the spatial structure functions for the diurnal, semidiurnal, and
281 ! low-frequency tidal components.
282 allocate(cs%sin_struct(isd:ied,jsd:jed,3), source=0.0)
283 allocate(cs%cos_struct(isd:ied,jsd:jed,3), source=0.0)
284 deg_to_rad = 4.0*atan(1.0)/180.0
285 do j=js-1,je+1 ; do i=is-1,ie+1
286 lat_rad(i,j) = g%geoLatT(i,j)*deg_to_rad
287 lon_rad(i,j) = g%geoLonT(i,j)*deg_to_rad
288 enddo ; enddo
289 do j=js-1,je+1 ; do i=is-1,ie+1
290 cs%sin_struct(i,j,1) = -sin(2.0*lat_rad(i,j)) * sin(lon_rad(i,j))
291 cs%cos_struct(i,j,1) = sin(2.0*lat_rad(i,j)) * cos(lon_rad(i,j))
292 cs%sin_struct(i,j,2) = -cos(lat_rad(i,j))**2 * sin(2.0*lon_rad(i,j))
293 cs%cos_struct(i,j,2) = cos(lat_rad(i,j))**2 * cos(2.0*lon_rad(i,j))
294 cs%sin_struct(i,j,3) = 0.0
295 cs%cos_struct(i,j,3) = (0.5-1.5*sin(lat_rad(i,j))**2)
296 enddo ; enddo
297
298 call get_param(param_file, mdl, "TIDE_M2", use_m2, &
299 "If true, apply tidal momentum forcing at the M2 "//&
300 "frequency. This is only used if TIDES is true.", &
301 default=.false.)
302 call get_param(param_file, mdl, "TIDE_S2", use_s2, &
303 "If true, apply tidal momentum forcing at the S2 "//&
304 "frequency. This is only used if TIDES is true.", &
305 default=.false.)
306 call get_param(param_file, mdl, "TIDE_N2", use_n2, &
307 "If true, apply tidal momentum forcing at the N2 "//&
308 "frequency. This is only used if TIDES is true.", &
309 default=.false.)
310 call get_param(param_file, mdl, "TIDE_K2", use_k2, &
311 "If true, apply tidal momentum forcing at the K2 "//&
312 "frequency. This is only used if TIDES is true.", &
313 default=.false.)
314 call get_param(param_file, mdl, "TIDE_K1", use_k1, &
315 "If true, apply tidal momentum forcing at the K1 "//&
316 "frequency. This is only used if TIDES is true.", &
317 default=.false.)
318 call get_param(param_file, mdl, "TIDE_O1", use_o1, &
319 "If true, apply tidal momentum forcing at the O1 "//&
320 "frequency. This is only used if TIDES is true.", &
321 default=.false.)
322 call get_param(param_file, mdl, "TIDE_P1", use_p1, &
323 "If true, apply tidal momentum forcing at the P1 "//&
324 "frequency. This is only used if TIDES is true.", &
325 default=.false.)
326 call get_param(param_file, mdl, "TIDE_Q1", use_q1, &
327 "If true, apply tidal momentum forcing at the Q1 "//&
328 "frequency. This is only used if TIDES is true.", &
329 default=.false.)
330 call get_param(param_file, mdl, "TIDE_MF", use_mf, &
331 "If true, apply tidal momentum forcing at the MF "//&
332 "frequency. This is only used if TIDES is true.", &
333 default=.false.)
334 call get_param(param_file, mdl, "TIDE_MM", use_mm, &
335 "If true, apply tidal momentum forcing at the MM "//&
336 "frequency. This is only used if TIDES is true.", &
337 default=.false.)
338
339 ! Determine how many tidal components are to be used.
340 nc = 0
341 if (use_m2) nc=nc+1 ; if (use_s2) nc=nc+1
342 if (use_n2) nc=nc+1 ; if (use_k2) nc=nc+1
343 if (use_k1) nc=nc+1 ; if (use_o1) nc=nc+1
344 if (use_p1) nc=nc+1 ; if (use_q1) nc=nc+1
345 if (use_mf) nc=nc+1 ; if (use_mm) nc=nc+1
346 cs%nc = nc
347
348 if (nc == 0) then
349 call mom_error(fatal, "tidal_forcing_init: "// &
350 "TIDES are defined, but no tidal constituents are used.")
351 return
352 endif
353
354 call get_param(param_file, mdl, "TIDAL_SAL_FROM_FILE", cs%use_tidal_sal_file, &
355 "If true, read the tidal self-attraction and loading "//&
356 "from input files, specified by TIDAL_INPUT_FILE. "//&
357 "This is only used if TIDES is true.", default=.false.)
358 call get_param(param_file, mdl, "USE_PREVIOUS_TIDES", cs%use_tidal_sal_prev, &
359 "If true, use the SAL from the previous iteration of the "//&
360 "tides to facilitate convergent iteration. "//&
361 "This is only used if TIDES is true.", default=.false.)
362 if (cs%use_tidal_sal_prev) &
363 call get_param(param_file, mdl, "SAL_SCALAR_VALUE", cs%sal_scalar, "The constant of "//&
364 "proportionality between self-attraction and loading (SAL) geopotential "//&
365 "anomaly and barotropic geopotential anomalies. This is only used if "//&
366 "SAL_SCALAR_APPROX is true or USE_PREVIOUS_TIDES is true.", default=0.0, &
367 units="m m-1", do_not_log=(.not.cs%use_tidal_sal_prev), &
368 old_name='TIDE_SAL_SCALAR_VALUE')
369
370 if (nc > max_constituents) then
371 write(mesg,'("Increase MAX_CONSTITUENTS in MOM_tidal_forcing.F90 to at least ",I0, &
372 &" to accommodate all the registered tidal constituents.")') nc
373 call mom_error(fatal, "MOM_tidal_forcing"//mesg)
374 endif
375
376 do c=1,4*max_constituents ; tidal_input_files(c) = "" ; enddo
377
378 if (cs%use_tidal_sal_file .or. cs%use_tidal_sal_prev) then
379 call get_param(param_file, mdl, "TIDAL_INPUT_FILE", tidal_input_files, &
380 "A list of input files for tidal information.", &
381 default="", fail_if_missing=.true.)
382 endif
383
384 call get_param(param_file, mdl, "TIDE_REF_DATE", tide_ref_date, &
385 "Year,month,day to use as reference date for tidal forcing. "//&
386 "If not specified, defaults to 0.", &
387 old_name="OBC_TIDE_REF_DATE", defaults=(/0, 0, 0/))
388
389 call get_param(param_file, mdl, "TIDE_USE_EQ_PHASE", cs%use_eq_phase, &
390 "Correct phases by calculating equilibrium phase arguments for TIDE_REF_DATE. ", &
391 old_name="OBC_TIDE_ADD_EQ_PHASE", default=.false., fail_if_missing=.false.)
392
393 if (sum(tide_ref_date) == 0) then ! tide_ref_date defaults to 0.
394 cs%time_ref = set_date(1, 1, 1, 0, 0, 0)
395 else
396 if (.not. cs%use_eq_phase) then
397 ! Using a reference date but not using phase relative to equilibrium.
398 ! This makes sense as long as either phases are overridden, or
399 ! correctly simulating tidal phases is not desired.
400 call mom_mesg('Tidal phases will *not* be corrected with equilibrium arguments.')
401 endif
402 cs%time_ref = set_date(tide_ref_date(1), tide_ref_date(2), tide_ref_date(3), 0, 0, 0)
403 endif
404
405 ! Initialize reference time for tides and find relevant lunar and solar
406 ! longitudes at the reference time.
407 if (cs%use_eq_phase) call astro_longitudes_init(cs%time_ref, cs%tidal_longitudes)
408
409 ! Set the parameters for all components that are in use.
410 c=0
411 if (use_m2) then
412 c=c+1 ; cs%const_name(c) = "M2" ; cs%struct(c) = 2
413 cs%love_no(c) = 0.693 ; amp_def(c) = 0.242334 ! Default amplitude in m.
414 endif
415
416 if (use_s2) then
417 c=c+1 ; cs%const_name(c) = "S2" ; cs%struct(c) = 2
418 cs%love_no(c) = 0.693 ; amp_def(c) = 0.112743 ! Default amplitude in m.
419 endif
420
421 if (use_n2) then
422 c=c+1 ; cs%const_name(c) = "N2" ; cs%struct(c) = 2
423 cs%love_no(c) = 0.693 ; amp_def(c) = 0.046397 ! Default amplitude in m.
424 endif
425
426 if (use_k2) then
427 c=c+1 ; cs%const_name(c) = "K2" ; cs%struct(c) = 2
428 cs%love_no(c) = 0.693 ; amp_def(c) = 0.030684 ! Default amplitude in m.
429 endif
430
431 if (use_k1) then
432 c=c+1 ; cs%const_name(c) = "K1" ; cs%struct(c) = 1
433 cs%love_no(c) = 0.736 ; amp_def(c) = 0.141565 ! Default amplitude in m.
434 endif
435
436 if (use_o1) then
437 c=c+1 ; cs%const_name(c) = "O1" ; cs%struct(c) = 1
438 cs%love_no(c) = 0.695 ; amp_def(c) = 0.100661 ! Default amplitude in m.
439 endif
440
441 if (use_p1) then
442 c=c+1 ; cs%const_name(c) = "P1" ; cs%struct(c) = 1
443 cs%love_no(c) = 0.706 ; amp_def(c) = 0.046848 ! Default amplitude in m.
444 endif
445
446 if (use_q1) then
447 c=c+1 ; cs%const_name(c) = "Q1" ; cs%struct(c) = 1
448 cs%love_no(c) = 0.695 ; amp_def(c) = 0.019273 ! Default amplitude in m.
449 endif
450
451 if (use_mf) then
452 c=c+1 ; cs%const_name(c) = "MF" ; cs%struct(c) = 3
453 cs%love_no(c) = 0.693 ; amp_def(c) = 0.042041 ! Default amplitude in m.
454 endif
455
456 if (use_mm) then
457 c=c+1 ; cs%const_name(c) = "MM" ; cs%struct(c) = 3
458 cs%love_no(c) = 0.693 ; amp_def(c) = 0.022191 ! Default amplitude in m.
459 endif
460
461 ! Set defaults for all included constituents
462 ! and things that can be set by functions
463 do c=1,nc
464 freq_def(c) = tidal_frequency(cs%const_name(c))
465 love_def(c) = cs%love_no(c)
466 cs%phase0(c) = 0.0
467 if (cs%use_eq_phase) then
468 phase0_def(c) = eq_phase(cs%const_name(c), cs%tidal_longitudes)
469 else
470 phase0_def(c) = 0.0
471 endif
472 enddo
473
474 ! Parse the input file to potentially override the default values for the
475 ! frequency, amplitude and initial phase of each constituent, and log the
476 ! values that are actually used.
477 do c=1,nc
478 call get_param(param_file, mdl, "TIDE_"//trim(cs%const_name(c))//"_FREQ", cs%freq(c), &
479 "Frequency of the "//trim(cs%const_name(c))//" tidal constituent. "//&
480 "This is only used if TIDES and TIDE_"//trim(cs%const_name(c))// &
481 " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(cs%const_name(c))// &
482 " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=freq_def(c), &
483 scale=us%T_to_s)
484 call get_param(param_file, mdl, "TIDE_"//trim(cs%const_name(c))//"_AMP", cs%amp(c), &
485 "Amplitude of the "//trim(cs%const_name(c))//" tidal constituent. "//&
486 "This is only used if TIDES and TIDE_"//trim(cs%const_name(c))// &
487 " are true.", units="m", default=amp_def(c), scale=us%m_to_Z)
488 call get_param(param_file, mdl, "TIDE_"//trim(cs%const_name(c))//"_PHASE_T0", cs%phase0(c), &
489 "Phase of the "//trim(cs%const_name(c))//" tidal constituent at time 0. "//&
490 "This is only used if TIDES and TIDE_"//trim(cs%const_name(c))// &
491 " are true.", units="radians", default=phase0_def(c))
492 enddo
493
494 if (cs%use_tidal_sal_file) then
495 allocate(cs%cosphasesal(isd:ied,jsd:jed,nc))
496 allocate(cs%sinphasesal(isd:ied,jsd:jed,nc))
497 allocate(cs%ampsal(isd:ied,jsd:jed,nc))
498 do c=1,nc
499 ! Read variables with names like PHASE_SAL_M2 and AMP_SAL_M2.
500 call find_in_files(tidal_input_files, "PHASE_SAL_"//trim(cs%const_name(c)), phase, g)
501 call find_in_files(tidal_input_files, "AMP_SAL_"//trim(cs%const_name(c)), cs%ampsal(:,:,c), &
502 g, scale=us%m_to_Z)
503 call pass_var(phase, g%domain,complete=.false.)
504 call pass_var(cs%ampsal(:,:,c),g%domain,complete=.true.)
505 do j=js-1,je+1 ; do i=is-1,ie+1
506 cs%cosphasesal(i,j,c) = cos(phase(i,j)*deg_to_rad)
507 cs%sinphasesal(i,j,c) = sin(phase(i,j)*deg_to_rad)
508 enddo ; enddo
509 enddo
510 endif
511
512 if (cs%use_tidal_sal_prev) then
513 allocate(cs%cosphase_prev(isd:ied,jsd:jed,nc))
514 allocate(cs%sinphase_prev(isd:ied,jsd:jed,nc))
515 allocate(cs%amp_prev(isd:ied,jsd:jed,nc))
516 do c=1,nc
517 ! Read variables with names like PHASE_PREV_M2 and AMP_PREV_M2.
518 call find_in_files(tidal_input_files, "PHASE_PREV_"//trim(cs%const_name(c)), phase, g)
519 call find_in_files(tidal_input_files, "AMP_PREV_"//trim(cs%const_name(c)), cs%amp_prev(:,:,c), &
520 g, scale=us%m_to_Z)
521 call pass_var(phase, g%domain,complete=.false.)
522 call pass_var(cs%amp_prev(:,:,c),g%domain,complete=.true.)
523 do j=js-1,je+1 ; do i=is-1,ie+1
524 cs%cosphase_prev(i,j,c) = cos(phase(i,j)*deg_to_rad)
525 cs%sinphase_prev(i,j,c) = sin(phase(i,j)*deg_to_rad)
526 enddo ; enddo
527 enddo
528 endif
529
530 call get_param(param_file, mdl, "TIDE_ADD_NODAL", add_nodal_terms, &
531 "If true, include 18.6 year nodal modulation in the astronomical tidal forcing.", &
532 old_name="OBC_TIDE_ADD_NODAL", default=.false.)
533 call get_param(param_file, mdl, "TIDE_NODAL_REF_DATE", nodal_ref_date, &
534 "Fixed reference date to use for nodal modulation of astronomical tidal forcing.", &
535 old_name="OBC_TIDE_REF_DATE", fail_if_missing=.false., defaults=(/0, 0, 0/))
536
537 ! If the nodal correction is based on a different time, initialize that.
538 ! Otherwise, it can use N from the time reference.
539 if (add_nodal_terms) then
540 if (sum(nodal_ref_date) /= 0) then
541 ! A reference date was provided for the nodal correction
542 nodal_time = set_date(nodal_ref_date(1), nodal_ref_date(2), nodal_ref_date(3))
543 call astro_longitudes_init(nodal_time, nodal_longitudes)
544 elseif (cs%use_eq_phase) then
545 ! Astronomical longitudes were already calculated for use in equilibrium phases,
546 ! so use nodal longitude from that.
547 nodal_longitudes = cs%tidal_longitudes
548 else
549 ! Tidal reference time is a required parameter, so calculate the longitudes from that.
550 call astro_longitudes_init(cs%time_ref, nodal_longitudes)
551 endif
552 endif
553
554 allocate(cs%tide_fn(nc))
555 allocate(cs%tide_un(nc))
556
557 do c=1,nc
558 ! Find nodal corrections if needed
559 if (add_nodal_terms) then
560 call nodal_fu(trim(cs%const_name(c)), nodal_longitudes%N, cs%tide_fn(c), cs%tide_un(c))
561 else
562 cs%tide_fn(c) = 1.0
563 cs%tide_un(c) = 0.0
564 endif
565 enddo
566
567 id_clock_tides = cpu_clock_id('(Ocean tides)', grain=clock_module)
568
569end subroutine tidal_forcing_init
570
571!> This subroutine finds a named variable in a list of files and reads its
572!! values into a domain-decomposed 2-d array
573subroutine find_in_files(filenames, varname, array, G, scale)
574 character(len=*), dimension(:), intent(in) :: filenames !< The names of the files to search for the named variable
575 character(len=*), intent(in) :: varname !< The name of the variable to read
576 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
577 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: array !< The array to fill with the data [arbitrary]
578 real, optional, intent(in) :: scale !< A factor by which to rescale the array to translate it
579 !! into its desired units [arbitrary]
580 ! Local variables
581 integer :: nf
582
583 do nf=1,size(filenames)
584 if (len_trim(filenames(nf)) == 0) cycle
585 if (field_exists(filenames(nf), varname, mom_domain=g%Domain)) then
586 call mom_read_data(filenames(nf), varname, array, g%Domain, scale=scale)
587 return
588 endif
589 enddo
590
591 do nf=size(filenames),1,-1
592 if (file_exists(filenames(nf), g%Domain)) then
593 call mom_error(fatal, "MOM_tidal_forcing.F90: Unable to find "// &
594 trim(varname)//" in any of the tidal input files, last tried "// &
595 trim(filenames(nf)))
596 endif
597 enddo
598
599 call mom_error(fatal, "MOM_tidal_forcing.F90: Unable to find any of the "// &
600 "tidal input files, including "//trim(filenames(1)))
601
602end subroutine find_in_files
603
604!> This subroutine calculates the geopotential anomalies that drive the tides,
605!! including tidal self-attraction and loading from previous solutions.
606subroutine calc_tidal_forcing(Time, e_tide_eq, e_tide_sal, G, US, CS)
607 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
608 type(time_type), intent(in) :: time !< The time for the caluculation.
609 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: e_tide_eq !< The geopotential height anomalies
610 !! due to the equilibrium tides [Z ~> m].
611 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: e_tide_sal !< The geopotential height anomalies
612 !! due to the tidal SAL [Z ~> m].
613 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
614 type(tidal_forcing_cs), intent(in) :: cs !< The control structure returned by a
615 !! previous call to tidal_forcing_init.
616
617 ! Local variables
618 real :: now ! The relative time compared with the tidal reference [T ~> s]
619 real :: amp_cosomegat, amp_sinomegat ! The tidal amplitudes times the components of phase [Z ~> m]
620 real :: cosomegat, sinomegat ! The components of the phase [nondim]
621 integer :: i, j, c, m, is, ie, js, je, isq, ieq, jsq, jeq
622 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
623 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
624
625 call cpu_clock_begin(id_clock_tides)
626
627 do j=jsq,jeq+1 ; do i=isq,ieq+1
628 e_tide_eq(i,j) = 0.0
629 e_tide_sal(i,j) = 0.0
630 enddo ; enddo
631
632 if (cs%nc == 0) then
633 return
634 endif
635
636 now = time_minus_signed(time, cs%time_ref, scale=us%s_to_T)
637
638 do c=1,cs%nc
639 m = cs%struct(c)
640 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)))
641 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)))
642 do j=jsq,jeq+1 ; do i=isq,ieq+1
643 e_tide_eq(i,j) = e_tide_eq(i,j) + (amp_cosomegat*cs%cos_struct(i,j,m) + &
644 amp_sinomegat*cs%sin_struct(i,j,m))
645 enddo ; enddo
646 enddo
647
648 if (cs%use_tidal_sal_file) then ; do c=1,cs%nc
649 cosomegat = cs%tide_fn(c) * cos(cs%freq(c)*now + (cs%phase0(c) + cs%tide_un(c)))
650 sinomegat = cs%tide_fn(c) * sin(cs%freq(c)*now + (cs%phase0(c) + cs%tide_un(c)))
651 do j=jsq,jeq+1 ; do i=isq,ieq+1
652 e_tide_sal(i,j) = e_tide_sal(i,j) + cs%ampsal(i,j,c) * &
653 (cosomegat*cs%cosphasesal(i,j,c) + sinomegat*cs%sinphasesal(i,j,c))
654 enddo ; enddo
655 enddo ; endif
656
657 if (cs%use_tidal_sal_prev) then ; do c=1,cs%nc
658 cosomegat = cs%tide_fn(c) * cos(cs%freq(c)*now + (cs%phase0(c) + cs%tide_un(c)))
659 sinomegat = cs%tide_fn(c) * sin(cs%freq(c)*now + (cs%phase0(c) + cs%tide_un(c)))
660 do j=jsq,jeq+1 ; do i=isq,ieq+1
661 e_tide_sal(i,j) = e_tide_sal(i,j) - cs%sal_scalar * cs%amp_prev(i,j,c) * &
662 (cosomegat*cs%cosphase_prev(i,j,c) + sinomegat*cs%sinphase_prev(i,j,c))
663 enddo ; enddo
664 enddo ; endif
665
666 call cpu_clock_end(id_clock_tides)
667
668end subroutine calc_tidal_forcing
669
670!> This subroutine functions the same as calc_tidal_forcing but outputs a field that combines
671!! previously calculated self-attraction and loading (SAL) and tidal forcings, so that old answers
672!! can be preserved bitwise before SAL is separated out as an individual module.
673subroutine calc_tidal_forcing_legacy(Time, e_sal, e_sal_tide, e_tide_eq, e_tide_sal, G, US, CS)
674 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
675 type(time_type), intent(in) :: time !< The time for the caluculation.
676 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: e_sal !< The self-attraction and loading fields
677 !! calculated previously used to
678 !! initialized e_sal_tide [Z ~> m].
679 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: e_sal_tide !< The total geopotential height anomalies
680 !! due to both SAL and tidal forcings [Z ~> m].
681 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: e_tide_eq !< The geopotential height anomalies
682 !! due to the equilibrium tides [Z ~> m].
683 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: e_tide_sal !< The geopotential height anomalies
684 !! due to the tidal SAL [Z ~> m].
685 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
686 type(tidal_forcing_cs), intent(in) :: cs !< The control structure returned by a
687 !! previous call to tidal_forcing_init.
688
689 ! Local variables
690 real :: now ! The relative time compared with the tidal reference [T ~> s]
691 real :: amp_cosomegat, amp_sinomegat ! The tidal amplitudes times the components of phase [Z ~> m]
692 real :: cosomegat, sinomegat ! The components of the phase [nondim]
693 real :: amp_cossin ! A temporary field that adds cosines and sines [nondim]
694 integer :: i, j, c, m, is, ie, js, je, isq, ieq, jsq, jeq
695 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
696 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
697
698 call cpu_clock_begin(id_clock_tides)
699
700 do j=jsq,jeq+1 ; do i=isq,ieq+1
701 e_sal_tide(i,j) = 0.0
702 e_tide_eq(i,j) = 0.0
703 e_tide_sal(i,j) = 0.0
704 enddo ; enddo
705
706 if (cs%nc == 0) then
707 return
708 endif
709
710 now = time_minus_signed(time, cs%time_ref, scale=us%s_to_T)
711
712 do j=jsq,jeq+1 ; do i=isq,ieq+1
713 e_sal_tide(i,j) = e_sal(i,j)
714 enddo ; enddo
715
716 do c=1,cs%nc
717 m = cs%struct(c)
718 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)))
719 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)))
720 do j=jsq,jeq+1 ; do i=isq,ieq+1
721 amp_cossin = (amp_cosomegat*cs%cos_struct(i,j,m) + amp_sinomegat*cs%sin_struct(i,j,m))
722 e_sal_tide(i,j) = e_sal_tide(i,j) + amp_cossin
723 e_tide_eq(i,j) = e_tide_eq(i,j) + amp_cossin
724 enddo ; enddo
725 enddo
726
727 if (cs%use_tidal_sal_file) then ; do c=1,cs%nc
728 cosomegat = cs%tide_fn(c) * cos(cs%freq(c)*now + (cs%phase0(c) + cs%tide_un(c)))
729 sinomegat = cs%tide_fn(c) * sin(cs%freq(c)*now + (cs%phase0(c) + cs%tide_un(c)))
730 do j=jsq,jeq+1 ; do i=isq,ieq+1
731 amp_cossin = cs%ampsal(i,j,c) &
732 * (cosomegat*cs%cosphasesal(i,j,c) + sinomegat*cs%sinphasesal(i,j,c))
733 e_sal_tide(i,j) = e_sal_tide(i,j) + amp_cossin
734 e_tide_sal(i,j) = e_tide_sal(i,j) + amp_cossin
735 enddo ; enddo
736 enddo ; endif
737
738 if (cs%use_tidal_sal_prev) then ; do c=1,cs%nc
739 cosomegat = cs%tide_fn(c) * cos(cs%freq(c)*now + (cs%phase0(c) + cs%tide_un(c)))
740 sinomegat = cs%tide_fn(c) * sin(cs%freq(c)*now + (cs%phase0(c) + cs%tide_un(c)))
741 do j=jsq,jeq+1 ; do i=isq,ieq+1
742 amp_cossin = -cs%sal_scalar * cs%amp_prev(i,j,c) &
743 * (cosomegat*cs%cosphase_prev(i,j,c) + sinomegat*cs%sinphase_prev(i,j,c))
744 e_sal_tide(i,j) = e_sal_tide(i,j) + amp_cossin
745 e_tide_sal(i,j) = e_tide_sal(i,j) + amp_cossin
746 enddo ; enddo
747 enddo ; endif
748 call cpu_clock_end(id_clock_tides)
749
750end subroutine calc_tidal_forcing_legacy
751
752!> This subroutine deallocates memory associated with the tidal forcing module.
753subroutine tidal_forcing_end(CS)
754 type(tidal_forcing_cs), intent(inout) :: cs !< The control structure returned by a previous call
755 !! to tidal_forcing_init; it is deallocated here.
756
757 if (allocated(cs%sin_struct)) deallocate(cs%sin_struct)
758 if (allocated(cs%cos_struct)) deallocate(cs%cos_struct)
759
760 if (allocated(cs%cosphasesal)) deallocate(cs%cosphasesal)
761 if (allocated(cs%sinphasesal)) deallocate(cs%sinphasesal)
762 if (allocated(cs%ampsal)) deallocate(cs%ampsal)
763
764 if (allocated(cs%cosphase_prev)) deallocate(cs%cosphase_prev)
765 if (allocated(cs%sinphase_prev)) deallocate(cs%sinphase_prev)
766 if (allocated(cs%amp_prev)) deallocate(cs%amp_prev)
767end subroutine tidal_forcing_end
768
769!> \namespace tidal_forcing
770!!
771!! \section section_tides Tidal forcing
772!!
773!! Code by Robert Hallberg, August 2005, based on C-code by Harper
774!! Simmons, February, 2003, in turn based on code by Brian Arbic.
775!!
776!! The main subroutine in this file calculates the total tidal
777!! contribution to the geopotential, including self-attraction and
778!! loading terms and the astronomical contributions. All options
779!! are selected with entries in a file that is parsed at run-time.
780!! Overall tides are enabled with the run-time parameter 'TIDES=True'.
781!! Tidal constituents must be individually enabled with lines like
782!! 'TIDE_M2=True'. This file has default values of amplitude,
783!! frequency, Love number, and phase at time 0 for the Earth's M2,
784!! S2, N2, K2, K1, O1, P1, Q1, MF, and MM tidal constituents, but
785!! the frequency, amplitude and phase ant time 0 for each constituent
786!! can be changed at run time by setting variables like TIDE_M2_FREQ,
787!! TIDE_M2_AMP and TIDE_M2_PHASE_T0 (for M2).
788!!
789!! In addition, approaches to calculate self-attraction and loading
790!! due to tides (harmonics of astronomical forcing frequencies)
791!! are provided. <code>TIDAL_SAL_FROM_FILE</code> can be set to read the phase and
792!! amplitude of the tidal SAL. <code>USE_PREVIOUS_TIDES</code> may be useful in
793!! combination with the scalar approximation to iterate the SAL to
794!! convergence (for details, see \cite Arbic2004). With
795!! <code>TIDAL_SAL_FROM_FILE</code> or <code>USE_PREVIOUS_TIDES</code>, a list of input
796!! files must be provided to describe each constituent's properties from
797!! a previous solution. The online SAL calculations that are functions
798!! of SSH (rather should be bottom pressure anmoaly), either a scalar
799!! approximation or with spherical harmonic transforms, are located in
800!! <code>MOM_self_attr_load</code>.
801end module mom_tidal_forcing