@@ -1736,12 +1736,11 @@ contains
17361736 ! Input parameters
17371737 integer , intent (in ) :: wave_speeds
17381738 real (wp), intent (in ) :: rho_L, rho_R
1739- real (wp), dimension (:), intent (in ) :: vel_L, vel_R
1739+ real (wp), dimension (:), intent (in ) :: vel_L, vel_R, tau_e_L, tau_e_R
17401740 real (wp), intent (in ) :: pres_L, pres_R, c_L, c_R
17411741 real (wp), intent (in ) :: gamma_L, gamma_R, pi_inf_L, pi_inf_R
17421742 real (wp), intent (in ) :: rho_avg, c_avg, c_fast_L, c_fast_R
17431743 real (wp), intent (in ) :: G_L, G_R
1744- real (wp), dimension (:), intent (in ) :: tau_e_L, tau_e_R
17451744
17461745 ! Local variables
17471746 real (wp) :: pres_SL, pres_SR, Ms_L, Ms_R
@@ -1751,64 +1750,53 @@ contains
17511750
17521751 if (wave_speeds == 1 ) then
17531752 if (mhd) then
1754- s_L = min (vel_L - c_fast_L, vel_R - c_fast_R)
1755- s_R = max (vel_R + c_fast_R, vel_L + c_fast_L)
1753+ s_L = min (vel_L(dir_idx( 1 )) - c_fast_L, vel_R(dir_idx( 1 )) - c_fast_R)
1754+ s_R = max (vel_R(dir_idx( 1 )) + c_fast_R, vel_L(dir_idx( 1 )) + c_fast_L)
17561755 elseif (hypoelasticity .or. elasticity) then
1757- s_L = min (vel_L - sqrt (c_L* c_L + &
1758- (((4._wp * G_L)/ 3._wp ) + &
1759- tau_e_L)/ rho_L) &
1760- , vel_R - sqrt (c_R* c_R + &
1761- (((4._wp * G_R)/ 3._wp ) + &
1762- tau_e_R)/ rho_R))
1763- s_R = max (vel_R + sqrt (c_R* c_R + &
1764- (((4._wp * G_R)/ 3._wp ) + &
1765- tau_e_R)/ rho_R) &
1766- , vel_L + sqrt (c_L* c_L + &
1767- (((4._wp * G_L)/ 3._wp ) + &
1768- tau_e_L)/ rho_L))
1756+ s_L = min (vel_L(dir_idx(1 )) - sqrt (c_L* c_L + (((4._wp * G_L)/ 3._wp ) + &
1757+ tau_e_L(dir_idx_tau(1 )))/ rho_L) &
1758+ , vel_R(dir_idx(1 )) - sqrt (c_R* c_R + (((4._wp * G_R)/ 3._wp ) + &
1759+ tau_e_R(dir_idx_tau(1 )))/ rho_R))
1760+ s_R = max (vel_R(dir_idx(1 )) + sqrt (c_R* c_R + (((4._wp * G_R)/ 3._wp ) + &
1761+ tau_e_R(dir_idx_tau(1 )))/ rho_R) &
1762+ , vel_L(dir_idx(1 )) + sqrt (c_L* c_L + (((4._wp * G_L)/ 3._wp ) + &
1763+ tau_e_L(dir_idx_tau(1 )))/ rho_L))
17691764 else if (hyperelasticity) then
1770- s_L = min (vel_L - sqrt (c_L* c_L + (4._wp * G_L/ 3._wp )/ rho_L) &
1771- , vel_R - sqrt (c_R* c_R + (4._wp * G_R/ 3._wp )/ rho_R))
1772- s_R = max (vel_R + sqrt (c_R* c_R + (4._wp * G_R/ 3._wp )/ rho_R) &
1773- , vel_L + sqrt (c_L* c_L + (4._wp * G_L/ 3._wp )/ rho_L))
1765+ s_L = min (vel_L(dir_idx( 1 )) - sqrt (c_L* c_L + (4._wp * G_L/ 3._wp )/ rho_L) &
1766+ , vel_R(dir_idx( 1 )) - sqrt (c_R* c_R + (4._wp * G_R/ 3._wp )/ rho_R))
1767+ s_R = max (vel_R(dir_idx( 1 )) + sqrt (c_R* c_R + (4._wp * G_R/ 3._wp )/ rho_R) &
1768+ , vel_L(dir_idx( 1 )) + sqrt (c_L* c_L + (4._wp * G_L/ 3._wp )/ rho_L))
17741769 else
1775- s_L = min (vel_L - c_L, vel_R - c_R)
1776- s_R = max (vel_R + c_R, vel_L + c_L)
1770+ s_L = min (vel_L(dir_idx( 1 )) - c_L, vel_R(dir_idx( 1 )) - c_R)
1771+ s_R = max (vel_R(dir_idx( 1 )) + c_R, vel_L(dir_idx( 1 )) + c_L)
17771772 end if
1778- s_S = (pres_R - pres_L + rho_L* vel_L* &
1779- (s_L - vel_L) - &
1780- rho_R* vel_R* &
1781- (s_R - vel_R)) &
1782- / (rho_L* (s_L - vel_L) - &
1783- rho_R* (s_R - vel_R))
1773+ s_S = (pres_R - pres_L + rho_L* vel_L(dir_idx(1 ))* &
1774+ (s_L - vel_L(dir_idx(1 ))) - rho_R* vel_R(dir_idx(1 ))* (s_R - vel_R(dir_idx(1 )))) &
1775+ / (rho_L* (s_L - vel_L(dir_idx(1 ))) - rho_R* (s_R - vel_R(dir_idx(1 ))))
17841776 elseif (wave_speeds == 2 ) then
1785- pres_SL = 5e-1_wp * (pres_L + pres_R + rho_avg* c_avg* &
1786- (vel_L - &
1787- vel_R))
1777+ pres_SL = 5e-1_wp * (pres_L + pres_R + rho_avg* c_avg* (vel_L(dir_idx(1 )) - vel_R(dir_idx(1 ))))
17881778 pres_SR = pres_SL
17891779 Ms_L = max (1._wp , sqrt (1._wp + ((5e-1_wp + gamma_L)/ (1._wp + gamma_L))* &
17901780 (pres_SL/ pres_L - 1._wp )* pres_L/ &
17911781 ((pres_L + pi_inf_L/ (1._wp + gamma_L)))))
17921782 Ms_R = max (1._wp , sqrt (1._wp + ((5e-1_wp + gamma_R)/ (1._wp + gamma_R))* &
17931783 (pres_SR/ pres_R - 1._wp )* pres_R/ &
17941784 ((pres_R + pi_inf_R/ (1._wp + gamma_R)))))
1795- s_L = vel_L - c_L* Ms_L
1796- s_R = vel_R + c_R* Ms_R
1797- s_S = 5e-1_wp * ((vel_L + vel_R) + &
1798- (pres_L - pres_R)/ &
1799- (rho_avg* c_avg))
1785+ s_L = vel_L(dir_idx(1 )) - c_L* Ms_L
1786+ s_R = vel_R(dir_idx(1 )) + c_R* Ms_R
1787+ s_S = 5e-1_wp * ((vel_L(dir_idx(1 )) + vel_R(dir_idx(1 ))) + (pres_L - pres_R)/ (rho_avg* c_avg))
18001788 end if
18011789
18021790 ! follows Einfeldt et al.
18031791 ! s_M/ P = min/ max (0 .,s_L/ R)
1804- s_M = min (0._wp , s_L); s_P = max (0._wp , s_R)
1792+ s_M = min (0._wp , s_L), s_P = max (0._wp , s_R)
18051793
18061794#ifdef DEBUG
18071795 ! Check for potential issues in wave speed calculation
18081796 if (s_R <= s_L) then
18091797 print * , ' WARNING: Wave speed issue detected in s_compute_wave_speed'
18101798 print * , ' Left wave speed >= Right wave speed:' , s_L, s_R
1811- print * , ' Input velocities :' , vel_L, vel_R
1799+ print * , ' Input velocities :' , vel_L(dir_idx( 1 )) , vel_R(dir_idx( 1 ))
18121800 print * , ' Sound speeds:' , c_L, c_R
18131801 print * , ' Densities:' , rho_L, rho_R
18141802 print * , ' Pressures:' , pres_L, pres_R
0 commit comments