Total Pressure Riemann boundary condition
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| character(len=*), | intent(in) | :: | face | |||
| type(facetype), | intent(in), | dimension(-2:dims%imx+3,-2:dims%jmx+2,-2:dims%kmx+2) | :: | Ifaces | Input varaible which stores I faces' area and unit normal |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+3,-2:dims%kmx+2) | :: | Jfaces | Input varaible which stores J faces' area and unit normal |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+3) | :: | Kfaces | Input varaible which stores K faces' area and unit normal |
|
| type(boundarytype), | intent(in) | :: | bc | |||
| type(extent), | intent(in) | :: | dims |
subroutine total_pressure(face, Ifaces, Jfaces, Kfaces, bc, dims)
!< Total Pressure Riemann boundary condition
implicit none
type(extent), intent(in) :: dims
character(len=*), intent(in) :: face
type(facetype), dimension(-2:dims%imx+3,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: Ifaces
!< Input varaible which stores I faces' area and unit normal
type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+3,-2:dims%kmx+2), intent(in) :: Jfaces
!< Input varaible which stores J faces' area and unit normal
type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+3), intent(in) :: Kfaces
!< Input varaible which stores K faces' area and unit normal
type(boundarytype), intent(in) :: bc
real(wp) :: cinf, cexp ! speed of sound
real(wp) :: Rinf, Rexp ! Riemann invarient
real(wp) :: Uninf, Unexp ! face normal speed
real(wp) :: Unb ! normal velocity boundary
real(wp) :: Cb ! speed of sound boundary
real(wp) :: vel_diff
real(wp) :: u,v,w
real(wp) :: uf, vf, wf
real(wp) :: Mb
integer :: i,j,k
select case(face)
case("imin")
do k = 1,kmx-1
do j = 1,jmx-1
do i = 1,1
! interior cell
u = x_speed(i,j,k)
v = y_speed(i,j,k)
w = z_speed(i,j,k)
! ghost cell
uf = x_speed_inf!x_speed(i-1,j,k)
vf = y_speed_inf!y_speed(i-1,j,k)
wf = z_speed_inf!z_speed(i-1,j,k)
cexp = sqrt(gm*pressure(i,j,k)/density(i,j,k))
!cinf = sqrt(gm*pressure(i-1,j,k)/density(i-1,j,k))
cinf = sqrt(gm*pressure_inf/density_inf)
Unexp = u *(-Ifaces(i,j,k)%nx) + v *(-Ifaces(i,j,k)%ny) + w *(-Ifaces(i,j,k)%nz)
Uninf = uf*(-Ifaces(i,j,k)%nx) + vf*(-Ifaces(i,j,k)%ny) + wf*(-Ifaces(i,j,k)%nz)
Rinf = Uninf - 2*cinf/(gm-1.)
Rexp = Unexp + 2*cexp/(gm-1.)
Unb = 0.5*(Rexp + Rinf)
Cb = 0.25*(gm-1.)*(Rexp - Rinf)
if(Unb > 0.)then
vel_diff = Unb - Unexp
x_speed(i-1,j,k) = x_speed(i,j,k) + vel_diff*(-Ifaces(i,j,k)%nx)
y_speed(i-1,j,k) = y_speed(i,j,k) + vel_diff*(-Ifaces(i,j,k)%ny)
z_speed(i-1,j,k) = z_speed(i,j,k) + vel_diff*(-Ifaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call copy3(tv, "flat", face, bc, dims)
case('sst', 'sst2003')
call copy3(tk, "flat", face, bc, dims)
call copy3(tw, "flat", face, bc, dims)
case('kkl')
call copy3(tk, "flat", face, bc, dims)
call copy3(tkl, "flat", face, bc, dims)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
call copy3(tgm, "flat", face, bc, dims)
case DEFAULT
continue
end select
else
vel_diff = Unb - Uninf
x_speed(i-1,j,k) = x_speed_inf + vel_diff*(-Ifaces(i,j,k)%nx)
y_speed(i-1,j,k) = y_speed_inf + vel_diff*(-Ifaces(i,j,k)%ny)
z_speed(i-1,j,k) = z_speed_inf + vel_diff*(-Ifaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call fix(tv, bc%fixed_tv, face)
case('sst', 'sst2003')
!call check_if_value_fixed(bc, "sst")
call fix(tk, bc%fixed_tk, face)
call fix(tw, bc%fixed_tw, face)
case('kkl')
!call check_if_value_fixed(bc, "kkl")
call fix(tk, bc%fixed_tk, face)
call fix(tkl, bc%fixed_tkl, face)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
!call check_if_value_fixed(bc, "lctm2015")
call fix(tgm, bc%fixed_tgm, face)
case DEFAULT
continue
end select
end if
Mb = sqrt(x_speed(i-1,j,k)**2+y_speed(i-1,j,k)**2+z_speed(i-1,j,k)**2)/Cb
pressure(i-1,j,k) = bc%fixed_Tpressure(1)/(((1+0.5*(gm-1.)*Mb*Mb))**(gm/(gm-1.)))
density(i-1,j,k) = gm*pressure(i-1,j,k)/(Cb*Cb)
end do
end do
end do
qp(-1,:,:,:) = qp(0,:,:,:)
qp(-2,:,:,:) = qp(0,:,:,:)
case("imax")
do k = 1,kmx-1
do j = 1,jmx-1
do i = imx,imx
! interior cell
u = x_speed(i-1,j,k)
v = y_speed(i-1,j,k)
w = z_speed(i-1,j,k)
! ghost cell
uf = x_speed_inf!x_speed(i,j,k)
vf = y_speed_inf!y_speed(i,j,k)
wf = z_speed_inf!z_speed(i,j,k)
cexp = sqrt(gm*pressure(i-1,j,k)/density(i-1,j,k))
!cinf = sqrt(gm*pressure(i,j,k)/density(i,j,k))
cinf = sqrt(gm*pressure_inf/density_inf)
Unexp = u *(Ifaces(i,j,k)%nx) + v *(Ifaces(i,j,k)%ny) + w *(Ifaces(i,j,k)%nz)
Uninf = uf*(Ifaces(i,j,k)%nx) + vf*(Ifaces(i,j,k)%ny) + wf*(Ifaces(i,j,k)%nz)
Rinf = Uninf - 2*cinf/(gm-1.)
Rexp = Unexp + 2*cexp/(gm-1.)
Unb = 0.5*(Rexp + Rinf)
Cb = 0.25*(gm-1.)*(Rexp - Rinf)
if(Unb > 0.)then
vel_diff = Unb - Unexp
x_speed(i,j,k) = x_speed(i-1,j,k) + vel_diff*(Ifaces(i,j,k)%nx)
y_speed(i,j,k) = y_speed(i-1,j,k) + vel_diff*(Ifaces(i,j,k)%ny)
z_speed(i,j,k) = z_speed(i-1,j,k) + vel_diff*(Ifaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call copy3(tv, "flat", face, bc, dims)
case('sst', 'sst2003')
call copy3(tk, "flat", face, bc, dims)
call copy3(tw, "flat", face, bc, dims)
case('kkl')
call copy3(tk, "flat", face, bc, dims)
call copy3(tkl, "flat", face, bc, dims)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
call copy3(tgm, "flat", face, bc, dims)
case DEFAULT
continue
end select
else
vel_diff = Unb - Uninf
x_speed(i,j,k) = x_speed_inf + vel_diff*(Ifaces(i,j,k)%nx)
y_speed(i,j,k) = y_speed_inf + vel_diff*(Ifaces(i,j,k)%ny)
z_speed(i,j,k) = z_speed_inf + vel_diff*(Ifaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call fix(tv, bc%fixed_tv, face)
case('sst', 'sst2003')
!call check_if_value_fixed(bc, "sst")
call fix(tk, bc%fixed_tk, face)
call fix(tw, bc%fixed_tw, face)
case('kkl')
!call check_if_value_fixed(bc, "kkl")
call fix(tk, bc%fixed_tk, face)
call fix(tkl, bc%fixed_tkl, face)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
!call check_if_value_fixed(bc, "lctm2015")
call fix(tgm, bc%fixed_tgm, face)
case DEFAULT
continue
end select
end if
Mb = sqrt(x_speed(i,j,k)**2+y_speed(i,j,k)**2+z_speed(i,j,k)**2)/Cb
pressure(i,j,k) = bc%fixed_Tpressure(2)/(((1+0.5*(gm-1.)*Mb*Mb))**(gm/(gm-1.)))
density(i,j,k) = gm*pressure(i,j,k)/(Cb*Cb)
end do
end do
end do
qp(imx+1,:,:,:) = qp(imx,:,:,:)
qp(imx+2,:,:,:) = qp(imx,:,:,:)
case("jmin")
do k = 1,kmx-1
do j = 1,1
do i = 1,imx-1
! interior cell
u = x_speed(i,j,k)
v = y_speed(i,j,k)
w = z_speed(i,j,k)
! ghost cell
uf = x_speed_inf!x_speed(i,j-1,k)
vf = y_speed_inf!y_speed(i,j-1,k)
wf = z_speed_inf!z_speed(i,j-1,k)
cexp = sqrt(gm*pressure(i,j,k)/density(i,j,k))
!cinf = sqrt(gm*pressure(i,j-1,k)/density(i,j-1,k))
cinf = sqrt(gm*pressure_inf/density_inf)
Unexp = u *(-Jfaces(i,j,k)%nx) + v *(-Jfaces(i,j,k)%ny) + w *(-Jfaces(i,j,k)%nz)
Uninf = uf*(-Jfaces(i,j,k)%nx) + vf*(-Jfaces(i,j,k)%ny) + wf*(-Jfaces(i,j,k)%nz)
Rinf = Uninf - 2*cinf/(gm-1.)
Rexp = Unexp + 2*cexp/(gm-1.)
Unb = 0.5*(Rexp + Rinf)
Cb = 0.25*(gm-1.)*(Rexp - Rinf)
if(Unb > 0.)then
vel_diff = Unb - Unexp
x_speed(i,j-1,k) = x_speed(i,j,k) + vel_diff*(-Jfaces(i,j,k)%nx)
y_speed(i,j-1,k) = y_speed(i,j,k) + vel_diff*(-Jfaces(i,j,k)%ny)
z_speed(i,j-1,k) = z_speed(i,j,k) + vel_diff*(-Jfaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call copy3(tv, "flat", face, bc, dims)
case('sst', 'sst2003')
call copy3(tk, "flat", face, bc, dims)
call copy3(tw, "flat", face, bc, dims)
case('kkl')
call copy3(tk, "flat", face, bc, dims)
call copy3(tkl, "flat", face, bc, dims)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
call copy3(tgm, "flat", face, bc, dims)
case DEFAULT
continue
end select
else
vel_diff = Unb - Uninf
x_speed(i,j-1,k) = x_speed_inf + vel_diff*(-Jfaces(i,j,k)%nx)
y_speed(i,j-1,k) = y_speed_inf + vel_diff*(-Jfaces(i,j,k)%ny)
z_speed(i,j-1,k) = z_speed_inf + vel_diff*(-Jfaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call fix(tv, bc%fixed_tv, face)
case('sst', 'sst2003')
!call check_if_value_fixed(bc, "sst")
call fix(tk, bc%fixed_tk, face)
call fix(tw, bc%fixed_tw, face)
case('kkl')
!call check_if_value_fixed(bc, "kkl")
call fix(tk, bc%fixed_tk, face)
call fix(tkl, bc%fixed_tkl, face)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
!call check_if_value_fixed(bc, "lctm2015")
call fix(tgm, bc%fixed_tgm, face)
case DEFAULT
continue
end select
end if
Mb = sqrt(x_speed(i,j-1,k)**2+y_speed(i,j-1,k)**2+z_speed(i,j-1,k)**2)/Cb
pressure(i,j-1,k) = bc%fixed_Tpressure(3)/(((1+0.5*(gm-1.)*Mb*Mb))**(gm/(gm-1.)))
density(i,j-1,k) = gm*pressure(i,j-1,k)/(Cb*Cb)
end do
end do
end do
qp(:,-1,:,:) = qp(:,0,:,:)
qp(:,-2,:,:) = qp(:,0,:,:)
case("jmax")
do k = 1,kmx-1
do j = jmx,jmx
do i = 1,imx-1
! interior cell
u = x_speed(i,j-1,k)
v = y_speed(i,j-1,k)
w = z_speed(i,j-1,k)
! ghost cell
uf = x_speed_inf!x_speed(i,j,k)
vf = y_speed_inf!y_speed(i,j,k)
wf = z_speed_inf!z_speed(i,j,k)
cexp = sqrt(gm*pressure(i,j-1,k)/density(i,j-1,k))
!cinf = sqrt(gm*pressure(i,j,k)/density(i,j,k))
cinf = sqrt(gm*pressure_inf/density_inf)
Unexp = u *(Jfaces(i,j,k)%nx) + v *(Jfaces(i,j,k)%ny) + w *(Jfaces(i,j,k)%nz)
Uninf = uf*(Jfaces(i,j,k)%nx) + vf*(Jfaces(i,j,k)%ny) + wf*(Jfaces(i,j,k)%nz)
Rinf = Uninf - 2*cinf/(gm-1.)
Rexp = Unexp + 2*cexp/(gm-1.)
Unb = 0.5*(Rexp + Rinf)
Cb = 0.25*(gm-1.)*(Rexp - Rinf)
if(Unb > 0.)then
vel_diff = Unb - Unexp
x_speed(i,j,k) = x_speed(i,j-1,k) + vel_diff*(Jfaces(i,j,k)%nx)
y_speed(i,j,k) = y_speed(i,j-1,k) + vel_diff*(Jfaces(i,j,k)%ny)
z_speed(i,j,k) = z_speed(i,j-1,k) + vel_diff*(Jfaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call copy3(tv, "flat", face, bc, dims)
case('sst', 'sst2003')
call copy3(tk, "flat", face, bc, dims)
call copy3(tw, "flat", face, bc, dims)
case('kkl')
call copy3(tk, "flat", face, bc, dims)
call copy3(tkl, "flat", face, bc, dims)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
call copy3(tgm, "flat", face, bc, dims)
case DEFAULT
continue
end select
else
vel_diff = Unb - Uninf
x_speed(i,j,k) = x_speed_inf + vel_diff*(Jfaces(i,j,k)%nx)
y_speed(i,j,k) = y_speed_inf + vel_diff*(Jfaces(i,j,k)%ny)
z_speed(i,j,k) = z_speed_inf + vel_diff*(Jfaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call fix(tv, bc%fixed_tv, face)
case('sst', 'sst2003')
!call check_if_value_fixed(bc, "sst")
call fix(tk, bc%fixed_tk, face)
call fix(tw, bc%fixed_tw, face)
case('kkl')
!call check_if_value_fixed(bc, "kkl")
call fix(tk, bc%fixed_tk, face)
call fix(tkl, bc%fixed_tkl, face)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
!call check_if_value_fixed(bc, "lctm2015")
call fix(tgm, bc%fixed_tgm, face)
case DEFAULT
continue
end select
end if
Mb = sqrt(x_speed(i,j,k)**2+y_speed(i,j,k)**2+z_speed(i,j,k)**2)/Cb
pressure(i,j,k) = bc%fixed_Tpressure(4)/(((1+0.5*(gm-1.)*Mb*Mb))**(gm/(gm-1.)))
density(i,j,k) = gm*pressure(i,j,k)/(Cb*Cb)
end do
end do
end do
qp(:,jmx+1,:,:) = qp(:,jmx,:,:)
qp(:,jmx+2,:,:) = qp(:,jmx,:,:)
case("kmin")
do k = 1,1
do j = 1,jmx-1
do i = 1,imx-1
! interior cell
u = x_speed(i,j,k)
v = y_speed(i,j,k)
w = z_speed(i,j,k)
! ghost cell
uf = x_speed_inf!x_speed(i,j,k-1)
vf = y_speed_inf!y_speed(i,j,k-1)
wf = z_speed_inf!z_speed(i,j,k-1)
cexp = sqrt(gm*pressure(i,j,k)/density(i,j,k))
!cinf = sqrt(gm*pressure(i,j,k-1)/density(i,j,k-1))
cinf = sqrt(gm*pressure_inf/density_inf)
Unexp = u *(-Kfaces(i,j,k)%nx) + v *(-Kfaces(i,j,k)%ny) + w *(-Kfaces(i,j,k)%nz)
Uninf = uf*(-Kfaces(i,j,k)%nx) + vf*(-Kfaces(i,j,k)%ny) + wf*(-Kfaces(i,j,k)%nz)
Rinf = Uninf - 2*cinf/(gm-1.)
Rexp = Unexp + 2*cexp/(gm-1.)
Unb = 0.5*(Rexp + Rinf)
Cb = 0.25*(gm-1.)*(Rexp - Rinf)
if(Unb > 0.)then
vel_diff = Unb - Unexp
x_speed(i,j,k-1) = x_speed(i,j,k) + vel_diff*(-Kfaces(i,j,k)%nx)
y_speed(i,j,k-1) = y_speed(i,j,k) + vel_diff*(-Kfaces(i,j,k)%ny)
z_speed(i,j,k-1) = z_speed(i,j,k) + vel_diff*(-Kfaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call copy3(tv, "flat", face, bc, dims)
case('sst', 'sst2003')
call copy3(tk, "flat", face, bc, dims)
call copy3(tw, "flat", face, bc, dims)
case('kkl')
call copy3(tk, "flat", face, bc, dims)
call copy3(tkl, "flat", face, bc, dims)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
call copy3(tgm, "flat", face, bc, dims)
case DEFAULT
continue
end select
else
vel_diff = Unb - Uninf
x_speed(i,j,k-1) = x_speed_inf + vel_diff*(-Kfaces(i,j,k)%nx)
y_speed(i,j,k-1) = y_speed_inf + vel_diff*(-Kfaces(i,j,k)%ny)
z_speed(i,j,k-1) = z_speed_inf + vel_diff*(-Kfaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call fix(tv, bc%fixed_tv, face)
case('sst', 'sst2003')
!call check_if_value_fixed(bc, "sst")
call fix(tk, bc%fixed_tk, face)
call fix(tw, bc%fixed_tw, face)
case('kkl')
!call check_if_value_fixed(bc, "kkl")
call fix(tk, bc%fixed_tk, face)
call fix(tkl, bc%fixed_tkl, face)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
!call check_if_value_fixed(bc, "lctm2015")
call fix(tgm, bc%fixed_tgm, face)
case DEFAULT
continue
end select
end if
Mb = sqrt(x_speed(i,j,k)**2+y_speed(i,j,k)**2+z_speed(i,j,k)**2)/Cb
pressure(i,j,k-1) = bc%fixed_Tpressure(5)/(((1+0.5*(gm-1.)*Mb*Mb))**(gm/(gm-1.)))
density(i,j,k-1) = gm*pressure(i,j,k-1)/(Cb*Cb)
end do
end do
end do
qp(:,:,-1,:) = qp(:,:,0,:)
qp(:,:,-2,:) = qp(:,:,0,:)
case("kmax")
do k = kmx,kmx
do j = 1,jmx-1
do i = 1,imx-1
! interior cell
u = x_speed(i,j,k-1)
v = y_speed(i,j,k-1)
w = z_speed(i,j,k-1)
! ghost cell
uf = x_speed_inf!x_speed(i,j,k)
vf = y_speed_inf!y_speed(i,j,k)
wf = z_speed_inf!z_speed(i,j,k)
cexp = sqrt(gm*pressure(i,j,k-1)/density(i,j,k-1))
!cinf = sqrt(gm*pressure(i,j,k)/density(i,j,k))
cinf = sqrt(gm*pressure_inf/density_inf)
Unexp = u *(Kfaces(i,j,k)%nx) + v *(Kfaces(i,j,k)%ny) + w *(Kfaces(i,j,k)%nz)
Uninf = uf*(Kfaces(i,j,k)%nx) + vf*(Kfaces(i,j,k)%ny) + wf*(Kfaces(i,j,k)%nz)
Rinf = Uninf - 2*cinf/(gm-1.)
Rexp = Unexp + 2*cexp/(gm-1.)
Unb = 0.5*(Rexp + Rinf)
Cb = 0.25*(gm-1.)*(Rexp - Rinf)
if(Unb > 0.)then
vel_diff = Unb - Unexp
x_speed(i,j,k) = x_speed(i,j,k-1) + vel_diff*(Kfaces(i,j,k)%nx)
y_speed(i,j,k) = y_speed(i,j,k-1) + vel_diff*(Kfaces(i,j,k)%ny)
z_speed(i,j,k) = z_speed(i,j,k-1) + vel_diff*(Kfaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call copy3(tv, "flat", face, bc, dims)
case('sst', 'sst2003')
call copy3(tk, "flat", face, bc, dims)
call copy3(tw, "flat", face, bc, dims)
case('kkl')
call copy3(tk, "flat", face, bc, dims)
call copy3(tkl, "flat", face, bc, dims)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
call copy3(tgm, "flat", face, bc, dims)
case DEFAULT
continue
end select
else
vel_diff = Unb - Uninf
x_speed(i,j,k) = x_speed_inf + vel_diff*(Kfaces(i,j,k)%nx)
y_speed(i,j,k) = y_speed_inf + vel_diff*(Kfaces(i,j,k)%ny)
z_speed(i,j,k) = z_speed_inf + vel_diff*(Kfaces(i,j,k)%nz)
select case (turbulence)
case('none')
!do nothing
continue
case('sa', 'saBC')
call fix(tv, bc%fixed_tv, face)
case('sst', 'sst2003')
!call check_if_value_fixed(bc, "sst")
call fix(tk, bc%fixed_tk, face)
call fix(tw, bc%fixed_tw, face)
case('kkl')
!call check_if_value_fixed(bc, "kkl")
call fix(tk, bc%fixed_tk, face)
call fix(tkl, bc%fixed_tkl, face)
case DEFAULT
Fatal_error
end select
select case(trim(transition))
case('lctm2015')
!call check_if_value_fixed(bc, "lctm2015")
call fix(tgm, bc%fixed_tgm, face)
case DEFAULT
continue
end select
end if
Mb = sqrt(x_speed(i,j,k)**2+y_speed(i,j,k)**2+z_speed(i,j,k)**2)/Cb
pressure(i,j,k) = bc%fixed_Tpressure(6)/(((1+0.5*(gm-1.)*Mb*Mb))**(gm/(gm-1.)))
density(i,j,k) = gm*pressure(i,j,k)/(Cb*Cb)
end do
end do
end do
qp(:,:,kmx+1,:) = qp(:,:,kmx,:)
qp(:,:,kmx+2,:) = qp(:,:,kmx,:)
case DEFAULT
!print*, "ERROR: wrong face for boundary condition"
Fatal_error
end select
end subroutine total_pressure