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