Zero accel

src/core/MOM_barotropic.F90

@@ -2040,6 +2040,19 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
 !$OMP parallel do default(none) shared(is,ie,js,je,nz,accel_layer_u,u_accel_bt,pbce,gtot_W, &
 !$OMP                                  e_anom,gtot_E,CS,accel_layer_v,v_accel_bt,      &
 !$OMP                                  gtot_S,gtot_N)
+  if (apply_OBCs) then
+!   call open_boundary_set_bt_accel(OBC, G, u_accel_bt, v_accel_bt)
+    if (CS%BT_OBC%apply_u_OBCs) then ; do j=js,je ; do I=is-1,ie
+      if (OBC%segnum_u(I,j) /= OBC_NONE) then
+        u_accel_bt(I,j) = (ubt_wtd(I,j) - ubt_Cor(I,j)) / dt
+      endif
+    enddo ; enddo ; endif
+    if (CS%BT_OBC%apply_v_OBCs) then ; do J=js-1,je ; do i=is,ie
+      if (OBC%segnum_v(i,J) /= OBC_NONE) then
+        v_accel_bt(i,J) = (vbt_wtd(i,J) - vbt_Cor(i,J)) / dt
+      endif
+    enddo ; enddo ; endif
+  endif
   do k=1,nz
     do j=js,je ; do I=is-1,ie
       accel_layer_u(I,j,k) = u_accel_bt(I,j) - &

src/core/MOM_dynamics_legacy_split.F90

@@ -111,6 +111,7 @@ module MOM_dynamics_legacy_split
 use MOM_MEKE_types, only : MEKE_type
 use MOM_open_boundary, only : ocean_OBC_type
 use MOM_open_boundary, only : radiation_open_bdry_conds
+use MOM_open_boundary, only : open_boundary_zero_normal_flow
 use MOM_boundary_update, only : update_OBC_data, update_OBC_CS
 use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS
 use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS
@@ -511,6 +512,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
   do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
     v_bc_accel(i,J,k) = (CS%Cav(i,J,k) + CS%PFv(i,J,k)) + CS%diffv(i,J,k)
   enddo ; enddo ; enddo
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, u_bc_accel, v_bc_accel)
+  endif
   call cpu_clock_end(id_clock_btforce)
 
   if (CS%debug) then
@@ -828,6 +832,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
   do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
     v_bc_accel(i,J,k) = (CS%Cav(i,J,k) + CS%PFv(i,J,k)) + CS%diffv(i,J,k)
   enddo ; enddo ; enddo
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, u_bc_accel, v_bc_accel)
+  endif
   call cpu_clock_end(id_clock_btforce)
 
   if (CS%debug) then

src/core/MOM_dynamics_split_RK2.F90

@@ -48,8 +48,8 @@ module MOM_dynamics_split_RK2
 use MOM_interface_heights,     only : find_eta
 use MOM_lateral_mixing_coeffs, only : VarMix_CS
 use MOM_MEKE_types,            only : MEKE_type
-use MOM_open_boundary,         only : ocean_OBC_type
-use MOM_open_boundary,         only : radiation_open_bdry_conds
+use MOM_open_boundary,         only : ocean_OBC_type, radiation_open_bdry_conds
+use MOM_open_boundary,         only : open_boundary_zero_normal_flow
 use MOM_PressureForce,         only : PressureForce, PressureForce_init, PressureForce_CS
 use MOM_set_visc,              only : set_viscous_BBL, set_viscous_ML, set_visc_CS
 use MOM_tidal_forcing,         only : tidal_forcing_init, tidal_forcing_CS
@@ -451,6 +451,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
       v_bc_accel(i,J,k) = (CS%Cav(i,J,k) + CS%PFv(i,J,k)) + CS%diffv(i,J,k)
     enddo ; enddo
   enddo
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, u_bc_accel, v_bc_accel)
+  endif
   call cpu_clock_end(id_clock_btforce)
 
   if (CS%debug) then
@@ -726,6 +729,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
       v_bc_accel(i,J,k) = (CS%Cav(i,J,k) + CS%PFv(i,J,k)) + CS%diffv(i,J,k)
     enddo ; enddo
   enddo
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, u_bc_accel, v_bc_accel)
+  endif
   call cpu_clock_end(id_clock_btforce)
 
   if (CS%debug) then

src/core/MOM_dynamics_unsplit.F90

