! --- This holds generic interpolations from cell centers
! --- to cell walls.
! Optionally uses the correct one-sided differences
! --- to impose boundary conditions
MODULE Generic_FDInterp
CONTAINS
SUBROUTINE FDInterpolateX(n,fx,fxi,lol,hil,np)
USE Interp_coeffs
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
INTEGER, INTENT(IN) :: lol ! ---- local low bound
INTEGER, INTENT(IN) :: hil ! ---- local high bound
INTEGER, INTENT(IN) :: np ! ---- local number of points
! ---- to be returned in interpolated flux
DOUBLE PRECISION, INTENT(IN) :: fx(n,lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(OUT) :: fxi(n,1:np+1) ! --- interpolated flux
CALL FDInterpolate(n,fx,fxi,lol,hil,np,bc_ixlow,bc_ixup)
RETURN
END SUBROUTINE FDInterpolateX
SUBROUTINE FDInterpolateY(n,fx,fxi,lol,hil,np)
USE Interp_coeffs
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
INTEGER, INTENT(IN) :: lol ! ---- local low bound
INTEGER, INTENT(IN) :: hil ! ---- local high bound
INTEGER, INTENT(IN) :: np ! ---- local number of points
! ---- to be returned in interpolated flux
DOUBLE PRECISION, INTENT(IN) :: fx(n,lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(OUT) :: fxi(n,1:np+1) ! --- interpolated flux
CALL FDInterpolate(n,fx,fxi,lol,hil,np,bc_iylow,bc_iyup)
RETURN
END SUBROUTINE FDInterpolateY
SUBROUTINE FDInterpolateZ(n,fx,fxi,lol,hil,np)
USE Interp_coeffs
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
INTEGER, INTENT(IN) :: lol ! ---- local low bound
INTEGER, INTENT(IN) :: hil ! ---- local high bound
INTEGER, INTENT(IN) :: np ! ---- local number of points
! ---- to be returned in interpolated flux
DOUBLE PRECISION, INTENT(IN) :: fx(n,lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(OUT) :: fxi(n,1:np+1) ! --- interpolated flux
CALL FDInterpolate(n,fx,fxi,lol,hil,np,bc_izlow,bc_izup)
RETURN
END SUBROUTINE FDInterpolateZ
SUBROUTINE FDInterpolate(n, fx,fxi,lol,hil,np,ixlow,ixup)
! ---- Shared Variables
USE Interp_coeffs
! ----
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: n
INTEGER, INTENT(IN) :: lol ! ---- local low bound
INTEGER, INTENT(IN) :: hil ! ---- local high bound
INTEGER, INTENT(IN) :: np ! ---- local number of points
! ---- to be returned in interpolated flux
INTEGER, INTENT(IN) :: ixlow
INTEGER, INTENT(IN) :: ixup ! Which biased stensil to impose
DOUBLE PRECISION, INTENT(IN) :: fx(n,lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(OUT) :: fxi(n,1:np+1) ! --- interpolated flux
INTEGER :: ip, is, is_loc, ip_loc
DOUBLE PRECISION :: fxac(n)
! ---- fxi(:,i) hold the half point interpolation flux(:,i-1/2)
! ---- This does the interpolations consistant with a
! ---- 3-point, 5-Point or 7-point center difference.
! ---- the values for cd1, cd2, cd3 are stored in 'Interp_coeffs'
! ---- and are initialized when all the weno coeffs are.
fxi(:,1:np+1) = tcd_flux(1) * (fx(:,0:np) + fx(:,1:np+1))
if ( stencil .gt. 3 ) then
fxi(:,1:np+1) = fxi(:,1:np+1) + &
tcd_flux(2) * (fx(:,2:np+2) + fx(:,-1:np-1))
endif
if (stencil.gt.5) then
fxi(:,1:np+1) = fxi(:,1:np+1) + &
tcd_flux(3) * (fx(:,3:np+3)+fx(:,-2:np-2) )
endif
! Now correct fluxes for boundary conditions
if ( ixlow .ne. CLES_PATCH_CORE ) then
! store the last centered flux term
is_loc = tcd_bndry_width+2-ixlow
if ( is_loc .le. np+1 ) then
fxac(:) = fxi(:,is_loc)
else
fxac(:) = tcd_flux(1) * (fx(:,is_loc-1) + fx(:,is_loc))
if ( stencil .gt. 3 ) then
fxac(:) = fxac(:) + tcd_flux(2) * (fx(:,is_loc+1) + fx(:,is_loc-2))
endif
if (stencil.gt.5) then
fxac(:) = fxac(:) + tcd_flux(3) * (fx(:,is_loc+2)+fx(:,is_loc-3) )
endif
endif
! telescope
do is=tcd_bndry_width, ixlow, -1
is_loc = is+1-ixlow
do ip=1, tcd_bndry_points
ip_loc = ip+1-ixlow
fxac(:) = fxac(:) - tcd_bndry(ip,is)*fx(:,ip_loc)
enddo
if ( is_loc .ge. 1 .and. is_loc .le. np+1 ) then
fxi(:,is_loc) = fxac(:)
endif
enddo
endif
if ( ixup .ne. CLES_PATCH_CORE ) then
! store the last centered flux term
is_loc = np-tcd_bndry_width+ixup
if ( is_loc .ge. 1 ) then
fxac(:) = fxi(:,is_loc)
else
fxac(:) = tcd_flux(1) * (fx(:,is_loc-1) + fx(:,is_loc))
if ( stencil .gt. 3 ) then
fxac(:) = fxac(:) + tcd_flux(2) * (fx(:,is_loc+1) + fx(:,is_loc-2))
endif
if (stencil.gt.5) then
fxac(:) = fxac(:) + tcd_flux(3) * (fx(:,is_loc+2)+fx(:,is_loc-3) )
endif
endif
! telescope
do is=tcd_bndry_width, ixup, -1
is_loc = np-is+ixup
do ip=1, tcd_bndry_points
ip_loc = np-ip+ixup
fxac(:) = fxac(:) - tcd_bndry(ip,is)*fx(:,ip_loc)
enddo
if ( is_loc .ge. 0 .and. is_loc .le. np ) then
fxi(:,is_loc+1) = fxac(:)
endif
enddo
endif
END SUBROUTINE FDInterpolate
subroutine FDInterpolateTest(lol, hil, np, ixlow, ixup)
! ---- Shared Variables
USE Interp_coeffs
implicit none
include 'cles.i'
integer :: lol, hil, np, ixlow, ixup
integer is_loc, width
sel: select case (stencil)
case (3)
if ( lol .gt. 0 .or. hil .lt. np + 1 ) then
print *, 'Core stencil ghost size error'
stop
endif
if ( ixlow .ne. CLES_PATCH_CORE ) then
is_loc = tcd_bndry_width+2-ixlow
if ( is_loc .le. np+1 ) then
if ( is_loc .lt. 1 ) then
print *, 'Boundary stencil ghost-size error (flux)'
stop
endif
else
if ( is_loc-1 .lt. lol .or. is_loc .gt. hil ) then
print *, 'Boundary stencil ghost-size error (der)'
stop
endif
endif
endif
if ( ixup .ne. CLES_PATCH_CORE ) then
is_loc = np-tcd_bndry_width+ixup
if ( is_loc .ge. 1 ) then
if ( is_loc .gt. np+1 ) then
print *, 'Boundary stencil ghost-size error (flux)'
stop
endif
else
if ( is_loc-1 .lt. lol .or. is_loc .gt. hil ) then
print *, 'Boundary stencil ghost-size error (der)'
stop
endif
endif
endif
case (5)
if ( lol .gt. -1 .or. hil .lt. np + 2 ) then
print *, 'Core stencil ghost size error'
stop
endif
if ( ixlow .ne. CLES_PATCH_CORE ) then
is_loc = tcd_bndry_width+2-ixlow
if ( is_loc .le. np+1 ) then
if ( is_loc .lt. 1 ) then
print *, 'Boundary stencil ghost-size error (flux)'
stop
endif
else
if ( is_loc-2 .lt. lol .or. is_loc+1 .gt. hil ) then
print *, 'Boundary stencil ghost-size error (der)'
stop
endif
endif
endif
if ( ixup .ne. CLES_PATCH_CORE ) then
is_loc = np-tcd_bndry_width+ixup
if ( is_loc .ge. 1 ) then
if ( is_loc .gt. np+1 ) then
print *, 'Boundary stencil ghost-size error (flux)'
stop
endif
else
if ( is_loc-2 .lt. lol .or. is_loc+1 .gt. hil ) then
print *, 'Boundary stencil ghost-size error (der)'
stop
endif
endif
endif
case (7)
if ( lol .gt. -2 .or. hil .lt. np + 3 ) then
print *, 'Core stencil ghost size error'
stop
endif
if ( ixlow .ne. CLES_PATCH_CORE ) then
is_loc = tcd_bndry_width+2-ixlow
if ( is_loc .le. np+1 ) then
if ( is_loc .lt. 1 ) then
print *, 'Boundary stencil ghost-size error (flux)'
stop
endif
else
if ( is_loc-3 .lt. lol .or. is_loc+2 .gt. hil ) then
print *, 'Boundary stencil ghost-size error (der)'
stop
endif
endif
endif
if ( ixup .ne. CLES_PATCH_CORE ) then
is_loc = np-tcd_bndry_width+ixup
if ( is_loc .ge. 1 ) then
if ( is_loc .gt. np+1 ) then
print *, 'Boundary stencil ghost-size error (flux)'
stop
endif
else
if ( is_loc-3 .lt. lol .or. is_loc+2 .gt. hil ) then
print *, 'Boundary stencil ghost-size error (der)'
stop
endif
endif
endif
end select sel
width = max(tcd_bndry_points,tcd_bndry_width)
if ( ixlow .ne. CLES_PATCH_CORE ) then
if ( 2-ixlow .lt. lol .or. width+1-ixlow .gt. hil ) then
print *, 'Boundary stencil ghost-size error (biased)'
stop
endif
endif
if ( ixup .ne. CLES_PATCH_CORE ) then
if ( np-width+ixup .lt. lol .or. np-1+ixup .gt. hil ) then
print *, 'Boundary stencil ghost-size error (biased)'
stop
endif
endif
end subroutine FDInterpolateTest
END Module Generic_FDInterp