Calculate gradient of temperature
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real, | intent(out), | dimension( 0:imx , 0:jmx , 0:kmx ) | :: | grad | Output gradient of termperature |
|
| character(len=*), | intent(in) | :: | dir | Direction with respect to which gradients are calculated |
subroutine compute_gradient_T(grad, dir)
!< Calculate gradient of temperature
implicit none
real, dimension( 0:imx , 0:jmx , 0:kmx ), intent(out) :: grad
!< Output gradient of termperature
character(len=*) , intent(in) :: dir
!< Direction with respect to which gradients are calculated
real, dimension(6) :: T
real :: cell_T
real, dimension(:,:,:), pointer :: nx
real, dimension(:,:,:), pointer :: ny
real, dimension(:,:,:), pointer :: nz
integer :: i
integer :: j
integer :: k
select case(dir)
case('x')
nx(-2:imx+3,-2:jmx+2,-2:kmx+2) => xn(:,:,:,1)
ny(-2:imx+2,-2:jmx+3,-2:kmx+2) => yn(:,:,:,1)
nz(-2:imx+2,-2:jmx+2,-2:kmx+3) => zn(:,:,:,1)
case('y')
nx(-2:imx+3,-2:jmx+2,-2:kmx+2) => xn(:,:,:,2)
ny(-2:imx+2,-2:jmx+3,-2:kmx+2) => yn(:,:,:,2)
nz(-2:imx+2,-2:jmx+2,-2:kmx+3) => zn(:,:,:,2)
case('z')
nx(-2:imx+3,-2:jmx+2,-2:kmx+2) => xn(:,:,:,3)
ny(-2:imx+2,-2:jmx+3,-2:kmx+2) => yn(:,:,:,3)
nz(-2:imx+2,-2:jmx+2,-2:kmx+3) => zn(:,:,:,3)
case DEFAULT
nx(-2:imx+3,-2:jmx+2,-2:kmx+2) => xn(:,:,:,1)
ny(-2:imx+2,-2:jmx+3,-2:kmx+2) => yn(:,:,:,1)
nz(-2:imx+2,-2:jmx+2,-2:kmx+3) => zn(:,:,:,1)
print*, "ERROR: gradient direction error"
end select
grad = 0.0
do k=0,kmx
do j=0,jmx
do i=0,imx
cell_T = (pressure(i,j,k)/density(i,j,k))/R_gas
T(1) = (pressure(i-1,j,k)/density(i-1,j,k))/R_gas + cell_T
T(2) = (pressure(i,j-1,k)/density(i,j-1,k))/R_gas + cell_T
T(3) = (pressure(i,j,k-1)/density(i,j,k-1))/R_gas + cell_T
T(4) = (pressure(i+1,j,k)/density(i+1,j,k))/R_gas + cell_T
T(5) = (pressure(i,j+1,k)/density(i,j+1,k))/R_gas + cell_T
T(6) = (pressure(i,j,k+1)/density(i,j,k+1))/R_gas + cell_T
grad(i,j,k) =(-T(1)*nx(i,j,k)*xA(i,j,k) &
-T(2)*ny(i,j,k)*yA(i,j,k) &
-T(3)*nz(i,j,k)*zA(i,j,k) &
+T(4)*nx(i+1,j ,k )*xA(i+1,j ,k ) &
+T(5)*ny(i ,j+1,k )*yA(i ,j+1,k ) &
+T(6)*nz(i ,j ,k+1)*zA(i ,j ,k+1) &
)/(2*volume(i,j,k))
end do
end do
end do
if(any(isnan(grad)))then
Fatal_error
end if
end subroutine compute_gradient_T