!********************************************************************************
!* A collection of three-dimensional Euler numerical fluxes, Version 1 (2011),
!*
!*        written by Dr. Katate Masatsuka (info[at]cfdbooks.com),
!*
!* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com).
!*
!*------------------------
!* List of Flux Functions:
!*
!*   Roe
!*   Rotated-RHLL
!*
!*------------------------
!*
!* These F90 routines were written and made available for an educational purpose.
!* Detailed descripstion of each numerical flux can be found in the original paper,
!* or in popular textbooks, or in the (will-be-available) second volume of
!* "I do like CFD". 
!*
!* Note: These subroutines have been prepared for an educational purpose.
!*       It is not at all efficient. Think about how you can optimize it.
!*
!* Note: These subroutines have not been tested in a 3D code.
!*       Please let me know if you find bugs. I'll greatly appreciate it and
!*       fix the bugs.
!*
!* 05-06-11: Minor bugs fixed.
!*
!* This file may be updated in future. (to add more flux functions.)
!********************************************************************************


!********************************************************************************
!* -- 3D Roe's Flux Function with an entropy fix --
!*
!* This subroutine computes the Roe flux for the Euler equations
!* in the direction, njk=[nx,ny].
!*
!* P. L. Roe, Approximate Riemann Solvers, Parameter Vectors and Difference
!* Schemes, Journal of Computational Physics, 43, pp. 357-372.
!*
!* Conservative form of the Euler equations:
!*
!*     dU/dt + dF/dx + dG/dy = 0
!*
!* The physical flux is given by
!*
!*     Fn = F*nx + G*ny = | rho*qn          |
!*                        | rho*qn*u + p*nx |
!*                        | rho*qn*v + p*ny |
!*                        | rho*qn*w + p*nz |
!*                        | rho*qn*H        |    (qn = u*nx + v*ny + w*nz)
!*
!* The Roe flux is defined as
!*
!*   Numerical flux = 1/2 [ Fn(UR) + Fn(UL) - |An|dU ], 
!*
!*  where
!*
!*    An = dFn/dU,  |An| = R|Lambda|L, dU = UR - UL.
!*
!* The dissipation term, |An|dU, is actually computed as
!*
!*     sum_{k=1,5} |lambda_k| * (LdU)_k * r_k,
!*
!* where lambda_k is the k-th eigenvalue, (LdU)_k is the k-th wave strength,
!* and r_k is the k-th right-eigenvector.
!*
!* ------------------------------------------------------------------------------
!*  Input: primL(1:5) =  Left state (rhoL, uL, vL, wR, pL)
!*         primR(1:5) = Right state (rhoR, uR, vR, wR, pR)
!*             njk(2) = Unit normal vector (nx, ny), pointing from Left to Right.
!*
!*           njk
!*  Face normal ^   o Right data point
!*              |  .
!*              | .
!*              |. 
!*       -------x-------- Face
!*             .                 Left and right states are
!*            .                   1. Values at data points for 1st-order accuracy
!*           .                    2. Extrapolated values at the face midpoint 'x'
!*          o Left data point        for 2nd-order accuracy.
!*
!*
!* Output:  flux(1:5) = The Roe flux with an entropy fix
!* ------------------------------------------------------------------------------
!*
!* Note: This subroutine has been prepared for an educational purpose.
!*       It is not at all efficient. Think about how you can optimize it.
!*
!* Note: This subroutine has not been tested in a 3D code.
!*       Please let me know if you find bugs. I'll greatly appreciate it and
!*       fix the bugs.
!*
!* Katate Masatsuka, March 2011. http://www.cfdbooks.com
!********************************************************************************
 subroutine inviscid_roe(primL, primR, njk,  num_flux)

 implicit none
 integer , parameter :: p2 = selected_real_kind(15)

!Input
 real(p2), intent( in) :: primL(5), primR(5)!Input: primitive variables
 real(p2), intent( in) :: njk(3)            !Input: face normal vector

