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