pure function flux3D_LinearEuler3D_t(this,s,dsdx) result(flux)
class(LinearEuler3D_t),intent(in) :: this
real(prec),intent(in) :: s(1:this%nvar)
real(prec),intent(in) :: dsdx(1:this%nvar,1:3)
real(prec) :: flux(1:this%nvar,1:3)
flux(1,1) = s(4)/s(6) ! x-velocity, x flux; p/rho0
flux(1,2) = 0.0_prec ! x-velocity, y flux; 0
flux(1,3) = 0.0_prec ! x-velocity, z flux; 0
flux(2,1) = 0.0_prec ! y-velocity, x flux; 0
flux(2,2) = s(4)/s(6) ! y-velocity, y flux; p/rho0
flux(2,3) = 0.0_prec ! y-velocity, z flux; 0
flux(3,1) = 0.0_prec ! z-velocity, x flux; 0
flux(3,2) = 0.0_prec ! z-velocity, y flux; 0
flux(3,3) = s(4)/s(6) ! z-velocity, z flux; p/rho0
flux(4,1) = s(5)*s(5)*s(6)*s(1) ! pressure, x flux : rho0*c^2*u
flux(4,2) = s(5)*s(5)*s(6)*s(2) ! pressure, y flux : rho0*c^2*v
flux(4,3) = s(5)*s(5)*s(6)*s(3) ! pressure, z flux : rho0*c^2*w
flux(5,1) = 0.0_prec ! sound speed, x flux; 0 (c held fixed in time)
flux(5,2) = 0.0_prec ! sound speed, y flux; 0
flux(5,3) = 0.0_prec ! sound speed, z flux; 0
flux(6,1) = 0.0_prec ! background density, x flux; 0 (rho0 held fixed in time)
flux(6,2) = 0.0_prec ! background density, y flux; 0
flux(6,3) = 0.0_prec ! background density, z flux; 0
if(.false.) flux(1,1) = flux(1,1)+dsdx(1,1) ! suppress unused-dummy-argument warning
endfunction flux3D_LinearEuler3D_t