!Output
 real(p2), intent(out) :: num_flux(5)       !Output (numerical flux)

!Some constants
 real(p2) ::  zero = 0.0_p2
 real(p2) ::   one = 1.0_p2
 real(p2) ::   two = 2.0_p2
 real(p2) ::  half = 0.5_p2
 real(p2) :: fifth = 0.2_p2

!Local variables
 real(p2) :: nx, ny, nz                   ! Normal vector
 real(p2) :: mx, my, mz                   ! Orthogonal tangent vector
 real(p2) :: lx, ly, lz                   ! Another orthogonal tangent vector
 real(p2) :: abs_n_cross_l                ! Magnitude of n x l
 real(p2) :: uL, uR, vL, vR, wL, wR       ! Velocity components.
 real(p2) :: rhoL, rhoR, pL, pR           ! Primitive variables.
 real(p2) :: qnL, qnR, qmL, qmR, qlL, qlR ! Normal and tangent velocities
 real(p2) :: aL, aR, HL, HR               ! Speeds of sound.
 real(p2) :: RT,rho,u,v,w,H,a,qn, ql, qm  ! Roe-averages
 real(p2) :: drho,dqn,dql,dqm,dp,LdU(5)   ! Wave strengths
 real(p2) :: ws(5), R(5,5)                ! Wave speeds and right-eigenvectors
 real(p2) :: dws(5)                       !
 real(p2) :: fL(5), fR(5), diss(5)        ! Fluxes ad dissipation term
 integer  :: i, j

 real(p2) :: gamma = 1.4_p2               ! Ratio of specific heats

  nx = njk(1)
  ny = njk(2)
  nz = njk(3)

!Tangent vectors, m and l, such that the three vectors, n, m, and l
!are mutually orthogonal.

! Let's choose the one in the xy-plane.

  lx = -ny
  ly =  nx
  lz =  zero

! The other one is chosen as a vector orthogonal to both n and l
! defined by (mx,my,mz) = n x l / |n x l|

  mx = ny*lz - nz*ly
  my = nz*lx - nx*lz
  mz = nx*ly - ny*lx

  abs_n_cross_l = sqrt(mx**2 + my**2 + mz**2)
  mx = mx / abs_n_cross_l
  my = my / abs_n_cross_l
  mz = mz / abs_n_cross_l

!(Do you like such ambiguous tangent vectors? Actually, Roe flux can
! be implemented without any tangent vector. See "I do like CFD, VOL.1"
! or the subroutine "inviscid_rotated_rhll" for details. Hence, the resulting
! flux is independent of the choice of these tangent vectors.)

!Primitive and other variables.

!  Left state
    rhoL = primL(1)
      uL = primL(2)
      vL = primL(3)
      wL = primL(4)
     qnL = uL*nx + vL*ny + wL*nz
     qlL = uL*lx + vL*ly + wL*lz
     qmL = uL*mx + vL*my + wL*mz
      pL = primL(5)
      aL = sqrt(gamma*pL/rhoL)
      HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL)
!  Right state
    rhoR = primR(1)
      uR = primR(2)
      vR = primR(3)
      wR = primR(4)
     qnR = uR*nx + vR*ny + wR*nz
     qlR = uR*lx + vR*ly + wR*lz
     qmR = uR*mx + vR*my + wR*mz
      pR = primR(5)
      aR = sqrt(gamma*pR/rhoR)
      HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR)

!First compute the Roe Averages

    RT = sqrt(rhoR/rhoL)
   rho = RT*rhoL
     u = (uL + RT*uR)/(one + RT)
     v = (vL + RT*vR)/(one + RT)
     w = (wL + RT*wR)/(one + RT)
     H = (HL + RT*HR)/(one + RT)
     a = sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) )
    qn = u*nx + v*ny + w*nz
    ql = u*lx + v*ly + w*lz
    qm = u*mx + v*my + w*mz

