diff --git a/src/equations.F90 b/src/equations.F90 index bee1b63..b136eef 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -308,6 +308,95 @@ module equations !------------------------------------------------------------------------------- ! end subroutine fluxspeed +! +!=============================================================================== +! +! subroutine UPDATE_PRIMITIVE_VARIABLES: +! ------------------------------------- +! +! Subroutine updates primitive variables from their conservative +! representation. This process is done ones after advance of the conserved +! variables due to their evolution in time. +! +!=============================================================================== +! + subroutine update_primitive_variables(uu, qq) + +! include external procedures and variables +! + use coordinates, only : im, jm, km +#ifdef DEBUG + use variables , only : idn +#ifdef ADI + use variables , only : ipr +#endif /* ADI */ +#endif /* DEBUG */ + use variables , only : nt + +! local variables are not implicit by default +! + implicit none + +! input/output arguments +! + real, dimension(nt,im,jm,km), intent(in) :: uu + real, dimension(nt,im,jm,km), intent(inout) :: qq + +! temporary variables +! +#ifdef DEBUG + integer :: i +#endif /* DEBUG */ + integer :: j, k + +! temporary array to store conserved variable vector +! + real, dimension(nt,im) :: u +! +!------------------------------------------------------------------------------- +! +! update primitive variables +! + do k = 1, km + do j = 1, jm + +! copy variables to temporary array of conserved variables +! + u(1:nt,1:im) = uu(1:nt,1:im,j,k) + +! convert conserved variables to primitive ones +! + call cons2prim(im, u(1:nt,1:im), qq(1:nt,1:im,j,k)) + +#ifdef DEBUG +! check the positivity of density and pressure +! + do i = 1, im + if (qq(idn,i,j,k) <= 0.0d0) then + write(*,*) + write(*,*) 'Unphysical state: ρ ≤ 0', i, j, k + write(*,"('Q = ',4(1pe24.16))") qq(1:4 ,i,j,k) + write(*,"(' ',4(1pe24.16))") qq(5:nt,i,j,k) + stop + end if +#ifdef ADI + if (qq(ipr,i,j,k) <= 0.0d0) then + write(*,*) + write(*,*) 'Unphysical state: p ≤ 0', i, j, k + write(*,"('Q = ',4(1pe24.16))") qq(1:4 ,i,j,k) + write(*,"(' ',4(1pe24.16))") qq(5:nt,i,j,k) + stop + end if +#endif /* ADI */ + end do +#endif /* DEBUG */ + + end do + end do + +!------------------------------------------------------------------------------- +! + end subroutine update_primitive_variables #endif /* HYDRO */ #ifdef MHD ! diff --git a/src/makefile b/src/makefile index 6d92f25..7b5f534 100644 --- a/src/makefile +++ b/src/makefile @@ -212,7 +212,7 @@ driver.o : driver.F90 blocks.o config.o coordinates.o equations.o \ evolution.o forcing.o integrals.o interpolations.o io.o \ mesh.o mpitools.o parameters.o problems.o random.o \ refinement.o -equations.o : equations.F90 parameters.o variables.o +equations.o : equations.F90 coordinates.o parameters.o variables.o error.o : error.F90 evolution.o : evolution.F90 blocks.o boundaries.o config.o coordinates.o \ forcing.o interpolations.o mesh.o mpitools.o problems.o \