@@ -105,6 +105,7 @@ module MOM_dynamics_unsplit
 use MOM_MEKE_types, only : MEKE_type
 use MOM_open_boundary, only : ocean_OBC_type
 use MOM_open_boundary, only : radiation_open_bdry_conds
+use MOM_open_boundary, only : open_boundary_zero_normal_flow
 use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS
 use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS
 use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS
@@ -330,6 +331,10 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, &
   if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then
     call update_OBC_data(CS%OBC, G, GV, tv, h, CS%update_OBC_CSp, Time_local)
   endif; endif
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%PFu, CS%PFv)
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%CAu, CS%CAv)
+  endif
 
 ! up = u + dt_pred * (PFu + CAu)
   call cpu_clock_begin(id_clock_mom_update)
@@ -418,6 +423,10 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, &
   if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then
     call update_OBC_data(CS%OBC, G, GV, tv, h, CS%update_OBC_CSp, Time_local)
   endif; endif
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%PFu, CS%PFv)
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%CAu, CS%CAv)
+  endif
 
 ! upp = u + dt/2 * ( PFu + CAu )
   call cpu_clock_begin(id_clock_mom_update)
@@ -499,6 +508,10 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, &
   endif; endif
 
 ! u = u + dt * ( PFu + CAu )
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%PFu, CS%PFv)
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%CAu, CS%CAv)
+  endif
   do k=1,nz ; do j=js,je ; do I=Isq,Ieq
     u(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * &
             (CS%PFu(I,j,k) + CS%CAu(I,j,k)))

src/core/MOM_dynamics_unsplit_RK2.F90

@@ -102,6 +102,7 @@ module MOM_dynamics_unsplit_RK2
 use MOM_MEKE_types, only : MEKE_type
 use MOM_open_boundary, only : ocean_OBC_type
 use MOM_open_boundary, only : radiation_open_bdry_conds
+use MOM_open_boundary, only : open_boundary_zero_normal_flow
 use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS
 use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS
 use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS
@@ -324,6 +325,11 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
   if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then
     call update_OBC_data(CS%OBC, G, GV, tv, h_in, CS%update_OBC_CSp, Time_local)
   endif; endif
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%PFu, CS%PFv)
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%CAu, CS%CAv)
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%diffu, CS%diffv)
+  endif
 
 ! up+[n-1/2] = u[n-1] + dt_pred * (PFu + CAu)
   call cpu_clock_begin(id_clock_mom_update)
@@ -400,6 +406,9 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
   call CorAdCalc(up, vp, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, &
                  G, GV, CS%CoriolisAdv_CSp)
   call cpu_clock_end(id_clock_Cor)
+  if (associated(CS%OBC)) then
+    call open_boundary_zero_normal_flow(CS%OBC, G, CS%CAu, CS%CAv)
+  endif
 
 ! call enable_averaging(dt,Time_local, CS%diag)  ?????????????????????/
 

src/user/Kelvin_initialization.F90

@@ -43,6 +43,15 @@ module Kelvin_initialization
   integer :: mode = 0         !< Vertical mode
   real    :: coast_angle = 0  !< Angle of coastline
   real    :: coast_offset = 0 !< Longshore distance to coastal angle
+  real    :: N0 = 0           !< Brunt-Vaisala frequency
+  real    :: H0 = 0           !< Bottom depth
+  real    :: F_0              !< Coriolis parameter
+  real    :: plx = 0          !< Longshore wave parameter
+  real    :: pmz = 0          !< Vertical wave parameter
+  real    :: lambda = 0       !< Vertical wave parameter
+  real    :: omega            !< Frequency
+  real    :: rho_range        !< Density range
+  real    :: rho_0            !< Mean density
 end type Kelvin_OBC_CS
 
 ! This include declares and sets the variable "version".
@@ -71,6 +80,8 @@ function register_Kelvin_OBC(param_file, CS, OBC_Reg)
   call get_param(param_file, mod, "KELVIN_WAVE_MODE", CS%mode, &
                  "Vertical Kelvin wave mode imposed at upstream open boundary.", &
                  default=0)
+  call get_param(param_file, mod, "F_0", CS%F_0, &
+                 default=0.0, do_not_log=.true.)
   call get_param(param_file, mod, "TOPO_CONFIG", config, do_not_log=.true.)
   if (trim(config) == "Kelvin") then
     call get_param(param_file, mod, "KELVIN_COAST_OFFSET", CS%coast_offset, &
@@ -81,6 +92,15 @@ function register_Kelvin_OBC(param_file, CS, OBC_Reg)
                    "The angle of the southern bondary beyond X=KELVIN_COAST_OFFSET.", &
                    units="degrees", default=11.3)
     CS%coast_angle = CS%coast_angle * (atan(1.0)/45.) ! Convert to radians
