Module Generic_FDInterpFlux
! ---- This module computes the interpolation of the flux vector fx
! ---- from the grid points to the 1/2 way points, eg dx*(j + 1/2 )
! ---- values are then stored in fxi(j).
INTERFACE FDInterpFluxSkew
MODULE PROCEDURE OneDFDInterpFluxSkew
MODULE PROCEDURE TwoDFDInterpFluxSkew
MODULE PROCEDURE ThreeDFDInterpFluxSkew
END INTERFACE
CONTAINS
SUBROUTINE OneDFDInterpFluxSkew(ux, vx, fx, fxi, direction)
! ---- Shared Variables
USE array_bounds
USE method_parms
USE Interp_coeffs
! ----
! ---- Shared Procedures
USE Generic_FDInterp
USE Generic_FDSkewInterp
! ----
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi)
! the cell centered flux
DOUBLE PRECISION, INTENT(IN) :: fx(nvars,ixlo:ixhi)
! the 1/2 cell intrps
DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,1)
INTEGER, INTENT(IN) :: direction
DOUBLE PRECISION, ALLOCATABLE :: rayfx(:,:) ! fx used in sweep
DOUBLE PRECISION, ALLOCATABLE :: rayfxi(:,:) ! fxi computed in sweep
DOUBLE PRECISION, ALLOCATABLE :: skw(:), rhou(:), phi(:) ! auxiliary
INTEGER :: is
INTEGER :: lol ! local lo bound
INTEGER :: hil ! local hi bound
INTEGER :: iret
whichdirection: SELECT CASE(direction)
CASE (1) ! ---- the x-direction
ALLOCATE ( rayfx(nvars,ixlo:ixhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'OneDFDInterpFluxSkew/rayfx', iret)
endif
ALLOCATE ( rayfxi(nvars,1:nx+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'OneDFDInterpFluxSkew/rayfxi', iret)
endif
ALLOCATE ( skw(1:nx+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'OneDFDInterpFluxSkew/skw', iret)
endif
ALLOCATE ( rhou(ixlo:ixhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'OneDFDInterpFluxSkew/rhou', iret)
endif
ALLOCATE ( phi(ixlo:ixhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'OneDFDInterpFluxSkew/phi', iret)
endif
! ---- holds (horizontal, ie x) a slice of the flux
! ---- at a given y location including ghostcells
rayfx(:,:) = fx(:,:)
lol = ixlo
hil = ixhi
CALL FDInterpolateX(nvars, rayfx,rayfxi,lol,hil,nx)
fxi(1,1:nx+1,direction) = rayfxi(1,1:nx+1)
fxi(2:nvars,1:nx+1,direction) = rayfxi(2:nvars,1:nx+1)/2.0d0
! Add other pressure half
phi(:) = vx(5,:)
CALL FDInterpolateX(1, phi,skw,lol,hil,nx)
fxi(2,1:nx+1,direction) = fxi(2,1:nx+1,direction) + skw(1:nx+1)/2.0d0
! Add nonlinear part
rhou(:) = ux(2,:)
phi = vx(2,:)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(2,1:nx+1,direction) = fxi(2,1:nx+1,direction) + skw(1:nx+1)/2.0d0
fxi(3,1:nx+1,direction) = 0.0d0
fxi(4,1:nx+1,direction) = 0.0d0
! Scalars
do is=6, nvars
phi(:) = vx(is,:)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(is,1:nx+1,direction) = fxi(is,1:nx+1,direction) &
+ skw(1:nx+1)/2.0d0
enddo
! Handling of the energy equation
! compute internal energy
phi(:) = (ux(5,:)-0.5d0*(ux(2,:)**2)/ux(1,:))/ux(1,:)
! skew half
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,direction) = skw(1:nx+1)/2.0d0
! conservative half
phi(:) = phi(:)*rhou(:)
CALL FDInterpolateX(1, phi,skw,lol,hil,nx)
fxi(5,1:nx+1,direction) = fxi(5,1:nx+1,direction) &
+ skw(1:nx+1)/2.0d0
! Quadratic term
rhou(:) = vx(2,:)*ux(2,:)
phi(:) = vx(2,:)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,direction) = fxi(5,1:nx+1,direction) &
+ skw(1:nx+1)/2.0d0
! Add pressure-velocity non-conservative part
phi(:) = vx(2,:)
rhou(:) = vx(5,:)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,direction) = fxi(5,1:nx+1,direction) &
+ skw(1:nx+1)
DEALLOCATE ( rayfx , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'OneDFDInterpFluxSkew/rayfx', iret)
endif
DEALLOCATE ( rayfxi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'OneDFDInterpFluxSkew/rayfxi', iret)
endif
DEALLOCATE ( skw, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'OneDFDInterpFluxSkew/skw', iret)
endif
DEALLOCATE ( rhou, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'OneDFDInterpFluxSkew/rhou', iret)
endif
DEALLOCATE ( phi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'OneDFDInterpFluxSkew/phi', iret)
endif
END SELECT whichdirection
END SUBROUTINE OneDFDInterpFluxSkew
SUBROUTINE TwoDFDInterpFluxSkew(ux, vx, fx, fxi, direction)
! ---- Shared Variables
USE array_bounds
USE method_parms
USE Interp_coeffs
! ----
! ---- Shared Procedures
USE Generic_FDInterp
USE Generic_FDSkewInterp
! ----
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi)
! the cell centered flux
DOUBLE PRECISION, INTENT(IN) :: fx(nvars,ixlo:ixhi,iylo:iyhi)
! the 1/2 cell intrps
DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,2)
INTEGER, INTENT(IN) :: direction
DOUBLE PRECISION, ALLOCATABLE :: rayfx(:,:) ! fx used in sweep
DOUBLE PRECISION, ALLOCATABLE :: rayfxi(:,:) ! fxi computed in sweep
DOUBLE PRECISION, ALLOCATABLE :: skw(:), phi(:) ! auxiliary
DOUBLE PRECISION, ALLOCATABLE :: rhou(:), rhov(:)
INTEGER :: i,j, is
INTEGER :: lol ! local lo bound
INTEGER :: hil ! local hi bound
INTEGER :: iret
whichdirection: SELECT CASE(direction)
CASE (1) ! ---- the x-direction
ALLOCATE ( rayfx(nvars,ixlo:ixhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/rayfx', iret)
endif
ALLOCATE ( rayfxi(nvars,1:nx+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/rayfxi', iret)
endif
ALLOCATE ( skw(1:nx+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/skw', iret)
endif
ALLOCATE ( rhou(ixlo:ixhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/rhou', iret)
endif
ALLOCATE ( phi(ixlo:ixhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/phi', iret)
endif
DO j = 1, ny, 1
! ---- holds (horizontal, ie x) a slice of the flux
! ---- at a given y location including ghostcells
rayfx(:,:) = fx(:,:,j)
lol = ixlo
hil = ixhi
CALL FDInterpolateX(nvars, rayfx,rayfxi,lol,hil,nx)
fxi(1,1:nx+1,j,direction) = rayfxi(1,1:nx+1)
fxi(2:nvars,1:nx+1,j,direction) = rayfxi(2:nvars,1:nx+1)/2.0d0
! Add other pressure half
phi(:) = vx(5,:,j)
CALL FDInterpolateX(1, phi,skw,lol,hil,nx)
fxi(2,1:nx+1,j,direction) = fxi(2,1:nx+1,j,direction) &
+ skw(1:nx+1)/2.0d0
! Add nonlinear part
rhou(:) = ux(2,:,j)
phi(:) = vx(2,:,j)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(2,1:nx+1,j,direction) = fxi(2,1:nx+1,j,direction) + skw(1:nx+1)/2.0d0
! Add nonlinear part
phi(:) = vx(3,:,j)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(3,1:nx+1,j,direction) = fxi(3,1:nx+1,j,direction) + skw(1:nx+1)/2.0d0
fxi(4,1:nx+1,j,direction) = 0.0d0
! scalars
do is=6, nvars
phi(:) = vx(is,:,j)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(is,1:nx+1,j,direction) = fxi(is,1:nx+1,j,direction) &
+ skw(1:nx+1)/2.0d0
enddo
! Handling of the energy equation
! compute internal energy
phi(:) = (ux(5,:,j)-0.5d0*(ux(2,:,j)**2+ux(3,:,j)**2 &
)/ux(1,:,j))/ux(1,:,j)
! skew half
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,direction) = skw(1:nx+1)/2.0d0
! conservative half
phi(:) = phi(:)*rhou(:)
CALL FDInterpolateX(1, phi,skw,lol,hil,nx)
fxi(5,1:nx+1,j,direction) = fxi(5,1:nx+1,j,direction) &
+ skw(1:nx+1)/2.0d0
! Quadratic term
rhou(:) = vx(2,:,j)*ux(2,:,j)
phi(:) = vx(2,:,j)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,direction) = fxi(5,1:nx+1,j,direction) &
+ skw(1:nx+1)/2.0d0
rhou(:) = vx(3,:,j)*ux(2,:,j)
phi(:) = vx(3,:,j)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,direction) = fxi(5,1:nx+1,j,direction) &
+ skw(1:nx+1)/2.0d0
! Add pressure-velocity non-conservative part
phi(:) = vx(2,:,j)
rhou(:) = vx(5,:,j)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,direction) = fxi(5,1:nx+1,j,direction) &
+ skw(1:nx+1)
END DO
DEALLOCATE ( rayfx , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/rayfx', iret)
endif
DEALLOCATE ( rayfxi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/rayfxi', iret)
endif
DEALLOCATE ( skw, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/skw', iret)
endif
DEALLOCATE ( rhou, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/rhou', iret)
endif
DEALLOCATE ( phi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/phi', iret)
endif
CASE (2) ! ---- the y-direction
ALLOCATE ( rayfx(nvars,iylo:iyhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/rayfx', iret)
endif
ALLOCATE ( rayfxi(nvars,1:ny+1), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/rayfxi', iret)
endif
ALLOCATE ( skw(1:ny+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/skw', iret)
endif
ALLOCATE ( rhov(iylo:iyhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/rhov', iret)
endif
ALLOCATE ( phi(iylo:iyhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'TwoDFDInterpFluxSkew/phi', iret)
endif
DO i = 1, nx, 1
! ---- holds (vertical, ie y) a slice of the flux
! ---- at a given x location including ghostcells
rayfx(:,:) = fx(:,i,:)
lol = iylo
hil = iyhi
CALL FDInterpolateY(nvars, rayfx,rayfxi,lol,hil,ny)
fxi(1,i,1:ny+1,direction) = rayfxi(1,1:ny+1)
fxi(2:nvars,i,1:ny+1,direction) = rayfxi(2:nvars,1:ny+1)/2.0d0
! Add nonlinear part
rhov(:) = ux(3,i,:)
phi(:) = vx(2,i,:)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(2,i,1:ny+1,direction) = fxi(2,i,1:ny+1,direction) + skw(1:ny+1)/2.0d0
! Add other pressure half
phi(:) = vx(5,i,:)
CALL FDInterpolateY(1, phi,skw,lol,hil,ny)
fxi(3,i,1:ny+1,direction) = fxi(3,i,1:ny+1,direction) &
+ skw(1:ny+1)/2.0d0
! Add nonlinear part
phi(:) = vx(3,i,:)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(3,i,1:ny+1,direction) = fxi(3,i,1:ny+1,direction) + skw(1:ny+1)/2.0d0
fxi(4,i,1:ny+1,direction) = 0.0d0
! Scalars
do is=6, nvars
phi(:) = vx(is,i,:)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(is,i,1:ny+1,direction) = fxi(is,i,1:ny+1,direction) &
+ skw(1:ny+1)/2.0d0
enddo
! Handling of the energy equation
! compute internal energy
phi(:) = (ux(5,i,:)-0.5d0*(ux(2,i,:)**2+ux(3,i,:)**2 &
)/ux(1,i,:))/ux(1,i,:)
! skew half
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,direction) = skw(1:ny+1)/2.0d0
! conservative half
phi(:) = phi(:)*rhov(:)
CALL FDInterpolateY(1, phi,skw,lol,hil,ny)
fxi(5,i,1:ny+1,direction) = fxi(5,i,1:ny+1,direction) &
+ skw(1:ny+1)/2.0d0
! Quadratic term
rhov(:) = vx(2,i,:)*ux(3,i,:)
phi(:) = vx(2,i,:)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,direction) = fxi(5,i,1:ny+1,direction) &
+ skw(1:ny+1)/2.0d0
rhov(:) = vx(3,i,:)*ux(3,i,:)
phi(:) = vx(3,i,:)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,direction) = fxi(5,i,1:ny+1,direction) &
+ skw(1:ny+1)/2.0d0
! Add pressure-velocity non-conservative part
phi(:) = vx(3,i,:)
rhov(:) = vx(5,i,:)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,direction) = fxi(5,i,1:ny+1,direction) &
+ skw(1:ny+1)
END DO
DEALLOCATE ( rayfx, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/rayfx', iret)
endif
DEALLOCATE ( rayfxi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/rayfxi', iret)
endif
DEALLOCATE ( skw, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/skw', iret)
endif
DEALLOCATE ( rhov, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/rhov', iret)
endif
DEALLOCATE ( phi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'TwoDFDInterpFluxSkew/phi', iret)
endif
END SELECT whichdirection
END SUBROUTINE TwoDFDInterpFluxSkew
SUBROUTINE ThreeDFDInterpFluxSkew(ux, vx, fx, fxi, direction)
! ---- Shared Variables
USE array_bounds
USE method_parms
USE Interp_coeffs
! ----
! ---- Shared Procedures
USE Generic_FDInterp
USE Generic_FDSkewInterp
! ----
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
! the cell centered flux
DOUBLE PRECISION, INTENT(IN) :: fx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
! the 1/2 cell intrps
DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi,3)
INTEGER, INTENT(IN) :: direction
DOUBLE PRECISION, ALLOCATABLE :: rayfx(:,:) ! fx used in sweep
DOUBLE PRECISION, ALLOCATABLE :: rayfxi(:,:) ! fxi computed in sweep
DOUBLE PRECISION, ALLOCATABLE :: skw(:), phi(:) ! auxiliary
DOUBLE PRECISION, ALLOCATABLE :: rhou(:), rhov(:), rhow(:)
INTEGER :: i,j,k,is
INTEGER :: lol ! local lo bound
INTEGER :: hil ! local hi bound
INTEGER :: iret
call cleslog_log_enter('ThreeDFDInterpFluxSkew')
whichdirection: SELECT CASE(direction)
CASE (1) ! ---- the x-direction
ALLOCATE ( rayfx(nvars,ixlo:ixhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfx', iret)
endif
ALLOCATE ( rayfxi(nvars,1:nx+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfxi', iret)
endif
ALLOCATE ( skw(1:nx+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/skw', iret)
endif
ALLOCATE ( rhou(ixlo:ixhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rhou', iret)
endif
ALLOCATE ( phi(ixlo:ixhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/phi', iret)
endif
DO k=1, nz, 1
DO j = 1, ny, 1
! ---- holds (horizontal, ie x) a slice of the flux
! ---- at a given y location including ghostcells
rayfx(:,:) = fx(:,:,j,k)
lol = ixlo
hil = ixhi
CALL FDInterpolateX(nvars, rayfx,rayfxi,lol,hil,nx)
fxi(1,1:nx+1,j,k,direction) = rayfxi(1,1:nx+1)
fxi(2:nvars,1:nx+1,j,k,direction) = rayfxi(2:nvars,1:nx+1)/2.0d0
! Add other pressure half
phi(:) = vx(5,:,j,k)
CALL FDInterpolateX(1, phi,skw,lol,hil,nx)
fxi(2,1:nx+1,j,k,direction) = fxi(2,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
! Add nonlinear part
rhou(:) = ux(2,:,j,k)
phi(:) = vx(2,:,j,k)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(2,1:nx+1,j,k,direction) = fxi(2,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
! Add nonlinear part
phi(:) = vx(3,:,j,k)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(3,1:nx+1,j,k,direction) = fxi(3,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
! Add nonlinear part
phi(:) = vx(4,:,j,k)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(4,1:nx+1,j,k,direction) = fxi(4,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
! Scalars
do is=6, nvars
phi(:) = vx(is,:,j,k)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(is,1:nx+1,j,k,direction) = fxi(is,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
enddo
! Handling of the energy equation
! compute internal energy
phi(:) = (ux(5,:,j,k)-0.5d0*(ux(2,:,j,k)**2+ux(3,:,j,k)**2 &
+ux(4,:,j,k)**2)/ux(1,:,j,k))/ux(1,:,j,k)
! skew half
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,k,direction) = skw(1:nx+1)/2.0d0
! conservative half
phi(:) = phi(:)*rhou(:)
CALL FDInterpolateX(1, phi,skw,lol,hil,nx)
fxi(5,1:nx+1,j,k,direction) = fxi(5,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
! Quadratic term
rhou(:) = vx(2,:,j,k)*ux(2,:,j,k)
phi(:) = vx(2,:,j,k)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,k,direction) = fxi(5,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
rhou(:) = vx(3,:,j,k)*ux(2,:,j,k)
phi(:) = vx(3,:,j,k)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,k,direction) = fxi(5,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
rhou(:) = vx(4,:,j,k)*ux(2,:,j,k)
phi(:) = vx(4,:,j,k)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,k,direction) = fxi(5,1:nx+1,j,k,direction) &
+ skw(1:nx+1)/2.0d0
! Add pressure-velocity non-conservative part
phi(:) = vx(2,:,j,k)
rhou(:) = vx(5,:,j,k)
CALL FDSkewInterpolateX(rhou, phi, skw,lol,hil,nx)
fxi(5,1:nx+1,j,k,direction) = fxi(5,1:nx+1,j,k,direction) &
+ skw(1:nx+1)
END DO
END DO
DEALLOCATE ( rayfx , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfx', iret)
endif
DEALLOCATE ( rayfxi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfxi', iret)
endif
DEALLOCATE ( skw, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/skw', iret)
endif
DEALLOCATE ( rhou, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rhou', iret)
endif
DEALLOCATE ( phi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/phi', iret)
endif
CASE (2) ! ---- the y-direction
ALLOCATE ( rayfx(nvars,iylo:iyhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfx', iret)
endif
ALLOCATE ( rayfxi(nvars,1:ny+1), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfxi', iret)
endif
ALLOCATE ( skw(1:ny+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/skw', iret)
endif
ALLOCATE ( rhov(iylo:iyhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rhov', iret)
endif
ALLOCATE ( phi(iylo:iyhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/phi', iret)
endif
DO k = 1, nz, 1
DO i = 1, nx, 1
! ---- holds (vertical, ie y) a slice of the flux
! ---- at a given x location including ghostcells
rayfx(:,:) = fx(:,i,:,k)
lol = iylo
hil = iyhi
CALL FDInterpolateY(nvars, rayfx,rayfxi,lol,hil,ny)
fxi(1,i,1:ny+1,k,direction) = rayfxi(1,1:ny+1)
fxi(2:nvars,i,1:ny+1,k,direction) = rayfxi(2:nvars,1:ny+1)/2.0d0
! Add nonlinear part
rhov(:) = ux(3,i,:,k)
phi(:) = vx(2,i,:,k)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(2,i,1:ny+1,k,direction) = fxi(2,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
! Add other pressure half
phi(:) = vx(5,i,:,k)
CALL FDInterpolateY(1, phi,skw,lol,hil,ny)
fxi(3,i,1:ny+1,k,direction) = fxi(3,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
! Add nonlinear part
phi(:) = vx(3,i,:,k)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(3,i,1:ny+1,k,direction) = fxi(3,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
! Add nonlinear part
phi(:) = vx(4,i,:,k)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(4,i,1:ny+1,k,direction) = fxi(4,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
! Scalars
do is=6, nvars
phi(:) = vx(is,i,:,k)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(is,i,1:ny+1,k,direction) = fxi(is,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
enddo
! Handling of the energy equation
! compute internal energy
phi(:) = (ux(5,i,:,k)-0.5d0*(ux(2,i,:,k)**2+ux(3,i,:,k)**2 &
+ux(4,i,:,k)**2)/ux(1,i,:,k))/ux(1,i,:,k)
! skew half
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,k,direction) = skw(1:ny+1)/2.0d0
! conservative half
phi(:) = phi(:)*rhov(:)
CALL FDInterpolateY(1, phi,skw,lol,hil,ny)
fxi(5,i,1:ny+1,k,direction) = fxi(5,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
! Quadratic term
rhov(:) = vx(2,i,:,k)*ux(3,i,:,k)
phi(:) = vx(2,i,:,k)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,k,direction) = fxi(5,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
rhov(:) = vx(3,i,:,k)*ux(3,i,:,k)
phi(:) = vx(3,i,:,k)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,k,direction) = fxi(5,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
rhov(:) = vx(4,i,:,k)*ux(3,i,:,k)
phi(:) = vx(4,i,:,k)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,k,direction) = fxi(5,i,1:ny+1,k,direction) &
+ skw(1:ny+1)/2.0d0
! Add pressure-velocity non-conservative part
phi(:) = vx(3,i,:,k)
rhov(:) = vx(5,i,:,k)
CALL FDSkewInterpolateY(rhov, phi, skw,lol,hil,ny)
fxi(5,i,1:ny+1,k,direction) = fxi(5,i,1:ny+1,k,direction) &
+ skw(1:ny+1)
END DO
END DO
DEALLOCATE ( rayfx, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfx', iret)
endif
DEALLOCATE ( rayfxi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfxi', iret)
endif
DEALLOCATE ( skw, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/skw', iret)
endif
DEALLOCATE ( rhov, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rhov', iret)
endif
DEALLOCATE ( phi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/phi', iret)
endif
CASE (3) ! ---- the z-direction
ALLOCATE ( rayfx(nvars,izlo:izhi), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfx', iret)
endif
ALLOCATE ( rayfxi(nvars,1:nz+1), stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfxi', iret)
endif
ALLOCATE ( skw(1:nz+1) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/skw', iret)
endif
ALLOCATE ( rhow(izlo:izhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/rhow', iret)
endif
ALLOCATE ( phi(izlo:izhi) , stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'ThreeDFDInterpFluxSkew/phi', iret)
endif
DO j = 1, ny, 1
DO i = 1, nx, 1
! ---- holds (vertical, ie y) a slice of the flux
! ---- at a given x location including ghostcells
rayfx(:,:) = fx(:,i,j,:)
lol = izlo
hil = izhi
CALL FDInterpolateZ(nvars, rayfx,rayfxi,lol,hil,nz)
fxi(1,i,j,1:nz+1,direction) = rayfxi(1,1:nz+1)
fxi(2:nvars,i,j,1:nz+1,direction) = rayfxi(2:nvars,1:nz+1)/2.0d0
! Add nonlinear part
rhow(:) = ux(4,i,j,:)
phi(:) = vx(2,i,j,:)
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(2,i,j,1:nz+1,direction) = fxi(2,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
! Add nonlinear part
phi(:) = vx(3,i,j,:)
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(3,i,j,1:nz+1,direction) = fxi(3,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
! Add other pressure half
phi(:) = vx(5,i,j,:)
CALL FDInterpolateZ(1, phi,skw,lol,hil,nz)
fxi(4,i,j,1:nz+1,direction) = fxi(4,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
! Add nonlinear part
phi(:) = vx(4,i,j,:)
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(4,i,j,1:nz+1,direction) = fxi(4,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
! Scalars
do is=6, nvars
phi(:) = vx(is,i,j,:)
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(is,i,j,1:nz+1,direction) = fxi(is,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
enddo
! Handling of the energy equation
! compute internal energy
phi(:) = (ux(5,i,j,:)-0.5d0*(ux(2,i,j,:)**2+ux(3,i,j,:)**2 &
+ux(4,i,j,:)**2)/ux(1,i,j,:))/ux(1,i,j,:)
! skew half
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(5,i,j,1:nz+1,direction) = skw(1:nz+1)/2.0d0
! conservative half
phi(:) = phi(:)*rhow(:)
CALL FDInterpolateZ(1, phi,skw,lol,hil,nz)
fxi(5,i,j,1:nz+1,direction) = fxi(5,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
! Quadratic term
rhow(:) = vx(2,i,j,:)*ux(4,i,j,:)
phi(:) = vx(2,i,j,:)
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(5,i,j,1:nz+1,direction) = fxi(5,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
rhow(:) = vx(3,i,j,:)*ux(4,i,j,:)
phi(:) = vx(3,i,j,:)
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(5,i,j,1:nz+1,direction) = fxi(5,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
rhow(:) = vx(4,i,j,:)*ux(4,i,j,:)
phi(:) = vx(4,i,j,:)
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(5,i,j,1:nz+1,direction) = fxi(5,i,j,1:nz+1,direction) &
+ skw(1:nz+1)/2.0d0
! Add pressure-velocity non-conservative part
phi(:) = vx(4,i,j,:)
rhow(:) = vx(5,i,j,:)
CALL FDSkewInterpolateZ(rhow, phi, skw,lol,hil,nz)
fxi(5,i,j,1:nz+1,direction) = fxi(5,i,j,1:nz+1,direction) &
+ skw(1:nz+1)
END DO
END DO
DEALLOCATE ( rayfx, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfx', iret)
endif
DEALLOCATE ( rayfxi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rayfxi', iret)
endif
DEALLOCATE ( skw, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/skw', iret)
endif
DEALLOCATE ( rhow, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/rhow', iret)
endif
DEALLOCATE ( phi, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'ThreeDFDInterpFluxSkew/phi', iret)
endif
END SELECT whichdirection
call cleslog_log_exit('ThreeDFDInterpFluxSkew')
END SUBROUTINE ThreeDFDInterpFluxSkew
END Module Generic_FDInterpFlux