weno9_face_swept Function

public pure function weno9_face_swept(qm4, qm3, qm2, qm1, q0, qp1, qp2, qp3, qp4, sigma) result(face)

WENO9-Z swept-average face value (u > 0, upwind cell = q0).

Five quartic candidates with optimal weights d = (1/126, 10/63, 10/21, 20/63, 5/126). Z-weights: tau9 = |beta0 - beta4|. Smoothness indicators use the simplified 3-term form.

Reference: Balsara & Shu (2000).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: qm4

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: qm3

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: qm2

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: qm1

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: q0

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: qp1

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: qp2

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: qp3

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: qp4

Cell averages at i-4 through i+4.
real(kind=wp), intent(in) :: sigma

CFL of upwind cell: |u_eff|*dt/dx.

Return Value real(kind=wp)


Called by

proc~~weno9_face_swept~~CalledByGraph proc~weno9_face_swept weno9_face_swept proc~weno_face_conc_x weno_face_conc_x proc~weno_face_conc_x->proc~weno9_face_swept proc~weno_face_conc_y weno_face_conc_y proc~weno_face_conc_y->proc~weno9_face_swept proc~drain_swept_flux_x_weno drain_swept_flux_x_weno proc~drain_swept_flux_x_weno->proc~weno_face_conc_x proc~drain_swept_flux_y_weno drain_swept_flux_y_weno proc~drain_swept_flux_y_weno->proc~weno_face_conc_y proc~continuity_tracer_drain continuity_tracer_drain proc~continuity_tracer_drain->proc~drain_swept_flux_x_weno proc~continuity_tracer_drain->proc~drain_swept_flux_y_weno proc~ocean_dyn_flush_tracer_window ocean_dyn_flush_tracer_window proc~ocean_dyn_flush_tracer_window->proc~continuity_tracer_drain proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~continuity_tracer_drain proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~continuity_tracer_drain proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~ocean_dyn_flush_tracer_window proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~engine_step->proc~ocean_dyn_step_split proc~rdb_ocean_set_tracer rdb_ocean_set_tracer proc~rdb_ocean_set_tracer->proc~ocean_dyn_flush_tracer_window

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: A0
real(kind=wp), private :: A1
real(kind=wp), private :: A2
real(kind=wp), private :: A3
real(kind=wp), private :: A4
real(kind=wp), private :: B0
real(kind=wp), private :: B1
real(kind=wp), private :: B2
real(kind=wp), private :: B3
real(kind=wp), private :: B4
real(kind=wp), private :: C0
real(kind=wp), private :: C1
real(kind=wp), private :: C2
real(kind=wp), private :: C3
real(kind=wp), private :: C4
real(kind=wp), private :: D0
real(kind=wp), private :: D1
real(kind=wp), private :: D2
real(kind=wp), private :: D3
real(kind=wp), private :: D4
real(kind=wp), private :: E0
real(kind=wp), private :: E1
real(kind=wp), private :: E2
real(kind=wp), private :: E3
real(kind=wp), private :: E4
real(kind=wp), private :: alpha0
real(kind=wp), private :: alpha1
real(kind=wp), private :: alpha2
real(kind=wp), private :: alpha3
real(kind=wp), private :: alpha4
real(kind=wp), private :: alpha_sum
real(kind=wp), private :: beta0
real(kind=wp), private :: beta1
real(kind=wp), private :: beta2
real(kind=wp), private :: beta3
real(kind=wp), private :: beta4
real(kind=wp), private, parameter :: c13_12 = 13.0_wp/12.0_wp
real(kind=wp), private, parameter :: c1_4 = 0.25_wp
real(kind=wp), private, parameter :: c1_80 = 1.0_wp/80.0_wp
real(kind=wp), private, parameter :: d0_w = 1.0_wp/126.0_wp
real(kind=wp), private, parameter :: d1_w = 10.0_wp/63.0_wp
real(kind=wp), private, parameter :: d2_w = 10.0_wp/21.0_wp
real(kind=wp), private, parameter :: d3_w = 20.0_wp/63.0_wp
real(kind=wp), private, parameter :: d4_w = 5.0_wp/126.0_wp
real(kind=wp), private, parameter :: eps = 1.0e-36_wp
real(kind=wp), private :: face0
real(kind=wp), private :: face1
real(kind=wp), private :: face2
real(kind=wp), private :: face3
real(kind=wp), private :: face4
real(kind=wp), private :: omega0
real(kind=wp), private :: omega1
real(kind=wp), private :: omega2
real(kind=wp), private :: omega3
real(kind=wp), private :: omega4
real(kind=wp), private :: p_half
real(kind=wp), private :: pp_half
real(kind=wp), private :: ppp_half
real(kind=wp), private :: pppp_half
real(kind=wp), private :: ppppp_half
real(kind=wp), private :: tau9