+    CS%coast_offset = CS%coast_offset * 1.e3          ! Convert to m
+  endif
+  if (CS%mode /= 0) then
+    call get_param(param_file, mod, "DENSITY_RANGE", CS%rho_range, &
+                   default=2.0, do_not_log=.true.)
+    call get_param(param_file, mod, "RHO_0", CS%rho_0, &
+                   default=1035.0, do_not_log=.true.)
+    call get_param(param_file, mod, "MAXIMUM_DEPTH", CS%H0, &
+                   default=1000.0, do_not_log=.true.)
   endif
 
   ! Register the Kelvin open boundary.
@@ -164,7 +184,7 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, h, Time)
   real :: PI
   integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz
   integer :: IsdB, IedB, JsdB, JedB
-  real    :: fac, omega, x, y, x1, y1
+  real    :: fac, x, y, x1, y1
   real    :: val1, val2, sina, cosa
   type(OBC_segment_type), pointer :: segment
 
@@ -177,13 +197,18 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, h, Time)
 
   time_sec = time_type_to_real(Time)
   PI = 4.0*atan(1.0)
+  fac = 1.0
 
   if (CS%mode == 0) then
-    fac = 1.0
-    omega = 2.0 * PI / (12.42 * 3600.0)      ! M2 Tide period
-    val1 = sin(omega * time_sec)
+    CS%omega = 2.0 * PI / (12.42 * 3600.0)      ! M2 Tide period
+    val1 = sin(CS%omega * time_sec)
   else
-    omega = 2.0 * PI / (12.42 * 3600.0)      ! M2 Tide period
+    CS%N0 = sqrt(CS%rho_range / CS%rho_0 * G%g_Earth * CS%H0)
+    ! Two wavelengths in domain
+    CS%plx = 4.0 * PI / G%len_lon
+    CS%pmz = PI * CS%mode / CS%H0
+    CS%lambda = CS%pmz * CS%F_0 / CS%N0
+    CS%omega = CS%F_0 * CS%plx / CS%lambda
   endif
 
   sina = sin(CS%coast_angle)
@@ -205,11 +230,18 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, h, Time)
         y = - (x1 - CS%coast_offset) * sina + y1 * cosa
         if (CS%mode == 0) then
           cff = sqrt(G%g_Earth * 0.5 * (G%bathyT(i+1,j) + G%bathyT(i,j)))
-          val2 = fac * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1)) * y / cff)
-          segment%eta(I,j) = val2 * cos(omega * time_sec)
+          val2 = fac * exp(- CS%F_0 * y / cff)
+          segment%eta(I,j) = val2 * cos(CS%omega * time_sec)
           segment%normal_vel_bt(I,j) = val1 * cff * cosa /         &
                  (0.5 * (G%bathyT(i+1,j) + G%bathyT(i,j))) * val2
         else
+          segment%eta(I,j) = 0.0
+          segment%normal_vel_bt(I,j) = 0.0
+          do k=1,nz
+            segment%nudged_normal_vel(I,j,k) = fac * CS%lambda / CS%F_0 * &
+                   exp(- CS%lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * &
+                   cos(CS%omega * time_sec)
+          enddo
         endif
       enddo ; enddo
     else
@@ -223,10 +255,16 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, h, Time)
         if (CS%mode == 0) then
           cff = sqrt(G%g_Earth * 0.5 * (G%bathyT(i,j+1) + G%bathyT(i,j)))
           val2 = fac * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * y / cff)
-          segment%eta(I,j) = val2 * cos(omega * time_sec)
+          segment%eta(I,j) = val2 * cos(CS%omega * time_sec)
           segment%normal_vel_bt(I,j) = val1 * cff * sina /       &
                  (0.5 * (G%bathyT(i+1,j) + G%bathyT(i,j))) * val2
         else
+          segment%eta(i,J) = 0.0
+          segment%normal_vel_bt(i,J) = 0.0
+          do k=1,nz
+            segment%nudged_normal_vel(i,J,k) = fac * CS%lambda / CS%F_0 * &
+                   exp(- CS%lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * cosa
+          enddo
         endif
       enddo ; enddo
     endif