apply_maxvel_clamp Subroutine

private pure subroutine apply_maxvel_clamp(ms, maxvel)

Truncate face velocities to |u| ≤ maxvel. MOM6’s MAXVEL analogue. Called once at the end of each outer step (after the RK2 average so the clipped state is what gets carried into the next stage). No-op when maxvel <= 0.

Clip preserves the sign and the velocity direction — it’s a per-component clip, not a magnitude clip. Matches MOM6’s if (abs(u) > maxvel) u = sign(maxvel, u) form so the WBC jet that pushes through 6 m/s gets clipped to ±6 m/s without introducing direction-reversal artifacts.

Not conservative: clipping a face velocity from 10 m/s to 6 m/s loses momentum. That’s acceptable as a safety net — when the clamp is active, momentum conservation is already broken by whatever generated the runaway velocity. In production runs the clamp should fire rarely or never.

Arguments

Type IntentOptional Attributes Name
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: maxvel

Called by

proc~~apply_maxvel_clamp~~CalledByGraph proc~apply_maxvel_clamp apply_maxvel_clamp proc~apply_velocity_truncation apply_velocity_truncation proc~apply_velocity_truncation->proc~apply_maxvel_clamp proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~apply_velocity_truncation proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~apply_velocity_truncation proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~engine_step->proc~ocean_dyn_step_split proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step proc~driver_run driver_run proc~driver_run->proc~driver_run_ocean

Variables

Type Visibility Attributes Name Initial
integer, private :: i
integer, private :: j
integer, private :: k
integer, private :: nx_face
integer, private :: nx_vface
integer, private :: ny_face
integer, private :: ny_uface
integer, private :: nz

Source Code

   pure subroutine apply_maxvel_clamp(ms, maxvel)
      !! Truncate face velocities to `|u| ≤ maxvel`.  MOM6's MAXVEL
      !! analogue.  Called once at the end of each outer step (after
      !! the RK2 average so the clipped state is what gets carried
      !! into the next stage).  No-op when `maxvel <= 0`.
      !!
      !! Clip preserves the sign and the velocity direction — it's
      !! a per-component clip, not a magnitude clip.  Matches MOM6's
      !! `if (abs(u) > maxvel) u = sign(maxvel, u)` form so the WBC
      !! jet that pushes through 6 m/s gets clipped to ±6 m/s without
      !! introducing direction-reversal artifacts.
      !!
      !! Not conservative: clipping a face velocity from 10 m/s to
      !! 6 m/s loses momentum.  That's acceptable as a safety net —
      !! when the clamp is active, momentum conservation is already
      !! broken by whatever generated the runaway velocity.  In
      !! production runs the clamp should fire rarely or never.
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: maxvel
      integer :: i, j, k, nx_face, ny_uface, nx_vface, ny_face, nz

      if (maxvel <= 0.0_wp) return

      nx_face = size(ms%u_face_x_layer, 1)
      ny_uface = size(ms%u_face_x_layer, 2)
      nx_vface = size(ms%v_face_y_layer, 1)
      ny_face = size(ms%v_face_y_layer, 2)
      nz = ms%nz_ml

      ! NaN-safe (pdc 8c2fd674 + 8675b0d4): the ieee_is_finite guard both
      ! (a) leaves a non-finite value untouched (NEVER laundering NaN →
      ! ±maxvel) and (b) blocks nvfortran -fast from lowering the
      ! `if(u>hi)…else if(u<lo)` pair to a NaN-blind min/max clamp on GPU.
      ! Non-finite velocities are zeroed upstream by the truncation's
      ! NaN-catch; this is the second line of defence.
      do concurrent(k=1:nz, j=1:ny_uface, i=1:nx_face)
         if (ieee_is_finite(ms%u_face_x_layer(i, j, k))) then
            if (ms%u_face_x_layer(i, j, k) > maxvel) then
               ms%u_face_x_layer(i, j, k) = maxvel
            else if (ms%u_face_x_layer(i, j, k) < -maxvel) then
               ms%u_face_x_layer(i, j, k) = -maxvel
            end if
         end if
      end do
      do concurrent(k=1:nz, j=1:ny_face, i=1:nx_vface)
         if (ieee_is_finite(ms%v_face_y_layer(i, j, k))) then
            if (ms%v_face_y_layer(i, j, k) > maxvel) then
               ms%v_face_y_layer(i, j, k) = maxvel
            else if (ms%v_face_y_layer(i, j, k) < -maxvel) then
               ms%v_face_y_layer(i, j, k) = -maxvel
            end if
         end if
      end do
   end subroutine apply_maxvel_clamp