update_SST_variables Subroutine

public subroutine update_SST_variables(qp, residue, delta_t, cells, Ifaces, Jfaces, Kfaces, dims)

Update the RANS (SST) equation with LU-SGS

Arguments

Type IntentOptional AttributesName
real(kind=wp), intent(inout), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var):: qp

Store primitive variable at cell center

real(kind=wp), intent(in), dimension(:, :, :, :):: residue

Store residue at each cell-center

real(kind=wp), intent(in), dimension(1:dims%imx-1, 1:dims%jmx-1, 1:dims%kmx-1):: delta_t

Local time increment value at each cell center

type(celltype), intent(in), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2):: cells

Input cell quantities: volume

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(extent), intent(in) :: dims

Extent of the domain:imx,jmx,kmx


Calls

proc~~update_sst_variables~~CallsGraph proc~update_sst_variables update_SST_variables proc~sstflux SSTFlux proc~update_sst_variables->proc~sstflux proc~spectralradius SpectralRadius proc~update_sst_variables->proc~spectralradius

Called by

proc~~update_sst_variables~~CalledByGraph proc~update_sst_variables update_SST_variables proc~update_with_plusgs update_with_plusgs proc~update_with_plusgs->proc~update_sst_variables proc~get_next_solution get_next_solution proc~get_next_solution->proc~update_with_plusgs proc~iterate_one_more_time_step iterate_one_more_time_step proc~iterate_one_more_time_step->proc~get_next_solution program~main main program~main->proc~iterate_one_more_time_step

Contents

Source Code


