+Steps to correct Kelvin_initialization issues

  This commit takes the first steps toward correcting problems with the Kelvin
initialization code, including making the OBC nudging timescale an explicit
parameter and adding a number of comments identifying concerns and suggesting
possible improvements with the code that is used when KELVIN_WAVE_MODE > 0.
This commit includes the introduction of the new runtime parameter
KELVIN_WAVE_VEL_NUDGING_TIMESCALE, with a default value set to recover the
current hard-coded answers.  A new warning message cautions about the suspected
problems with the internal wave version of the code.  By default, all answers
are bitwise identical in the 4 Kelvin wave test cases that can be found at
github.com/ESMG/ESMG-configs, but there is a new parameter in the
MOM_parameter_doc.all files these cases.

src/user/Kelvin_initialization.F90

@@ -45,6 +45,9 @@ module Kelvin_initialization
   real :: wave_period       !< Period of the mode-0 waves [T ~> s]
   real :: ssh_amp           !< Amplitude of the sea surface height forcing for mode-0 waves [Z ~> m]
   real :: inflow_amp        !< Amplitude of the boundary velocity forcing for internal waves [L T-1 ~> m s-1]
+  real :: OBC_nudging_time  !< The timescale with which the inflowing open boundary velocities are nudged toward
+                            !! their intended values with the Kelvin wave test case [T ~> s], or a negative
+                            !! value to retain the value that is set when the OBC segments are initialized.
 end type Kelvin_OBC_CS
 
 ! This include declares and sets the variable "version".
@@ -114,10 +117,21 @@ function register_Kelvin_OBC(param_file, CS, US, OBC_Reg)
                    "at the open boundaries.", units="m s-1", default=1.0, scale=US%m_s_to_L_T)
   endif
 
+  call get_param(param_file, mdl, "KELVIN_WAVE_VEL_NUDGING_TIMESCALE", CS%OBC_nudging_time, &
+                 "The timescale with which the inflowing open boundary velocities are nudged toward "//&
+                 "their intended values with the Kelvin wave test case, or a negative value to keep "//&
+                 "the value that is set when the OBC segments are initialized.", &
+                 units="s", default=1.0/(0.3*86400.), scale=US%s_to_T)
+                 !### Change the default nudging timescale to -1. or another value?
+
   ! Register the Kelvin open boundary.
   call register_OBC(casename, param_file, OBC_Reg)
   register_Kelvin_OBC = .true.
 
+  if (CS%mode > 0) call MOM_error(WARNING, &
+      "register_Kelvin_OBC: The Kelvin_initialization code is not yet working properly unless KELVIN_WAVE_MODE = 0.")
+  ! TODO: Revisit and correct the internal Kelvin wave test case.
+
 end function register_Kelvin_OBC
 
 !> Clean up the Kelvin wave OBC from registry.
@@ -193,7 +207,7 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
   real :: time_sec     ! The time in the run [T ~> s]
   real :: cff          ! The wave speed [L T-1 ~> m s-1]
   real :: N0           ! Brunt-Vaisala frequency times a rescaling of slopes [L Z-1 T-1 ~> s-1]
-  real :: lambda       ! Offshore decay scale [L-1 ~> m-1]
+  real :: lambda       ! Offshore decay scale, i.e. the inverse of the deformation radius of a mode [L-1 ~> m-1]
   real :: omega        ! Wave frequency [T-1 ~> s-1]
   real :: PI           ! The ratio of the circumference of a circle to its diameter [nondim]
   real :: depth_tot(SZI_(G),SZJ_(G))  ! The total depth of the ocean [Z ~> m]
@@ -233,13 +247,17 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
     omega = 2.0 * PI / CS%wave_period
     val1 = sin(omega * time_sec)
   else
+    ! This is supposed to be a linear internal Kelvin wave case, so I do not understand the purpose
+    ! of the squared velocity here. -RWH
     mag_int = CS%inflow_amp**2
     N0 = sqrt((CS%rho_range / CS%rho_0) * (GV%g_Earth / CS%H0))
     lambda = PI * CS%mode * CS%F_0 / (CS%H0 * N0)
     ! Two wavelengths in domain
-    omega = (4.0 * CS%H0 * N0)  / (CS%mode * US%m_to_L*G%len_lon)
-    !### There is a bug here when len_lon is in km.  This should be
-    ! omega = (4.0 * CS%H0 * N0)  / (CS%mode * G%grid_unit_to_L*G%len_lon)
+    ! The reason for the factor of 0.001 is unclear, but it is needed to recreate the previous answers. -RWH
+    omega = (4.0 * CS%H0 * N0)  / (CS%mode * (0.001*G%grid_unit_to_L)*G%len_lon)
+    ! If the modal wave speed were calculated via wave_speeds(), we should have
+    !   lambda = CS%F_0 / CS%cg_mode
+    !   omega = (4.0 * PI / (G%grid_unit_to_L*G%len_lon)) * CS%cg_mode
   endif
 
   sina = sin(CS%coast_angle)
@@ -251,9 +269,9 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
     if (segment%direction == OBC_DIRECTION_E) cycle
     if (segment%direction == OBC_DIRECTION_N) cycle
 
-    ! This should be somewhere else...
-    !### This is supposed to be a timescale [T ~> s] but appears to be a rate in [s-1].
-    segment%Velocity_nudging_timescale_in = US%s_to_T * 1.0/(0.3*86400)
+    ! If OBC_nudging_time is negative, the value of Velocity_nudging_timescale_in that was set
+    ! when the segments are initialized is retained.
+    if (CS%OBC_nudging_time >= 0.0) segment%Velocity_nudging_timescale_in = CS%OBC_nudging_time
 
     if (segment%direction == OBC_DIRECTION_W) then
       IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB
@@ -281,11 +299,22 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
             enddo
           endif
         else
-          ! Baroclinic, not rotated yet
+          ! Baroclinic, not rotated yet (and apparently not working as intended yet).
           segment%SSH(I,j) = 0.0
           segment%normal_vel_bt(I,j) = 0.0
+          ! I suspect that the velocities in both of the following loops should instead be
+          !   normal_vel(I,j,k) = CS%inflow_amp * CS%u_struct(k) * exp(-lambda * y) * cos(omega * time_sec)
+          ! In addition, there should be a specification of the interface-height anomalies at the
+          ! open boundaries that are specified as something like
+          !   eta_anom(I,j,K) = (CS%inflow_amp*depth_tot/CS%cg_mode) * CS%w_struct(K) * &
+          !                     exp(-lambda * y) * cos(omega * time_sec)
+          ! In these expressions CS%u_struct and CS%w_struct could be returned from the subroutine wave_speeds
+          ! in MOM_wave_speed() based on the horizontally uniform initial state.
           if (segment%nudged) then
             do k=1,nz
+              ! Note that mag_int is the square of the specified inflow amplitude, in [L2 T-2 ~> m2 s-2].
+              ! The following expression is dimensionally correct, but I do not understand why the
+              ! normal velocities should scale with the square of their amplitude. -RWH
               segment%nudged_normal_vel(I,j,k) = mag_int * lambda / CS%F_0 * &
                    exp(-lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * &
                    cos(omega * time_sec)
@@ -339,7 +368,7 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
             enddo
           endif
         else
-          ! Not rotated yet
+          ! Not rotated yet (also see the notes above on how this case might be improved)
           segment%SSH(i,J) = 0.0
           segment%normal_vel_bt(i,J) = 0.0
           if (segment%nudged) then