!Wave Strengths

   drho = rhoR - rhoL 
     dp =   pR - pL
    dqn =  qnR - qnL
    dql =  qlR - qlL
    dqm =  qmR - qmL

  LdU(1) = (dp - rho*a*dqn )/(two*a*a)
  LdU(2) =  drho - dp/(a*a)
  LdU(3) = (dp + rho*a*dqn )/(two*a*a)
  LdU(4) = rho*dql
  LdU(5) = rho*dqm

!Absolute values of the wave Speeds

  ws(1) = abs(qn-a)
  ws(2) = abs(qn)
  ws(3) = abs(qn+a)
  ws(4) = abs(qn)
  ws(5) = abs(qn)

!Harten's Entropy Fix JCP(1983), 49, pp357-393: only for the nonlinear fields.

  dws(1) = fifth
   if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) )
  dws(3) = fifth
   if ( ws(3) < dws(3) ) ws(3) = half * ( ws(3)*ws(3)/dws(3)+dws(3) )

!Right Eigenvectors

  R(1,1) = one    
  R(2,1) = u - a*nx
  R(3,1) = v - a*ny
  R(4,1) = w - a*nz
  R(5,1) = H - a*qn

  R(1,2) = one
  R(2,2) = u
  R(3,2) = v 
  R(4,2) = w
  R(5,2) = half*(u*u + v*v + w*w)

  R(1,3) = one
  R(2,3) = u + a*nx
  R(3,3) = v + a*ny
  R(4,3) = w + a*nz
  R(5,3) = H + a*qn

  R(1,4) = zero
  R(2,4) = lx
  R(3,4) = ly
  R(4,4) = lz
  R(5,4) = ql

  R(1,5) = zero
  R(2,5) = mx
  R(3,5) = my
  R(4,5) = mz
  R(5,5) = qm

!Dissipation Term: |An|(UR-UL)

  diss = zero
  do i=1,5
   do j=1,5
    diss(i) = diss(i) + ws(j)*LdU(j)*R(i,j)
   end do
  end do

!Compute the physical flux: fL = Fn(UL) and fR = Fn(UR)

  fL(1) = rhoL*qnL
  fL(2) = rhoL*qnL * uL + pL*nx
  fL(3) = rhoL*qnL * vL + pL*ny
  fL(4) = rhoL*qnL * wL + pL*nz
  fL(5) = rhoL*qnL * HL

  fR(1) = rhoR*qnR
  fR(2) = rhoR*qnR * uR + pR*nx
  fR(3) = rhoR*qnR * vR + pR*ny
  fR(4) = rhoR*qnR * wR + pR*nz
  fR(5) = rhoR*qnR * HR

! This is the numerical flux: Roe flux = 1/2 *[  Fn(UL)+Fn(UR) - |An|(UR-UL) ]

  num_flux = half * (fL + fR - diss)

!Normal max wave speed in the normal direction.
!  wsn = abs(qn) + a

 end subroutine inviscid_roe
!--------------------------------------------------------------------------------


