diff --git a/sources/schemes.F90 b/sources/schemes.F90
index 0a1b8bc..3e7e391 100644
--- a/sources/schemes.F90
+++ b/sources/schemes.F90
@@ -689,6 +689,1178 @@ module schemes
!
!===============================================================================
!
+! subroutine UPDATE_FLUX_HD_ISO:
+! -----------------------------
+!
+! Subroutine solves the Riemann problem along each direction and calculates
+! the numerical fluxes, which are used later to calculate the conserved
+! variable increment.
+!
+! Arguments:
+!
+! dx - the spatial step;
+! q - the array of primitive variables;
+! f - the array of numerical fluxes;
+!
+!===============================================================================
+!
+ subroutine update_flux_hd_iso(dx, q, f)
+
+! include external variables
+!
+ use coordinates, only : nn => bcells, nbl, neu
+ use equations , only : nf
+ use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
+
+! local variables are not implicit by default
+!
+ implicit none
+
+! input arguments
+!
+ real(kind=8), dimension(:) , intent(in) :: dx
+ real(kind=8), dimension(:,:,:,:) , intent(in) :: q
+ real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
+
+! local variables
+!
+ integer :: i, j, k = 1
+
+! local temporary arrays
+!
+#if NDIMS == 3
+ real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
+#else /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
+#endif /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,2) :: qi
+ real(kind=8), dimension(nf,nn) :: fi
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for flux update
+!
+ call start_timer(imf)
+#endif /* PROFILE */
+
+! initialize fluxes
+!
+ f(:,:,:,:,:) = 0.0d+00
+
+! reconstruct interfaces
+!
+ call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
+
+! calculate the flux along the X-direction
+!
+#if NDIMS == 3
+ do k = nbl, neu
+#endif /* NDIMS == 3 */
+ do j = nbl, neu
+
+! copy states to directional lines with proper vector component ordering
+!
+ qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
+ qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
+ qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
+ qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(1,idn,:,j,k) = fi(idn,:)
+ f(1,imx,:,j,k) = fi(imx,:)
+ f(1,imy,:,j,k) = fi(imy,:)
+ f(1,imz,:,j,k) = fi(imz,:)
+
+ end do ! j = nbl, neu
+
+! calculate the flux along the Y direction
+!
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
+ qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
+ qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
+ qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(2,idn,i,:,k) = fi(idn,:)
+ f(2,imx,i,:,k) = fi(imz,:)
+ f(2,imy,i,:,k) = fi(imx,:)
+ f(2,imz,i,:,k) = fi(imy,:)
+
+ end do ! i = nbl, neu
+#if NDIMS == 3
+ end do ! k = nbl, neu
+#endif /* NDIMS == 3 */
+
+#if NDIMS == 3
+! calculate the flux along the Z direction
+!
+ do j = nbl, neu
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
+ qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
+ qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
+ qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(3,idn,i,j,:) = fi(idn,:)
+ f(3,imx,i,j,:) = fi(imy,:)
+ f(3,imy,i,j,:) = fi(imz,:)
+ f(3,imz,i,j,:) = fi(imx,:)
+
+ end do ! i = nbl, neu
+ end do ! j = nbl, neu
+#endif /* NDIMS == 3 */
+
+#ifdef PROFILE
+! stop accounting time for flux update
+!
+ call stop_timer(imf)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine update_flux_hd_iso
+!
+!===============================================================================
+!
+! subroutine UPDATE_FLUX_HD_ADI:
+! -----------------------------
+!
+! Subroutine solves the Riemann problem along each direction and calculates
+! the numerical fluxes, which are used later to calculate the conserved
+! variable increment.
+!
+! Arguments:
+!
+! dx - the spatial step;
+! q - the array of primitive variables;
+! f - the array of numerical fluxes;
+!
+!===============================================================================
+!
+ subroutine update_flux_hd_adi(dx, q, f)
+
+! include external variables
+!
+ use coordinates, only : nn => bcells, nbl, neu
+ use equations , only : nf
+ use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
+
+! local variables are not implicit by default
+!
+ implicit none
+
+! input arguments
+!
+ real(kind=8), dimension(:) , intent(in) :: dx
+ real(kind=8), dimension(:,:,:,:) , intent(in) :: q
+ real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
+
+! local variables
+!
+ integer :: i, j, k = 1
+
+! local temporary arrays
+!
+#if NDIMS == 3
+ real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
+#else /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
+#endif /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,2) :: qi
+ real(kind=8), dimension(nf,nn) :: fi
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for flux update
+!
+ call start_timer(imf)
+#endif /* PROFILE */
+
+! initialize fluxes
+!
+ f(:,:,:,:,:) = 0.0d+00
+
+! reconstruct interfaces
+!
+ call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
+
+! calculate the flux along the X-direction
+!
+#if NDIMS == 3
+ do k = nbl, neu
+#endif /* NDIMS == 3 */
+ do j = nbl, neu
+
+! copy states to directional lines with proper vector component ordering
+!
+ qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
+ qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
+ qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
+ qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
+ qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(1,idn,:,j,k) = fi(idn,:)
+ f(1,imx,:,j,k) = fi(imx,:)
+ f(1,imy,:,j,k) = fi(imy,:)
+ f(1,imz,:,j,k) = fi(imz,:)
+ f(1,ien,:,j,k) = fi(ien,:)
+
+ end do ! j = nbl, neu
+
+! calculate the flux along the Y direction
+!
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
+ qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
+ qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
+ qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
+ qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(2,idn,i,:,k) = fi(idn,:)
+ f(2,imx,i,:,k) = fi(imz,:)
+ f(2,imy,i,:,k) = fi(imx,:)
+ f(2,imz,i,:,k) = fi(imy,:)
+ f(2,ien,i,:,k) = fi(ien,:)
+
+ end do ! i = nbl, neu
+#if NDIMS == 3
+ end do ! k = nbl, neu
+#endif /* NDIMS == 3 */
+
+#if NDIMS == 3
+! calculate the flux along the Z direction
+!
+ do j = nbl, neu
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
+ qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
+ qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
+ qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
+ qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(3,idn,i,j,:) = fi(idn,:)
+ f(3,imx,i,j,:) = fi(imy,:)
+ f(3,imy,i,j,:) = fi(imz,:)
+ f(3,imz,i,j,:) = fi(imx,:)
+ f(3,ien,i,j,:) = fi(ien,:)
+
+ end do ! i = nbl, neu
+ end do ! j = nbl, neu
+#endif /* NDIMS == 3 */
+
+#ifdef PROFILE
+! stop accounting time for flux update
+!
+ call stop_timer(imf)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine update_flux_hd_adi
+!
+!===============================================================================
+!
+! subroutine UPDATE_FLUX_MHD_ISO:
+! ------------------------------
+!
+! Subroutine solves the Riemann problem along each direction and calculates
+! the numerical fluxes, which are used later to calculate the conserved
+! variable increment.
+!
+! Arguments:
+!
+! dx - the spatial step;
+! q - the array of primitive variables;
+! f - the array of numerical fluxes;
+!
+!===============================================================================
+!
+ subroutine update_flux_mhd_iso(dx, q, f)
+
+! include external variables
+!
+ use coordinates, only : nn => bcells, nbl, neu
+ use equations , only : nf
+ use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
+ use equations , only : ibx, iby, ibz, ibp
+
+! local variables are not implicit by default
+!
+ implicit none
+
+! input arguments
+!
+ real(kind=8), dimension(:) , intent(in) :: dx
+ real(kind=8), dimension(:,:,:,:) , intent(in) :: q
+ real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
+
+! local variables
+!
+ integer :: i, j, k = 1
+
+! local temporary arrays
+!
+#if NDIMS == 3
+ real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
+#else /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
+#endif /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,2) :: qi
+ real(kind=8), dimension(nf,nn) :: fi
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for flux update
+!
+ call start_timer(imf)
+#endif /* PROFILE */
+
+! initialize fluxes
+!
+ f(:,:,:,:,:) = 0.0d+00
+
+! reconstruct interfaces
+!
+ call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
+
+! calculate the flux along the X-direction
+!
+#if NDIMS == 3
+ do k = nbl, neu
+#endif /* NDIMS == 3 */
+ do j = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
+ qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
+ qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
+ qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
+ qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
+ qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
+ qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
+ qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(1,idn,:,j,k) = fi(idn,:)
+ f(1,imx,:,j,k) = fi(imx,:)
+ f(1,imy,:,j,k) = fi(imy,:)
+ f(1,imz,:,j,k) = fi(imz,:)
+ f(1,ibx,:,j,k) = fi(ibx,:)
+ f(1,iby,:,j,k) = fi(iby,:)
+ f(1,ibz,:,j,k) = fi(ibz,:)
+ f(1,ibp,:,j,k) = fi(ibp,:)
+
+ end do ! j = nbl, neu
+
+! calculate the flux along the Y direction
+!
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
+ qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
+ qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
+ qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
+ qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
+ qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
+ qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
+ qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(2,idn,i,:,k) = fi(idn,:)
+ f(2,imx,i,:,k) = fi(imz,:)
+ f(2,imy,i,:,k) = fi(imx,:)
+ f(2,imz,i,:,k) = fi(imy,:)
+ f(2,ibx,i,:,k) = fi(ibz,:)
+ f(2,iby,i,:,k) = fi(ibx,:)
+ f(2,ibz,i,:,k) = fi(iby,:)
+ f(2,ibp,i,:,k) = fi(ibp,:)
+
+ end do ! i = nbl, neu
+#if NDIMS == 3
+ end do ! k = nbl, neu
+#endif /* NDIMS == 3 */
+
+#if NDIMS == 3
+! calculate the flux along the Z direction
+!
+ do j = nbl, neu
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
+ qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
+ qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
+ qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
+ qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
+ qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
+ qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
+ qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(3,idn,i,j,:) = fi(idn,:)
+ f(3,imx,i,j,:) = fi(imy,:)
+ f(3,imy,i,j,:) = fi(imz,:)
+ f(3,imz,i,j,:) = fi(imx,:)
+ f(3,ibx,i,j,:) = fi(iby,:)
+ f(3,iby,i,j,:) = fi(ibz,:)
+ f(3,ibz,i,j,:) = fi(ibx,:)
+ f(3,ibp,i,j,:) = fi(ibp,:)
+
+ end do ! i = nbl, neu
+ end do ! j = nbl, neu
+#endif /* NDIMS == 3 */
+
+#ifdef PROFILE
+! stop accounting time for flux update
+!
+ call stop_timer(imf)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine update_flux_mhd_iso
+!
+!===============================================================================
+!
+! subroutine UPDATE_FLUX_MHD_ADI:
+! ------------------------------
+!
+! Subroutine solves the Riemann problem along each direction and calculates
+! the numerical fluxes, which are used later to calculate the conserved
+! variable increment.
+!
+! Arguments:
+!
+! dx - the spatial step;
+! q - the array of primitive variables;
+! f - the array of numerical fluxes;
+!
+!===============================================================================
+!
+ subroutine update_flux_mhd_adi(dx, q, f)
+
+! include external variables
+!
+ use coordinates, only : nn => bcells, nbl, neu
+ use equations , only : nf
+ use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
+ use equations , only : ibx, iby, ibz, ibp
+
+! local variables are not implicit by default
+!
+ implicit none
+
+! input arguments
+!
+ real(kind=8), dimension(:) , intent(in) :: dx
+ real(kind=8), dimension(:,:,:,:) , intent(in) :: q
+ real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
+
+! local variables
+!
+ integer :: i, j, k = 1
+
+! local temporary arrays
+!
+#if NDIMS == 3
+ real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
+#else /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
+#endif /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,2) :: qi
+ real(kind=8), dimension(nf,nn) :: fi
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for flux update
+!
+ call start_timer(imf)
+#endif /* PROFILE */
+
+! initialize fluxes
+!
+ f(:,:,:,:,:) = 0.0d+00
+
+! reconstruct interfaces
+!
+ call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
+
+! calculate the flux along the X-direction
+!
+#if NDIMS == 3
+ do k = nbl, neu
+#endif /* NDIMS == 3 */
+ do j = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
+ qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
+ qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
+ qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
+ qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
+ qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
+ qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
+ qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
+ qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(1,idn,:,j,k) = fi(idn,:)
+ f(1,imx,:,j,k) = fi(imx,:)
+ f(1,imy,:,j,k) = fi(imy,:)
+ f(1,imz,:,j,k) = fi(imz,:)
+ f(1,ibx,:,j,k) = fi(ibx,:)
+ f(1,iby,:,j,k) = fi(iby,:)
+ f(1,ibz,:,j,k) = fi(ibz,:)
+ f(1,ibp,:,j,k) = fi(ibp,:)
+ f(1,ien,:,j,k) = fi(ien,:)
+
+ end do ! j = nbl, neu
+
+! calculate the flux along the Y direction
+!
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
+ qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
+ qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
+ qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
+ qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
+ qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
+ qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
+ qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
+ qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(2,idn,i,:,k) = fi(idn,:)
+ f(2,imx,i,:,k) = fi(imz,:)
+ f(2,imy,i,:,k) = fi(imx,:)
+ f(2,imz,i,:,k) = fi(imy,:)
+ f(2,ibx,i,:,k) = fi(ibz,:)
+ f(2,iby,i,:,k) = fi(ibx,:)
+ f(2,ibz,i,:,k) = fi(iby,:)
+ f(2,ibp,i,:,k) = fi(ibp,:)
+ f(2,ien,i,:,k) = fi(ien,:)
+
+ end do ! i = nbl, neu
+#if NDIMS == 3
+ end do ! k = nbl, neu
+#endif /* NDIMS == 3 */
+
+#if NDIMS == 3
+! calculate the flux along the Z direction
+!
+ do j = nbl, neu
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
+ qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
+ qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
+ qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
+ qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
+ qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
+ qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
+ qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
+ qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(3,idn,i,j,:) = fi(idn,:)
+ f(3,imx,i,j,:) = fi(imy,:)
+ f(3,imy,i,j,:) = fi(imz,:)
+ f(3,imz,i,j,:) = fi(imx,:)
+ f(3,ibx,i,j,:) = fi(iby,:)
+ f(3,iby,i,j,:) = fi(ibz,:)
+ f(3,ibz,i,j,:) = fi(ibx,:)
+ f(3,ibp,i,j,:) = fi(ibp,:)
+ f(3,ien,i,j,:) = fi(ien,:)
+
+ end do ! i = nbl, neu
+ end do ! j = nbl, neu
+#endif /* NDIMS == 3 */
+
+#ifdef PROFILE
+! stop accounting time for flux update
+!
+ call stop_timer(imf)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine update_flux_mhd_adi
+!
+!===============================================================================
+!
+! subroutine UPDATE_FLUX_SRHD_ADI:
+! -------------------------------
+!
+! Subroutine solves the Riemann problem along each direction and calculates
+! the numerical fluxes, which are used later to calculate the conserved
+! variable increment.
+!
+! Arguments:
+!
+! dx - the spatial step;
+! q - the array of primitive variables;
+! f - the array of numerical fluxes;
+!
+!===============================================================================
+!
+ subroutine update_flux_srhd_adi(dx, q, f)
+
+! include external variables
+!
+ use coordinates, only : nn => bcells, nbl, neu
+ use equations , only : nf
+ use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
+
+! local variables are not implicit by default
+!
+ implicit none
+
+! input arguments
+!
+ real(kind=8), dimension(:) , intent(in) :: dx
+ real(kind=8), dimension(:,:,:,:) , intent(in) :: q
+ real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
+
+! local variables
+!
+ integer :: i, j, k = 1, l, p
+ real(kind=8) :: vm
+
+! local temporary arrays
+!
+#if NDIMS == 3
+ real(kind=8), dimension(nf,nn,nn,nn) :: qq
+ real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
+#else /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,nn, 1) :: qq
+ real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
+#endif /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,2) :: qi
+ real(kind=8), dimension(nf,nn) :: fi
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for flux update
+!
+ call start_timer(imf)
+#endif /* PROFILE */
+
+! initialize fluxes
+!
+ f(:,:,:,:,:) = 0.0d+00
+
+! apply reconstruction on variables using 4-vector or 3-vector
+!
+ if (states_4vec) then
+
+! convert velocities to four-velocities for physical reconstruction
+!
+#if NDIMS == 3
+ do k = 1, nn
+#endif /* NDIMS == 3 */
+ do j = 1, nn
+ do i = 1, nn
+
+ vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))
+
+ qq(idn,i,j,k) = q(idn,i,j,k) / vm
+ qq(ivx,i,j,k) = q(ivx,i,j,k) / vm
+ qq(ivy,i,j,k) = q(ivy,i,j,k) / vm
+ qq(ivz,i,j,k) = q(ivz,i,j,k) / vm
+ qq(ipr,i,j,k) = q(ipr,i,j,k)
+
+ end do ! i = 1, nn
+ end do ! j = 1, nn
+#if NDIMS == 3
+ end do ! k = 1, nn
+#endif /* NDIMS == 3 */
+
+! reconstruct interfaces
+!
+ call reconstruct_interfaces(dx(:), qq(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
+
+! convert state four-velocities back to velocities
+!
+#if NDIMS == 3
+ do k = 1, nn
+#endif /* NDIMS == 3 */
+ do j = 1, nn
+ do i = 1, nn
+
+ do l = 1, 2
+ do p = 1, NDIMS
+
+ vm = sqrt(1.0d+00 + sum(qs(ivx:ivz,i,j,k,l,p)**2))
+
+ qs(idn,i,j,k,l,p) = qs(idn,i,j,k,l,p) / vm
+ qs(ivx,i,j,k,l,p) = qs(ivx,i,j,k,l,p) / vm
+ qs(ivy,i,j,k,l,p) = qs(ivy,i,j,k,l,p) / vm
+ qs(ivz,i,j,k,l,p) = qs(ivz,i,j,k,l,p) / vm
+
+ end do ! p = 1, ndims
+ end do ! l = 1, 2
+ end do ! i = 1, nn
+ end do ! j = 1, nn
+#if NDIMS == 3
+ end do ! k = 1, nn
+#endif /* NDIMS == 3 */
+
+ else
+
+! reconstruct interfaces
+!
+ call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
+
+ end if
+
+! calculate the flux along the X-direction
+!
+#if NDIMS == 3
+ do k = nbl, neu
+#endif /* NDIMS == 3 */
+ do j = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
+ qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
+ qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
+ qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
+ qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(1,idn,:,j,k) = fi(idn,:)
+ f(1,imx,:,j,k) = fi(imx,:)
+ f(1,imy,:,j,k) = fi(imy,:)
+ f(1,imz,:,j,k) = fi(imz,:)
+ f(1,ien,:,j,k) = fi(ien,:)
+
+ end do ! j = nbl, neu
+
+! calculate the flux along the Y direction
+!
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
+ qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
+ qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
+ qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
+ qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(2,idn,i,:,k) = fi(idn,:)
+ f(2,imx,i,:,k) = fi(imz,:)
+ f(2,imy,i,:,k) = fi(imx,:)
+ f(2,imz,i,:,k) = fi(imy,:)
+ f(2,ien,i,:,k) = fi(ien,:)
+
+ end do ! i = nbl, neu
+#if NDIMS == 3
+ end do ! k = nbl, neu
+#endif /* NDIMS == 3 */
+
+#if NDIMS == 3
+! calculate the flux along the Z direction
+!
+ do j = nbl, neu
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
+ qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
+ qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
+ qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
+ qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(3,idn,i,j,:) = fi(idn,:)
+ f(3,imx,i,j,:) = fi(imy,:)
+ f(3,imy,i,j,:) = fi(imz,:)
+ f(3,imz,i,j,:) = fi(imx,:)
+ f(3,ien,i,j,:) = fi(ien,:)
+
+ end do ! i = nbl, neu
+ end do ! j = nbl, neu
+#endif /* NDIMS == 3 */
+
+#ifdef PROFILE
+! stop accounting time for flux update
+!
+ call stop_timer(imf)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine update_flux_srhd_adi
+!
+!===============================================================================
+!
+! subroutine UPDATE_FLUX_SRMHD_ADI:
+! --------------------------------
+!
+! Subroutine solves the Riemann problem along each direction and calculates
+! the numerical fluxes, which are used later to calculate the conserved
+! variable increment.
+!
+! Arguments:
+!
+! dx - the spatial step;
+! q - the array of primitive variables;
+! f - the array of numerical fluxes;
+!
+!===============================================================================
+!
+ subroutine update_flux_srmhd_adi(dx, q, f)
+
+! include external variables
+!
+ use coordinates, only : nn => bcells, nbl, neu
+ use equations , only : nf
+ use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
+ use equations , only : ibx, iby, ibz, ibp
+
+! local variables are not implicit by default
+!
+ implicit none
+
+! input arguments
+!
+ real(kind=8), dimension(:) , intent(in) :: dx
+ real(kind=8), dimension(:,:,:,:) , intent(in) :: q
+ real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
+
+! local variables
+!
+ integer :: i, j, k = 1, l, p
+ real(kind=8) :: vm
+
+! local temporary arrays
+!
+#if NDIMS == 3
+ real(kind=8), dimension(nf,nn,nn,nn) :: qq
+ real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
+#else /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,nn, 1) :: qq
+ real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
+#endif /* NDIMS == 3 */
+ real(kind=8), dimension(nf,nn,2) :: qi
+ real(kind=8), dimension(nf,nn) :: fi
+!
+!-------------------------------------------------------------------------------
+!
+#ifdef PROFILE
+! start accounting time for flux update
+!
+ call start_timer(imf)
+#endif /* PROFILE */
+
+! initialize fluxes
+!
+ f(:,:,:,:,:) = 0.0d+00
+
+! apply reconstruction on variables using 4-vector or 3-vector
+!
+ if (states_4vec) then
+
+! convert velocities to four-velocities for physical reconstruction
+!
+#if NDIMS == 3
+ do k = 1, nn
+#endif /* NDIMS == 3 */
+ do j = 1, nn
+ do i = 1, nn
+
+ vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))
+
+ qq(idn,i,j,k) = q(idn,i,j,k) / vm
+ qq(ivx,i,j,k) = q(ivx,i,j,k) / vm
+ qq(ivy,i,j,k) = q(ivy,i,j,k) / vm
+ qq(ivz,i,j,k) = q(ivz,i,j,k) / vm
+ qq(ipr,i,j,k) = q(ipr,i,j,k)
+ qq(ibx,i,j,k) = q(ibx,i,j,k)
+ qq(iby,i,j,k) = q(iby,i,j,k)
+ qq(ibz,i,j,k) = q(ibz,i,j,k)
+ qq(ibp,i,j,k) = q(ibp,i,j,k)
+ qq(ipr,i,j,k) = q(ipr,i,j,k)
+
+ end do ! i = 1, nn
+ end do ! j = 1, nn
+#if NDIMS == 3
+ end do ! k = 1, nn
+#endif /* NDIMS == 3 */
+
+! reconstruct interfaces
+!
+ call reconstruct_interfaces(dx(:), qq(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
+
+! convert state four-velocities back to velocities
+!
+#if NDIMS == 3
+ do k = 1, nn
+#endif /* NDIMS == 3 */
+ do j = 1, nn
+ do i = 1, nn
+
+ do l = 1, 2
+ do p = 1, NDIMS
+
+ vm = sqrt(1.0d+00 + sum(qs(ivx:ivz,i,j,k,l,p)**2))
+
+ qs(idn,i,j,k,l,p) = qs(idn,i,j,k,l,p) / vm
+ qs(ivx,i,j,k,l,p) = qs(ivx,i,j,k,l,p) / vm
+ qs(ivy,i,j,k,l,p) = qs(ivy,i,j,k,l,p) / vm
+ qs(ivz,i,j,k,l,p) = qs(ivz,i,j,k,l,p) / vm
+
+ end do ! p = 1, ndims
+ end do ! l = 1, 2
+ end do ! i = 1, nn
+ end do ! j = 1, nn
+#if NDIMS == 3
+ end do ! k = 1, nn
+#endif /* NDIMS == 3 */
+
+ else
+
+! reconstruct interfaces
+!
+ call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
+
+ end if
+
+! calculate the flux along the X-direction
+!
+#if NDIMS == 3
+ do k = nbl, neu
+#endif /* NDIMS == 3 */
+ do j = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
+ qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
+ qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
+ qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
+ qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
+ qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
+ qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
+ qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
+ qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(1,idn,:,j,k) = fi(idn,:)
+ f(1,imx,:,j,k) = fi(imx,:)
+ f(1,imy,:,j,k) = fi(imy,:)
+ f(1,imz,:,j,k) = fi(imz,:)
+ f(1,ibx,:,j,k) = fi(ibx,:)
+ f(1,iby,:,j,k) = fi(iby,:)
+ f(1,ibz,:,j,k) = fi(ibz,:)
+ f(1,ibp,:,j,k) = fi(ibp,:)
+ f(1,ien,:,j,k) = fi(ien,:)
+
+ end do ! j = nbl, neu
+
+! calculate the flux along the Y direction
+!
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
+ qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
+ qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
+ qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
+ qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
+ qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
+ qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
+ qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
+ qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(2,idn,i,:,k) = fi(idn,:)
+ f(2,imx,i,:,k) = fi(imz,:)
+ f(2,imy,i,:,k) = fi(imx,:)
+ f(2,imz,i,:,k) = fi(imy,:)
+ f(2,ibx,i,:,k) = fi(ibz,:)
+ f(2,iby,i,:,k) = fi(ibx,:)
+ f(2,ibz,i,:,k) = fi(iby,:)
+ f(2,ibp,i,:,k) = fi(ibp,:)
+ f(2,ien,i,:,k) = fi(ien,:)
+
+ end do ! i = nbl, neu
+#if NDIMS == 3
+ end do ! k = nbl, neu
+#endif /* NDIMS == 3 */
+
+#if NDIMS == 3
+! calculate the flux along the Z direction
+!
+ do j = nbl, neu
+ do i = nbl, neu
+
+! copy directional variable vectors to pass to the one dimensional solver
+!
+ qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
+ qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
+ qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
+ qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
+ qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
+ qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
+ qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
+ qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
+ qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
+
+! call one dimensional Riemann solver in order to obtain numerical fluxes
+!
+ call numerical_flux(qi(:,:,:), fi(:,:))
+
+! update the array of fluxes
+!
+ f(3,idn,i,j,:) = fi(idn,:)
+ f(3,imx,i,j,:) = fi(imy,:)
+ f(3,imy,i,j,:) = fi(imz,:)
+ f(3,imz,i,j,:) = fi(imx,:)
+ f(3,ibx,i,j,:) = fi(iby,:)
+ f(3,iby,i,j,:) = fi(ibz,:)
+ f(3,ibz,i,j,:) = fi(ibx,:)
+ f(3,ibp,i,j,:) = fi(ibp,:)
+ f(3,ien,i,j,:) = fi(ien,:)
+
+ end do ! i = nbl, neu
+ end do ! j = nbl, neu
+#endif /* NDIMS == 3 */
+
+#ifdef PROFILE
+! stop accounting time for flux update
+!
+ call stop_timer(imf)
+#endif /* PROFILE */
+
+!-------------------------------------------------------------------------------
+!
+ end subroutine update_flux_srmhd_adi
+!
+!===============================================================================
+!
! subroutine RECONSTRUCT_INTERFACES:
! ---------------------------------
!
@@ -1019,7 +2191,7 @@ module schemes
ql(ibp,i) = bp
qr(ibp,i) = bp
- end do ! i = 1, n
+ end do
! calculate corresponding conserved variables of the left and right states
!
@@ -1621,167 +2793,6 @@ module schemes
!
!===============================================================================
!
-!***** ISOTHERMAL HYDRODYNAMICS *****
-!
-!===============================================================================
-!
-! subroutine UPDATE_FLUX_HD_ISO:
-! -----------------------------
-!
-! Subroutine solves the Riemann problem along each direction and calculates
-! the numerical fluxes, which are used later to calculate the conserved
-! variable increment.
-!
-! Arguments:
-!
-! dx - the spatial step;
-! q - the array of primitive variables;
-! f - the array of numerical fluxes;
-!
-!===============================================================================
-!
- subroutine update_flux_hd_iso(dx, q, f)
-
-! include external variables
-!
- use coordinates, only : nn => bcells, nbl, neu
- use equations , only : nf
- use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
-
-! local variables are not implicit by default
-!
- implicit none
-
-! input arguments
-!
- real(kind=8), dimension(:) , intent(in) :: dx
- real(kind=8), dimension(:,:,:,:) , intent(in) :: q
- real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
-
-! local variables
-!
- integer :: i, j, k = 1
-
-! local temporary arrays
-!
-#if NDIMS == 3
- real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
-#else /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
-#endif /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,2) :: qi
- real(kind=8), dimension(nf,nn) :: fi
-!
-!-------------------------------------------------------------------------------
-!
-#ifdef PROFILE
-! start accounting time for flux update
-!
- call start_timer(imf)
-#endif /* PROFILE */
-
-! initialize fluxes
-!
- f(:,:,:,:,:) = 0.0d+00
-
-! reconstruct interfaces
-!
- call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
-
-! calculate the flux along the X-direction
-!
-#if NDIMS == 3
- do k = nbl, neu
-#endif /* NDIMS == 3 */
- do j = nbl, neu
-
-! copy states to directional lines with proper vector component ordering
-!
- qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
- qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
- qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
- qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(1,idn,:,j,k) = fi(idn,:)
- f(1,imx,:,j,k) = fi(imx,:)
- f(1,imy,:,j,k) = fi(imy,:)
- f(1,imz,:,j,k) = fi(imz,:)
-
- end do ! j = nbl, neu
-
-! calculate the flux along the Y direction
-!
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
- qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
- qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
- qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(2,idn,i,:,k) = fi(idn,:)
- f(2,imx,i,:,k) = fi(imz,:)
- f(2,imy,i,:,k) = fi(imx,:)
- f(2,imz,i,:,k) = fi(imy,:)
-
- end do ! i = nbl, neu
-#if NDIMS == 3
- end do ! k = nbl, neu
-#endif /* NDIMS == 3 */
-
-#if NDIMS == 3
-! calculate the flux along the Z direction
-!
- do j = nbl, neu
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
- qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
- qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
- qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(3,idn,i,j,:) = fi(idn,:)
- f(3,imx,i,j,:) = fi(imy,:)
- f(3,imy,i,j,:) = fi(imz,:)
- f(3,imz,i,j,:) = fi(imx,:)
-
- end do ! i = nbl, neu
- end do ! j = nbl, neu
-#endif /* NDIMS == 3 */
-
-#ifdef PROFILE
-! stop accounting time for flux update
-!
- call stop_timer(imf)
-#endif /* PROFILE */
-
-!-------------------------------------------------------------------------------
-!
- end subroutine update_flux_hd_iso
-!
-!===============================================================================
-!
! subroutine RIEMANN_HD_ISO_ROE:
! -----------------------------
!
@@ -1936,173 +2947,6 @@ module schemes
!
!===============================================================================
!
-!***** ADIABATIC HYDRODYNAMICS *****
-!
-!===============================================================================
-!
-! subroutine UPDATE_FLUX_HD_ADI:
-! -----------------------------
-!
-! Subroutine solves the Riemann problem along each direction and calculates
-! the numerical fluxes, which are used later to calculate the conserved
-! variable increment.
-!
-! Arguments:
-!
-! dx - the spatial step;
-! q - the array of primitive variables;
-! f - the array of numerical fluxes;
-!
-!===============================================================================
-!
- subroutine update_flux_hd_adi(dx, q, f)
-
-! include external variables
-!
- use coordinates, only : nn => bcells, nbl, neu
- use equations , only : nf
- use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
-
-! local variables are not implicit by default
-!
- implicit none
-
-! input arguments
-!
- real(kind=8), dimension(:) , intent(in) :: dx
- real(kind=8), dimension(:,:,:,:) , intent(in) :: q
- real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
-
-! local variables
-!
- integer :: i, j, k = 1
-
-! local temporary arrays
-!
-#if NDIMS == 3
- real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
-#else /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
-#endif /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,2) :: qi
- real(kind=8), dimension(nf,nn) :: fi
-!
-!-------------------------------------------------------------------------------
-!
-#ifdef PROFILE
-! start accounting time for flux update
-!
- call start_timer(imf)
-#endif /* PROFILE */
-
-! initialize fluxes
-!
- f(:,:,:,:,:) = 0.0d+00
-
-! reconstruct interfaces
-!
- call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
-
-! calculate the flux along the X-direction
-!
-#if NDIMS == 3
- do k = nbl, neu
-#endif /* NDIMS == 3 */
- do j = nbl, neu
-
-! copy states to directional lines with proper vector component ordering
-!
- qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
- qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
- qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
- qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
- qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(1,idn,:,j,k) = fi(idn,:)
- f(1,imx,:,j,k) = fi(imx,:)
- f(1,imy,:,j,k) = fi(imy,:)
- f(1,imz,:,j,k) = fi(imz,:)
- f(1,ien,:,j,k) = fi(ien,:)
-
- end do ! j = nbl, neu
-
-! calculate the flux along the Y direction
-!
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
- qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
- qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
- qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
- qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(2,idn,i,:,k) = fi(idn,:)
- f(2,imx,i,:,k) = fi(imz,:)
- f(2,imy,i,:,k) = fi(imx,:)
- f(2,imz,i,:,k) = fi(imy,:)
- f(2,ien,i,:,k) = fi(ien,:)
-
- end do ! i = nbl, neu
-#if NDIMS == 3
- end do ! k = nbl, neu
-#endif /* NDIMS == 3 */
-
-#if NDIMS == 3
-! calculate the flux along the Z direction
-!
- do j = nbl, neu
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
- qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
- qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
- qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
- qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(3,idn,i,j,:) = fi(idn,:)
- f(3,imx,i,j,:) = fi(imy,:)
- f(3,imy,i,j,:) = fi(imz,:)
- f(3,imz,i,j,:) = fi(imx,:)
- f(3,ien,i,j,:) = fi(ien,:)
-
- end do ! i = nbl, neu
- end do ! j = nbl, neu
-#endif /* NDIMS == 3 */
-
-#ifdef PROFILE
-! stop accounting time for flux update
-!
- call stop_timer(imf)
-#endif /* PROFILE */
-
-!-------------------------------------------------------------------------------
-!
- end subroutine update_flux_hd_adi
-!
-!===============================================================================
-!
! subroutine RIEMANN_HD_ADI_ROE:
! -----------------------------
!
@@ -2259,192 +3103,6 @@ module schemes
!
!===============================================================================
!
-!***** ISOTHERMAL MAGNETOHYDRODYNAMICS *****
-!
-!===============================================================================
-!
-! subroutine UPDATE_FLUX_MHD_ISO:
-! ------------------------------
-!
-! Subroutine solves the Riemann problem along each direction and calculates
-! the numerical fluxes, which are used later to calculate the conserved
-! variable increment.
-!
-! Arguments:
-!
-! dx - the spatial step;
-! q - the array of primitive variables;
-! f - the array of numerical fluxes;
-!
-!===============================================================================
-!
- subroutine update_flux_mhd_iso(dx, q, f)
-
-! include external variables
-!
- use coordinates, only : nn => bcells, nbl, neu
- use equations , only : nf
- use equations , only : idn, ivx, ivy, ivz, imx, imy, imz
- use equations , only : ibx, iby, ibz, ibp
-
-! local variables are not implicit by default
-!
- implicit none
-
-! input arguments
-!
- real(kind=8), dimension(:) , intent(in) :: dx
- real(kind=8), dimension(:,:,:,:) , intent(in) :: q
- real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
-
-! local variables
-!
- integer :: i, j, k = 1
-
-! local temporary arrays
-!
-#if NDIMS == 3
- real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
-#else /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
-#endif /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,2) :: qi
- real(kind=8), dimension(nf,nn) :: fi
-!
-!-------------------------------------------------------------------------------
-!
-#ifdef PROFILE
-! start accounting time for flux update
-!
- call start_timer(imf)
-#endif /* PROFILE */
-
-! initialize fluxes
-!
- f(:,:,:,:,:) = 0.0d+00
-
-! reconstruct interfaces
-!
- call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
-
-! calculate the flux along the X-direction
-!
-#if NDIMS == 3
- do k = nbl, neu
-#endif /* NDIMS == 3 */
- do j = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
- qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
- qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
- qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
- qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
- qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
- qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
- qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(1,idn,:,j,k) = fi(idn,:)
- f(1,imx,:,j,k) = fi(imx,:)
- f(1,imy,:,j,k) = fi(imy,:)
- f(1,imz,:,j,k) = fi(imz,:)
- f(1,ibx,:,j,k) = fi(ibx,:)
- f(1,iby,:,j,k) = fi(iby,:)
- f(1,ibz,:,j,k) = fi(ibz,:)
- f(1,ibp,:,j,k) = fi(ibp,:)
-
- end do ! j = nbl, neu
-
-! calculate the flux along the Y direction
-!
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
- qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
- qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
- qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
- qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
- qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
- qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
- qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(2,idn,i,:,k) = fi(idn,:)
- f(2,imx,i,:,k) = fi(imz,:)
- f(2,imy,i,:,k) = fi(imx,:)
- f(2,imz,i,:,k) = fi(imy,:)
- f(2,ibx,i,:,k) = fi(ibz,:)
- f(2,iby,i,:,k) = fi(ibx,:)
- f(2,ibz,i,:,k) = fi(iby,:)
- f(2,ibp,i,:,k) = fi(ibp,:)
-
- end do ! i = nbl, neu
-#if NDIMS == 3
- end do ! k = nbl, neu
-#endif /* NDIMS == 3 */
-
-#if NDIMS == 3
-! calculate the flux along the Z direction
-!
- do j = nbl, neu
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
- qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
- qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
- qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
- qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
- qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
- qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
- qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(3,idn,i,j,:) = fi(idn,:)
- f(3,imx,i,j,:) = fi(imy,:)
- f(3,imy,i,j,:) = fi(imz,:)
- f(3,imz,i,j,:) = fi(imx,:)
- f(3,ibx,i,j,:) = fi(iby,:)
- f(3,iby,i,j,:) = fi(ibz,:)
- f(3,ibz,i,j,:) = fi(ibx,:)
- f(3,ibp,i,j,:) = fi(ibp,:)
-
- end do ! i = nbl, neu
- end do ! j = nbl, neu
-#endif /* NDIMS == 3 */
-
-#ifdef PROFILE
-! stop accounting time for flux update
-!
- call stop_timer(imf)
-#endif /* PROFILE */
-
-!-------------------------------------------------------------------------------
-!
- end subroutine update_flux_mhd_iso
-!
-!===============================================================================
-!
! subroutine RIEMANN_MHD_ISO_HLLD:
! -------------------------------
!
@@ -3443,198 +4101,6 @@ module schemes
!
!===============================================================================
!
-!***** ADIABATIC MAGNETOHYDRODYNAMICS *****
-!
-!===============================================================================
-!
-! subroutine UPDATE_FLUX_MHD_ADI:
-! ------------------------------
-!
-! Subroutine solves the Riemann problem along each direction and calculates
-! the numerical fluxes, which are used later to calculate the conserved
-! variable increment.
-!
-! Arguments:
-!
-! dx - the spatial step;
-! q - the array of primitive variables;
-! f - the array of numerical fluxes;
-!
-!===============================================================================
-!
- subroutine update_flux_mhd_adi(dx, q, f)
-
-! include external variables
-!
- use coordinates, only : nn => bcells, nbl, neu
- use equations , only : nf
- use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
- use equations , only : ibx, iby, ibz, ibp
-
-! local variables are not implicit by default
-!
- implicit none
-
-! input arguments
-!
- real(kind=8), dimension(:) , intent(in) :: dx
- real(kind=8), dimension(:,:,:,:) , intent(in) :: q
- real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
-
-! local variables
-!
- integer :: i, j, k = 1
-
-! local temporary arrays
-!
-#if NDIMS == 3
- real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
-#else /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
-#endif /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,2) :: qi
- real(kind=8), dimension(nf,nn) :: fi
-!
-!-------------------------------------------------------------------------------
-!
-#ifdef PROFILE
-! start accounting time for flux update
-!
- call start_timer(imf)
-#endif /* PROFILE */
-
-! initialize fluxes
-!
- f(:,:,:,:,:) = 0.0d+00
-
-! reconstruct interfaces
-!
- call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
-
-! calculate the flux along the X-direction
-!
-#if NDIMS == 3
- do k = nbl, neu
-#endif /* NDIMS == 3 */
- do j = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
- qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
- qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
- qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
- qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
- qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
- qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
- qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
- qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(1,idn,:,j,k) = fi(idn,:)
- f(1,imx,:,j,k) = fi(imx,:)
- f(1,imy,:,j,k) = fi(imy,:)
- f(1,imz,:,j,k) = fi(imz,:)
- f(1,ibx,:,j,k) = fi(ibx,:)
- f(1,iby,:,j,k) = fi(iby,:)
- f(1,ibz,:,j,k) = fi(ibz,:)
- f(1,ibp,:,j,k) = fi(ibp,:)
- f(1,ien,:,j,k) = fi(ien,:)
-
- end do ! j = nbl, neu
-
-! calculate the flux along the Y direction
-!
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
- qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
- qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
- qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
- qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
- qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
- qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
- qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
- qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(2,idn,i,:,k) = fi(idn,:)
- f(2,imx,i,:,k) = fi(imz,:)
- f(2,imy,i,:,k) = fi(imx,:)
- f(2,imz,i,:,k) = fi(imy,:)
- f(2,ibx,i,:,k) = fi(ibz,:)
- f(2,iby,i,:,k) = fi(ibx,:)
- f(2,ibz,i,:,k) = fi(iby,:)
- f(2,ibp,i,:,k) = fi(ibp,:)
- f(2,ien,i,:,k) = fi(ien,:)
-
- end do ! i = nbl, neu
-#if NDIMS == 3
- end do ! k = nbl, neu
-#endif /* NDIMS == 3 */
-
-#if NDIMS == 3
-! calculate the flux along the Z direction
-!
- do j = nbl, neu
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
- qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
- qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
- qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
- qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
- qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
- qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
- qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
- qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(3,idn,i,j,:) = fi(idn,:)
- f(3,imx,i,j,:) = fi(imy,:)
- f(3,imy,i,j,:) = fi(imz,:)
- f(3,imz,i,j,:) = fi(imx,:)
- f(3,ibx,i,j,:) = fi(iby,:)
- f(3,iby,i,j,:) = fi(ibz,:)
- f(3,ibz,i,j,:) = fi(ibx,:)
- f(3,ibp,i,j,:) = fi(ibp,:)
- f(3,ien,i,j,:) = fi(ien,:)
-
- end do ! i = nbl, neu
- end do ! j = nbl, neu
-#endif /* NDIMS == 3 */
-
-#ifdef PROFILE
-! stop accounting time for flux update
-!
- call stop_timer(imf)
-#endif /* PROFILE */
-
-!-------------------------------------------------------------------------------
-!
- end subroutine update_flux_mhd_adi
-!
-!===============================================================================
-!
! subroutine RIEMANN_MHD_ADI_HLLD:
! -------------------------------
!
@@ -4496,236 +4962,6 @@ module schemes
!
!===============================================================================
!
-!***** ADIABATIC SPECIAL RELATIVITY HYDRODYNAMICS *****
-!
-!===============================================================================
-!
-! subroutine UPDATE_FLUX_SRHD_ADI:
-! -------------------------------
-!
-! Subroutine solves the Riemann problem along each direction and calculates
-! the numerical fluxes, which are used later to calculate the conserved
-! variable increment.
-!
-! Arguments:
-!
-! dx - the spatial step;
-! q - the array of primitive variables;
-! f - the array of numerical fluxes;
-!
-!===============================================================================
-!
- subroutine update_flux_srhd_adi(dx, q, f)
-
-! include external variables
-!
- use coordinates, only : nn => bcells, nbl, neu
- use equations , only : nf
- use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
-
-! local variables are not implicit by default
-!
- implicit none
-
-! input arguments
-!
- real(kind=8), dimension(:) , intent(in) :: dx
- real(kind=8), dimension(:,:,:,:) , intent(in) :: q
- real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
-
-! local variables
-!
- integer :: i, j, k = 1, l, p
- real(kind=8) :: vm
-
-! local temporary arrays
-!
-#if NDIMS == 3
- real(kind=8), dimension(nf,nn,nn,nn) :: qq
- real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
-#else /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,nn, 1) :: qq
- real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
-#endif /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,2) :: qi
- real(kind=8), dimension(nf,nn) :: fi
-!
-!-------------------------------------------------------------------------------
-!
-#ifdef PROFILE
-! start accounting time for flux update
-!
- call start_timer(imf)
-#endif /* PROFILE */
-
-! initialize fluxes
-!
- f(:,:,:,:,:) = 0.0d+00
-
-! apply reconstruction on variables using 4-vector or 3-vector
-!
- if (states_4vec) then
-
-! convert velocities to four-velocities for physical reconstruction
-!
-#if NDIMS == 3
- do k = 1, nn
-#endif /* NDIMS == 3 */
- do j = 1, nn
- do i = 1, nn
-
- vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))
-
- qq(idn,i,j,k) = q(idn,i,j,k) / vm
- qq(ivx,i,j,k) = q(ivx,i,j,k) / vm
- qq(ivy,i,j,k) = q(ivy,i,j,k) / vm
- qq(ivz,i,j,k) = q(ivz,i,j,k) / vm
- qq(ipr,i,j,k) = q(ipr,i,j,k)
-
- end do ! i = 1, nn
- end do ! j = 1, nn
-#if NDIMS == 3
- end do ! k = 1, nn
-#endif /* NDIMS == 3 */
-
-! reconstruct interfaces
-!
- call reconstruct_interfaces(dx(:), qq(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
-
-! convert state four-velocities back to velocities
-!
-#if NDIMS == 3
- do k = 1, nn
-#endif /* NDIMS == 3 */
- do j = 1, nn
- do i = 1, nn
-
- do l = 1, 2
- do p = 1, NDIMS
-
- vm = sqrt(1.0d+00 + sum(qs(ivx:ivz,i,j,k,l,p)**2))
-
- qs(idn,i,j,k,l,p) = qs(idn,i,j,k,l,p) / vm
- qs(ivx,i,j,k,l,p) = qs(ivx,i,j,k,l,p) / vm
- qs(ivy,i,j,k,l,p) = qs(ivy,i,j,k,l,p) / vm
- qs(ivz,i,j,k,l,p) = qs(ivz,i,j,k,l,p) / vm
-
- end do ! p = 1, ndims
- end do ! l = 1, 2
- end do ! i = 1, nn
- end do ! j = 1, nn
-#if NDIMS == 3
- end do ! k = 1, nn
-#endif /* NDIMS == 3 */
-
- else
-
-! reconstruct interfaces
-!
- call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
-
- end if
-
-! calculate the flux along the X-direction
-!
-#if NDIMS == 3
- do k = nbl, neu
-#endif /* NDIMS == 3 */
- do j = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
- qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
- qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
- qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
- qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(1,idn,:,j,k) = fi(idn,:)
- f(1,imx,:,j,k) = fi(imx,:)
- f(1,imy,:,j,k) = fi(imy,:)
- f(1,imz,:,j,k) = fi(imz,:)
- f(1,ien,:,j,k) = fi(ien,:)
-
- end do ! j = nbl, neu
-
-! calculate the flux along the Y direction
-!
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
- qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
- qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
- qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
- qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(2,idn,i,:,k) = fi(idn,:)
- f(2,imx,i,:,k) = fi(imz,:)
- f(2,imy,i,:,k) = fi(imx,:)
- f(2,imz,i,:,k) = fi(imy,:)
- f(2,ien,i,:,k) = fi(ien,:)
-
- end do ! i = nbl, neu
-#if NDIMS == 3
- end do ! k = nbl, neu
-#endif /* NDIMS == 3 */
-
-#if NDIMS == 3
-! calculate the flux along the Z direction
-!
- do j = nbl, neu
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
- qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
- qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
- qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
- qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(3,idn,i,j,:) = fi(idn,:)
- f(3,imx,i,j,:) = fi(imy,:)
- f(3,imy,i,j,:) = fi(imz,:)
- f(3,imz,i,j,:) = fi(imx,:)
- f(3,ien,i,j,:) = fi(ien,:)
-
- end do ! i = nbl, neu
- end do ! j = nbl, neu
-#endif /* NDIMS == 3 */
-
-#ifdef PROFILE
-! stop accounting time for flux update
-!
- call stop_timer(imf)
-#endif /* PROFILE */
-
-!-------------------------------------------------------------------------------
-!
- end subroutine update_flux_srhd_adi
-!
-!===============================================================================
-!
! subroutine RIEMANN_SRHD_ADI_HLLC:
! --------------------------------
!
@@ -4956,266 +5192,6 @@ module schemes
!
!===============================================================================
!
-!***** ADIABATIC SPECIAL RELATIVITY MAGNETOHYDRODYNAMICS *****
-!
-!===============================================================================
-!
-! subroutine UPDATE_FLUX_SRMHD_ADI:
-! --------------------------------
-!
-! Subroutine solves the Riemann problem along each direction and calculates
-! the numerical fluxes, which are used later to calculate the conserved
-! variable increment.
-!
-! Arguments:
-!
-! dx - the spatial step;
-! q - the array of primitive variables;
-! f - the array of numerical fluxes;
-!
-!===============================================================================
-!
- subroutine update_flux_srmhd_adi(dx, q, f)
-
-! include external variables
-!
- use coordinates, only : nn => bcells, nbl, neu
- use equations , only : nf
- use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
- use equations , only : ibx, iby, ibz, ibp
-
-! local variables are not implicit by default
-!
- implicit none
-
-! input arguments
-!
- real(kind=8), dimension(:) , intent(in) :: dx
- real(kind=8), dimension(:,:,:,:) , intent(in) :: q
- real(kind=8), dimension(:,:,:,:,:), intent(out) :: f
-
-! local variables
-!
- integer :: i, j, k = 1, l, p
- real(kind=8) :: vm
-
-! local temporary arrays
-!
-#if NDIMS == 3
- real(kind=8), dimension(nf,nn,nn,nn) :: qq
- real(kind=8), dimension(nf,nn,nn,nn,2,NDIMS) :: qs
-#else /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,nn, 1) :: qq
- real(kind=8), dimension(nf,nn,nn, 1,2,NDIMS) :: qs
-#endif /* NDIMS == 3 */
- real(kind=8), dimension(nf,nn,2) :: qi
- real(kind=8), dimension(nf,nn) :: fi
-!
-!-------------------------------------------------------------------------------
-!
-#ifdef PROFILE
-! start accounting time for flux update
-!
- call start_timer(imf)
-#endif /* PROFILE */
-
-! initialize fluxes
-!
- f(:,:,:,:,:) = 0.0d+00
-
-! apply reconstruction on variables using 4-vector or 3-vector
-!
- if (states_4vec) then
-
-! convert velocities to four-velocities for physical reconstruction
-!
-#if NDIMS == 3
- do k = 1, nn
-#endif /* NDIMS == 3 */
- do j = 1, nn
- do i = 1, nn
-
- vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2))
-
- qq(idn,i,j,k) = q(idn,i,j,k) / vm
- qq(ivx,i,j,k) = q(ivx,i,j,k) / vm
- qq(ivy,i,j,k) = q(ivy,i,j,k) / vm
- qq(ivz,i,j,k) = q(ivz,i,j,k) / vm
- qq(ipr,i,j,k) = q(ipr,i,j,k)
- qq(ibx,i,j,k) = q(ibx,i,j,k)
- qq(iby,i,j,k) = q(iby,i,j,k)
- qq(ibz,i,j,k) = q(ibz,i,j,k)
- qq(ibp,i,j,k) = q(ibp,i,j,k)
- qq(ipr,i,j,k) = q(ipr,i,j,k)
-
- end do ! i = 1, nn
- end do ! j = 1, nn
-#if NDIMS == 3
- end do ! k = 1, nn
-#endif /* NDIMS == 3 */
-
-! reconstruct interfaces
-!
- call reconstruct_interfaces(dx(:), qq(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
-
-! convert state four-velocities back to velocities
-!
-#if NDIMS == 3
- do k = 1, nn
-#endif /* NDIMS == 3 */
- do j = 1, nn
- do i = 1, nn
-
- do l = 1, 2
- do p = 1, NDIMS
-
- vm = sqrt(1.0d+00 + sum(qs(ivx:ivz,i,j,k,l,p)**2))
-
- qs(idn,i,j,k,l,p) = qs(idn,i,j,k,l,p) / vm
- qs(ivx,i,j,k,l,p) = qs(ivx,i,j,k,l,p) / vm
- qs(ivy,i,j,k,l,p) = qs(ivy,i,j,k,l,p) / vm
- qs(ivz,i,j,k,l,p) = qs(ivz,i,j,k,l,p) / vm
-
- end do ! p = 1, ndims
- end do ! l = 1, 2
- end do ! i = 1, nn
- end do ! j = 1, nn
-#if NDIMS == 3
- end do ! k = 1, nn
-#endif /* NDIMS == 3 */
-
- else
-
-! reconstruct interfaces
-!
- call reconstruct_interfaces(dx(:), q(:,:,:,:), qs(:,:,:,:,1:2,1:NDIMS))
-
- end if
-
-! calculate the flux along the X-direction
-!
-#if NDIMS == 3
- do k = nbl, neu
-#endif /* NDIMS == 3 */
- do j = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,:,j,k,1:2,1)
- qi(ivx,:,1:2) = qs(ivx,:,j,k,1:2,1)
- qi(ivy,:,1:2) = qs(ivy,:,j,k,1:2,1)
- qi(ivz,:,1:2) = qs(ivz,:,j,k,1:2,1)
- qi(ibx,:,1:2) = qs(ibx,:,j,k,1:2,1)
- qi(iby,:,1:2) = qs(iby,:,j,k,1:2,1)
- qi(ibz,:,1:2) = qs(ibz,:,j,k,1:2,1)
- qi(ibp,:,1:2) = qs(ibp,:,j,k,1:2,1)
- qi(ipr,:,1:2) = qs(ipr,:,j,k,1:2,1)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(1,idn,:,j,k) = fi(idn,:)
- f(1,imx,:,j,k) = fi(imx,:)
- f(1,imy,:,j,k) = fi(imy,:)
- f(1,imz,:,j,k) = fi(imz,:)
- f(1,ibx,:,j,k) = fi(ibx,:)
- f(1,iby,:,j,k) = fi(iby,:)
- f(1,ibz,:,j,k) = fi(ibz,:)
- f(1,ibp,:,j,k) = fi(ibp,:)
- f(1,ien,:,j,k) = fi(ien,:)
-
- end do ! j = nbl, neu
-
-! calculate the flux along the Y direction
-!
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,:,k,1:2,2)
- qi(ivx,:,1:2) = qs(ivy,i,:,k,1:2,2)
- qi(ivy,:,1:2) = qs(ivz,i,:,k,1:2,2)
- qi(ivz,:,1:2) = qs(ivx,i,:,k,1:2,2)
- qi(ibx,:,1:2) = qs(iby,i,:,k,1:2,2)
- qi(iby,:,1:2) = qs(ibz,i,:,k,1:2,2)
- qi(ibz,:,1:2) = qs(ibx,i,:,k,1:2,2)
- qi(ibp,:,1:2) = qs(ibp,i,:,k,1:2,2)
- qi(ipr,:,1:2) = qs(ipr,i,:,k,1:2,2)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(2,idn,i,:,k) = fi(idn,:)
- f(2,imx,i,:,k) = fi(imz,:)
- f(2,imy,i,:,k) = fi(imx,:)
- f(2,imz,i,:,k) = fi(imy,:)
- f(2,ibx,i,:,k) = fi(ibz,:)
- f(2,iby,i,:,k) = fi(ibx,:)
- f(2,ibz,i,:,k) = fi(iby,:)
- f(2,ibp,i,:,k) = fi(ibp,:)
- f(2,ien,i,:,k) = fi(ien,:)
-
- end do ! i = nbl, neu
-#if NDIMS == 3
- end do ! k = nbl, neu
-#endif /* NDIMS == 3 */
-
-#if NDIMS == 3
-! calculate the flux along the Z direction
-!
- do j = nbl, neu
- do i = nbl, neu
-
-! copy directional variable vectors to pass to the one dimensional solver
-!
- qi(idn,:,1:2) = qs(idn,i,j,:,1:2,3)
- qi(ivx,:,1:2) = qs(ivz,i,j,:,1:2,3)
- qi(ivy,:,1:2) = qs(ivx,i,j,:,1:2,3)
- qi(ivz,:,1:2) = qs(ivy,i,j,:,1:2,3)
- qi(ibx,:,1:2) = qs(ibz,i,j,:,1:2,3)
- qi(iby,:,1:2) = qs(ibx,i,j,:,1:2,3)
- qi(ibz,:,1:2) = qs(iby,i,j,:,1:2,3)
- qi(ibp,:,1:2) = qs(ibp,i,j,:,1:2,3)
- qi(ipr,:,1:2) = qs(ipr,i,j,:,1:2,3)
-
-! call one dimensional Riemann solver in order to obtain numerical fluxes
-!
- call numerical_flux(qi(:,:,:), fi(:,:))
-
-! update the array of fluxes
-!
- f(3,idn,i,j,:) = fi(idn,:)
- f(3,imx,i,j,:) = fi(imy,:)
- f(3,imy,i,j,:) = fi(imz,:)
- f(3,imz,i,j,:) = fi(imx,:)
- f(3,ibx,i,j,:) = fi(iby,:)
- f(3,iby,i,j,:) = fi(ibz,:)
- f(3,ibz,i,j,:) = fi(ibx,:)
- f(3,ibp,i,j,:) = fi(ibp,:)
- f(3,ien,i,j,:) = fi(ien,:)
-
- end do ! i = nbl, neu
- end do ! j = nbl, neu
-#endif /* NDIMS == 3 */
-
-#ifdef PROFILE
-! stop accounting time for flux update
-!
- call stop_timer(imf)
-#endif /* PROFILE */
-
-!-------------------------------------------------------------------------------
-!
- end subroutine update_flux_srmhd_adi
-!
-!===============================================================================
-!
! subroutine RIEMANN_SRMHD_ADI_HLLC:
! ---------------------------------
!