MOM_harmonic_analysis.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!> Inline harmonic analysis (conventional)
7
8use mom_time_manager, only : time_type, real_to_time, time_to_real, time_minus_signed
9use mom_time_manager, only : set_date, get_date, increment_date
10use mom_time_manager, only : operator(+), operator(-), operator(<), operator(>), operator(>=)
11use mom_grid, only : ocean_grid_type
13use mom_file_parser, only : param_file_type, get_param
14use mom_io, only : file_exists, open_ascii_file, readonly_file, close_file
15use mom_io, only : mom_infra_file, vardesc, mom_field
16use mom_io, only : var_desc, create_mom_file, single_file, mom_write_field
17use mom_error_handler, only : mom_mesg, mom_error, note
19
20implicit none ; private
21
22public ha_init, ha_accum
23
24#include <MOM_memory.h>
25
26!> The private control structure for storing the HA info of a particular field
27type, private :: ha_type
28 character(len=16) :: key = "none" !< Name of the field of which harmonic analysis is to be performed
29 character(len=1) :: grid !< The grid on which the field is defined ('h', 'q', 'u', or 'v')
30 real :: old_time = -1.0 !< The time of the previous accumulating step [T ~> s]
31 real, allocatable :: ref(:,:) !< The initial field in arbitrary units [A]
32 real, allocatable :: ftf(:,:) !< Accumulator of (F' * F) [nondim]
33 real, allocatable :: ftssh(:,:,:) !< Accumulator of (F' * SSH_in) in arbitrary units [A]
34 !>@{ Lower and upper bounds of input data
35 integer :: is, ie, js, je
36 !>@}
37end type ha_type
38
39!> A linked list of control structures that store the HA info of different fields
40type, private :: ha_node
41 type(ha_type) :: this !< Control structure of the current field in the list
42 type(ha_node), pointer :: next !< The list of other fields
43end type ha_node
44
45!> The public control structure of the MOM_harmonic_analysis module
46type, public :: harmonic_analysis_cs ; private
47 logical :: haready = .false. !< If true, perform harmonic analysis
48 type(time_type) :: &
49 time_start, & !< Start time of harmonic analysis
50 time_end, & !< End time of harmonic analysis
51 time_ref !< Reference time (t = 0) used to calculate tidal forcing
52 real, allocatable, dimension(:) :: &
53 freq, & !< The frequency of a tidal constituent [T-1 ~> s-1]
54 phase0, & !< The phase of a tidal constituent at time 0 [rad]
55 tide_fn, & !< Amplitude modulation of tides by nodal cycle [nondim].
56 tide_un !< Phase modulation of tides by nodal cycle [rad].
57 integer :: nc !< The number of tidal constituents in use
58 integer :: length !< Number of fields of which harmonic analysis is to be performed
59 character(len=4), allocatable, dimension(:) :: const_name !< The name of each constituent
60 character(len=255) :: path !< Path to directory where output will be written
61 type(unit_scale_type) :: us !< A dimensional unit scaling type
62 type(ha_node), pointer :: list => null() !< A linked list for storing the HA info of different fields
64
65contains
66
67!> This subroutine sets static variables used by this module and initializes CS%list.
68!! THIS MUST BE CALLED AT THE END OF tidal_forcing_init.
69subroutine ha_init(Time, US, param_file, nc, CS)
70 type(time_type), intent(in) :: time !< The current model time
71 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
72 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
73 integer, intent(in) :: nc !< The number of tidal constituents in use
74 type(harmonic_analysis_cs), intent(out) :: cs !< Control structure of the MOM_harmonic_analysis module
75
76 ! Local variables
77 logical :: tides !< True if tidal forcing module is enabled
78 logical :: use_eq_phase !< If true, tidal forcing is phase-shifted to match
79 !! equilibrium tide. Set to false if providing tidal phases
80 !! that have already been shifted by the
81 !! astronomical/equilibrium argument
82 logical :: add_nodal_terms !< If true, insert terms for the 18.6 year modulation when
83 !! calculating tidal forcing.
84 integer, dimension(3) :: tide_ref_date !< Reference date (t = 0) for tidal forcing (year, month, day)
85 integer, dimension(3) :: nodal_ref_date !< Date to calculate nodal modulation for (year, month, day)
86 type(time_type) :: nodal_time !< Model time to calculate nodal modulation for.
87 type(astro_longitudes) :: tidal_longitudes !< Astronomical longitudes used to calculate
88 !! tidal phases at t = 0.
89 type(astro_longitudes) :: nodal_longitudes !< Solar and lunar longitudes for tidal forcing
90 character(len=50) :: const_name !< Names of all tidal constituents to be harmonically analyzed
91 integer :: c
92
93 type(ha_type) :: ha1 !< A temporary, null field used for initializing CS%list
94 real :: ha_start_time !< Start time of harmonic analysis [T ~> s]
95 real :: ha_end_time !< End time of harmonic analysis [T ~> s]
96 logical :: ha_ssh, ha_ubt, ha_vbt
97 character(len=40) :: mdl="MOM_harmonic_analysis" !< This module's name
98 character(len=255) :: mesg
99 integer :: year, month, day, hour, minute, second
100
101 call get_param(param_file, mdl, "TIDES", tides, &
102 "If true, apply tidal momentum forcing.", default=.false., do_not_log=.true.)
103 call get_param(param_file, mdl, "TIDE_USE_EQ_PHASE", use_eq_phase, &
104 "If true, add the equilibrium phase argument to the specified tidal phases.", &
105 old_name="OBC_TIDE_ADD_EQ_PHASE", default=.false., do_not_log=tides)
106 call get_param(param_file, mdl, "TIDE_ADD_NODAL", add_nodal_terms, &
107 "If true, include 18.6 year nodal modulation in the boundary tidal forcing.", &
108 old_name="OBC_TIDE_ADD_NODAL", default=.false., do_not_log=tides)
109 call get_param(param_file, mdl, "TIDE_REF_DATE", tide_ref_date, &
110 "Reference date to use for tidal calculations and equilibrium phase.", &
111 old_name="OBC_TIDE_REF_DATE", defaults=(/0, 0, 0/), do_not_log=tides)
112 call get_param(param_file, mdl, "TIDE_NODAL_REF_DATE", nodal_ref_date, &
113 "Fixed reference date to use for nodal modulation.", &
114 old_name="OBC_TIDE_NODAL_REF_DATE", defaults=(/0, 0, 0/), do_not_log=tides)
115 call get_param(param_file, mdl, "HA_CONSTITUENTS", const_name, &
116 "Names of tidal constituents to be harmonically analyzed. "//&
117 "They don't have to be the same as those used in MOM_tidal_forcing.", &
118 fail_if_missing=.true.)
119
120 if (sum(tide_ref_date) == 0) then ! tide_ref_date defaults to 0.
121 cs%time_ref = set_date(1, 1, 1, 0, 0, 0)
122 else
123 if (.not. use_eq_phase) then
124 ! Using a reference date but not using phase relative to equilibrium.
125 ! This makes sense as long as either phases are overridden, or
126 ! correctly simulating tidal phases is not desired.
127 call mom_mesg('Tidal phases will *not* be corrected with equilibrium arguments.')
128 endif
129 cs%time_ref = set_date(tide_ref_date(1), tide_ref_date(2), tide_ref_date(3), 0, 0, 0)
130 endif
131
132 ! Initialize reference time for tides and find relevant lunar and solar
133 ! longitudes at the reference time.
134 if (use_eq_phase) call astro_longitudes_init(cs%time_ref, tidal_longitudes)
135
136 ! If the nodal correction is based on a different time, initialize that.
137 ! Otherwise, it can use N from the time reference.
138 if (add_nodal_terms) then
139 if (sum(nodal_ref_date) /= 0) then
140 ! A reference date was provided for the nodal correction
141 nodal_time = set_date(nodal_ref_date(1), nodal_ref_date(2), nodal_ref_date(3))
142 call astro_longitudes_init(nodal_time, nodal_longitudes)
143 elseif (use_eq_phase) then
144 ! Astronomical longitudes were already calculated for use in equilibrium phases,
145 ! so use nodal longitude from that.
146 nodal_longitudes = tidal_longitudes
147 else
148 ! Tidal reference time is a required parameter, so calculate the longitudes from that.
149 call astro_longitudes_init(cs%time_ref, nodal_longitudes)
150 endif
151 endif
152
153 allocate(cs%const_name(nc))
154 allocate(cs%freq(nc))
155 allocate(cs%phase0(nc))
156 allocate(cs%tide_fn(nc))
157 allocate(cs%tide_un(nc))
158
159 ! Tidal constituents for harmonic analysis can be different from those defined in MOM_tidal_forcing
160 read(const_name, *) cs%const_name
161
162 ! For major tidal constituents, tidal parameters defined in MOM_tidal_forcing will be used.
163 ! For those not available in MOM_tidal_forcing, parameters needs to be defined in MOM_input.
164 do c=1,nc
165 call get_param(param_file, mdl, "HA_"//trim(cs%const_name(c))//"_FREQ", &
166 cs%freq(c), "Frequency of the "//trim(cs%const_name(c))//&
167 " constituent. This is used if USE_HA is true and "//trim(cs%const_name(c))//&
168 " is in HA_CONSTITUENTS.", units="rad s-1", scale=us%T_to_s, default=0.0)
169 if (cs%freq(c)<=0.0) then
170 select case (trim(cs%const_name(c)))
171 case ('M4')
172 cs%freq(c) = tidal_frequency('M2') * 2
173 case ('M6')
174 cs%freq(c) = tidal_frequency('M2') * 3
175 case ('M8')
176 cs%freq(c) = tidal_frequency('M2') * 4
177 case ('S4')
178 cs%freq(c) = tidal_frequency('S2') * 2
179 case ('S6')
180 cs%freq(c) = tidal_frequency('S2') * 3
181 case ('MK3')
182 cs%freq(c) = tidal_frequency('M2') + tidal_frequency('K1')
183 case ('MS4')
184 cs%freq(c) = tidal_frequency('M2') + tidal_frequency('S2')
185 case ('MN4')
186 cs%freq(c) = tidal_frequency('M2') + tidal_frequency('N2')
187 case default
188 cs%freq(c) = tidal_frequency(trim(cs%const_name(c)))
189 end select
190 endif
191
192 call get_param(param_file, mdl, "HA_"//trim(cs%const_name(c))//"_PHASE_T0", cs%phase0(c), &
193 "Phase of the "//trim(cs%const_name(c))//" tidal constituent at time 0. "//&
194 "This is only used if USE_HA is true and "//trim(cs%const_name(c))// &
195 " is in HA_CONSTITUENTS.", units="radians", default=0.0)
196 if (use_eq_phase) cs%phase0(c) = eq_phase(trim(cs%const_name(c)), tidal_longitudes)
197
198 ! Nodal modulation should be turned off for tidal constituents not available in MOM_tidal_forcing
199 if (add_nodal_terms) then
200 call nodal_fu(trim(trim(cs%const_name(c))), nodal_longitudes%N, cs%tide_fn(c), cs%tide_un(c))
201 else
202 cs%tide_fn(c) = 1.0
203 cs%tide_un(c) = 0.0
204 endif
205 enddo
206
207 ! Determine CS%time_start and CS%time_end
208 call get_param(param_file, mdl, "HA_START_TIME", ha_start_time, &
209 "Start time of harmonic analysis, in units of days after "//&
210 "the start of the current run segment. Must be smaller than "//&
211 "HA_END_TIME, otherwise harmonic analysis will not be performed. "//&
212 "If negative, |HA_START_TIME| determines the length of harmonic analysis, "//&
213 "and harmonic analysis will start |HA_START_TIME| days before HA_END_TIME, "//&
214 "or at the beginning of the run segment, whichever occurs later.", &
215 units="days", default=0.0, scale=86400.0*us%s_to_T)
216 call get_param(param_file, mdl, "HA_END_TIME", ha_end_time, &
217 "End time of harmonic analysis, in units of days after "//&
218 "the start of the current run segment. Must be positive "//&
219 "and smaller than the length of the currnet run segment, "//&
220 "otherwise harmonic analysis will not be performed.", &
221 units="days", default=0.0, scale=86400.0*us%s_to_T)
222
223 if (ha_end_time <= 0.0) then
224 call mom_mesg('MOM_harmonic_analysis: HA_END_TIME is zero or negative. '//&
225 'Harmonic analysis will not be performed.')
226 cs%HAready = .false. ; return
227 endif
228
229 if (ha_end_time <= ha_start_time) then
230 call mom_mesg('MOM_harmonic_analysis: HA_END_TIME is smaller than or equal to HA_START_TIME. '//&
231 'Harmonic analysis will not be performed.')
232 cs%HAready = .false. ; return
233 endif
234
235 cs%HAready = .true.
236
237 if (ha_start_time < 0.0) then
238 ha_start_time = ha_end_time + ha_start_time
239 if (ha_start_time <= 0.0) ha_start_time = 0.0
240 endif
241
242 cs%time_start = time + real_to_time(ha_start_time, unscale=us%T_to_s)
243 cs%time_end = time + real_to_time(ha_end_time, unscale=us%T_to_s)
244
245 call get_date(time, year, month, day, hour, minute, second)
246 write(mesg,*) "MOM_harmonic_analysis: run segment starts on ", year, month, day, hour, minute, second
247 call mom_error(note, trim(mesg))
248 call get_date(cs%time_start, year, month, day, hour, minute, second)
249 write(mesg,*) "MOM_harmonic_analysis: harmonic analysis starts on ", year, month, day, hour, minute, second
250 call mom_error(note, trim(mesg))
251 call get_date(cs%time_end, year, month, day, hour, minute, second)
252 write(mesg,*) "MOM_harmonic_analysis: harmonic analysis ends on ", year, month, day, hour, minute, second
253 call mom_error(note, trim(mesg))
254
255 ! Set path to directory where output will be written
256 call get_param(param_file, mdl, "HA_PATH", cs%path, &
257 "Path to output files for runtime harmonic analysis.", default="./")
258
259 ! Populate some parameters of the control structure
260 cs%nc = nc
261 cs%length = 0
262 cs%US = us
263
264 ! Initialize CS%list
265 allocate(cs%list)
266 cs%list%this = ha1
267 nullify(cs%list%next)
268
269 ! Register variables/fields to be analyzed
270 call get_param(param_file, mdl, "HA_SSH", ha_ssh, &
271 "If true, perform harmonic analysis of sea serface height.", default=.false.)
272 if (ha_ssh) call ha_register('ssh', 'h', cs)
273 call get_param(param_file, mdl, "HA_UBT", ha_ubt, &
274 "If true, perform harmonic analysis of zonal barotropic velocity.", default=.false.)
275 if (ha_ubt) call ha_register('ubt', 'u', cs)
276 call get_param(param_file, mdl, "HA_VBT", ha_vbt, &
277 "If true, perform harmonic analysis of meridional barotropic velocity.", default=.false.)
278 if (ha_vbt) call ha_register('vbt', 'v', cs)
279
280end subroutine ha_init
281
282!> This subroutine registers each of the fields on which HA is to be performed.
283subroutine ha_register(key, grid, CS)
284 character(len=*), intent(in) :: key !< Name of the current field
285 character(len=1), intent(in) :: grid !< The grid on which the key field is defined
286 type(harmonic_analysis_cs), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module
287
288 ! Local variables
289 type(ha_type) :: ha1 !< Control structure for the current field
290 type(ha_node), pointer :: tmp !< A temporary list to hold the current field
291
292 if (.not. cs%HAready) return
293
294 allocate(tmp)
295 ha1%key = trim(key)
296 ha1%grid = trim(grid)
297 tmp%this = ha1
298 tmp%next => cs%list
299 cs%list => tmp
300 cs%length = cs%length + 1
301
302end subroutine ha_register
303
304!> This subroutine accumulates the temporal basis functions in FtF and FtSSH and then calls HA_write to compute
305!! harmonic constants and write results. The tidal constituents are those used in MOM_tidal_forcing, plus the
306!! mean (of zero frequency). For FtF, only the main diagonal and entries below it are calculated, which are needed
307!! for Cholesky decomposition.
308subroutine ha_accum(key, data, Time, G, CS)
309 character(len=*), intent(in) :: key !< Name of the current field
310 real, dimension(:,:), intent(in) :: data !< Input data of which harmonic analysis is to be performed [A]
311 type(time_type), intent(in) :: time !< The current model time
312 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
313 type(harmonic_analysis_cs), intent(inout) :: cs !< Control structure of the MOM_harmonic_analysis module
314
315 ! Local variables
316 type(ha_type), pointer :: ha1
317 type(ha_node), pointer :: tmp
318 real :: now !< The relative time compared with the tidal reference [T ~> s]
319 real :: dt !< The current time step size of the accumulator [T ~> s]
320 real :: cosomegat, sinomegat, ccosomegat, ssinomegat !< The components of the phase [nondim]
321 integer :: nc, i, j, k, c, cc, icos, isin, iccos, issin, is, ie, js, je
322 character(len=128) :: mesg
323
324 ! Exit the accumulator in the following cases
325 if (.not. cs%HAready) return
326 if (cs%length == 0) return
327 if (time < cs%time_start) return
328 if (time > cs%time_end) return
329
330 ! Loop through the full list to find the current field
331 tmp => cs%list
332 do k=1,cs%length
333 ha1 => tmp%this
334 if (trim(key) == trim(ha1%key)) exit
335 tmp => tmp%next
336 if (k == cs%length) return !< Do not perform harmonic analysis of a field that is not registered
337 enddo
338
339 nc = cs%nc
340 now = time_minus_signed(time, cs%time_ref, scale=cs%US%s_to_T)
341
342 !!! Additional processing at the initial accumulating step !!!
343 if (ha1%old_time < 0.0) then
344 ha1%old_time = now
345
346 write(mesg,*) "MOM_harmonic_analysis: initializing accumulator, key = ", trim(ha1%key)
347 call mom_error(note, trim(mesg))
348
349 ! Get the lower and upper bounds of input data
350 ha1%is = lbound(data,1) ; is = ha1%is
351 ha1%ie = ubound(data,1) ; ie = ha1%ie
352 ha1%js = lbound(data,2) ; js = ha1%js
353 ha1%je = ubound(data,2) ; je = ha1%je
354
355 allocate(ha1%ref(is:ie,js:je), source=0.0)
356 allocate(ha1%FtF(2*nc+1,2*nc+1), source=0.0)
357 allocate(ha1%FtSSH(is:ie,js:je,2*nc+1), source=0.0)
358 ha1%ref(:,:) = data(:,:)
359 endif
360
361 dt = now - ha1%old_time
362 ha1%old_time = now !< Keep track of time so we know when Time approaches CS%time_end
363
364 is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je
365
366 !!! Accumulator of FtF !!!
367 !< First entry, corresponding to the zero frequency constituent (mean)
368 ha1%FtF(1,1) = ha1%FtF(1,1) + 1.0
369
370 do c=1,nc
371 icos = 2*c
372 isin = 2*c+1
373 cosomegat = cs%tide_fn(c) * cos(cs%freq(c) * now + (cs%phase0(c) + cs%tide_un(c)))
374 sinomegat = cs%tide_fn(c) * sin(cs%freq(c) * now + (cs%phase0(c) + cs%tide_un(c)))
375
376 ! First column, corresponding to the zero frequency constituent (mean)
377 ha1%FtF(icos,1) = ha1%FtF(icos,1) + cosomegat
378 ha1%FtF(isin,1) = ha1%FtF(isin,1) + sinomegat
379
380 do cc=1,c
381 iccos = 2*cc
382 issin = 2*cc+1
383 ccosomegat = cs%tide_fn(cc) * cos(cs%freq(cc) * now + (cs%phase0(cc) + cs%tide_un(cc)))
384 ssinomegat = cs%tide_fn(cc) * sin(cs%freq(cc) * now + (cs%phase0(cc) + cs%tide_un(cc)))
385
386 ! Interior of the matrix, corresponding to the products of cosine and sine terms
387 ha1%FtF(icos,iccos) = ha1%FtF(icos,iccos) + cosomegat * ccosomegat
388 ha1%FtF(icos,issin) = ha1%FtF(icos,issin) + cosomegat * ssinomegat
389 ha1%FtF(isin,iccos) = ha1%FtF(isin,iccos) + sinomegat * ccosomegat
390 ha1%FtF(isin,issin) = ha1%FtF(isin,issin) + sinomegat * ssinomegat
391 enddo ! cc=1,c
392 enddo ! c=1,nc
393
394 !!! Accumulator of FtSSH !!!
395 !< First entry, corresponding to the zero frequency constituent (mean)
396 do j=js,je ; do i=is,ie
397 ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j))
398 enddo ; enddo
399
400 !< The remaining entries
401 do c=1,nc
402 icos = 2*c
403 isin = 2*c+1
404 cosomegat = cs%tide_fn(c) * cos(cs%freq(c) * now + (cs%phase0(c) + cs%tide_un(c)))
405 sinomegat = cs%tide_fn(c) * sin(cs%freq(c) * now + (cs%phase0(c) + cs%tide_un(c)))
406 do j=js,je ; do i=is,ie
407 ha1%FtSSH(i,j,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat
408 ha1%FtSSH(i,j,isin) = ha1%FtSSH(i,j,isin) + (data(i,j) - ha1%ref(i,j)) * sinomegat
409 enddo ; enddo
410 enddo ! c=1,nc
411
412 !!! Compute harmonic constants and write output as Time approaches CS%time_end !!!
413 ! This guarantees that HA_write will be called before Time becomes larger than CS%time_end.
414 ! Result of subtracting time types is always >= 0, which is acceptable here.
415 if (time_to_real(cs%time_end - time, scale=cs%US%s_to_T) <= dt) then
416 call ha_write(ha1, time, g, cs)
417
418 write(mesg,*) "MOM_harmonic_analysis: harmonic analysis done, key = ", trim(ha1%key)
419 call mom_error(note, trim(mesg))
420
421 ! De-register the current field and deallocate memory
422 ha1%key = 'none'
423 deallocate(ha1%ref)
424 deallocate(ha1%FtSSH)
425 endif
426
427end subroutine ha_accum
428
429!> This subroutine computes the harmonic constants and write output for the current field
430subroutine ha_write(ha1, Time, G, CS)
431 type(ha_type), pointer, intent(in) :: ha1 !< Control structure for the current field
432 type(time_type), intent(in) :: Time !< The current model time
433 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
434 type(harmonic_analysis_cs), intent(in) :: CS !< Control structure of the MOM_harmonic_analysis module
435
436 ! Local variables
437 real, dimension(:,:,:), allocatable :: FtSSHw !< An array containing the harmonic constants [A]
438 integer :: year, month, day, hour, minute, second
439 integer :: nc, i, j, k, is, ie, js, je
440
441 character(len=255) :: filename !< Output file name
442 type(mom_infra_file) :: cdf !< The file handle for output harmonic constants
443 type(vardesc), allocatable :: cdf_vars(:) !< Output variable names
444 type(mom_field), allocatable :: cdf_fields(:) !< Field type variables for the output fields
445
446 nc = cs%nc ; is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je
447
448 allocate(ftsshw(is:ie,js:je,2*nc+1), source=0.0)
449
450 ! Compute the harmonic coefficients
451 call ha_solver(ha1, nc, ha1%FtF, ftsshw)
452
453 ! Output file name
454 call get_date(time, year, month, day, hour, minute, second)
455 write(filename, '(a,"HA_",a,i0.4,i0.2,i0.2,".nc")') &
456 trim(cs%path), trim(ha1%key), year, month, day
457
458 allocate(cdf_vars(2*nc+1))
459 allocate(cdf_fields(2*nc+1))
460
461 ! Variable names
462 cdf_vars(1) = var_desc("z0", "m" ,"mean value", ha1%grid, '1', '1')
463 do k=1,nc
464 cdf_vars(2*k ) = var_desc(trim(cs%const_name(k))//"cos", "m", "cosine coefficient", ha1%grid, '1', '1')
465 cdf_vars(2*k+1) = var_desc(trim(cs%const_name(k))//"sin", "m", "sine coefficient", ha1%grid, '1', '1')
466 enddo
467
468 ! Create output file
469 call create_mom_file(cdf, trim(filename), cdf_vars, &
470 2*nc+1, cdf_fields, single_file, 86400.0, g=g)
471
472 ! Add the initial field back to the mean state
473 do j=js,je ; do i=is,ie
474 ftsshw(i,j,1) = ftsshw(i,j,1) + ha1%ref(i,j)
475 enddo ; enddo
476
477 ! Write data
478 call mom_write_field(cdf, cdf_fields(1), g%domain, ftsshw(:,:,1), 0.0)
479 do k=1,nc
480 call mom_write_field(cdf, cdf_fields(2*k ), g%domain, ftsshw(:,:,2*k ), 0.0)
481 call mom_write_field(cdf, cdf_fields(2*k+1), g%domain, ftsshw(:,:,2*k+1), 0.0)
482 enddo
483
484 call cdf%flush()
485 deallocate(cdf_vars)
486 deallocate(cdf_fields)
487 deallocate(ftsshw)
488
489end subroutine ha_write
490
491!> This subroutine computes the harmonic constants (stored in x) using the dot products of the temporal
492!! basis functions accumulated in FtF, and the dot products of the SSH (or other fields) with the temporal basis
493!! functions accumulated in FtSSH. The system is solved by Cholesky decomposition,
494!!
495!! FtF * x = FtSSH, => L * (L' * x) = FtSSH, => L * y = FtSSH,
496!!
497!! where L is the lower triangular matrix, y = L' * x, and x is the solution vector.
498!!
499subroutine ha_solver(ha1, nc, FtF, x)
500 type(ha_type), pointer, intent(in) :: ha1 !< Control structure for the current field
501 integer, intent(in) :: nc !< Number of harmonic constituents
502 real, dimension(:,:), intent(in) :: FtF !< Accumulator of (F' * F) for all fields [nondim]
503 real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je,2*nc+1), &
504 intent(out) :: x !< Solution vector of harmonic constants [A]
505
506 ! Local variables
507 real :: tmp0 !< Temporary variable for Cholesky decomposition [nondim]
508 real, dimension(2*nc+1,2*nc+1) :: L !< Lower triangular matrix of Cholesky decomposition [nondim]
509 real, dimension(2*nc+1) :: tmp1 !< Inverse of the diagonal entries of L [nondim]
510 real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je) :: tmp2 !< 2D temporary array involving FtSSH [A]
511 real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je,2*nc+1) :: y !< 3D temporary array, i.e., L' * x [A]
512 integer :: k, m, n
513
514 ! Cholesky decomposition
515 do m=1,2*nc+1
516
517 ! First, calculate the diagonal entries
518 tmp0 = 0.0
519 do k=1,m-1 ! This loop operates along the m-th row
520 tmp0 = tmp0 + l(m,k) * l(m,k)
521 enddo
522 l(m,m) = sqrt(ftf(m,m) - tmp0) ! This is the m-th diagonal entry
523
524 ! Now calculate the off-diagonal entries
525 tmp1(m) = 1 / l(m,m)
526 do k=m+1,2*nc+1 ! This loop operates along the column below the m-th diagonal entry
527 tmp0 = 0.0
528 do n=1,m-1
529 tmp0 = tmp0 + l(k,n) * l(m,n)
530 enddo
531 l(k,m) = (ftf(k,m) - tmp0) * tmp1(m) ! This is the k-th off-diagonal entry below the m-th diagonal entry
532 enddo
533 enddo
534
535 ! Solve for y from L * y = FtSSH
536 do k=1,2*nc+1
537 tmp2(:,:) = 0.0
538 do m=1,k-1
539 tmp2(:,:) = tmp2(:,:) + l(k,m) * y(:,:,m)
540 enddo
541 y(:,:,k) = (ha1%FtSSH(:,:,k) - tmp2(:,:)) * tmp1(k)
542 enddo
543
544 ! Solve for x from L' * x = y
545 do k=2*nc+1,1,-1
546 tmp2(:,:) = 0.0
547 do m=k+1,2*nc+1
548 tmp2(:,:) = tmp2(:,:) + l(m,k) * x(:,:,m)
549 enddo
550 x(:,:,k) = (y(:,:,k) - tmp2(:,:)) * tmp1(k)
551 enddo
552
553end subroutine ha_solver
554
555!> \namespace harmonic_analysis
556!!
557!! Major revision (August, 2025)
558!!
559!! This module is now independent of MOM_tidal_forcing, providing more flexibility for performing harmonic analyses
560!! on tidal constituents not available in MOM_tidal_forcing (e.g., MK3, M4), with the following conditions:
561!! 1) For tidal constituents not available in MOM_tidal_forcing, the frequencies and equilibrium phases (if used)
562!! must be specified manually in MOM_input.
563!! 2) If any tidal constituents not available in MOM_tidal_forcing are used, the nodal modulation cannot be added.
564!! Or, if nodal modulation is added, then harmonic analysis can only be performed on major tidal constituents.
565!!
566!! Original version (April, 2024)
567!!
568!! This module computes the harmonic constants which can be used to reconstruct the tidal elevation (or other
569!! fields) through SSH = F * x, where F is an nt-by-2*nc matrix (nt is the number of time steps and nc is the
570!! number of tidal constituents) containing the cosine/sine functions for each frequency evaluated at each time
571!! step, and x is a 2*nc-by-1 vector containing the constant coefficients of the sine/cosine for each constituent
572!! (i.e., the harmonic constants). At each grid point, the harmonic constants are computed using least squares,
573!!
574!! (F' * F) * x = F' * SSH_in, => FtF * x = FtSSH,
575!!
576!! where the prime denotes matrix transpose, and SSH_in is the sea surface height (or other fields) determined by
577!! the model. The dot products (F' * F) and (F' * SSH_in) are computed by accumulating the sums as the model is
578!! running and stored in the arrays FtF and FtSSH, respectively. The FtF matrix is inverted as needed before
579!! computing and writing out the harmonic constants.
580!!
581!! Ed Zaron and William Xu (chengzhu.xu@oregonstate.edu)
582
583end module mom_harmonic_analysis
584