Inline Harmonic Analysis

Update 6: The frequencies of 8 overtides/compound tides (MK3, MN4,
M4, MS4, S4, M6, S6, M8) have been added for the harmonic analysis.

Author: c2xu

src/core/MOM.F90

@@ -1064,7 +1064,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('ssh', ssh, Time_local, CS%HA_CSp)
+    if (associated(CS%HA_CSp)) call HA_accum('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)

src/core/MOM_barotropic.F90

@@ -1935,8 +1935,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
   ! 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('ubt', ubt, CS%Time, CS%HA_CSp)
-    call HA_accum('vbt', vbt, CS%Time, CS%HA_CSp)
+    call HA_accum('ubt', ubt, CS%Time, G, CS%HA_CSp)
+    call HA_accum('vbt', vbt, CS%Time, G, CS%HA_CSp)
   endif
 
   if (id_clock_calc_post > 0) call cpu_clock_end(id_clock_calc_post)

src/core/MOM_dynamics_split_RK2.F90

@@ -1573,7 +1573,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p
   if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp, restart_CS)
   if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   if (CS%use_HA) then
-    call HA_init(Time, G, US, param_file, nc, CS%HA_CSp)
+    call HA_init(Time, US, param_file, nc, CS%HA_CSp)
     HA_CSp => CS%HA_CSp
   else
     HA_CSp => NULL()

src/core/MOM_dynamics_split_RK2b.F90

@@ -1456,7 +1456,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US,
   if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp, restart_CS)
   if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   if (CS%use_HA) then
-    call HA_init(Time, G, US, param_file, nc, CS%HA_CSp)
+    call HA_init(Time, US, param_file, nc, CS%HA_CSp)
     HA_CSp => CS%HA_CSp
   else
     HA_CSp => NULL()

src/diagnostics/MOM_harmonic_analysis.F90

@@ -1,15 +1,15 @@
 !> Inline harmonic analysis (conventional)
 module MOM_harmonic_analysis
 
-use MOM_time_manager,  only : time_type, real_to_time, time_type_to_real, &
-                              set_date, get_date, increment_date, &
-                              operator(+), operator(-), operator(<), operator(>), operator(>=)
+use MOM_time_manager,  only : time_type, real_to_time, time_type_to_real
+use MOM_time_manager,  only : set_date, get_date, increment_date
+use MOM_time_manager,  only : 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_io,            only : file_exists, open_ASCII_file, READONLY_FILE, close_file
+use MOM_io,            only : MOM_infra_file, vardesc, MOM_field
+use MOM_io,            only : var_desc, create_MOM_file, SINGLE_FILE, MOM_write_field
 use MOM_error_handler, only : MOM_mesg, MOM_error, NOTE
 use MOM_tidal_forcing, only : astro_longitudes, astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency
 
@@ -54,7 +54,6 @@ module MOM_harmonic_analysis
   integer :: length                          !< Number of fields of which harmonic analysis is to be performed
   character(len=4), allocatable, dimension(:) :: const_name !< The name of each constituent
   character(len=255) :: path                 !< Path to directory where output will be written
-  type(ocean_grid_type)  :: G                !< The ocean's grid structure
   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
@@ -63,9 +62,8 @@ module MOM_harmonic_analysis
 
 !> This subroutine sets static variables used by this module and initializes CS%list.
 !! THIS MUST BE CALLED AT THE END OF tidal_forcing_init.
-subroutine HA_init(Time, G, US, param_file, nc, CS)
+subroutine HA_init(Time, US, param_file, nc, CS)
   type(time_type),       intent(in)  :: Time        !< The current model time
-  type(ocean_grid_type), intent(in)  :: 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
   integer,               intent(in)  :: nc          !< The number of tidal constituents in use
@@ -164,7 +162,28 @@ subroutine HA_init(Time, G, US, param_file, nc, CS)
                    CS%freq(c), "Frequency of the "//trim(CS%const_name(c))//&
                    " constituent. This is used if USE_HA is true and "//trim(CS%const_name(c))//&
                    " is in HA_CONSTITUENTS.", units="rad s-1", scale=US%T_to_s, default=0.0)
-    if (CS%freq(c)<=0.0) CS%freq(c) = tidal_frequency(trim(CS%const_name(c)))
+    if (CS%freq(c)<=0.0) then
+      select case (trim(CS%const_name(c)))
+        case ('M4')
+          CS%freq(c) = tidal_frequency('M2') * 2
+        case ('M6')
+          CS%freq(c) = tidal_frequency('M2') * 3
+        case ('M8')
+          CS%freq(c) = tidal_frequency('M2') * 4
+        case ('S4')
+          CS%freq(c) = tidal_frequency('S2') * 2
+        case ('S6')
+          CS%freq(c) = tidal_frequency('S2') * 3
+        case ('MK3')
+          CS%freq(c) = tidal_frequency('M2') + tidal_frequency('K1')
+        case ('MS4')
+          CS%freq(c) = tidal_frequency('M2') + tidal_frequency('S2')
+        case ('MN4')
+          CS%freq(c) = tidal_frequency('M2') + tidal_frequency('N2')
+        case default
+          CS%freq(c) = tidal_frequency(trim(CS%const_name(c)))
+      end select
+    endif
 
     call get_param(param_file, mdl, "HA_"//trim(CS%const_name(c))//"_PHASE_T0", CS%phase0(c), &
                    "Phase of the "//trim(CS%const_name(c))//" tidal constituent at time 0. "//&
@@ -236,7 +255,6 @@ subroutine HA_init(Time, G, US, param_file, nc, CS)
   ! Populate some parameters of the control structure
   CS%nc         =  nc
   CS%length     =  0
-  CS%G          =  G
   CS%US         =  US
 
   ! Initialize CS%list
@@ -283,10 +301,11 @@ end subroutine HA_register
 !! harmonic constants and write results. The tidal constituents are those used in MOM_tidal_forcing, plus the
 !! mean (of zero frequency). For FtF, only the main diagonal and entries below it are calculated, which are needed
 !! for Cholesky decomposition.
-subroutine HA_accum(key, data, Time, CS)
+subroutine HA_accum(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
@@ -389,7 +408,7 @@ subroutine HA_accum(key, data, Time, CS)
   !!! 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, CS)
+    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))
@@ -403,9 +422,10 @@ subroutine HA_accum(key, data, Time, CS)
 end subroutine HA_accum
 
 !> This subroutine computes the harmonic constants and write output for the current field
-subroutine HA_write(ha1, Time, CS)
+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
@@ -442,18 +462,18 @@ subroutine HA_write(ha1, Time, CS)
 
   ! Create output file
   call create_MOM_file(cdf, trim(filename), cdf_vars, &
-                       2*nc+1, cdf_fields, SINGLE_FILE, 86400.0, G=CS%G)
+                       2*nc+1, cdf_fields, SINGLE_FILE, 86400.0, G=G)
 
   ! Add the initial field back to the mean state
   do j=js,je ; do i=is,ie
     FtSSHw(i,j,1) = FtSSHw(i,j,1) + ha1%ref(i,j)
   enddo ; enddo
 
   ! Write data
-  call MOM_write_field(cdf, cdf_fields(1), CS%G%domain, FtSSHw(:,:,1), 0.0)
+  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  ), CS%G%domain, FtSSHw(:,:,2*k  ), 0.0)
-    call MOM_write_field(cdf, cdf_fields(2*k+1), CS%G%domain, FtSSHw(:,:,2*k+1), 0.0)
+    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()