!********************************************************************************
!* -- 3D Rotated-Roe-HLL (Rotated-RHLL) Flux Function ---
!*
!* This subroutine computes the Rotated-RHLL flux for the Euler equations
!* in the direction, njk=[nx,ny].
!*
!* H. Nishikawa and K. Kitamura, Very Simple, Carbuncle-Free, Boundary-Layer
!* Resolving, Rotated-Hybrid Riemann Solvers,
!* Journal of Computational Physics, 227, pp. 2560-2581, 2008.
!*
!* Very robust for nonlinear instability (carbuncle).
!* Recommended for high-speed flows involving strong shocks.
!*
!*
!* Conservative form of the Euler equations:
!*
!*     dU/dt + dF/dx + dG/dy = 0
!*
!* This subroutine computes the numerical flux for the flux in the direction,
!* njk=[nx,ny]:
!*
!*     Fn = F*nx + G*ny = | rho*qn          |
!*                        | rho*qn*u + p*nx |
!*                        | rho*qn*v + p*ny |
!*                        | rho*qn*w + p*nz |
!*                        | rho*qn*H        |    (qn = u*nx + v*ny + w*nz)
!*
!* The Rotated-RHLL flux is defined as
!*
!*   Numerical flux = alpha1*HLL(n1) + alpha2*Roe(n2),  if |dq| > eps
!*                    Roe(n)                         ,  if |dq| < eps
!*
!* where n1 is taken as the velocity difference vector (normal to shock or
!* tangent to shear waves), n2 is a vector perpendicular to n1,
!* alpha1 = n*n1, alpha2 = n*n2. That is, we decompose the normal vector, n,
!* in the two directions, n1 and n2, and apply the HLL in n1 and the Roe in n2.
!*
!* Note: HLL is introduced only near discontinuities and only by a fraction, alpha1.
!* Note: For full convergence to machine zero for staedy state computations,
!*       the vectors, n1 and n2, may have to be frozen.
!* Note: Here, the Roe flux is computed without tangent vectors.
!*
!* ------------------------------------------------------------------------------
!*  Input: primL(1:5) =  left state (rhoL, uL, vL, wR, pL)
!*         primR(1:5) = right state (rhoR, uR, vR, wR, pR)
!*             njk(2) = unit normal vector (nx, ny), pointing from Left to Right
!*
!*           njk
!*  Face normal ^   o Right data point
!*              |  .
!*              | .
!*              |. 
!*       -------x-------- Face
!*             .                 Left and right states are
!*            .                   1. Values at data points for 1st-order accuracy
!*           .                    2. Extrapolated values at the face midpoint 'x'
!*          o Left data point        for 2nd-order accuracy.
!*
!* Output:  flux(1:5) = the Roe flux
!* ------------------------------------------------------------------------------
!*
!* Note: This subroutine has been prepared for an educational purpose.
!*       It is not at all efficient. Think about how you can optimize it.
!*
!* Note: This subroutine has not been tested in a 3D code.
!*       Please let me know if you find bugs. I'll greatly appreciate it and
!*       fix the bugs.
!* 
!* Katate Masatsuka, March 2011. http://www.cfdbooks.com
!********************************************************************************
 subroutine inviscid_rotated_rhll(primL, primR, njk,  num_flux)

 implicit none
 integer , parameter :: p2 = selected_real_kind(15) !Double precision

!Input
 real(p2), intent( in) :: primL(5), primR(5)!Input: primitive variables
 real(p2), intent( in) :: njk(3)            !Input: face normal vector

!Output
 real(p2), intent(out) :: num_flux(5)       !Output (numerical flux)

!Some constants
 real(p2) ::  zero = 0.0_p2
 real(p2) ::   one = 1.0_p2
 real(p2) ::   two = 2.0_p2
 real(p2) ::  half = 0.5_p2
 real(p2) :: fifth = 0.2_p2

!Local variables
 real(p2) :: nx, ny, nz                     !Face normal vector
 real(p2) :: uL, uR, vL, vR, wL, wR         !Velocity components.
 real(p2) :: rhoL, rhoR, pL, pR             !Primitive variables.
 real(p2) :: qnL, qnR                       !Normal velocity
 real(p2) :: aL, aR, HL, HR                 !Speeds of sound.
 real(p2) :: RT,rho,u,v,w,H,a,qn            !Roe-averages
 real(p2) :: drho,dqn,dp,LdU(4)             !Wave strengths
 real(p2) :: du, dv, dw                     !Differences
 real(p2) :: eig(4)                         !Eigenvalues
 real(p2) :: ws(4), R(5,4)                  !Absolute Wave speeds and right-eigenvectors
 real(p2) :: dws(4)                         !
 real(p2) :: fL(5), fR(5), diss(5)          !Fluxes ad dissipation term

 real(p2) :: gamma = 1.4_p2 !Ratio of specific heats

 real(p2) :: SRp,SLm                        !Wave speeds for the HLL part
 real(p2) :: nx1, ny1, nz1                  !Vector along which HLL is applied
 real(p2) :: nx2, ny2, nz2                  !Vector along which Roe is applied
 real(p2) :: alpha1, alpha2                 !Projections of the new normals
 real(p2) :: abs_dq                         !Magnitude of the velocity difference
 real(p2) :: temp

