riemannflux3D_LinearEuler3D_t Function

  pure function riemannflux3D_LinearEuler3D_t(this,sL,sR,dsdx,nhat) result(flux)
    !! Uses a local lax-friedrich's upwind flux
    !! The max eigenvalue is taken as the sound speed
    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) :: fL(1:this%nvar)
    real(prec) :: fR(1:this%nvar)
    real(prec) :: u,v,w,p,c,rho0

    u = sL(2)
    v = sL(3)
    w = sL(4)
    p = sL(5)
    rho0 = this%rho0
    c = this%c
    fL(1) = rho0*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! density
    fL(2) = p*nhat(1)/rho0 ! u
    fL(3) = p*nhat(2)/rho0 ! v
    fL(4) = p*nhat(3)/rho0 ! w
    fL(5) = rho0*c*c*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! pressure

    u = sR(2)
    v = sR(3)
    w = sR(4)
    p = sR(5)
    fR(1) = rho0*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! density
    fR(2) = p*nhat(1)/rho0 ! u
    fR(3) = p*nhat(2)/rho0 ! v'
    fR(4) = p*nhat(3)/rho0 ! w
    fR(5) = rho0*c*c*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! pressure

    flux(1:5) = 0.5_prec*(fL(1:5)+fR(1:5))+c*(sL(1:5)-sR(1:5))

  endfunction riemannflux3D_LinearEuler3D_t