Source Code

   pure function weno9_face_swept(qm4, qm3, qm2, qm1, q0, qp1, qp2, qp3, qp4, sigma) &
      result(face)
      !! WENO9-Z swept-average face value (u > 0, upwind cell = q0).
      !!
      !! Five quartic candidates with optimal weights
      !! d = (1/126, 10/63, 10/21, 20/63, 5/126).
      !! Z-weights: tau9 = |beta0 - beta4|.
      !! Smoothness indicators use the simplified 3-term form.
      !!
      !! Reference: Balsara & Shu (2000).
      !$acc routine seq
      real(wp), intent(in) :: qm4, qm3, qm2, qm1, q0, qp1, qp2, qp3, qp4
         !! Cell averages at i-4 through i+4.
      real(wp), intent(in) :: sigma
         !! CFL of upwind cell: |u_eff|*dt/dx.
      real(wp) :: face

      ! Optimal weights (BS00): d = (1/126, 10/63, 10/21, 20/63, 5/126)
      real(wp), parameter :: d0_w = 1.0_wp/126.0_wp
      real(wp), parameter :: d1_w = 10.0_wp/63.0_wp
      real(wp), parameter :: d2_w = 10.0_wp/21.0_wp
      real(wp), parameter :: d3_w = 20.0_wp/63.0_wp
      real(wp), parameter :: d4_w = 5.0_wp/126.0_wp
      real(wp), parameter :: eps = 1.0e-36_wp
      real(wp), parameter :: c13_12 = 13.0_wp/12.0_wp
      real(wp), parameter :: c1_4 = 0.25_wp
      real(wp), parameter :: c1_80 = 1.0_wp/80.0_wp

      ! Quartic polynomial A + B*xi + C*xi^2 + D*xi^3 + E*xi^4 per stencil
      real(wp) :: A0, B0, C0, D0, E0
      real(wp) :: A1, B1, C1, D1, E1
      real(wp) :: A2, B2, C2, D2, E2
      real(wp) :: A3, B3, C3, D3, E3
      real(wp) :: A4, B4, C4, D4, E4

      ! Simplified 3-term beta (centred on stencil centre c):
      !   beta = (13/12)*(q_{c-1}-2*q_c+q_{c+1})^2 + (1/4)*(q_{c-1}-q_{c+1})^2
      !          + (1/80)*(q_{c-2}-4*q_{c-1}+6*q_c-4*q_{c+1}+q_{c+2})^2
      real(wp) :: beta0, beta1, beta2, beta3, beta4, tau9

      ! Z-weights
      real(wp) :: alpha0, alpha1, alpha2, alpha3, alpha4, alpha_sum
      real(wp) :: omega0, omega1, omega2, omega3, omega4

      ! Swept-average face values per candidate
      real(wp) :: face0, face1, face2, face3, face4

      ! Quartic swept-average helpers
      real(wp) :: p_half, pp_half, ppp_half, pppp_half, ppppp_half

      ! -----------------------------------------------------------------------
      ! Quartic coefficients (exact rational from BS00)
      ! -----------------------------------------------------------------------

      ! Stencil 0: offsets (-4,-3,-2,-1,0) → qm4,qm3,qm2,qm1,q0
      A0 = (-71.0_wp/1920.0_wp)*qm4 + (91.0_wp/480.0_wp)*qm3 &
           + (-373.0_wp/960.0_wp)*qm2 + (57.0_wp/160.0_wp)*qm1 &
           + (563.0_wp/640.0_wp)*q0
      B0 = (3.0_wp/16.0_wp)*qm4 + (-25.0_wp/24.0_wp)*qm3 &
           + (5.0_wp/2.0_wp)*qm2 + (-29.0_wp/8.0_wp)*qm1 &
           + (95.0_wp/48.0_wp)*q0
      C0 = (7.0_wp/16.0_wp)*qm4 + (-9.0_wp/4.0_wp)*qm3 &
           + (37.0_wp/8.0_wp)*qm2 + (-17.0_wp/4.0_wp)*qm1 &
           + (23.0_wp/16.0_wp)*q0
      D0 = (1.0_wp/4.0_wp)*qm4 + (-7.0_wp/6.0_wp)*qm3 &
           + 2.0_wp*qm2 + (-3.0_wp/2.0_wp)*qm1 &
           + (5.0_wp/12.0_wp)*q0
      E0 = (1.0_wp/24.0_wp)*qm4 + (-1.0_wp/6.0_wp)*qm3 &
           + (1.0_wp/4.0_wp)*qm2 + (-1.0_wp/6.0_wp)*qm1 &
           + (1.0_wp/24.0_wp)*q0

      ! Stencil 1: offsets (-3,-2,-1,0,+1) → qm3,qm2,qm1,q0,qp1
      A1 = (3.0_wp/640.0_wp)*qm3 + (-3.0_wp/160.0_wp)*qm2 &
           + (-13.0_wp/960.0_wp)*qm1 + (511.0_wp/480.0_wp)*q0 &
           + (-71.0_wp/1920.0_wp)*qp1
      B1 = (-5.0_wp/48.0_wp)*qm3 + (5.0_wp/8.0_wp)*qm2 &
           + (-7.0_wp/4.0_wp)*qm1 + (25.0_wp/24.0_wp)*q0 &
           + (3.0_wp/16.0_wp)*qp1
      C1 = (-1.0_wp/16.0_wp)*qm3 + (1.0_wp/4.0_wp)*qm2 &
           + (1.0_wp/8.0_wp)*qm1 + (-3.0_wp/4.0_wp)*q0 &
           + (7.0_wp/16.0_wp)*qp1
      D1 = (1.0_wp/12.0_wp)*qm3 + (-1.0_wp/2.0_wp)*qm2 &
           + 1.0_wp*qm1 + (-5.0_wp/6.0_wp)*q0 &
           + (1.0_wp/4.0_wp)*qp1
      E1 = (1.0_wp/24.0_wp)*qm3 + (-1.0_wp/6.0_wp)*qm2 &
           + (1.0_wp/4.0_wp)*qm1 + (-1.0_wp/6.0_wp)*q0 &
           + (1.0_wp/24.0_wp)*qp1

      ! Stencil 2: offsets (-2,-1,0,+1,+2) → qm2,qm1,q0,qp1,qp2
      A2 = (3.0_wp/640.0_wp)*qm2 + (-29.0_wp/480.0_wp)*qm1 &
           + (1067.0_wp/960.0_wp)*q0 + (-29.0_wp/480.0_wp)*qp1 &
           + (3.0_wp/640.0_wp)*qp2
      B2 = (5.0_wp/48.0_wp)*qm2 + (-17.0_wp/24.0_wp)*qm1 &
           + 0.0_wp*q0 + (17.0_wp/24.0_wp)*qp1 &
           + (-5.0_wp/48.0_wp)*qp2
      C2 = (-1.0_wp/16.0_wp)*qm2 + (3.0_wp/4.0_wp)*qm1 &
           + (-11.0_wp/8.0_wp)*q0 + (3.0_wp/4.0_wp)*qp1 &
           + (-1.0_wp/16.0_wp)*qp2
      D2 = (-1.0_wp/12.0_wp)*qm2 + (1.0_wp/6.0_wp)*qm1 &
           + 0.0_wp*q0 + (-1.0_wp/6.0_wp)*qp1 &
           + (1.0_wp/12.0_wp)*qp2
      E2 = (1.0_wp/24.0_wp)*qm2 + (-1.0_wp/6.0_wp)*qm1 &
           + (1.0_wp/4.0_wp)*q0 + (-1.0_wp/6.0_wp)*qp1 &
           + (1.0_wp/24.0_wp)*qp2

      ! Stencil 3: offsets (-1,0,+1,+2,+3) → qm1,q0,qp1,qp2,qp3
      A3 = (-71.0_wp/1920.0_wp)*qm1 + (511.0_wp/480.0_wp)*q0 &
           + (-13.0_wp/960.0_wp)*qp1 + (-3.0_wp/160.0_wp)*qp2 &
           + (3.0_wp/640.0_wp)*qp3
      B3 = (-3.0_wp/16.0_wp)*qm1 + (-25.0_wp/24.0_wp)*q0 &
           + (7.0_wp/4.0_wp)*qp1 + (-5.0_wp/8.0_wp)*qp2 &
           + (5.0_wp/48.0_wp)*qp3
      C3 = (7.0_wp/16.0_wp)*qm1 + (-3.0_wp/4.0_wp)*q0 &
           + (1.0_wp/8.0_wp)*qp1 + (1.0_wp/4.0_wp)*qp2 &
           + (-1.0_wp/16.0_wp)*qp3
      D3 = (-1.0_wp/4.0_wp)*qm1 + (5.0_wp/6.0_wp)*q0 &
           + (-1.0_wp)*qp1 + (1.0_wp/2.0_wp)*qp2 &
           + (-1.0_wp/12.0_wp)*qp3
      E3 = (1.0_wp/24.0_wp)*qm1 + (-1.0_wp/6.0_wp)*q0 &
           + (1.0_wp/4.0_wp)*qp1 + (-1.0_wp/6.0_wp)*qp2 &
           + (1.0_wp/24.0_wp)*qp3

      ! Stencil 4: offsets (0,+1,+2,+3,+4) → q0,qp1,qp2,qp3,qp4
      A4 = (563.0_wp/640.0_wp)*q0 + (57.0_wp/160.0_wp)*qp1 &
           + (-373.0_wp/960.0_wp)*qp2 + (91.0_wp/480.0_wp)*qp3 &
           + (-71.0_wp/1920.0_wp)*qp4
      B4 = (-95.0_wp/48.0_wp)*q0 + (29.0_wp/8.0_wp)*qp1 &
           + (-5.0_wp/2.0_wp)*qp2 + (25.0_wp/24.0_wp)*qp3 &
           + (-3.0_wp/16.0_wp)*qp4
      C4 = (23.0_wp/16.0_wp)*q0 + (-17.0_wp/4.0_wp)*qp1 &
           + (37.0_wp/8.0_wp)*qp2 + (-9.0_wp/4.0_wp)*qp3 &
           + (7.0_wp/16.0_wp)*qp4
      D4 = (-5.0_wp/12.0_wp)*q0 + (3.0_wp/2.0_wp)*qp1 &
           + (-2.0_wp)*qp2 + (7.0_wp/6.0_wp)*qp3 &
           + (-1.0_wp/4.0_wp)*qp4
      E4 = (1.0_wp/24.0_wp)*q0 + (-1.0_wp/6.0_wp)*qp1 &
           + (1.0_wp/4.0_wp)*qp2 + (-1.0_wp/6.0_wp)*qp3 &
           + (1.0_wp/24.0_wp)*qp4

      ! -----------------------------------------------------------------------
      ! Simplified 3-term smoothness indicators (centred on each stencil)
      ! Stencil r centre c:
      !   beta_r = (13/12)*(q_{c-1}-2*q_c+q_{c+1})^2
      !          + (1/4)  *(q_{c-1}-q_{c+1})^2
      !          + (1/80) *(q_{c-2}-4*q_{c-1}+6*q_c-4*q_{c+1}+q_{c+2})^2
      ! -----------------------------------------------------------------------
      ! Stencil 0: c = qm2 (index i-2), c±1 = qm3/qm1, c±2 = qm4/q0
      beta0 = c13_12*(qm3 - 2.0_wp*qm2 + qm1)**2 &
              + c1_4*(qm3 - qm1)**2 &
              + c1_80*(qm4 - 4.0_wp*qm3 + 6.0_wp*qm2 - 4.0_wp*qm1 + q0)**2

      ! Stencil 1: c = qm1 (index i-1), c±1 = qm2/q0, c±2 = qm3/qp1
      beta1 = c13_12*(qm2 - 2.0_wp*qm1 + q0)**2 &
              + c1_4*(qm2 - q0)**2 &
              + c1_80*(qm3 - 4.0_wp*qm2 + 6.0_wp*qm1 - 4.0_wp*q0 + qp1)**2

      ! Stencil 2: c = q0 (index i), c±1 = qm1/qp1, c±2 = qm2/qp2
      beta2 = c13_12*(qm1 - 2.0_wp*q0 + qp1)**2 &
              + c1_4*(qm1 - qp1)**2 &
              + c1_80*(qm2 - 4.0_wp*qm1 + 6.0_wp*q0 - 4.0_wp*qp1 + qp2)**2

      ! Stencil 3: c = qp1 (index i+1), c±1 = q0/qp2, c±2 = qm1/qp3
      beta3 = c13_12*(q0 - 2.0_wp*qp1 + qp2)**2 &
              + c1_4*(q0 - qp2)**2 &
              + c1_80*(qm1 - 4.0_wp*q0 + 6.0_wp*qp1 - 4.0_wp*qp2 + qp3)**2

      ! Stencil 4: c = qp2 (index i+2), c±1 = qp1/qp3, c±2 = q0/qp4
      beta4 = c13_12*(qp1 - 2.0_wp*qp2 + qp3)**2 &
              + c1_4*(qp1 - qp3)**2 &
              + c1_80*(q0 - 4.0_wp*qp1 + 6.0_wp*qp2 - 4.0_wp*qp3 + qp4)**2

      tau9 = abs(beta0 - beta4)
      alpha0 = d0_w*(1.0_wp + (tau9/(eps + beta0))**2)
      alpha1 = d1_w*(1.0_wp + (tau9/(eps + beta1))**2)
      alpha2 = d2_w*(1.0_wp + (tau9/(eps + beta2))**2)
      alpha3 = d3_w*(1.0_wp + (tau9/(eps + beta3))**2)
      alpha4 = d4_w*(1.0_wp + (tau9/(eps + beta4))**2)
      alpha_sum = alpha0 + alpha1 + alpha2 + alpha3 + alpha4
      omega0 = alpha0/alpha_sum
      omega1 = alpha1/alpha_sum
      omega2 = alpha2/alpha_sum
      omega3 = alpha3/alpha_sum
      omega4 = alpha4/alpha_sum

      ! Quartic swept-average over [1/2 - sigma, 1/2]:
      !   q_face = p(1/2) - sigma/2 * p'(1/2) + sigma^2/6 * p''(1/2)
      !            - sigma^3/24 * p'''(1/2) + sigma^4/120 * p''''(1/2)
      ! where for p = A + B*xi + C*xi^2 + D*xi^3 + E*xi^4:
      !   p(1/2)    = A + B/2 + C/4 + D/8 + E/16
      !   p'(1/2)   = B + C   + 3D/4 + E/2
      !   p''(1/2)  = 2C + 3D + 3E
      !   p'''(1/2) = 6D + 12E
      !   p''''(1/2)= 24E

      ! Stencil 0
      p_half = A0 + 0.5_wp*B0 + 0.25_wp*C0 + 0.125_wp*D0 + 0.0625_wp*E0
      pp_half = B0 + C0 + 0.75_wp*D0 + 0.5_wp*E0
      ppp_half = 2.0_wp*C0 + 3.0_wp*D0 + 3.0_wp*E0
      pppp_half = 6.0_wp*D0 + 12.0_wp*E0
      ppppp_half = 24.0_wp*E0
      face0 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      ! Stencil 1
      p_half = A1 + 0.5_wp*B1 + 0.25_wp*C1 + 0.125_wp*D1 + 0.0625_wp*E1
      pp_half = B1 + C1 + 0.75_wp*D1 + 0.5_wp*E1
      ppp_half = 2.0_wp*C1 + 3.0_wp*D1 + 3.0_wp*E1
      pppp_half = 6.0_wp*D1 + 12.0_wp*E1
      ppppp_half = 24.0_wp*E1
      face1 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      ! Stencil 2
      p_half = A2 + 0.5_wp*B2 + 0.25_wp*C2 + 0.125_wp*D2 + 0.0625_wp*E2
      pp_half = B2 + C2 + 0.75_wp*D2 + 0.5_wp*E2
      ppp_half = 2.0_wp*C2 + 3.0_wp*D2 + 3.0_wp*E2
      pppp_half = 6.0_wp*D2 + 12.0_wp*E2
      ppppp_half = 24.0_wp*E2
      face2 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      ! Stencil 3
      p_half = A3 + 0.5_wp*B3 + 0.25_wp*C3 + 0.125_wp*D3 + 0.0625_wp*E3
      pp_half = B3 + C3 + 0.75_wp*D3 + 0.5_wp*E3
      ppp_half = 2.0_wp*C3 + 3.0_wp*D3 + 3.0_wp*E3
      pppp_half = 6.0_wp*D3 + 12.0_wp*E3
      ppppp_half = 24.0_wp*E3
      face3 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      ! Stencil 4
      p_half = A4 + 0.5_wp*B4 + 0.25_wp*C4 + 0.125_wp*D4 + 0.0625_wp*E4
      pp_half = B4 + C4 + 0.75_wp*D4 + 0.5_wp*E4
      ppp_half = 2.0_wp*C4 + 3.0_wp*D4 + 3.0_wp*E4
      pppp_half = 6.0_wp*D4 + 12.0_wp*E4
      ppppp_half = 24.0_wp*E4
      face4 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      face = omega0*face0 + omega1*face1 + omega2*face2 + omega3*face3 + omega4*face4
   end function weno9_face_swept