! Normal vector

  nx = njk(1)
  ny = njk(2)
  nz = njk(3)

!Primitive and other variables.

!  Left state
    rhoL = primL(1)
      uL = primL(2)
      vL = primL(3)
      wL = primL(4)
     qnL = uL*nx + vL*ny + wL*nz
      pL = primL(5)
      aL = sqrt(gamma*pL/rhoL)
      HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL)

!  Right state
    rhoR = primR(1)
      uR = primR(2)
      vR = primR(3)
      wR = primR(4)
     qnR = uR*nx + vR*ny + wR*nz
      pR = primR(5)
      aR = sqrt(gamma*pR/rhoR)
      HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR)

!Compute the physical flux: fL = Fn(UL) and fR = Fn(UR)

  fL(1) = rhoL*qnL
  fL(2) = rhoL*qnL * uL + pL*nx
  fL(3) = rhoL*qnL * vL + pL*ny
  fL(4) = rhoL*qnL * wL + pL*nz
  fL(5) = rhoL*qnL * HL

  fR(1) = rhoR*qnR
  fR(2) = rhoR*qnR * uR + pR*nx
  fR(3) = rhoR*qnR * vR + pR*ny
  fR(4) = rhoR*qnR * wR + pR*nz
  fR(5) = rhoR*qnR * HR

!--------------------------------------------------------------------------------
!Define n1 and n2, and compute alpha1 and alpha2: (4.2) in the original paper.
! Note: n1 and n2 may need to be frozen at some point during 
!       a steady calculation to fully make it converge. For time-accurate 
!       calculation, it is not necessary.
! Note: For a boundary face, you may want toset (nx2,ny2)=(nx,ny),
!       (nx1,ny1)=(-ny,nx), i.e., use the Roe flux.

    abs_dq = sqrt( (uR-uL)**2 + (vR-vL)**2 + (wR-wL)**2 )

!n1 = Velocity difference vector: normal to shock or tangent to shear
  if ( abs_dq > 1.0e-12_p2) then

       nx1 = (uR-uL)/abs_dq
       ny1 = (vR-vL)/abs_dq
       nz1 = (wR-wL)/abs_dq

!n1 = Face tangent vector if abs_dq is too small.
!     Note: There are infinitely many choices for the tangent vector.
!           The best choice may be discovered in future.
  else
       nx1 = -ny
       ny1 =  nx
       nz1 =  zero
  endif

    alpha1 = nx*nx1 + ny*ny1 + nz*nz1

! Make alpha1 always positive.
      temp = sign(one,alpha1)
       nx1 = temp * nx1
       ny1 = temp * ny1
       nz1 = temp * nz1
    alpha1 = temp * alpha1

!n2 = direction perpendicular to n1.
!     Note: There are infinitely many choices for this vector.
!           The best choice may be discovered in future.
       nx2 = -ny1
       ny2 =  nx1
       nz2 =  zero
    alpha2 = nx*nx2 + ny*ny2 + nz*nz2

! Make alpha2 always positive.
      temp = sign(one,alpha2)
       nx2 = temp * nx2
       ny2 = temp * ny2
       nz2 = temp * nz2
    alpha2 = temp * alpha2

