Better Kelvin wave results in layer mode

Author: kshedstrom

src/core/MOM_dynamics_split_RK2.F90

@@ -798,7 +798,6 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
   call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
 
   if (associated(CS%OBC)) then
-    call update_segment_thickness_reservoirs(G, GV, uhtr, vhtr, hp, CS%OBC)
 
     if (CS%debug) &
       call uvchksum("Pre OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
@@ -1065,7 +1064,6 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
   endif
 
   if (associated(CS%OBC)) then
-    call update_segment_thickness_reservoirs(G, GV, uhtr, vhtr, h, CS%OBC)
     !### I suspect that there is a bug here when u_inst is compared with a previous value of u_av
     ! to estimate the dominant outward group velocity, but a fix is not available yet.
     call radiation_open_bdry_conds(CS%OBC, u_inst, u_old_rad_OBC, v_inst, v_old_rad_OBC, G, GV, US, dt)
@@ -1090,6 +1088,10 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
     enddo ; enddo
   enddo
 
+  if (associated(CS%OBC)) then
+    call update_segment_thickness_reservoirs(G, GV, uhtr, vhtr, h, CS%OBC)
+  endif
+
   if (CS%store_CAu) then
     ! Calculate a predictor-step estimate of the Coriolis and momentum advection terms
     ! for use in the next time step, possibly after it has been vertically remapped.

src/core/MOM_dynamics_split_RK2b.F90

@@ -59,6 +59,8 @@ module MOM_dynamics_split_RK2b
 use MOM_open_boundary,         only : ocean_OBC_type, radiation_open_bdry_conds
 use MOM_open_boundary,         only : open_boundary_zero_normal_flow, open_boundary_query
 use MOM_open_boundary,         only : open_boundary_test_extern_h, update_OBC_ramp
+use MOM_open_boundary,         only : copy_thickness_reservoirs
+use MOM_open_boundary,         only : update_segment_thickness_reservoirs
 use MOM_PressureForce,         only : PressureForce, PressureForce_CS
 use MOM_PressureForce,         only : PressureForce_init
 use MOM_set_visc,              only : set_viscous_ML, set_visc_CS
@@ -73,7 +75,7 @@ module MOM_dynamics_split_RK2b
 use MOM_vert_friction,         only : updateCFLtruncationValue, vertFPmix
 use MOM_verticalGrid,          only : verticalGrid_type, get_thickness_units
 use MOM_verticalGrid,          only : get_flux_units, get_tr_flux_units
-use MOM_wave_interface,        only: wave_parameters_CS, Stokes_PGF
+use MOM_wave_interface,        only : wave_parameters_CS, Stokes_PGF
 
 implicit none ; private
 
@@ -543,6 +545,9 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
   call disable_averaging(CS%diag)
   if (showCallTree) call callTree_wayPoint("done with PressureForce (step_MOM_dyn_split_RK2b)")
 
+  if (associated(CS%OBC)) then ; if (CS%OBC%update_OBC) then
+    call update_OBC_data(CS%OBC, G, GV, US, tv, h, CS%update_OBC_CSp, Time_local)
+  endif ; endif
   if (associated(CS%OBC) .and. CS%debug_OBC) &
     call open_boundary_zero_normal_flow(CS%OBC, G, GV, CS%PFu, CS%PFv)
 
@@ -667,6 +672,9 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
   call cpu_clock_end(id_clock_mom_update)
   call do_group_pass(CS%pass_uv_inst, G%Domain, clock=id_clock_pass)
 
+  if (associated(CS%OBC)) &
+    call copy_thickness_reservoirs(CS%OBC, G, GV)
+
 ! u_accel_bt = layer accelerations due to barotropic solver
   call cpu_clock_begin(id_clock_continuity)
   call continuity(u_inst, v_inst, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, &
@@ -798,6 +806,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
   call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
 
   if (associated(CS%OBC)) then
+
     if (CS%debug) &
       call uvchksum("Pre OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
 
@@ -1051,6 +1060,10 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
     enddo ; enddo
   enddo
 
+  if (associated(CS%OBC)) then
+    call update_segment_thickness_reservoirs(G, GV, uhtr, vhtr, h, CS%OBC)
+  endif
+
   !  if (CS%fpmix) then
   !    if (CS%id_uold > 0) call post_data(CS%id_uold, uold, CS%diag)
   !    if (CS%id_vold > 0) call post_data(CS%id_vold, vold, CS%diag)

src/core/MOM_open_boundary.F90

@@ -6002,9 +6002,9 @@ subroutine update_segment_thickness_reservoirs(G, GV, uhr, vhr, h, OBC)
           ! When InvLscale_in is 0 and inflow, only nudged data is applied to reservoirs
           a_out = b_out * max(0.0, sign(1.0, idir*uhr(I,j,k)))
           a_in  = b_in  * min(0.0, sign(1.0, idir*uhr(I,j,k)))
-          u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_out*resrv_lfac_out / &
+          u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Th_InvLscale_out*resrv_lfac_out / &
                     ((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j)))
-          u_L_in  = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_in*resrv_lfac_in  / &
+          u_L_in  = min(0.0, (idir*uhr(I,j,k))*segment%Th_InvLscale_in*resrv_lfac_in  / &
                     ((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j)))
           fac1 = (1.0 - (a_out - a_in)) + ((u_L_out + a_out) - (u_L_in + a_in))
           segment%h_Reg%h_res(I,j,k) = (1.0/fac1) * &
@@ -6039,9 +6039,9 @@ subroutine update_segment_thickness_reservoirs(G, GV, uhr, vhr, h, OBC)
         if (allocated(segment%h_Reg%h_res)) then ; do k=1,nz
           a_out = b_out * max(0.0, sign(1.0, jdir*vhr(i,J,k)))
           a_in  = b_in  * min(0.0, sign(1.0, jdir*vhr(i,J,k)))
-          v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_out*resrv_lfac_out / &
+          v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Th_InvLscale_out*resrv_lfac_out / &
                     ((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J)))
-          v_L_in  = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_in*resrv_lfac_in  / &
+          v_L_in  = min(0.0, (jdir*vhr(i,J,k))*segment%Th_InvLscale_in*resrv_lfac_in  / &
                     ((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J)))
           fac1 = (1.0 - (a_out - a_in)) + ((v_L_out + a_out) - (v_L_in + a_in))
           segment%h_Reg%h_res(i,J,k) = (1.0/fac1) * &

src/user/Kelvin_initialization.F90

@@ -283,8 +283,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
     segment => OBC%segment(n)
     if (.not. segment%on_pe) cycle
 
-   unrot_dir = segment%direction
-   if (turns /= 0) unrot_dir = rotate_OBC_segment_direction(segment%direction, -turns)
+    unrot_dir = segment%direction
+    if (turns /= 0) unrot_dir = rotate_OBC_segment_direction(segment%direction, -turns)
 
     ! Apply values to the inflow end only.
     if ((unrot_dir == OBC_DIRECTION_E) .or. (unrot_dir == OBC_DIRECTION_N)) cycle
@@ -346,9 +346,19 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
             enddo
           endif
         else
-          ! Baroclinic, not rotated yet (and apparently not working as intended yet).
+          ! Baroclinic, not rotated yet
           segment%SSH(I,j) = 0.0
           segment%normal_vel_bt(I,j) = 0.0
+          ! Use inside bathymetry
+          if (segment%direction == OBC_DIRECTION_W) then
+            depth_tot_vel = depth_tot(i+1,j)
+          elseif (segment%direction == OBC_DIRECTION_S) then
+            depth_tot_vel = depth_tot(i,j+1)
+          elseif (segment%direction == OBC_DIRECTION_E) then
+            depth_tot_vel = depth_tot(i,j)
+          elseif (segment%direction == OBC_DIRECTION_N) then
+            depth_tot_vel = depth_tot(i,j)
+          endif
           ! 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_wt
           ! In addition, there should be a specification of the interface-height anomalies at the
@@ -368,6 +378,16 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
                    exp(-lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * cos_wt
             enddo
           endif
+          if (associated(segment%h_Reg)) then
+            if (allocated(segment%h_Reg%h)) then
+              do k=1,nz
+                segment%h_Reg%h(I,j,k) = depth_tot_vel / nz + &
+                              ((CS%mode * PI) * CS%inflow_amp / (N0 * nz)) * &
+                                          cos(((PI * k) * CS%mode) / nz) * &
+                                          exp(-lambda * y) * cos_wt
+              enddo
+            endif
+          endif
         endif
       enddo ; enddo
     endif