Source Code

    subroutine update_SST_variables(qp, residue, delta_t, cells, Ifaces, Jfaces, Kfaces, dims)
      !< Update the RANS (SST) equation with LU-SGS
      implicit none
      type(extent), intent(in) :: dims
      !< Extent of the domain:imx,jmx,kmx
      real(wp), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var), intent(inout) :: qp
      !< Store primitive variable at cell center
      real(wp), dimension(:, :, :, :), intent(in)  :: residue
      !< Store residue at each cell-center
      real(wp) , dimension(1:dims%imx-1, 1:dims%jmx-1, 1:dims%kmx-1), intent(in) :: delta_t
      !< Local time increment value at each cell center
      type(celltype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: cells
      !< Input cell quantities: volume
      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
      integer :: i,j,k
        real(wp), dimension(1:7)     :: deltaU
        real(wp), dimension(1:7)     :: D
        real(wp), dimension(1:7)     :: conservativeQ
        real(wp), dimension(1:7)     :: OldIminusFlux
        real(wp), dimension(1:7)     :: OldJminusFlux
        real(wp), dimension(1:7)     :: OldKminusFlux
        real(wp), dimension(1:7)     :: NewIminusFlux
        real(wp), dimension(1:7)     :: NewJminusFlux
        real(wp), dimension(1:7)     :: NewKminusFlux
        real(wp), dimension(1:7)     :: DelIminusFlux
        real(wp), dimension(1:7)     :: DelJminusFlux
        real(wp), dimension(1:7)     :: DelKminusFlux
        real(wp), dimension(1:6)     :: LambdaTimesArea
        real(wp), dimension(1:7)     :: Q0 ! state at cell
        real(wp), dimension(1:7)     :: Q1 ! state at neighbours 
        real(wp), dimension(1:7)     :: Q2
        real(wp), dimension(1:7)     :: Q3
        real(wp), dimension(1:7)     :: Q4
        real(wp), dimension(1:7)     :: Q5
        real(wp), dimension(1:7)     :: Q6
        real(wp), dimension(1:7)     :: DQ0! change in state
        real(wp), dimension(1:7)     :: DQ1
        real(wp), dimension(1:7)     :: DQ2
        real(wp), dimension(1:7)     :: DQ3
        real(wp), dimension(1:7)     :: DQ4
        real(wp), dimension(1:7)     :: DQ5
        real(wp), dimension(1:7)     :: DQ6
        real(wp), dimension(1:8)     :: Flist1
        real(wp), dimension(1:8)     :: Flist2
        real(wp), dimension(1:8)     :: Flist3
        real(wp), dimension(1:8)     :: Flist4
        real(wp), dimension(1:8)     :: Flist5
        real(wp), dimension(1:8)     :: Flist6
        real(wp), dimension(1:3)     :: C0
        real(wp), dimension(1:3)     :: C1
        real(wp), dimension(1:3)     :: C2
        real(wp), dimension(1:3)     :: C3
        real(wp), dimension(1:3)     :: C4
        real(wp), dimension(1:3)     :: C5
        real(wp), dimension(1:3)     :: C6
        real(wp)                     :: beta
        real(wp)                     :: eps
        real(wp)                     :: M
        real(wp)                     :: VMag
        real(wp)                     :: SoundMag
        real(wp)                     :: u,v,w,r,p,kk,ww,H
        real(wp)                     :: factor
        real(wp), dimension(1:7,1:7) :: PrecondInv

        ! intermittency
        real(wp) :: De, Dp

        De = 0.0
        Dp = 0.0


        !intialize delQ
        delQstar = 0.0

        !forward sweep
        do k=1,dims%kmx-1
          do j=1,dims%jmx-1
            do i=1,dims%imx-1
              C0  = (/Cells(i  ,j  ,k  )%Centerx,Cells(i  ,j  ,k  )%Centery,Cells(i  ,j  ,k  )%Centerz/)
              C1  = (/Cells(i-1,j  ,k  )%Centerx,Cells(i-1,j  ,k  )%Centery,Cells(i-1,j  ,k  )%Centerz/)
              C2  = (/Cells(i  ,j-1,k  )%Centerx,Cells(i  ,j-1,k  )%Centery,Cells(i  ,j-1,k  )%Centerz/)
              C3  = (/Cells(i  ,j  ,k-1)%Centerx,Cells(i  ,j  ,k-1)%Centery,Cells(i  ,j  ,k-1)%Centerz/)
              C4  = (/Cells(i+1,j  ,k  )%Centerx,Cells(i+1,j  ,k  )%Centery,Cells(i+1,j  ,k  )%Centerz/)
              C5  = (/Cells(i  ,j+1,k  )%Centerx,Cells(i  ,j+1,k  )%Centery,Cells(i  ,j+1,k  )%Centerz/)
              C6  = (/Cells(i  ,j  ,k+1)%Centerx,Cells(i  ,j  ,k+1)%Centery,Cells(i  ,j  ,k+1)%Centerz/)

              Q0  = qp(i  , j  , k  , 1:7)
              Q1  = qp(i-1, j  , k  , 1:7)
              Q2  = qp(i  , j-1, k  , 1:7)
              Q3  = qp(i  , j  , k-1, 1:7)
              Q4  = qp(i+1, j  , k  , 1:7)
              Q5  = qp(i  , j+1, k  , 1:7)
              Q6  = qp(i  , j  , k+1, 1:7)

              DQ0 = 0.0
              DQ1 = delQstar(i-1, j  , k  , 1:7)
              DQ2 = delQstar(i  , j-1, k  , 1:7)
              DQ3 = delQstar(i  , j  , k-1, 1:7)

              Flist1(1) =  Ifaces(i,j,k)%A
              Flist1(2) = -Ifaces(i,j,k)%nx
              Flist1(3) = -Ifaces(i,j,k)%ny
              Flist1(4) = -Ifaces(i,j,k)%nz
              Flist1(5) = 0.5*(cells(i-1, j  , k  )%volume + cells(i,j,k)%volume)
              Flist1(6) = 0.5*(   mmu(i-1, j  , k  ) +    mmu(i,j,k))
              Flist1(7) = 0.5*(   tmu(i-1, j  , k  ) +    tmu(i,j,k))
              Flist1(8) = 0.5*(sst_F1(i-1, j  , k  ) + sst_F1(i,j,k))

              Flist2(1) =  Jfaces(i,j,k)%A
              Flist2(2) = -Jfaces(i,j,k)%nx
              Flist2(3) = -Jfaces(i,j,k)%ny
              Flist2(4) = -Jfaces(i,j,k)%nz
              Flist2(5) = 0.5*(cells(i  , j-1, k  )%volume + cells(i,j,k)%volume)
              Flist2(6) = 0.5*(   mmu(i  , j-1, k  ) +    mmu(i,j,k))
              Flist2(7) = 0.5*(   tmu(i  , j-1, k  ) +    tmu(i,j,k))
              Flist2(8) = 0.5*(sst_F1(i  , j-1, k  ) + sst_F1(i,j,k))

              Flist3(1) =  Kfaces(i,j,k)%A
              Flist3(2) = -Kfaces(i,j,k)%nx
              Flist3(3) = -Kfaces(i,j,k)%ny
              Flist3(4) = -Kfaces(i,j,k)%nz
              Flist3(5) = 0.5*(cells(i  , j  , k-1)%volume + cells(i,j,k)%volume)
              Flist3(6) = 0.5*(   mmu(i  , j  , k-1) +    mmu(i,j,k))
              Flist3(7) = 0.5*(   tmu(i  , j  , k-1) +    tmu(i,j,k))
              Flist3(8) = 0.5*(sst_F1(i  , j  , k-1) + sst_F1(i,j,k))

              Flist4(1) =  Ifaces(i+1,j,k)%A
              Flist4(2) = +Ifaces(i+1,j,k)%nx
              Flist4(3) = +Ifaces(i+1,j,k)%ny
              Flist4(4) = +Ifaces(i+1,j,k)%nz
              Flist4(5) = 0.5*(cells(i+1, j  , k  )%volume + cells(i,j,k)%volume)
              Flist4(6) = 0.5*(   mmu(i+1, j  , k  ) +    mmu(i,j,k))
              Flist4(7) = 0.5*(   tmu(i+1, j  , k  ) +    tmu(i,j,k))
              Flist4(8) = 0.5*(sst_F1(i+1, j  , k  ) + sst_F1(i,j,k))

              Flist5(1) =  Jfaces(i,j+1,k)%A
              Flist5(2) = +Jfaces(i,j+1,k)%nx
              Flist5(3) = +Jfaces(i,j+1,k)%ny
              Flist5(4) = +Jfaces(i,j+1,k)%nz
              Flist5(5) = 0.5*(cells(i  , j+1, k  )%volume + cells(i,j,k)%volume)
              Flist5(6) = 0.5*(   mmu(i  , j+1, k  ) +    mmu(i,j,k))
              Flist5(7) = 0.5*(   tmu(i  , j+1, k  ) +    tmu(i,j,k))
              Flist5(8) = 0.5*(sst_F1(i  , j+1, k  ) + sst_F1(i,j,k))

              Flist6(1) =  Kfaces(i,j,k+1)%A
              Flist6(2) = +Kfaces(i,j,k+1)%nx
              Flist6(3) = +Kfaces(i,j,k+1)%ny
              Flist6(4) = +Kfaces(i,j,k+1)%nz
              Flist6(5) = 0.5*(cells(i  , j  , k+1)%volume + cells(i,j,k)%volume)
              Flist6(6) = 0.5*(   mmu(i  , j  , k+1) +    mmu(i,j,k))
              Flist6(7) = 0.5*(   tmu(i  , j  , k+1) +    tmu(i,j,k))
              Flist6(8) = 0.5*(sst_F1(i  , j  , k+1) + sst_F1(i,j,k))

              NewIminusFlux     = SSTFlux(Q1, Q0, DQ1, Flist1)
              NewJminusFlux     = SSTFlux(Q2, Q0, DQ2, Flist2)
              NewKminusFlux     = SSTFlux(Q3, Q0, DQ3, Flist3)
              OldIminusFlux     = SSTFlux(Q1, Q0, DQ0, Flist1)
              OldJminusFlux     = SSTFlux(Q2, Q0, DQ0, Flist2)
              OldKminusFlux     = SSTFlux(Q3, Q0, DQ0, Flist3)

              !---preconditioning---
              r  = Q0(1)
              u  = Q0(2)
              v  = Q0(3)
              w  = Q0(4)
              p  = Q0(5)
              kk = Q0(6)
              ww = Q0(7)
              VMag     = sqrt(u*u + v*v + w*w)
              SoundMag = sqrt(gm*p/r)
              M        = VMag/SoundMag 
              H  = (gm*p/(r*(gm-1.0))) + 0.5*(VMag)
              eps = min(1.0, max(M*M, Minf*Minf))
              factor = (1.0-eps)*(gm-1.0)/(SoundMag*SoundMag)

              LambdaTimesArea(1)= SpectralRadius(Q1, Q0, Flist1, C1, C0, eps)
              LambdaTimesArea(2)= SpectralRadius(Q2, Q0, Flist2, C2, C0, eps)
              LambdaTimesArea(3)= SpectralRadius(Q3, Q0, Flist3, C3, C0, eps)
              LambdaTimesArea(4)= SpectralRadius(Q4, Q0, Flist4, C4, C0, eps)
              LambdaTimesArea(5)= SpectralRadius(Q5, Q0, Flist5, C5, C0, eps)
              LambdaTimesArea(6)= SpectralRadius(Q6, Q0, Flist6, C6, C0, eps)


              PrecondInv(1,1) = 1.0 - factor*1*VMag*VMag/2.0
              PrecondInv(2,1) = 0.0 - factor*u*VMag*VMag/2.0
              PrecondInv(3,1) = 0.0 - factor*v*VMag*VMag/2.0
              PrecondInv(4,1) = 0.0 - factor*w*VMag*VMag/2.0
              PrecondInv(5,1) = 0.0 - factor*H*VMag*VMag/2.0
              PrecondInv(6,1) = 0.0 - factor*kk*VMag*VMag/2.0
              PrecondInv(7,1) = 0.0 - factor*ww*VMag*VMag/2.0
              PrecondInv(1,2) = 0.0 - factor*1*(-u)
              PrecondInv(2,2) = 1.0 - factor*u*(-u)
              PrecondInv(3,2) = 0.0 - factor*v*(-u)
              PrecondInv(4,2) = 0.0 - factor*w*(-u)
              PrecondInv(5,2) = 0.0 - factor*H*(-u)
              PrecondInv(6,2) = 0.0 - factor*kk*(-u)
              PrecondInv(7,2) = 0.0 - factor*ww*(-u)
              PrecondInv(1,3) = 0.0 - factor*1*(-v)
              PrecondInv(2,3) = 0.0 - factor*u*(-v)
              PrecondInv(3,3) = 1.0 - factor*v*(-v)
              PrecondInv(4,3) = 0.0 - factor*w*(-v)
              PrecondInv(5,3) = 0.0 - factor*H*(-v)
              PrecondInv(6,3) = 0.0 - factor*kk*(-v)
              PrecondInv(7,3) = 0.0 - factor*ww*(-v)
              PrecondInv(1,4) = 0.0 - factor*1*(-w)
              PrecondInv(2,4) = 0.0 - factor*u*(-w)
              PrecondInv(3,4) = 0.0 - factor*v*(-w)
              PrecondInv(4,4) = 1.0 - factor*w*(-w)
              PrecondInv(5,4) = 0.0 - factor*H*(-w)
              PrecondInv(6,4) = 0.0 - factor*kk*(-w)
              PrecondInv(7,4) = 0.0 - factor*ww*(-w)
              PrecondInv(1,5) = 0.0 - factor*1*(1.)
              PrecondInv(2,5) = 0.0 - factor*u*(1.)
              PrecondInv(3,5) = 0.0 - factor*v*(1.)
              PrecondInv(4,5) = 0.0 - factor*w*(1.)
              PrecondInv(5,5) = 1.0 - factor*H*(1.)
              PrecondInv(6,5) = 0.0 - factor*kk*(1.)
              PrecondInv(7,5) = 0.0 - factor*ww*(1.)
              PrecondInv(1,6) = 0.0 - factor*1*(-1.)
              PrecondInv(2,6) = 0.0 - factor*u*(-1.)
              PrecondInv(3,6) = 0.0 - factor*v*(-1.)
              PrecondInv(4,6) = 0.0 - factor*w*(-1.)
              PrecondInv(5,6) = 0.0 - factor*H*(-1.)
              PrecondInv(6,6) = 1.0 - factor*kk*(-1.)
              PrecondInv(7,6) = 0.0 - factor*ww*(-1.)
              PrecondInv(1,7) = 0.0 - factor*1*(0.)
              PrecondInv(2,7) = 0.0 - factor*u*(0.)
              PrecondInv(3,7) = 0.0 - factor*v*(0.)
              PrecondInv(4,7) = 0.0 - factor*w*(0.)
              PrecondInv(5,7) = 0.0 - factor*H*(0.)
              PrecondInv(6,7) = 0.0 - factor*kk*(0.)
              PrecondInv(7,7) = 1.0 - factor*ww*(0.)
              !---end preconditioning


              ! multiply above flux with area to get correct values
              DelIminusFlux =  NewIminusFlux - OldIminusFlux
              DelJminusFlux =  NewJminusFlux - OldJminusFlux
              DelKminusFlux =  NewKminusFlux - OldKminusFlux


              D = (cells(i,j,k)%volume/delta_t(i,j,k)) + 0.5*SUM(LambdaTimesArea)
              beta = sst_F1(i,j,k)*beta1 + (1.0-sst_F1(i,j,k))*beta2
              D(6) = (D(6) + (bstar*qp(i,j,k,7))*cells(i,j,k)%volume)
              D(7) = (D(7) + 2.0*beta*qp(i,j,k,7)*cells(i,j,k)%volume)
              !storing D in Iflux array for backward sweep
              !F_p(i,j,k,1) = D
              

              deltaU(1:7) = -matmul(PrecondInv,residue(i,j,k,1:7)) &
                - 0.5*((matmul(PrecondInv,DelIminusFlux) - LambdaTimesArea(1)*delQstar(i-1,j,k,1:7)) &
                     + (matmul(PrecondInv,DelJminusFlux) - LambdaTimesArea(2)*delQstar(i,j-1,k,1:7)) &
                     + (matmul(PrecondInv,DelKminusFlux) - LambdaTimesArea(3)*delQstar(i,j,k-1,1:7)) )

              delQstar(i,j,k,1:7) = deltaU(1:7)/D
            end do
          end do
        end do

        delQ=0.0
        !backward sweep
            do i=dims%imx-1,1,-1
          do j=dims%jmx-1,1,-1
        do k=dims%kmx-1,1,-1
              C0  = (/Cells(i  ,j  ,k  )%Centerx,Cells(i  ,j  ,k  )%Centery,Cells(i  ,j  ,k  )%Centerz/)
              C1  = (/Cells(i-1,j  ,k  )%Centerx,Cells(i-1,j  ,k  )%Centery,Cells(i-1,j  ,k  )%Centerz/)
              C2  = (/Cells(i  ,j-1,k  )%Centerx,Cells(i  ,j-1,k  )%Centery,Cells(i  ,j-1,k  )%Centerz/)
              C3  = (/Cells(i  ,j  ,k-1)%Centerx,Cells(i  ,j  ,k-1)%Centery,Cells(i  ,j  ,k-1)%Centerz/)
              C4  = (/Cells(i+1,j  ,k  )%Centerx,Cells(i+1,j  ,k  )%Centery,Cells(i+1,j  ,k  )%Centerz/)
              C5  = (/Cells(i  ,j+1,k  )%Centerx,Cells(i  ,j+1,k  )%Centery,Cells(i  ,j+1,k  )%Centerz/)
              C6  = (/Cells(i  ,j  ,k+1)%Centerx,Cells(i  ,j  ,k+1)%Centery,Cells(i  ,j  ,k+1)%Centerz/)

              Q0  = qp(i  , j  , k  , 1:7)
              Q1  = qp(i-1, j  , k  , 1:7)
              Q2  = qp(i  , j-1, k  , 1:7)
              Q3  = qp(i  , j  , k-1, 1:7)
              Q4  = qp(i+1, j  , k  , 1:7)
              Q5  = qp(i  , j+1, k  , 1:7)
              Q6  = qp(i  , j  , k+1, 1:7)

              DQ0 = 0.0
              DQ4 = delQ(i+1, j  , k  , 1:7)
              DQ5 = delQ(i  , j+1, k  , 1:7)
              DQ6 = delQ(i  , j  , k+1, 1:7)

              Flist1(1) =  Ifaces(i,j,k)%A
              Flist1(2) = -Ifaces(i,j,k)%nx
              Flist1(3) = -Ifaces(i,j,k)%ny
              Flist1(4) = -Ifaces(i,j,k)%nz
              Flist1(5) = 0.5*(cells(i-1, j  , k  )%volume + cells(i,j,k)%volume)
              Flist1(6) = 0.5*(   mmu(i-1, j  , k  ) +    mmu(i,j,k))
              Flist1(7) = 0.5*(   tmu(i-1, j  , k  ) +    tmu(i,j,k))
              Flist1(8) = 0.5*(sst_F1(i-1, j  , k  ) + sst_F1(i,j,k))

              Flist2(1) =  Jfaces(i,j,k)%A
              Flist2(2) = -Jfaces(i,j,k)%nx
              Flist2(3) = -Jfaces(i,j,k)%ny
              Flist2(4) = -Jfaces(i,j,k)%nz
              Flist2(5) = 0.5*(cells(i  , j-1, k  )%volume + cells(i,j,k)%volume)
              Flist2(6) = 0.5*(   mmu(i  , j-1, k  ) +    mmu(i,j,k))
              Flist2(7) = 0.5*(   tmu(i  , j-1, k  ) +    tmu(i,j,k))
              Flist2(8) = 0.5*(sst_F1(i  , j-1, k  ) + sst_F1(i,j,k))

              Flist3(1) =  Kfaces(i,j,k)%A
              Flist3(2) = -Kfaces(i,j,k)%nx
              Flist3(3) = -Kfaces(i,j,k)%ny
              Flist3(4) = -Kfaces(i,j,k)%nz
              Flist3(5) = 0.5*(cells(i  , j  , k-1)%volume + cells(i,j,k)%volume)
              Flist3(6) = 0.5*(   mmu(i  , j  , k-1) +    mmu(i,j,k))
              Flist3(7) = 0.5*(   tmu(i  , j  , k-1) +    tmu(i,j,k))
              Flist3(8) = 0.5*(sst_F1(i  , j  , k-1) + sst_F1(i,j,k))

              Flist4(1) =  Ifaces(i+1,j,k)%A
              Flist4(2) = +Ifaces(i+1,j,k)%nx
              Flist4(3) = +Ifaces(i+1,j,k)%ny
              Flist4(4) = +Ifaces(i+1,j,k)%nz
              Flist4(5) = 0.5*(cells(i+1, j  , k  )%volume + cells(i,j,k)%volume)
              Flist4(6) = 0.5*(   mmu(i+1, j  , k  ) +    mmu(i,j,k))
              Flist4(7) = 0.5*(   tmu(i+1, j  , k  ) +    tmu(i,j,k))
              Flist4(8) = 0.5*(sst_F1(i+1, j  , k  ) + sst_F1(i,j,k))

              Flist5(1) =  Jfaces(i,j+1,k)%A
              Flist5(2) = +Jfaces(i,j+1,k)%nx
              Flist5(3) = +Jfaces(i,j+1,k)%ny
              Flist5(4) = +Jfaces(i,j+1,k)%nz
              Flist5(5) = 0.5*(cells(i  , j+1, k  )%volume + cells(i,j,k)%volume)
              Flist5(6) = 0.5*(   mmu(i  , j+1, k  ) +    mmu(i,j,k))
              Flist5(7) = 0.5*(   tmu(i  , j+1, k  ) +    tmu(i,j,k))
              Flist5(8) = 0.5*(sst_F1(i  , j+1, k  ) + sst_F1(i,j,k))

              Flist6(1) =  Kfaces(i,j,k+1)%A
              Flist6(2) = +Kfaces(i,j,k+1)%nx
              Flist6(3) = +Kfaces(i,j,k+1)%ny
              Flist6(4) = +Kfaces(i,j,k+1)%nz
              Flist6(5) = 0.5*(cells(i  , j  , k+1)%volume + cells(i,j,k)%volume)
              Flist6(6) = 0.5*(   mmu(i  , j  , k+1) +    mmu(i,j,k))
              Flist6(7) = 0.5*(   tmu(i  , j  , k+1) +    tmu(i,j,k))
              Flist6(8) = 0.5*(sst_F1(i  , j  , k+1) + sst_F1(i,j,k))

              NewIminusFlux     = SSTFlux(Q4, Q0, DQ4, Flist4)
              NewJminusFlux     = SSTFlux(Q5, Q0, DQ5, Flist5)
              NewKminusFlux     = SSTFlux(Q6, Q0, DQ6, Flist6)
              OldIminusFlux     = SSTFlux(Q4, Q0, DQ0, Flist4)
              OldJminusFlux     = SSTFlux(Q5, Q0, DQ0, Flist5)
              OldKminusFlux     = SSTFlux(Q6, Q0, DQ0, Flist6)

              !---preconditioning---
              r  = Q0(1)
              u  = Q0(2)
              v  = Q0(3)
              w  = Q0(4)
              p  = Q0(5)
              kk = Q0(6)
              ww = Q0(7)
              VMag     = sqrt(u*u + v*v + w*w)
              SoundMag = sqrt(gm*p/r)
              M        = VMag/SoundMag 
              H  = (gm*p/(r*(gm-1.0))) + 0.5*(VMag)
              eps = min(1.0, max(M*M, Minf*Minf))
              factor = (1.0-eps)*(gm-1.0)/(SoundMag*SoundMag)

              LambdaTimesArea(1)= SpectralRadius(Q1, Q0, Flist1, C1, C0, eps)
              LambdaTimesArea(2)= SpectralRadius(Q2, Q0, Flist2, C2, C0, eps)
              LambdaTimesArea(3)= SpectralRadius(Q3, Q0, Flist3, C3, C0, eps)
              LambdaTimesArea(4)= SpectralRadius(Q4, Q0, Flist4, C4, C0, eps)
              LambdaTimesArea(5)= SpectralRadius(Q5, Q0, Flist5, C5, C0, eps)
              LambdaTimesArea(6)= SpectralRadius(Q6, Q0, Flist6, C6, C0, eps)


              PrecondInv(1,1) = 1.0 - factor*1*VMag*VMag/2.0
              PrecondInv(2,1) = 0.0 - factor*u*VMag*VMag/2.0
              PrecondInv(3,1) = 0.0 - factor*v*VMag*VMag/2.0
              PrecondInv(4,1) = 0.0 - factor*w*VMag*VMag/2.0
              PrecondInv(5,1) = 0.0 - factor*H*VMag*VMag/2.0
              PrecondInv(6,1) = 0.0 - factor*kk*VMag*VMag/2.0
              PrecondInv(7,1) = 0.0 - factor*ww*VMag*VMag/2.0
              PrecondInv(1,2) = 0.0 - factor*1*(-u)
              PrecondInv(2,2) = 1.0 - factor*u*(-u)
              PrecondInv(3,2) = 0.0 - factor*v*(-u)
              PrecondInv(4,2) = 0.0 - factor*w*(-u)
              PrecondInv(5,2) = 0.0 - factor*H*(-u)
              PrecondInv(6,2) = 0.0 - factor*kk*(-u)
              PrecondInv(7,2) = 0.0 - factor*ww*(-u)
              PrecondInv(1,3) = 0.0 - factor*1*(-v)
              PrecondInv(2,3) = 0.0 - factor*u*(-v)
              PrecondInv(3,3) = 1.0 - factor*v*(-v)
              PrecondInv(4,3) = 0.0 - factor*w*(-v)
              PrecondInv(5,3) = 0.0 - factor*H*(-v)
              PrecondInv(6,3) = 0.0 - factor*kk*(-v)
              PrecondInv(7,3) = 0.0 - factor*ww*(-v)
              PrecondInv(1,4) = 0.0 - factor*1*(-w)
              PrecondInv(2,4) = 0.0 - factor*u*(-w)
              PrecondInv(3,4) = 0.0 - factor*v*(-w)
              PrecondInv(4,4) = 1.0 - factor*w*(-w)
              PrecondInv(5,4) = 0.0 - factor*H*(-w)
              PrecondInv(6,4) = 0.0 - factor*kk*(-w)
              PrecondInv(7,4) = 0.0 - factor*ww*(-w)
              PrecondInv(1,5) = 0.0 - factor*1*(1.)
              PrecondInv(2,5) = 0.0 - factor*u*(1.)
              PrecondInv(3,5) = 0.0 - factor*v*(1.)
              PrecondInv(4,5) = 0.0 - factor*w*(1.)
              PrecondInv(5,5) = 1.0 - factor*H*(1.)
              PrecondInv(6,5) = 0.0 - factor*kk*(1.)
              PrecondInv(7,5) = 0.0 - factor*ww*(1.)
              PrecondInv(1,6) = 0.0 - factor*1*(-1.)
              PrecondInv(2,6) = 0.0 - factor*u*(-1.)
              PrecondInv(3,6) = 0.0 - factor*v*(-1.)
              PrecondInv(4,6) = 0.0 - factor*w*(-1.)
              PrecondInv(5,6) = 0.0 - factor*H*(-1.)
              PrecondInv(6,6) = 1.0 - factor*kk*(-1.)
              PrecondInv(7,6) = 0.0 - factor*ww*(-1.)
              PrecondInv(1,7) = 0.0 - factor*1*(0.)
              PrecondInv(2,7) = 0.0 - factor*u*(0.)
              PrecondInv(3,7) = 0.0 - factor*v*(0.)
              PrecondInv(4,7) = 0.0 - factor*w*(0.)
              PrecondInv(5,7) = 0.0 - factor*H*(0.)
              PrecondInv(6,7) = 0.0 - factor*kk*(0.)
              PrecondInv(7,7) = 1.0 - factor*ww*(0.)
              !---end preconditioning


              ! multiply above flux with area to get correct values
              DelIminusFlux =  NewIminusFlux - OldIminusFlux
              DelJminusFlux =  NewJminusFlux - OldJminusFlux
              DelKminusFlux =  NewKminusFlux - OldKminusFlux

              D = (cells(i,j,k)%volume/delta_t(i,j,k)) + 0.5*SUM(LambdaTimesArea)
              beta = sst_F1(i,j,k)*beta1 + (1.0-sst_F1(i,j,k))*beta2
              D(6) = (D(6) + (bstar*qp(i,j,k,7))*cells(i,j,k)%volume)
              D(7) = (D(7) + 2.0*beta*qp(i,j,k,7)*cells(i,j,k)%volume)


              delQ(i,j,k,1:7) = delQstar(i,j,k,1:7) &
                - 0.5*((matmul(PrecondInv,DelIminusFlux) - LambdaTimesArea(4)*delQ(i+1,j,k,1:7)) &
                     + (matmul(PrecondInv,DelJminusFlux) - LambdaTimesArea(5)*delQ(i,j+1,k,1:7)) &
                     + (matmul(PrecondInv,DelKminusFlux) - LambdaTimesArea(6)*delQ(i,j,k+1,1:7)) )/D

            end do
          end do
        end do

        
        do k=1,dims%kmx-1
          do j = 1,dims%jmx-1
            do i = 1,dims%imx-1
              conservativeQ(1) = qp(i,j,k,1)
              conservativeQ(2) = qp(i,j,k,1) * qp(i,j,k,2)
              conservativeQ(3) = qp(i,j,k,1) * qp(i,j,k,3)
              conservativeQ(4) = qp(i,j,k,1) * qp(i,j,k,4)
              conservativeQ(5) = (qp(i,j,k,5) / (gm-1.0)) + ( 0.5 * qp(i,j,k,1) * sum( qp(i,j,k,2:4)**2) )
              conservativeQ(6) = qp(i,j,k,1) * qp(i,j,k,6)
              conservativeQ(7) = qp(i,j,k,1) * qp(i,j,k,7)
              
              ! add new change into conservative solution
              conservativeQ(1:7) = conservativeQ(1:7) + delQ(i,j,k,1:7)

              ! convert back conservative to primitive
              qp(i,j,k,1) = conservativeQ(1)
              qp(i,j,k,2) = conservativeQ(2) / conservativeQ(1)
              qp(i,j,k,3) = conservativeQ(3) / conservativeQ(1)
              qp(i,j,k,4) = conservativeQ(4) / conservativeQ(1)
              qp(i,j,k,5) = (gm-1.0) * ( conservativeQ(5) - (0.5 * sum(conservativeQ(2:4)**2) / conservativeQ(1)) )
              if(conservativeQ(6)>0)then
              qp(i,j,k,6) = conservativeQ(6) / conservativeQ(1)
              end if
              if(conservativeQ(7)>0)then
              qp(i,j,k,7) = conservativeQ(7) / conservativeQ(1)
              end if
            end do
          end do
        end do

    end subroutine update_SST_variables