! --- 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