ppm_limit_pos Subroutine

public pure subroutine ppm_limit_pos(h_centre, h_left, h_right, h_min)

Positivity-preserving limiter on the PPM reconstruction. Mirrors MOM6’s PPM_limit_pos: when the parabolic fit predicts a minimum interior to the cell that dips below h_min, shrink h_left / h_right toward h_centre so the minimum sits at exactly h_min. Pure scalar form per cell; runs after ppm_cell_limiter so the monotonic-limited reconstruction is the input.

Algorithm: curv = 3·(h_L + h_R − 2·h_in) ! +ve ⇒ interior min if curv > 0 and |dh| < curv: ! min inside cell if h_in ≤ h_min: flatten (h_L = h_R = h_in) elif 12·curv·(h_in − h_min) < curv² + 3·dh²: loc_scale = 12·curv·(h_in − h_min) / (curv² + 3·dh²) ∈ (0,1) h_L = h_in + loc_scale·(h_L − h_in) h_R = h_in + loc_scale·(h_R − h_in)

h_min = 0 ⇒ pure positivity (parabola can’t go negative inside the cell). Larger h_min ⇒ harder floor; matches MOM6’s GV%Angstrom_H for the reduced-gravity setup.

No-op when curv ≤ 0 (maximum interior, or linear / monotone profile) or |dh| ≥ curv (minimum outside the cell, edges already control).

Consumer: rdb_ice_transport (sea-ice PR 4b) — SIS2 runs this unconditionally on the PD continuity scheme (SIS_continuity.F90:1585, PPM_limit_pos), so the ice-mass reconstruction always applies it (not optional there, unlike the ocean’s use_ppm_limit_pos knob).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_centre
real(kind=wp), intent(inout) :: h_left
real(kind=wp), intent(inout) :: h_right
real(kind=wp), intent(in) :: h_min

Called by

proc~~ppm_limit_pos~~CalledByGraph proc~ppm_limit_pos ppm_limit_pos proc~continuity_compute_fluxes continuity_compute_fluxes proc~continuity_compute_fluxes->proc~ppm_limit_pos proc~continuity_compute_fluxes_barotropic continuity_compute_fluxes_barotropic proc~continuity_compute_fluxes_barotropic->proc~ppm_limit_pos proc~continuity_meridional_flux continuity_meridional_flux proc~continuity_meridional_flux->proc~ppm_limit_pos proc~continuity_zonal_flux continuity_zonal_flux proc~continuity_zonal_flux->proc~ppm_limit_pos proc~ice_cat_flux_x_impl ice_cat_flux_x_impl proc~ice_cat_flux_x_impl->proc~ppm_limit_pos proc~ice_cat_flux_y_impl ice_cat_flux_y_impl proc~ice_cat_flux_y_impl->proc~ppm_limit_pos proc~continuity_step_split continuity_step_split proc~continuity_step_split->proc~continuity_meridional_flux proc~continuity_step_split->proc~continuity_zonal_flux proc~continuity_tracer_step_split continuity_tracer_step_split proc~continuity_tracer_step_split->proc~continuity_meridional_flux proc~continuity_tracer_step_split->proc~continuity_zonal_flux proc~ice_pass_x ice_pass_x proc~ice_pass_x->proc~ice_cat_flux_x_impl proc~ice_pass_y ice_pass_y proc~ice_pass_y->proc~ice_cat_flux_y_impl proc~ocean_dyn_step_barotropic ocean_dyn_step_barotropic proc~ocean_dyn_step_barotropic->proc~continuity_compute_fluxes_barotropic proc~ice_transport_step ice_transport_step proc~ice_transport_step->proc~ice_pass_x proc~ice_transport_step->proc~ice_pass_y proc~run_continuity_chain run_continuity_chain proc~run_continuity_chain->proc~continuity_tracer_step_split proc~run_stage run_stage proc~run_stage->proc~continuity_tracer_step_split proc~engine_step_ice engine_step_ice proc~engine_step_ice->proc~ice_transport_step proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~run_stage proc~run_stage_split run_stage_split proc~run_stage_split->proc~run_continuity_chain proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step_ice proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step_ice

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: curv
real(kind=wp), private :: dh
real(kind=wp), private :: loc_scale

Source Code

   pure subroutine ppm_limit_pos(h_centre, h_left, h_right, h_min)
      !$acc routine seq
      !! Positivity-preserving limiter on the PPM reconstruction.
      !! Mirrors MOM6's `PPM_limit_pos`:
      !! when the parabolic fit predicts a minimum interior to the
      !! cell that dips below `h_min`, shrink h_left / h_right toward
      !! h_centre so the minimum sits at exactly `h_min`.  Pure
      !! scalar form per cell; runs after `ppm_cell_limiter` so the
      !! monotonic-limited reconstruction is the input.
      !!
      !! Algorithm:
      !!   curv = 3·(h_L + h_R − 2·h_in)         ! +ve ⇒ interior min
      !!   if curv > 0 and |dh| < curv:           ! min inside cell
      !!     if h_in ≤ h_min: flatten (h_L = h_R = h_in)
      !!     elif 12·curv·(h_in − h_min) < curv² + 3·dh²:
      !!        loc_scale = 12·curv·(h_in − h_min) / (curv² + 3·dh²) ∈ (0,1)
      !!        h_L = h_in + loc_scale·(h_L − h_in)
      !!        h_R = h_in + loc_scale·(h_R − h_in)
      !!
      !! `h_min = 0` ⇒ pure positivity (parabola can't go negative
      !! inside the cell).  Larger `h_min` ⇒ harder floor; matches
      !! MOM6's `GV%Angstrom_H` for the reduced-gravity setup.
      !!
      !! No-op when curv ≤ 0 (maximum interior, or linear / monotone
      !! profile) or |dh| ≥ curv (minimum outside the cell, edges
      !! already control).
      !!
      !! Consumer: `rdb_ice_transport` (sea-ice PR 4b) — SIS2 runs this
      !! unconditionally on the PD continuity scheme
      !! (`SIS_continuity.F90:1585`, `PPM_limit_pos`), so the ice-mass
      !! reconstruction always applies it (not optional there, unlike
      !! the ocean's `use_ppm_limit_pos` knob).
      real(wp), intent(in) :: h_centre, h_min
      real(wp), intent(inout) :: h_left, h_right
      real(wp) :: curv, dh, loc_scale

      curv = 3.0_wp*((h_left + h_right) - 2.0_wp*h_centre)
      if (curv > 0.0_wp) then
         dh = h_right - h_left
         if (abs(dh) < curv) then
            if (h_centre <= h_min) then
               h_left = h_centre
               h_right = h_centre
            else if (12.0_wp*curv*(h_centre - h_min) < (curv*curv + 3.0_wp*dh*dh)) then
               loc_scale = 12.0_wp*curv*(h_centre - h_min)/(curv*curv + 3.0_wp*dh*dh)
               h_left = h_centre + loc_scale*(h_left - h_centre)
               h_right = h_centre + loc_scale*(h_right - h_centre)
            end if
         end if
      end if
   end subroutine ppm_limit_pos