!--------------------------------------------------------------------------------
!Now we are going to compute the Roe flux with n2 as the normal with modified 
!wave speeds (5.12). NOTE: the Roe flux here is computed without tangent vectors.
!See "I do like CFD, VOL.1" for details: page 57, Equation (3.6.31).

!First compute the Roe Averages

      RT = sqrt(rhoR/rhoL)
     rho = RT*rhoL
       u = (uL + RT*uR)/(one + RT)
       v = (vL + RT*vR)/(one + RT)
       w = (wL + RT*wR)/(one + RT)
       H = (HL + RT*HR)/(one + RT)
       a = sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) )

!----------------------------------------------------
!Compute the wave speed estimates for the HLL part.
!Note: HLL is actually applied to n1, but this is
!      all we need to incorporate HLL. See JCP paper.

      qn = u*nx1 + v*ny1 + w*nz1
     SLm = min( zero, qn - a, qnL - aL )
     SRp = max( zero, qn + a, qnR + aR )
!----------------------------------------------------

!Wave Strengths

      qn = u*nx2 + v*ny2 + w*nz2

    drho = rhoR - rhoL 
      dp =   pR - pL
     dqn =  qnR - qnL

  LdU(1) = (dp - rho*a*dqn )/(two*a*a)
  LdU(2) =  drho - dp/(a*a)
  LdU(3) = (dp + rho*a*dqn )/(two*a*a)
  LdU(4) = rho 

!Wave Speed (Eigenvalues)

  eig(1) = qn-a
  eig(2) = qn
  eig(3) = qn+a
  eig(4) = qn

!Absolute values of the wave speeds (Eigenvalues)

   ws(1) = abs(qn-a)
   ws(2) = abs(qn)
   ws(3) = abs(qn+a)
   ws(4) = abs(qn)

!Harten's Entropy Fix JCP(1983), 49, pp357-393: only for the nonlinear fields.

  dws(1) = fifth
   if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) )
  dws(3) = fifth
   if ( ws(3) < dws(3) ) ws(3) = half * ( ws(3)*ws(3)/dws(3)+dws(3) )

!Combine the wave speeds for Rotated-RHLL: Eq.(5.12) in the original JCP paper.

      ws = alpha2*ws - (alpha1*two*SRp*SLm + alpha2*(SRp+SLm)*eig)/(SRp-SLm)

!Right Eigenvectors:
!Note: Two shear wave components are combined into one, so that tangent vectors
!      are not required. And that's why there are only 4 vectors here.

  R(1,1) = one    
  R(2,1) = u - a*nx2
  R(3,1) = v - a*ny2
  R(4,1) = w - a*nz2
  R(5,1) = H - a*qn

  R(1,2) = one
  R(2,2) = u
  R(3,2) = v 
  R(4,2) = w
  R(5,2) = half*(u*u + v*v + w*w)

  R(1,3) = one
  R(2,3) = u + a*nx2
  R(3,3) = v + a*ny2
  R(4,3) = w + a*nz2
  R(5,3) = H + a*qn

  du = uR - uL
  dv = vR - vL
  dw = wR - wL

  R(1,4) = zero
  R(2,4) = du - dqn*nx2
  R(3,4) = dv - dqn*ny2
  R(4,4) = dw - dqn*nz2
  R(5,4) = u*du + v*dv + w*dw - qn*dqn

!Dissipation Term: Roe dissipation with the modified wave speeds.

    diss = ws(1)*LdU(1)*R(:,1) + ws(2)*LdU(2)*R(:,2) &
         + ws(3)*LdU(3)*R(:,3) + ws(4)*LdU(4)*R(:,4)

!Compute the Rotated-RHLL flux.

  num_flux = (SRp*fL - SLm*fR)/(SRp-SLm) - half*diss

!Normal max wave speed in the normal direction.
!  wsn = abs(qn) + a

 end subroutine inviscid_rotated_rhll
!--------------------------------------------------------------------------------