diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 118d76ac93..2edd508ca8 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/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