riemannflux3D_LinearEuler3D_t Function

  pure function riemannflux3D_LinearEuler3D_t(this,sL,sR,dsdx,nhat) result(flux)
    !! Impedance-matched (characteristic/Godunov) Riemann flux, identical in
    !! form to the 2-D model's. The interface states are resolved with the
    !! per-side acoustic impedances Z = rho0*c (each side using its own
    !! background density rho0 and sound speed c), so material interfaces in a
    !! heterogeneous field are handled with the physically correct
    !! transmission/reflection:
    !!
    !!   un* = (ZL*unL + ZR*unR + (pL - pR)) / (ZL + ZR)
    !!   p*  = (ZR*pL + ZL*pR + ZL*ZR*(unL - unR)) / (ZL + ZR)
    !!
    !! The reconstructed momentum/pressure fluxes use the arithmetic
    !! averages rho0_avg and c2_avg at the face (see the 2-D model for the
    !! rationale). The sound speed and background density variables carry zero
    !! flux.
    class(LinearEuler3D_t),intent(in) :: this
    real(prec),intent(in) :: sL(1:this%nvar)
    real(prec),intent(in) :: sR(1:this%nvar)
    real(prec),intent(in) :: dsdx(1:this%nvar,1:3)
    real(prec),intent(in) :: nhat(1:3)
    real(prec) :: flux(1:this%nvar)
    ! Local
    real(prec) :: rho0L,rho0R,rho0_avg,cL,cR,ZL,ZR,unL,unR,pL,pR,un_star,p_star,c2_avg

    rho0L = sL(6)
    rho0R = sR(6)
    rho0_avg = 0.5_prec*(rho0L+rho0R)
    cL = sL(5)
    cR = sR(5)
    ZL = rho0L*cL
    ZR = rho0R*cR

    unL = sL(1)*nhat(1)+sL(2)*nhat(2)+sL(3)*nhat(3)
    unR = sR(1)*nhat(1)+sR(2)*nhat(2)+sR(3)*nhat(3)
    pL = sL(4)
    pR = sR(4)

    un_star = (ZL*unL+ZR*unR+(pL-pR))/(ZL+ZR)
    p_star = (ZR*pL+ZL*pR+ZL*ZR*(unL-unR))/(ZL+ZR)
    c2_avg = 0.5_prec*(cL*cL+cR*cR)

    flux(1) = p_star*nhat(1)/rho0_avg
    flux(2) = p_star*nhat(2)/rho0_avg
    flux(3) = p_star*nhat(3)/rho0_avg
    flux(4) = rho0_avg*c2_avg*un_star
    flux(5) = 0.0_prec
    flux(6) = 0.0_prec
    if(.false.) flux(1) = flux(1)+dsdx(1,1) ! suppress unused-dummy-argument warning

  endfunction riemannflux3D_LinearEuler3D_t