MODULE Generic_EvalsAndETAX
INTERFACE EvalsAndETAX
MODULE PROCEDURE OneDEValsAndETAX
MODULE PROCEDURE TwoDEValsAndETAX
MODULE PROCEDURE ThreeDEValsAndETAX
END INTERFACE
CONTAINS
SUBROUTINE OneDEValsAndETAX(ux,vx,i,lami,ETA)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_EvalGamma
use cles_interfaces
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi)
INTEGER, INTENT(IN) :: i
DOUBLE PRECISION, DIMENSION(2,3),INTENT(OUT) :: lami
DOUBLE PRECISION, INTENT(OUT) :: ETA
DOUBLE PRECISION :: gamma
DOUBLE PRECISION, DIMENSION(0:1) :: csnd
DOUBLE PRECISION :: csndsq
DOUBLE PRECISION, DIMENSION(9) :: etaH
INTEGER :: ii,l
double precision :: xc
! ---- Pack the local speed of sounds
DO l=0,1
ii = i+l
CALL GetGamma(ux(:,ii), vx(:,ii), gamma)
csndsq = gamma*vx(5,ii)/vx(1,ii)
IF(csndsq.le.0.00) THEN
print *, 'Warning: CarbFix-X imaginary speed of sound at x, p, r'
call cles_xLocation(i,xc)
WRITE(*,'(3(1X,E12.6))') xc, vx(5,ii), vx(1,ii)
! STOP
csnd(l) = 0.0
ELSE
csnd(l) = dsqrt(csndsq)
END IF
END DO
! ---- Compute the local speeds ( eigen values)
lami(1,1) = vx(2,i)
lami(1,2) = vx(2,i) + csnd(0)
lami(1,3) = vx(2,i) - csnd(0)
lami(2,1) = vx(2,i+1)
lami(2,2) = vx(2,i+1) + csnd(1)
lami(2,3) = vx(2,i+1) - csnd(1)
etaH(:) = 0.0d0
! ---- Compute ETA from the H-fix (Carbuncle)
etaH(1)=maxval(abs(lami(2,:)-lami(1,:)))
ETA =0.5*maxval(etaH)
END SUBROUTINE OneDEValsAndETAX
SUBROUTINE TwoDEValsAndETAX(ux,vx,i,j,lami,ETA)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_EvalGamma
use cles_interfaces
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi)
INTEGER, INTENT(IN) :: i,j
DOUBLE PRECISION, DIMENSION(2,3),INTENT(OUT) :: lami
DOUBLE PRECISION, INTENT(OUT) :: ETA
DOUBLE PRECISION, DIMENSION(2,3) :: lamj,lamjp,lamjm
DOUBLE PRECISION, DIMENSION(0:1,-1:1) :: csnd
DOUBLE PRECISION :: csndsq
DOUBLE PRECISION :: gamma
DOUBLE PRECISION, DIMENSION(9) :: etaH
INTEGER :: ii,jj,l,m
double precision :: xc, yc
! ---- Pack the local speed of sounds
DO l=0,1
DO m = -1,1
ii = i+l
jj = j+m
CALL GetGamma(ux(:,ii,jj), vx(:,ii,jj), gamma)
csndsq = gamma * vx(5,ii,jj)/vx(1,ii,jj)
IF (csndsq.le.0.00) THEN
print *, 'Warning: CarbFix-X imaginary speed of sound at x, y, p, r'
call cles_xLocation(i,xc)
call cles_yLocation(j,yc)
WRITE(*,'(4(1X,E12.6))') xc, yc, vx(5,ii,jj), vx(1,ii,jj)
! STOP
csnd(l,m) = 0.0
ELSE
csnd(l,m) = dsqrt(csndsq)
END IF
END DO
END DO
! ---- Compute the local speeds ( eigen values)
lami(1,1) = vx(2,i,j)
lami(1,2) = vx(2,i,j) + csnd(0,0)
lami(1,3) = vx(2,i,j) - csnd(0,0)
lami(2,1) = vx(2,i+1,j)
lami(2,2) = vx(2,i+1,j) + csnd(1,0)
lami(2,3) = vx(2,i+1,j) - csnd(1,0)
lamj(1,1)=vx(3,i,j)
lamj(1,2)=vx(3,i,j)+csnd(0,0)
lamj(1,3)=vx(3,i,j)-csnd(0,0)
lamj(2,1)=vx(3,i+1,j)
lamj(2,2)=vx(3,i+1,j)+csnd(1,0)
lamj(2,3)=vx(3,i+1,j)-csnd(1,0)
lamjp(1,1)=vx(3,i,j+1)
lamjp(1,2)=vx(3,i,j+1)+csnd(0,1)
lamjp(1,3)=vx(3,i,j+1)-csnd(0,1)
lamjp(2,1)=vx(3,i+1,j+1)
lamjp(2,2)=vx(3,i+1,j+1)+csnd(1,1)
lamjp(2,3)=vx(3,i+1,j+1)-csnd(1,1)
lamjm(1,1)=vx(3,i,j-1)
lamjm(1,2)=vx(3,i,j-1)+csnd(0,-1)
lamjm(1,3)=vx(3,i,j-1)-csnd(0,-1)
lamjm(2,1)=vx(3,i+1,j-1)
lamjm(2,2)=vx(3,i+1,j-1)+csnd(1,-1)
lamjm(2,3)=vx(3,i+1,j-1)-csnd(1,-1)
etaH(:) = 0.0d0
! ---- Compute ETA from the H-fix (Carbuncle)
etaH(1)=maxval(abs(lami(2,:)-lami(1,:)))
etaH(2)=maxval(abs(lamjp(1,:)-lamj(1,:)))
etaH(3)=maxval(abs(lamj(1,:)-lamjm(1,:)))
etaH(4)=maxval(abs(lamjp(2,:)-lamj(2,:)))
etaH(5)=maxval(abs(lamj(2,:)-lamjm(2,:)))
ETA =0.5*maxval(etaH)
! ---- defeat Carbuncle fix
! ETA =0.5*etaH(1)
END SUBROUTINE TwoDEValsAndETAX
SUBROUTINE ThreeDEValsAndETAX(ux,vx,i,j,k,lami,ETA)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_EvalGamma
use cles_interfaces
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
INTEGER, INTENT(IN) :: i,j,k
DOUBLE PRECISION, DIMENSION(2,3),INTENT(OUT) :: lami
DOUBLE PRECISION, INTENT(OUT) :: ETA
DOUBLE PRECISION, DIMENSION(2,3) :: lamj,lamk,lamjp,lamjm,lamkp,lamkm
DOUBLE PRECISION, DIMENSION(0:1,-1:1,-1:1) :: csnd
DOUBLE PRECISION :: csndsq
DOUBLE PRECISION :: gamma
DOUBLE PRECISION, DIMENSION(9) :: etaH
INTEGER :: ii,jj,kk,l,m,n
double precision :: xc, yc, zc
! call cleslog_log_enter('ThreeDEValsAndETAX')
! ---- Pack the local speed of sounds
DO l=0,1
DO m = -1,1
DO n = -1,1
ii = i+l
jj = j+m
kk = k+n
CALL GetGamma(ux(:,ii,jj,kk),vx(:,ii,jj,kk), gamma)
csndsq = gamma*vx(5,ii,jj,kk)/vx(1,ii,jj,kk)
IF (csndsq.le.0.00) THEN
print *, 'Warning: CarbFix-X imaginary speed of sound at x, y, z, p, r'
call cles_xLocation(i,xc)
call cles_yLocation(j,yc)
call cles_zLocation(k,zc)
WRITE(*,'(5(1X,E12.6))') xc, yc, zc, &
vx(5,ii,jj,kk), vx(1,ii,jj,kk)
! STOP
csnd(l,m,n) = 0.0
ELSE
csnd(l,m,n) = dsqrt(csndsq)
END IF
END DO
END DO
END DO
! ---- Compute the local speeds ( eigen values)
lami(1,1) = vx(2,i,j,k)
lami(1,2) = vx(2,i,j,k) + csnd(0,0,0)
lami(1,3) = vx(2,i,j,k) - csnd(0,0,0)
lami(2,1) = vx(2,i+1,j,k)
lami(2,2) = vx(2,i+1,j,k) + csnd(1,0,0)
lami(2,3) = vx(2,i+1,j,k) - csnd(1,0,0)
lamj(1,1)=vx(3,i,j,k)
lamj(1,2)=vx(3,i,j,k)+csnd(0,0,0)
lamj(1,3)=vx(3,i,j,k)-csnd(0,0,0)
lamj(2,1)=vx(3,i+1,j,k)
lamj(2,2)=vx(3,i+1,j,k)+csnd(1,0,0)
lamj(2,3)=vx(3,i+1,j,k)-csnd(1,0,0)
lamjp(1,1)=vx(3,i,j+1,k)
lamjp(1,2)=vx(3,i,j+1,k)+csnd(0,1,0)
lamjp(1,3)=vx(3,i,j+1,k)-csnd(0,1,0)
lamjp(2,1)=vx(3,i+1,j+1,k)
lamjp(2,2)=vx(3,i+1,j+1,k)+csnd(1,1,0)
lamjp(2,3)=vx(3,i+1,j+1,k)-csnd(1,1,0)
lamjm(1,1)=vx(3,i,j-1,k)
lamjm(1,2)=vx(3,i,j-1,k)+csnd(0,-1,0)
lamjm(1,3)=vx(3,i,j-1,k)-csnd(0,-1,0)
lamjm(2,1)=vx(3,i+1,j-1,k)
lamjm(2,2)=vx(3,i+1,j-1,k)+csnd(1,-1,0)
lamjm(2,3)=vx(3,i+1,j-1,k)-csnd(1,-1,0)
lamk(1,1) = vx(4,i,j,k)
lamk(1,2) = vx(4,i,j,k) + csnd(0,0,0)
lamk(1,3) = vx(4,i,j,k) - csnd(0,0,0)
lamk(2,1) = vx(4,i+1,j,k)
lamk(2,2) = vx(4,i+1,j,k) + csnd(1,0,0)
lamk(2,3) = vx(4,i+1,j,k) - csnd(1,0,0)
lamkp(1,1) = vx(4,i,j,k+1)
lamkp(1,2) = vx(4,i,j,k+1) + csnd(0,0,1)
lamkp(1,3) = vx(4,i,j,k+1) - csnd(0,0,1)
lamkp(2,1) = vx(4,i+1,j,k+1)
lamkp(2,2) = vx(4,i+1,j,k+1) + csnd(1,0,1)
lamkp(2,3) = vx(4,i+1,j,k+1) - csnd(1,0,1)
lamkm(1,1) = vx(4,i,j,k-1)
lamkm(1,2) = vx(4,i,j,k-1) + csnd(0,0,-1)
lamkm(1,3) = vx(4,i,j,k-1) - csnd(0,0,-1)
lamkm(2,1) = vx(4,i+1,j,k-1)
lamkm(2,2) = vx(4,i+1,j,k-1) + csnd(1,0,-1)
lamkm(2,3) = vx(4,i+1,j,k-1) - csnd(1,0,-1)
! ---- Compute ETA from the H-fix (Carbuncle)
etaH(1)=maxval(abs(lami(2,:)-lami(1,:)))
etaH(2)=maxval(abs(lamjp(1,:)-lamj(1,:)))
etaH(3)=maxval(abs(lamj(1,:)-lamjm(1,:)))
etaH(4)=maxval(abs(lamjp(2,:)-lamj(2,:)))
etaH(5)=maxval(abs(lamj(2,:)-lamjm(2,:)))
etaH(6)=maxval(abs(lamkp(1,:)-lamk(1,:)))
etaH(7)=maxval(abs(lamk(1,:)-lamkm(1,:)))
etaH(8)=maxval(abs(lamkp(2,:)-lamk(2,:)))
etaH(9)=maxval(abs(lamk(2,:)-lamkm(2,:)))
ETA =0.5*maxval(etaH)
! call cleslog_log_exit('ThreeDEValsAndETAX')
END SUBROUTINE ThreeDEValsAndETAX
END MODULE Generic_EvalsAndETAX
MODULE Generic_EvalsAndETAY
INTERFACE EvalsAndETAY
MODULE PROCEDURE TwoDEValsAndETAY
MODULE PROCEDURE ThreeDEValsAndETAY
END INTERFACE
CONTAINS
SUBROUTINE TwoDEValsAndETAY(ux,vx,i,j,lami,ETA)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_EvalGamma
use cles_interfaces
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi)
INTEGER, INTENT(IN) :: i,j
DOUBLE PRECISION, DIMENSION(2,3),INTENT(OUT) :: lami
DOUBLE PRECISION, INTENT(OUT) :: ETA
DOUBLE PRECISION, DIMENSION(2,3) :: lamj,lamjp,lamjm
DOUBLE PRECISION, DIMENSION(-1:1,-1:1) :: csnd
DOUBLE PRECISION :: csndsq
DOUBLE PRECISION :: gamma
DOUBLE PRECISION, DIMENSION(9) :: etaH
INTEGER :: ii,jj,kk,l,m,n
double precision :: xc, yc
! ---- Pack the local speed of sounds
DO l=-1,1
DO m = -1,1
ii = i+l
jj = j+m
CALL GetGamma(ux(:,ii,jj), vx(:,ii,jj), gamma )
csndsq = gamma *vx(5,ii,jj)/vx(1,ii,jj)
IF (csndsq.le.0.00) THEN
print *, 'Warning: CarbFix-Y imaginary speed of sound at x, y, p, r'
call cles_xLocation(i,xc)
call cles_yLocation(j,yc)
WRITE(*,'(4(1X,E12.6))') xc, yc, vx(5,ii,jj), vx(1,ii,jj)
! STOP
csnd(l,m) = 0.0
ELSE
csnd(l,m) = sqrt(csndsq)
END IF
END DO
END DO
! ---- Compute the local speeds ( eigen values)
lami(1,1) = vx(3,i,j)
lami(1,2) = vx(3,i,j) + csnd(0,0)
lami(1,3) = vx(3,i,j) - csnd(0,0)
lami(2,1) = vx(3,i,j+1)
lami(2,2) = vx(3,i,j+1) + csnd(0,1)
lami(2,3) = vx(3,i,j+1) - csnd(0,1)
lamj(1,1)=vx(2,i,j)
lamj(1,2)=vx(2,i,j)+csnd(0,0)
lamj(1,3)=vx(2,i,j)-csnd(0,0)
lamj(2,1)=vx(2,i,j+1)
lamj(2,2)=vx(2,i,j+1)+csnd(0,1)
lamj(2,3)=vx(2,i,j+1)-csnd(0,1)
lamjp(1,1)=vx(2,i+1,j)
lamjp(1,2)=vx(2,i+1,j)+csnd(1,0)
lamjp(1,3)=vx(2,i+1,j)-csnd(1,0)
lamjp(2,1)=vx(2,i+1,j+1)
lamjp(2,2)=vx(2,i+1,j+1)+csnd(1,1)
lamjp(2,3)=vx(2,i+1,j+1)-csnd(1,1)
lamjm(1,1)=vx(2,i-1,j)
lamjm(1,2)=vx(2,i-1,j)+csnd(-1,0)
lamjm(1,3)=vx(2,i-1,j)-csnd(-1,0)
lamjm(2,1)=vx(2,i-1,j+1)
lamjm(2,2)=vx(2,i-1,j+1)+csnd(-1,1)
lamjm(2,3)=vx(2,i-1,j+1)-csnd(-1,1)
etaH(:) = 0.0d0
! ---- Compute ETA from the H-fix (Carbuncle)
etaH(1)=maxval(abs(lami(2,:)-lami(1,:)))
etaH(2)=maxval(abs(lamjp(1,:)-lamj(1,:)))
etaH(3)=maxval(abs(lamj(1,:)-lamjm(1,:)))
etaH(4)=maxval(abs(lamjp(2,:)-lamj(2,:)))
etaH(5)=maxval(abs(lamj(2,:)-lamjm(2,:)))
ETA =0.5*maxval(etaH)
! ---- defeat carbuncle
! ETA =0.5*etaH(1)
END SUBROUTINE TwoDEValsAndETAY
SUBROUTINE ThreeDEValsAndETAY(ux,vx,i,j,k,lami,ETA)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_EvalGamma
use cles_interfaces
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
INTEGER, INTENT(IN) :: i,j,k
DOUBLE PRECISION, DIMENSION(2,3),INTENT(OUT) :: lami
DOUBLE PRECISION, INTENT(OUT) :: ETA
DOUBLE PRECISION, DIMENSION(2,3) :: lamj,lamk,lamjp,lamjm,lamkp,lamkm
DOUBLE PRECISION, DIMENSION(-1:1,0:1,-1:1) :: csnd
DOUBLE PRECISION :: csndsq
DOUBLE PRECISION :: gamma
DOUBLE PRECISION, DIMENSION(9) :: etaH
INTEGER :: ii,jj,kk,l,m,n
double precision :: xc, yc, zc
! call cleslog_log_enter('ThreeDEValsAndETAY')
! ---- Pack the local speed of sounds
DO l=-1,1
DO m = 0,1
DO n = -1,1
ii = i+l
jj = j+m
kk = k+n
CALL GetGamma(ux(:,ii,jj,kk),vx(:,ii,jj,kk), gamma)
csndsq = gamma &
& *vx(5,ii,jj,kk)/vx(1,ii,jj,kk)
IF (csndsq.le.0.00) THEN
print *, 'Warning: CarbFix-Y imaginary speed of sound at x, y, z, p, r'
call cles_xLocation(i,xc)
call cles_yLocation(j,yc)
call cles_zLocation(k,zc)
WRITE(*,'(5(1X,E12.6))') xc, yc, zc, &
vx(5,ii,jj,kk), vx(1,ii,jj,kk)
! STOP
csnd(l,m,n) = 0.0
ELSE
csnd(l,m,n) = sqrt(csndsq)
END IF
END DO
END DO
END DO
! ---- Compute the local speeds ( eigen values)
lami(1,1) = vx(3,i,j,k)
lami(1,2) = vx(3,i,j,k) + csnd(0,0,0)
lami(1,3) = vx(3,i,j,k) - csnd(0,0,0)
lami(2,1) = vx(3,i,j+1,k)
lami(2,2) = vx(3,i,j+1,k) + csnd(0,1,0)
lami(2,3) = vx(3,i,j+1,k) - csnd(0,1,0)
lamj(1,1)=vx(2,i,j,k)
lamj(1,2)=vx(2,i,j,k)+csnd(0,0,0)
lamj(1,3)=vx(2,i,j,k)-csnd(0,0,0)
lamj(2,1)=vx(2,i,j+1,k)
lamj(2,2)=vx(2,i,j+1,k)+csnd(0,1,0)
lamj(2,3)=vx(2,i,j+1,k)-csnd(0,1,0)
lamjp(1,1)=vx(2,i+1,j,k)
lamjp(1,2)=vx(2,i+1,j,k)+csnd(1,0,0)
lamjp(1,3)=vx(2,i+1,j,k)-csnd(1,0,0)
lamjp(2,1)=vx(2,i+1,j+1,k)
lamjp(2,2)=vx(2,i+1,j+1,k)+csnd(1,1,0)
lamjp(2,3)=vx(2,i+1,j+1,k)-csnd(1,1,0)
lamjm(1,1)=vx(2,i-1,j,k)
lamjm(1,2)=vx(2,i-1,j,k)+csnd(-1,0,0)
lamjm(1,3)=vx(2,i-1,j,k)-csnd(-1,0,0)
lamjm(2,1)=vx(2,i-1,j+1,k)
lamjm(2,2)=vx(2,i-1,j+1,k)+csnd(-1,1,0)
lamjm(2,3)=vx(2,i-1,j+1,k)-csnd(-1,1,0)
lamk(1,1) = vx(4,i,j,k)
lamk(1,2) = vx(4,i,j,k) + csnd(0,0,0)
lamk(1,3) = vx(4,i,j,k) - csnd(0,0,0)
lamk(2,1) = vx(4,i,j+1,k)
lamk(2,2) = vx(4,i,j+1,k) + csnd(0,1,0)
lamk(2,3) = vx(4,i,j+1,k) - csnd(0,1,0)
lamkp(1,1) = vx(4,i,j,k+1)
lamkp(1,2) = vx(4,i,j,k+1) + csnd(0,0,1)
lamkp(1,3) = vx(4,i,j,k+1) - csnd(0,0,1)
lamkp(2,1) = vx(4,i,j+1,k+1)
lamkp(2,2) = vx(4,i,j+1,k+1) + csnd(0,1,1)
lamkp(2,3) = vx(4,i,j+1,k+1) - csnd(0,1,1)
lamkm(1,1) = vx(4,i,j,k-1)
lamkm(1,2) = vx(4,i,j,k-1) + csnd(0,0,-1)
lamkm(1,3) = vx(4,i,j,k-1) - csnd(0,0,-1)
lamkm(2,1) = vx(4,i,j+1,k-1)
lamkm(2,2) = vx(4,i,j+1,k-1) + csnd(0,1,-1)
lamkm(2,3) = vx(4,i,j+1,k-1) - csnd(0,1,-1)
! ---- Compute ETA from the H-fix (Carbuncle)
etaH(1)=maxval(abs(lami(2,:)-lami(1,:)))
etaH(2)=maxval(abs(lamjp(1,:)-lamj(1,:)))
etaH(3)=maxval(abs(lamj(1,:)-lamjm(1,:)))
etaH(4)=maxval(abs(lamjp(2,:)-lamj(2,:)))
etaH(5)=maxval(abs(lamj(2,:)-lamjm(2,:)))
etaH(6)=maxval(abs(lamkp(1,:)-lamk(1,:)))
etaH(7)=maxval(abs(lamk(1,:)-lamkm(1,:)))
etaH(8)=maxval(abs(lamkp(2,:)-lamk(2,:)))
etaH(9)=maxval(abs(lamk(2,:)-lamkm(2,:)))
ETA =0.5*maxval(etaH)
! call cleslog_log_exit('ThreeDEValsAndETAY')
END SUBROUTINE ThreeDEValsAndETAY
END MODULE Generic_EvalsAndETAY
MODULE Generic_EvalsAndETAZ
INTERFACE EvalsAndETAZ
MODULE PROCEDURE ThreeDEValsAndETAZ
END INTERFACE
CONTAINS
SUBROUTINE ThreeDEValsAndETAZ(ux,vx,i,j,k,lami,ETA)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_EvalGamma
use cles_interfaces
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
INTEGER, INTENT(IN) :: i,j,k
DOUBLE PRECISION, DIMENSION(2,3),INTENT(OUT) :: lami
DOUBLE PRECISION, INTENT(OUT) :: ETA
DOUBLE PRECISION, DIMENSION(2,3) :: lamj,lamk,lamjp,lamjm,lamkp,lamkm
DOUBLE PRECISION, DIMENSION(-1:1,-1:1,0:1) :: csnd
DOUBLE PRECISION :: csndsq
DOUBLE PRECISION :: gamma
DOUBLE PRECISION, DIMENSION(9) :: etaH
INTEGER :: ii,jj,kk,l,m,n
double precision :: xc, yc, zc
! call cleslog_log_enter('ThreeDEValsAndETAZ')
! ---- Pack the local speed of sounds
DO l=-1,1
DO m = -1,1
DO n = 0,1
ii = i+l
jj = j+m
kk = k+n
CALL GetGamma(ux(:,ii,jj,kk), vx(:,ii,jj,kk), gamma)
csndsq = gamma * vx(5,ii,jj,kk)/vx(1,ii,jj,kk)
IF (csndsq.le.0.00) THEN
print *, 'Warning: CarbFix-Z imaginary speed of sound at x, y, z, p, r'
call cles_xLocation(i,xc)
call cles_yLocation(j,yc)
call cles_zLocation(k,zc)
WRITE(*,'(5(1X,E12.6))') xc, yc, zc, &
vx(5,ii,jj,kk), vx(1,ii,jj,kk)
! STOP
csnd(l,m,n) = 0.0
ELSE
csnd(l,m,n) = dsqrt(csndsq)
END IF
END DO
END DO
END DO
! ---- Compute the local speeds ( eigen values)
lami(1,1) = vx(4,i,j,k)
lami(1,2) = vx(4,i,j,k) + csnd(0,0,0)
lami(1,3) = vx(4,i,j,k) - csnd(0,0,0)
lami(2,1) = vx(4,i,j,k+1)
lami(2,2) = vx(4,i,j,k+1) + csnd(0,0,1)
lami(2,3) = vx(4,i,j,k+1) - csnd(0,0,1)
lamj(1,1)=vx(2,i,j,k)
lamj(1,2)=vx(2,i,j,k)+csnd(0,0,0)
lamj(1,3)=vx(2,i,j,k)-csnd(0,0,0)
lamj(2,1)=vx(2,i,j,k+1)
lamj(2,2)=vx(2,i,j,k+1)+csnd(0,0,1)
lamj(2,3)=vx(2,i,j,k+1)-csnd(0,0,1)
lamjp(1,1)=vx(2,i+1,j,k)
lamjp(1,2)=vx(2,i+1,j,k)+csnd(1,0,0)
lamjp(1,3)=vx(2,i+1,j,k)-csnd(1,0,0)
lamjp(2,1)=vx(2,i+1,j,k+1)
lamjp(2,2)=vx(2,i+1,j,k+1)+csnd(1,0,1)
lamjp(2,3)=vx(2,i+1,j,k+1)-csnd(1,0,1)
lamjm(1,1)=vx(2,i-1,j,k)
lamjm(1,2)=vx(2,i-1,j,k)+csnd(-1,0,0)
lamjm(1,3)=vx(2,i-1,j,k)-csnd(-1,0,0)
lamjm(2,1)=vx(2,i-1,j,k+1)
lamjm(2,2)=vx(2,i-1,j,k+1)+csnd(-1,0,1)
lamjm(2,3)=vx(2,i-1,j,k+1)-csnd(-1,0,1)
lamk(1,1) = vx(3,i,j,k)
lamk(1,2) = vx(3,i,j,k) + csnd(0,0,0)
lamk(1,3) = vx(3,i,j,k) - csnd(0,0,0)
lamk(2,1) = vx(3,i,j,k+1)
lamk(2,2) = vx(3,i,j,k+1) + csnd(0,0,1)
lamk(2,3) = vx(3,i,j,k+1) - csnd(0,0,1)
lamkp(1,1) = vx(3,i,j+1,k)
lamkp(1,2) = vx(3,i,j+1,k) + csnd(0,1,0)
lamkp(1,3) = vx(3,i,j+1,k) - csnd(0,1,0)
lamkp(2,1) = vx(3,i,j+1,k+1)
lamkp(2,2) = vx(3,i,j+1,k+1) + csnd(0,1,1)
lamkp(2,3) = vx(3,i,j+1,k+1) - csnd(0,1,1)
lamkm(1,1) = vx(3,i,j-1,k)
lamkm(1,2) = vx(3,i,j-1,k) + csnd(0,-1,0)
lamkm(1,3) = vx(3,i,j-1,k) - csnd(0,-1,0)
lamkm(2,1) = vx(3,i,j-1,k+1)
lamkm(2,2) = vx(3,i,j-1,k+1) + csnd(0,-1,1)
lamkm(2,3) = vx(3,i,j-1,k+1) - csnd(0,-1,1)
! ---- Compute ETA from the H-fix (Carbuncle)
etaH(1)=maxval(abs(lami(2,:)-lami(1,:)))
etaH(2)=maxval(abs(lamjp(1,:)-lamj(1,:)))
etaH(3)=maxval(abs(lamj(1,:)-lamjm(1,:)))
etaH(4)=maxval(abs(lamjp(2,:)-lamj(2,:)))
etaH(5)=maxval(abs(lamj(2,:)-lamjm(2,:)))
etaH(6)=maxval(abs(lamkp(1,:)-lamk(1,:)))
etaH(7)=maxval(abs(lamk(1,:)-lamkm(1,:)))
etaH(8)=maxval(abs(lamkp(2,:)-lamk(2,:)))
etaH(9)=maxval(abs(lamk(2,:)-lamkm(2,:)))
ETA =0.5*maxval(etaH)
! call cleslog_log_exit('ThreeDEValsAndETAZ')
END SUBROUTINE ThreeDEValsAndETAZ
END MODULE Generic_EvalsAndETAZ