From 85c72b8c0a0cca2341c836afcd379c03cb002136 Mon Sep 17 00:00:00 2001 From: William Xu Date: Wed, 12 Jun 2024 05:27:19 -0600 Subject: [PATCH] Inline harmonic analysis The module MOM_harmonic_analysis computes the constant coefficients of sine/cosine functions for the SSH and barotropic velocity fields of each tidal constituent. --- src/core/MOM.F90 | 9 +- src/core/MOM_barotropic.F90 | 16 +- src/core/MOM_dynamics_split_RK2.F90 | 17 +- src/core/MOM_dynamics_split_RK2b.F90 | 16 +- src/diagnostics/MOM_harmonic_analysis.F90 | 430 ++++++++++++++++++ .../lateral/MOM_tidal_forcing.F90 | 19 +- 6 files changed, 496 insertions(+), 11 deletions(-) create mode 100644 src/diagnostics/MOM_harmonic_analysis.F90 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 8b34b3e6a3..9e358ab17e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -95,6 +95,7 @@ module MOM use MOM_forcing_type, only : homogenize_forcing, homogenize_mech_forcing use MOM_grid, only : ocean_grid_type, MOM_grid_init, MOM_grid_end use MOM_grid, only : set_first_direction +use MOM_harmonic_analysis, only : HA_accum_FtF, HA_accum_FtSSH, harmonic_analysis_CS use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_hor_index, only : rotate_hor_index use MOM_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz @@ -389,6 +390,8 @@ module MOM !< Pointer to the control structure used for the mode-split RK2 dynamics type(MOM_dyn_split_RK2b_CS), pointer :: dyn_split_RK2b_CSp => NULL() !< Pointer to the control structure used for an alternate version of the mode-split RK2 dynamics + type(harmonic_analysis_CS), pointer :: HA_CSp => NULL() + !< Pointer to the control structure for harmonic analysis type(thickness_diffuse_CS) :: thickness_diffuse_CSp !< Pointer to the control structure used for the isopycnal height diffusive transport. !! This is also common referred to as Gent-McWilliams diffusion @@ -911,6 +914,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS enddo ; enddo endif + if (associated(CS%HA_CSp)) call HA_accum_FtF(Time_Local, CS%HA_CSp) call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, & dt_therm_here, bbl_time_int, CS, & @@ -1020,6 +1024,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ssh(i,j) = CS%ssh_rint(i,j)*I_wt_ssh CS%ave_ssh_ibc(i,j) = ssh(i,j) enddo ; enddo + if (associated(CS%HA_CSp)) call HA_accum_FtSSH('ssh', ssh, Time_local, G, CS%HA_CSp) if (do_dyn) then call adjust_ssh_for_p_atm(CS%tv, G, GV, US, CS%ave_ssh_ibc, forces%p_surf_SSH, & CS%calc_rho_for_sea_lev) @@ -3247,13 +3252,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & allocate(eta(SZI_(G),SZJ_(G)), source=0.0) if (CS%use_alt_split) then call initialize_dyn_split_RK2b(CS%u, CS%v, CS%h, CS%tv, CS%uh, CS%vh, eta, Time, & - G, GV, US, param_file, diag, CS%dyn_split_RK2b_CSp, restart_CSp, & + G, GV, US, param_file, diag, CS%dyn_split_RK2b_CSp, CS%HA_CSp, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%thickness_diffuse_CSp, CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & CS%visc, dirs, CS%ntrunc, CS%pbv, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil) else call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%tv, CS%uh, CS%vh, eta, Time, & - G, GV, US, param_file, diag, CS%dyn_split_RK2_CSp, restart_CSp, & + G, GV, US, param_file, diag, CS%dyn_split_RK2_CSp, CS%HA_CSp, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%thickness_diffuse_CSp, CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & CS%visc, dirs, CS%ntrunc, CS%pbv, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index a132200759..bc8fddbdde 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -16,6 +16,7 @@ module MOM_barotropic use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : mech_forcing use MOM_grid, only : ocean_grid_type +use MOM_harmonic_analysis, only : HA_accum_FtSSH, harmonic_analysis_CS use MOM_hor_index, only : hor_index_type use MOM_io, only : vardesc, var_desc, MOM_read_data, slasher use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, open_boundary_query @@ -289,6 +290,7 @@ module MOM_barotropic type(MOM_domain_type), pointer :: BT_Domain => NULL() !< Barotropic MOM domain type(hor_index_type), pointer :: debug_BT_HI => NULL() !< debugging copy of horizontal index_type type(SAL_CS), pointer :: SAL_CSp => NULL() !< Control structure for SAL + type(harmonic_analysis_CS), pointer :: HA_CSp => NULL() !< Control structure for harmonic analysis logical :: module_is_initialized = .false. !< If true, module has been initialized integer :: isdw !< The lower i-memory limit for the wide halo arrays. @@ -2551,6 +2553,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, eta_out(i,j) = eta_wtd(i,j) * I_sum_wt_eta enddo ; enddo + ! Accumulator is updated at the end of every baroclinic time step. + ! Harmonic analysis will not be performed of a field that is not registered. + if (associated(CS%HA_CSp) .and. find_etaav) then + call HA_accum_FtSSH('ubt', ubt, CS%Time, G, CS%HA_CSp) + call HA_accum_FtSSH('vbt', vbt, CS%Time, G, CS%HA_CSp) + endif + if (id_clock_calc_post > 0) call cpu_clock_end(id_clock_calc_post) if (id_clock_pass_post > 0) call cpu_clock_begin(id_clock_pass_post) if (G%nonblocking_updates) then @@ -4402,7 +4411,7 @@ end subroutine bt_mass_source !! barotropic calculation and initializes any barotropic fields that have not !! already been initialized. subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, & - restart_CS, calc_dtbt, BT_cont, SAL_CSp) + restart_CS, calc_dtbt, BT_cont, SAL_CSp, HA_CSp) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -4428,6 +4437,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, !! barotropic flow. type(SAL_CS), target, optional :: SAL_CSp !< A pointer to the control structure of the !! SAL module. + type(harmonic_analysis_CS), target, optional :: HA_CSp !< A pointer to the control structure of the + !! harmonic analysis module ! This include declares and sets the variable "version". # include "version_variable.h" @@ -4489,6 +4500,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, if (present(SAL_CSp)) then CS%SAL_CSp => SAL_CSp endif + if (present(HA_CSp)) then + CS%HA_CSp => HA_CSp + endif ! Read all relevant parameters and write them to the model log. call get_param(param_file, mdl, "SPLIT", CS%split, default=.true., do_not_log=.true.) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index b545d99daa..11c3ff1873 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -47,6 +47,7 @@ module MOM_dynamics_split_RK2 use MOM_CoriolisAdv, only : CoriolisAdv_init, CoriolisAdv_end use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type +use MOM_harmonic_analysis, only : harmonic_analysis_CS use MOM_hor_index, only : hor_index_type use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS use MOM_hor_visc, only : hor_visc_init, hor_visc_end @@ -242,6 +243,8 @@ module MOM_dynamics_split_RK2 type(SAL_CS) :: SAL_CSp !> A pointer to the tidal forcing control structure type(tidal_forcing_CS) :: tides_CSp + !> A pointer to the harmonic analysis control structure + type(harmonic_analysis_CS) :: HA_CSp !> A pointer to the ALE control structure. type(ALE_CS), pointer :: ALE_CSp => NULL() @@ -482,7 +485,6 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f call cpu_clock_end(id_clock_pass) !--- end set up for group halo pass - ! PFu = d/dx M(h,T,S) ! pbce = dM/deta if (CS%begw == 0.0) call enable_averages(dt, Time_local, CS%diag) @@ -1308,7 +1310,7 @@ end subroutine remap_dyn_split_RK2_aux_vars !> This subroutine initializes all of the variables that are used by this !! dynamic core, including diagnostics and the cpu clocks. subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, & - diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, & + diag, CS, HA_CSp, restart_CS, dt, Accel_diag, Cont_diag, MIS, & VarMix, MEKE, thickness_diffuse_CSp, & OBC, update_OBC_CSp, ALE_CSp, set_visc, & visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil) @@ -1331,6 +1333,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p type(param_file_type), intent(in) :: param_file !< parameter file for parsing type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure + type(harmonic_analysis_CS), pointer :: HA_CSp !< A pointer to the control structure of the + !! harmonic analysis module type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure real, intent(in) :: dt !< time step [T ~> s] type(accel_diag_ptrs), target, intent(inout) :: Accel_diag !< points to momentum equation terms for @@ -1504,7 +1508,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv) if (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp) - if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp) + if (CS%use_tides) then + call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp) + HA_CSp => CS%HA_CSp + else + HA_CSp => NULL() + endif call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & CS%SAL_CSp, CS%tides_CSp) call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp) @@ -1538,7 +1547,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo call barotropic_init(u, v, h, CS%eta, Time, G, GV, US, param_file, diag, & - CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp) + CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp, CS%HA_CSp) if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. & .not. query_initialized(CS%diffv, "diffv", restart_CS)) then diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index c6642506f9..0220db7993 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -49,6 +49,7 @@ module MOM_dynamics_split_RK2b use MOM_CoriolisAdv, only : CoriolisAdv_init, CoriolisAdv_end use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type +use MOM_harmonic_analysis, only : harmonic_analysis_CS use MOM_hor_index, only : hor_index_type use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS use MOM_hor_visc, only : hor_visc_init, hor_visc_end @@ -239,6 +240,8 @@ module MOM_dynamics_split_RK2b type(SAL_CS) :: SAL_CSp !> A pointer to the tidal forcing control structure type(tidal_forcing_CS) :: tides_CSp + !> A pointer to the harmonic analysis control structure + type(harmonic_analysis_CS) :: HA_CSp !> A pointer to the ALE control structure. type(ALE_CS), pointer :: ALE_CSp => NULL() @@ -1233,7 +1236,7 @@ end subroutine remap_dyn_split_RK2b_aux_vars !> This subroutine initializes all of the variables that are used by this !! dynamic core, including diagnostics and the cpu clocks. subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, & - diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, & + diag, CS, HA_CSp, restart_CS, dt, Accel_diag, Cont_diag, MIS, & VarMix, MEKE, thickness_diffuse_CSp, & OBC, update_OBC_CSp, ALE_CSp, set_visc, & visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil) @@ -1256,6 +1259,8 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, type(param_file_type), intent(in) :: param_file !< parameter file for parsing type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics type(MOM_dyn_split_RK2b_CS), pointer :: CS !< module control structure + type(harmonic_analysis_CS), pointer :: HA_CSp !< A pointer to the control structure of the + !! harmonic analysis module type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure real, intent(in) :: dt !< time step [T ~> s] type(accel_diag_ptrs), target, intent(inout) :: Accel_diag !< points to momentum equation terms for @@ -1415,7 +1420,12 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv) if (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp) - if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp) + if (CS%use_tides) then + call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp) + HA_CSp => CS%HA_CSp + else + HA_CSp => NULL() + endif call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & CS%SAL_CSp, CS%tides_CSp) call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp) @@ -1453,7 +1463,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, call barotropic_init(u, v, h, CS%eta, Time, G, GV, US, param_file, diag, & CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, & - CS%SAL_CSp) + CS%SAL_CSp, CS%HA_CSp) flux_units = get_flux_units(GV) thickness_units = get_thickness_units(GV) diff --git a/src/diagnostics/MOM_harmonic_analysis.F90 b/src/diagnostics/MOM_harmonic_analysis.F90 new file mode 100644 index 0000000000..4a4b4ccd1f --- /dev/null +++ b/src/diagnostics/MOM_harmonic_analysis.F90 @@ -0,0 +1,430 @@ +!> Inline harmonic analysis (conventional) +module MOM_harmonic_analysis + +use MOM_time_manager, only : time_type, real_to_time, time_type_to_real, get_date, increment_date, & + operator(+), operator(-), operator(<), operator(>), operator(>=) +use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_file_parser, only : param_file_type, get_param +use MOM_io, only : file_exists, open_ASCII_file, READONLY_FILE, close_file, & + MOM_infra_file, vardesc, MOM_field, & + var_desc, create_MOM_file, SINGLE_FILE, MOM_write_field +use MOM_error_handler, only : MOM_mesg, MOM_error, NOTE + +implicit none ; private + +public HA_init, HA_register, HA_accum_FtF, HA_accum_FtSSH + +#include + +integer, parameter :: MAX_CONSTITUENTS = 10 !< The maximum number of tidal constituents + +!> The private control structure for storing the HA info of a particular field +type, private :: HA_type + character(len=16) :: key = "none" !< Name of the field of which harmonic analysis is to be performed + character(len=1) :: grid !< The grid on which the field is defined ('h', 'q', 'u', or 'v') + real :: old_time = -1.0 !< The time of the previous accumulating step [T ~> s] + real, allocatable :: ref(:,:) !< The initial field in arbitrary units [A] + real, allocatable :: FtSSH(:,:,:) !< Accumulator of (F' * SSH_in) in arbitrary units [A] + !>@{ Lower and upper bounds of input data + integer :: is, ie, js, je + !>@} +end type HA_type + +!> A linked list of control structures that store the HA info of different fields +type, private :: HA_node + type(HA_type) :: this !< Control structure of the current field in the list + type(HA_node), pointer :: next !< The list of other fields +end type HA_node + +!> The public control structure of the MOM_harmonic_analysis module +type, public :: harmonic_analysis_CS ; private + logical :: HAready = .true. !< If true, perform harmonic analysis + type(time_type) :: & + time_start, & !< Start time of harmonic analysis + time_end, & !< End time of 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] + 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 + character(len=16) :: const_name(MAX_CONSTITUENTS) !< The name of each constituent + character(len=255) :: path !< Path to directory where output will be written + type(unit_scale_type) :: US !< A dimensional unit scaling type + type(HA_node), pointer :: list => NULL() !< A linked list for storing the HA info of different fields +end type harmonic_analysis_CS + +contains + +!> 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) + 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] + 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 + + ! Local variables + type(HA_type) :: ha1 !< A temporary, null field used for initializing CS%list + real :: HA_start_time !< Start time of harmonic analysis [T ~> s] + real :: HA_end_time !< End time of harmonic analysis [T ~> s] + character(len=40) :: mdl="MOM_harmonic_analysis" !< This module's name + character(len=255) :: mesg + integer :: year, month, day, hour, minute, second + + ! Determine CS%time_start and CS%time_end + call get_param(param_file, mdl, "HA_START_TIME", HA_start_time, & + "Start time of harmonic analysis, in units of days after "//& + "the start of the current run segment. Must be smaller than "//& + "HA_END_TIME, otherwise harmonic analysis will not be performed. "//& + "If negative, |HA_START_TIME| determines the length of harmonic analysis, "//& + "and harmonic analysis will start |HA_START_TIME| days before HA_END_TIME, "//& + "or at the beginning of the run segment, whichever occurs later.", & + units="days", default=0.0, scale=86400.0*US%s_to_T) + call get_param(param_file, mdl, "HA_END_TIME", HA_end_time, & + "End time of harmonic analysis, in units of days after "//& + "the start of the current run segment. Must be positive "//& + "and smaller than the length of the currnet run segment, "//& + "otherwise harmonic analysis will not be performed.", & + units="days", default=0.0, scale=86400.0*US%s_to_T) + + if (HA_end_time <= 0.0) then + call MOM_mesg('MOM_harmonic_analysis: HA_END_TIME is zero or negative. '//& + 'Harmonic analysis will not be performed.') + CS%HAready = .false. ; return + endif + + if (HA_end_time <= HA_start_time) then + call MOM_mesg('MOM_harmonic_analysis: HA_END_TIME is smaller than or equal to HA_START_TIME. '//& + 'Harmonic analysis will not be performed.') + CS%HAready = .false. ; return + endif + + if (HA_start_time < 0.0) then + HA_start_time = HA_end_time + HA_start_time + if (HA_start_time <= 0.0) HA_start_time = 0.0 + endif + + CS%time_start = Time + real_to_time(US%T_to_s * HA_start_time) + CS%time_end = Time + real_to_time(US%T_to_s * HA_end_time) + + call get_date(Time, year, month, day, hour, minute, second) + write(mesg,*) "MOM_harmonic_analysis: run segment starts on ", year, month, day, hour, minute, second + call MOM_error(NOTE, trim(mesg)) + call get_date(CS%time_start, year, month, day, hour, minute, second) + write(mesg,*) "MOM_harmonic_analysis: harmonic analysis starts on ", year, month, day, hour, minute, second + call MOM_error(NOTE, trim(mesg)) + call get_date(CS%time_end, year, month, day, hour, minute, second) + write(mesg,*) "MOM_harmonic_analysis: harmonic analysis ends on ", year, month, day, hour, minute, second + call MOM_error(NOTE, trim(mesg)) + + ! Set path to directory where output will be written + call get_param(param_file, mdl, "HA_PATH", CS%path, & + "Path to output files for runtime harmonic analysis.", default="./") + + ! Populate some parameters of the control structure + CS%time_ref = time_ref + CS%freq = freq + CS%phase0 = phase0 + CS%nc = nc + CS%const_name = const_name + CS%length = 0 + CS%US = US + + allocate(CS%FtF(2*nc+1,2*nc+1), source=0.0) + + ! Initialize CS%list + allocate(CS%list) + CS%list%this = ha1 + nullify(CS%list%next) + +end subroutine HA_init + +!> This subroutine registers each of the fields on which HA is to be performed. +subroutine HA_register(key, grid, CS) + character(len=*), intent(in) :: key !< Name of the current field + character(len=1), intent(in) :: grid !< The grid on which the key field is defined + type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + type(HA_type) :: ha1 !< Control structure for the current field + type(HA_node), pointer :: tmp !< A temporary list to hold the current field + + if (.not. CS%HAready) return + + allocate(tmp) + ha1%key = trim(key) + ha1%grid = trim(grid) + tmp%this = ha1 + tmp%next => CS%list + CS%list => tmp + CS%length = CS%length + 1 + +end subroutine HA_register + +!> This subroutine accumulates the temporal basis functions in FtF. +!! The tidal constituents are those used in MOM_tidal_forcing, plus the mean (of zero frequency). +subroutine HA_accum_FtF(Time, CS) + type(time_type), intent(in) :: Time !< The current model time + type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + real :: now !< The relative time compared with tidal reference [T ~> s] + real :: cosomegat, sinomegat, ccosomegat, ssinomegat !< The components of the phase [nondim] + integer :: nc, c, icos, isin, cc, iccos, issin + + ! Exit the accumulator in the following cases + if (.not. CS%HAready) return + if (CS%length == 0) return + if (Time < CS%time_start) return + if (Time > CS%time_end) return + + nc = CS%nc + now = CS%US%s_to_T * time_type_to_real(Time - CS%time_ref) + + ! Accumulate FtF + CS%FtF(1,1) = CS%FtF(1,1) + 1.0 !< For the zero frequency + 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)) + CS%FtF(icos,1) = CS%FtF(icos,1) + cosomegat + CS%FtF(isin,1) = CS%FtF(isin,1) + sinomegat + CS%FtF(1,icos) = CS%FtF(icos,1) + CS%FtF(1,isin) = CS%FtF(isin,1) + do cc=c,nc + iccos = 2*cc + issin = 2*cc+1 + ccosomegat = cos(CS%freq(cc) * now + CS%phase0(cc)) + ssinomegat = sin(CS%freq(cc) * now + CS%phase0(cc)) + CS%FtF(icos,iccos) = CS%FtF(icos,iccos) + cosomegat * ccosomegat + CS%FtF(icos,issin) = CS%FtF(icos,issin) + cosomegat * ssinomegat + CS%FtF(isin,iccos) = CS%FtF(isin,iccos) + sinomegat * ccosomegat + CS%FtF(isin,issin) = CS%FtF(isin,issin) + sinomegat * ssinomegat + enddo ! cc=c,nc + enddo ! c=1,nc + +end subroutine HA_accum_FtF + +!> This subroutine accumulates the temporal basis functions in FtSSH and then calls HA_write to compute +!! harmonic constants and write results. The tidal constituents are those used in MOM_tidal_forcing, plus the +!! mean (of zero frequency). +subroutine HA_accum_FtSSH(key, data, Time, G, CS) + character(len=*), intent(in) :: key !< Name of the current field + real, dimension(:,:), intent(in) :: data !< Input data of which harmonic analysis is to be performed [A] + type(time_type), intent(in) :: Time !< The current model time + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + type(HA_type), pointer :: ha1 + type(HA_node), pointer :: tmp + real :: now !< The relative time compared with the tidal reference [T ~> s] + real :: dt !< The current time step size of the accumulator [T ~> s] + real :: cosomegat, sinomegat !< The components of the phase [nondim] + integer :: nc, i, j, k, c, icos, isin, is, ie, js, je + character(len=128) :: mesg + + ! Exit the accumulator in the following cases + if (.not. CS%HAready) return + if (CS%length == 0) return + if (Time < CS%time_start) return + if (Time > CS%time_end) return + + ! Loop through the full list to find the current field + tmp => CS%list + do k=1,CS%length + ha1 => tmp%this + if (trim(key) == trim(ha1%key)) exit + tmp => tmp%next + if (k == CS%length) return !< Do not perform harmonic analysis of a field that is not registered + enddo + + nc = CS%nc + now = CS%US%s_to_T * time_type_to_real(Time - CS%time_ref) + + ! Additional processing at the initial accumulating step + if (ha1%old_time < 0.0) then + ha1%old_time = now + + write(mesg,*) "MOM_harmonic_analysis: initializing accumulator, key = ", trim(ha1%key) + call MOM_error(NOTE, trim(mesg)) + + ! Get the lower and upper bounds of input data + ha1%is = LBOUND(data,1) ; is = ha1%is + ha1%ie = UBOUND(data,1) ; ie = ha1%ie + ha1%js = LBOUND(data,2) ; js = ha1%js + ha1%je = UBOUND(data,2) ; je = ha1%je + + allocate(ha1%ref(is:ie,js:je), source=0.0) + allocate(ha1%FtSSH(is:ie,js:je,2*nc+1), source=0.0) + ha1%ref(:,:) = data(:,:) + endif + + dt = now - ha1%old_time + ha1%old_time = now !< Keep track of time so we know when Time approaches CS%time_end + + is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je + + ! Accumulate FtF and FtSSH + 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)) + do j=js,je ; do i=is,ie + ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j)) + ha1%FtSSH(i,j,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat + ha1%FtSSH(i,j,isin) = ha1%FtSSH(i,j,isin) + (data(i,j) - ha1%ref(i,j)) * sinomegat + enddo ; enddo + enddo ! c=1,nc + + ! Compute harmonic constants and write output as Time approaches CS%time_end + ! This guarantees that HA_write will be called before Time becomes larger than CS%time_end + if (time_type_to_real(CS%time_end - Time) <= dt) then + call HA_write(ha1, Time, G, CS) + + write(mesg,*) "MOM_harmonic_analysis: harmonic analysis done, key = ", trim(ha1%key) + call MOM_error(NOTE, trim(mesg)) + + ! De-register the current field and deallocate memory + ha1%key = 'none' + deallocate(ha1%ref) + deallocate(ha1%FtSSH) + endif + +end subroutine HA_accum_FtSSH + +!> This subroutine computes the harmonic constants and write output for the current field +subroutine HA_write(ha1, Time, G, CS) + type(HA_type), pointer, intent(in) :: ha1 !< Control structure for the current field + type(time_type), intent(in) :: Time !< The current model time + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(harmonic_analysis_CS), intent(in) :: CS !< Control structure of the MOM_harmonic_analysis module + + ! Local variables + real, dimension(:,:,:), allocatable :: FtSSHw !< An array containing the harmonic constants [A] + integer :: year, month, day, hour, minute, second + integer :: nc, k, is, ie, js, je + + character(len=255) :: filename !< Output file name + type(MOM_infra_file) :: cdf !< The file handle for output harmonic constants + type(vardesc), allocatable :: cdf_vars(:) !< Output variable names + type(MOM_field), allocatable :: cdf_fields(:) !< Field type variables for the output fields + + nc = CS%nc ; is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je + + allocate(FtSSHw(is:ie,js:je,2*nc+1), source=0.0) + + ! Compute the harmonic coefficients + call HA_solver(ha1, nc, CS%FtF, FtSSHw) + + ! Output file name + call get_date(Time, year, month, day, hour, minute, second) + write(filename, '(a,"HA_",a,i0.4,i0.2,i0.2,".nc")') & + trim(CS%path), trim(ha1%key), year, month, day + + allocate(cdf_vars(2*nc+1)) + allocate(cdf_fields(2*nc+1)) + + ! Variable names + cdf_vars(1) = var_desc("z0", "m" ,"mean value", ha1%grid, '1', '1') + do k=1,nc + cdf_vars(2*k ) = var_desc(trim(CS%const_name(k))//"cos", "m", "cosine coefficient", ha1%grid, '1', '1') + cdf_vars(2*k+1) = var_desc(trim(CS%const_name(k))//"sin", "m", "sine coefficient", ha1%grid, '1', '1') + enddo + + ! Create output file + call create_MOM_file(cdf, trim(filename), cdf_vars, & + 2*nc+1, cdf_fields, SINGLE_FILE, 86400.0, G=G) + + ! Write data + call MOM_write_field(cdf, cdf_fields(1), G%domain, FtSSHw(:,:,1), 0.0) + do k=1,nc + call MOM_write_field(cdf, cdf_fields(2*k ), G%domain, FtSSHw(:,:,2*k ), 0.0) + call MOM_write_field(cdf, cdf_fields(2*k+1), G%domain, FtSSHw(:,:,2*k+1), 0.0) + enddo + + call cdf%flush() + deallocate(cdf_vars) + deallocate(cdf_fields) + deallocate(FtSSHw) + +end subroutine HA_write + +!> This subroutine computes the harmonic constants (stored in FtSSHw) using the dot products of the temporal +!! basis functions accumulated in FtF, and the dot products of the SSH (or other fields) with the temporal basis +!! functions accumulated in FtSSH. The system is solved by Cholesky decomposition. +subroutine HA_solver(ha1, nc, FtF, FtSSHw) + type(HA_type), pointer, intent(in) :: ha1 !< Control structure for the current field + integer, intent(in) :: nc !< Number of harmonic constituents + real, dimension(:,:), intent(in) :: FtF !< Accumulator of (F' * F) for all fields [nondim] + real, dimension(:,:,:), allocatable, intent(out) :: FtSSHw !< Work array for Cholesky decomposition [A] + + ! Local variables + real, dimension(:,:), allocatable :: tmp !< Work array for Cholesky decomposition [A] + real, dimension(:,:), allocatable :: FtFw !< Work array for Cholesky decomposition [nondim] + integer :: k, l, is, ie, js, je + + is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je + + allocate(tmp(is:ie,js:je), source=0.0) + allocate(FtFw(1:2*nc+1,1:2*nc+1), source=0.0) + allocate(FtSSHw(is:ie,js:je,2*nc+1), source=0.0) + + FtFw(:,:) = 0.0 + do l=1,2*nc+1 + FtFw(l,l) = sqrt(FtF(l,l) - dot_product(FtFw(l,1:l-1), FtFw(l,1:l-1))) + do k=l+1,2*nc+1 + FtFw(k,l) = (FtF(k,l) - dot_product(FtFw(k,1:l-1), FtFw(l,1:l-1))) / FtFw(l,l) + enddo + enddo + + FtSSHw(:,:,:) = ha1%FtSSH(:,:,:) + do k=1,2*nc+1 + tmp(:,:) = 0.0 + do l=1,k-1 + tmp(:,:) = tmp(:,:) + FtFw(k,l) * FtSSHw(:,:,l) + enddo + FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp(:,:)) / FtFw(k,k) + enddo + do k=2*nc+1,1,-1 + tmp(:,:) = 0.0 + do l=k+1,2*nc+1 + tmp(:,:) = tmp(:,:) + FtSSHw(:,:,l) * FtFw(l,k) + enddo + FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp(:,:)) / FtFw(k,k) + enddo + + deallocate(tmp) + deallocate(FtFw) + +end subroutine HA_solver + +!> \namespace harmonic_analysis +!! +!! This module computes the harmonic constants which can be used to reconstruct the tidal elevation (or other +!! 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 +!! number of tidal constituents) containing the cosine/sine functions for each frequency evaluated at each time +!! step, and x is a 2*nc-by-1 vector containing the constant coefficients of the sine/cosine for each constituent +!! (i.e., the harmonic constants). At each grid point, the harmonic constants are computed using least squares, +!! +!! x = (F' * F)^{-1} * (F' * SSH_in), +!! +!! where the prime denotes matrix transpose, and SSH_in is the sea surface height (or other fields) determined by +!! the model. The dot products (F' * F) and (F' * SSH_in) are computed by accumulating the sums as the model is +!! running and stored in the arrays FtF and FtSSH, respectively. The FtF matrix is inverted as needed before +!! computing and writing out the harmonic constants. +!! +!! Ed Zaron and William Xu (chengzhu.xu@oregonstate.edu), April 2024. + +end module MOM_harmonic_analysis + diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 1cd8a45a78..d53c3321f9 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -9,6 +9,8 @@ module MOM_tidal_forcing use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type +use MOM_harmonic_analysis, & + only : HA_init, HA_register, harmonic_analysis_CS use MOM_io, only : field_exists, file_exists, MOM_read_data use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) use MOM_unit_scaling, only : unit_scale_type @@ -231,12 +233,13 @@ end subroutine nodal_fu !! while fields like the background viscosities are 2-D arrays. !! ALLOC is a macro defined in MOM_memory.h for allocate or nothing with !! static memory. -subroutine tidal_forcing_init(Time, G, US, param_file, CS) +subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) type(time_type), intent(in) :: Time !< The current model time. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. 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. type(tidal_forcing_CS), intent(inout) :: CS !< Tidal forcing control structure + type(harmonic_analysis_CS), optional, intent(out) :: HA_CS !< Control structure for harmonic analysis ! Local variables real, dimension(SZI_(G), SZJ_(G)) :: & @@ -251,6 +254,7 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS) 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 :: HA_ssh, HA_ubt, HA_vbt ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_tidal_forcing" ! This module's name. @@ -523,6 +527,19 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS) enddo endif + 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 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) + call get_param(param_file, mdl, "HA_UBT", HA_ubt, & + "If true, perform harmonic analysis of zonal barotropic velocity.", default=.false.) + if (HA_ubt) call HA_register('ubt', 'u', HA_CS) + call get_param(param_file, mdl, "HA_VBT", HA_vbt, & + "If true, perform harmonic analysis of meridional barotropic velocity.", default=.false.) + if (HA_vbt) call HA_register('vbt', 'v', HA_CS) + endif + id_clock_tides = cpu_clock_id('(Ocean tides)', grain=CLOCK_MODULE) end subroutine tidal_forcing_init