! --- This holds interpolations from cell centers
! --- to cell walls used in product rule
! ---- d(ab)/dx = a d(b)/dx + b d(a)/dx.
! Optionally uses the correct one-sided differences
! --- to impose boundary conditions
Module Generic_FDSkewInterp
CONTAINS
SUBROUTINE FDSkewInterpolateX(a, b, skw,lol,hil,np)
! ---- Shared Variables
USE Interp_coeffs
! ----
IMPLICIT NONE
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) :: a(lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(IN) :: b(lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(OUT) :: skw(1:np+1) ! --- interpolated flux
CALL FDSkewInterpolate(a, b, skw,lol,hil,np,bc_ixlow,bc_ixup)
END SUBROUTINE FDSkewInterpolateX
SUBROUTINE FDSkewInterpolateY(a, b, skw,lol,hil,np)
! ---- Shared Variables
USE Interp_coeffs
! ----
IMPLICIT NONE
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) :: a(lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(IN) :: b(lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(OUT) :: skw(1:np+1) ! --- interpolated flux
CALL FDSkewInterpolate(a, b, skw,lol,hil,np,bc_iylow,bc_iyup)
END SUBROUTINE FDSkewInterpolateY
SUBROUTINE FDSkewInterpolateZ(a, b, skw,lol,hil,np)
! ---- Shared Variables
USE Interp_coeffs
! ----
IMPLICIT NONE
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) :: a(lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(IN) :: b(lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(OUT) :: skw(1:np+1) ! --- interpolated flux
CALL FDSkewInterpolate(a, b, skw,lol,hil,np,bc_izlow,bc_izup)
END SUBROUTINE FDSkewInterpolateZ
SUBROUTINE FDSkewInterpolate(a, b, skw,lol,hil,np,ixlow,ixup)
! ---- Shared Variables
USE Interp_coeffs
! ----
IMPLICIT NONE
include 'cles.i'
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) :: a(lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(IN) :: b(lol:hil) ! --- incoming flux
DOUBLE PRECISION, INTENT(OUT) :: skw(1:np+1) ! --- interpolated flux
INTEGER :: ip, is, is_loc, ip_loc
DOUBLE PRECISION :: fxac
skw(1:np+1) = tcd_center(1)*(a(0:np)*b(1:np+1) + a(1:np+1)*b(0:np))
if ( stencil .gt. 3 ) then
skw(1:np+1) = skw(1:np+1) + tcd_center(2)*(a(0:np)*b(2:np+2) &
+ a(1:np+1)*b(-1:np-1) + a(2:np+2)*b(0:np) &
+ a(-1:np-1)*b(1:np+1))
endif
if (stencil .gt. 5 ) then
skw(1:np+1) = skw(1:np+1) + tcd_center(3)*(a(0:np)*b(3:np+3)&
+ a(-1:np-1)*b(2:np+2) + a(-2:np-2)*b(1:np+1) &
+ a(3:np+3)*b(0:np) + a(2:np+2)*b(-1:np-1) + a(1:np+1)*b(-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 = skw(is_loc)
else
fxac = tcd_center(1)*(a(is_loc-1)*b(is_loc) + a(is_loc)*b(is_loc-1))
if ( stencil .gt. 3 ) then
fxac = fxac + tcd_center(2)*(a(is_loc-1)*b(is_loc+1) &
+ a(is_loc)*b(is_loc-2) + a(is_loc+1)*b(is_loc-1) &
+ a(is_loc-2)*b(is_loc))
endif
if (stencil .gt. 5 ) then
fxac = fxac + tcd_center(3)*(a(is_loc-1)*b(is_loc+2)&
+ a(is_loc-2)*b(is_loc+1) + a(is_loc-3)*b(is_loc) &
+ a(is_loc+2)*b(is_loc-1) + a(is_loc+1)*b(is_loc-2) &
+ a(is_loc)*b(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)*(a(ip_loc)*b(is_loc)+b(ip_loc)*a(is_loc))
enddo
if ( is_loc .ge. 1 .and. is_loc .le. np+1 ) then
skw(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 = skw(is_loc)
else
fxac = tcd_center(1)*(a(is_loc-1)*b(is_loc) + a(is_loc)*b(is_loc-1))
if ( stencil .gt. 3 ) then
fxac = fxac + tcd_center(2)*(a(is_loc-1)*b(is_loc+1) &
+ a(is_loc)*b(is_loc-2) + a(is_loc+1)*b(is_loc-1) &
+ a(is_loc-2)*b(is_loc))
endif
if (stencil .gt. 5 ) then
fxac = fxac + tcd_center(3)*(a(is_loc-1)*b(is_loc+2)&
+ a(is_loc-2)*b(is_loc+1) + a(is_loc-3)*b(is_loc) &
+ a(is_loc+2)*b(is_loc-1) + a(is_loc+1)*b(is_loc-2) &
+ a(is_loc)*b(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)*(a(ip_loc)*b(is_loc)+b(ip_loc)*a(is_loc))
enddo
if ( is_loc .ge. 0 .and. is_loc .le. np ) then
skw(is_loc+1) = fxac
endif
enddo
endif
RETURN
END SUBROUTINE FDSkewInterpolate
END MODULE Generic_FDSkewInterp