add_sst_source Subroutine

    subroutine add_sst_source(qp, residue, cells,scheme,dims)
      !< Add residual due to source terms of the SST turbulence model
      implicit none
      type(extent), intent(in) :: dims
      !< Extent of the domain:imx,jmx,kmx
      type(schemetype), intent(in) :: scheme
      !< finite-volume Schemes
      real(wp), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var), intent(in) :: qp
      !< Store primitive variable at cell center
      real(wp), dimension(:, :, :, :), intent(inout)  :: residue
      !< Store residue 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
      !< Store residue at each cell-center
      integer :: i,j,k

      real(wp) :: CD
      !< cross diffusion term
      real(wp) :: F1
      !< single cell belding fuction 
      real(wp) :: vort
      !< vorticity magnitude
      real(wp) :: S_k
      !< Total source term of TKE equation
      real(wp) :: S_w
      !< Total source term of omega equation
      real(wp) :: D_k
      !< destruction term of TKE equation
      real(wp) :: D_w
      !< destruction term of omega equation
      real(wp) :: P_k
      !< production term of TKE equation
      real(wp) :: P_w
      !< production term of Omega equation
      real(wp) :: lamda
      !< additional source term in Omega equation
      integer :: limiter
      !< production term limiter
      real(wp) :: divergence
      !< del.V
      real(wp) :: density
      !< single cell density
      real(wp) :: tk
      !< single cell TKE
      real(wp) :: tw
      !< single cell Omega
      
      if(trim(scheme%turbulence) == 'sst2003')then
        limiter = 10
        gama1 = 5.0/9.0
        gama2 = 0.44
      else 
        limiter = 20
      end if


      do k = 1,dims%kmx-1
        do j = 1,dims%jmx-1
          do i = 1,dims%imx-1

            density = qp(i,j,k,1)
            tk      = qp(i,j,k,6)
            tw      = qp(i,j,k,7)

            ! __ vorticity __
            vort = sqrt(     ((gradw_y(i,j,k)- gradv_z(i,j,k))**2 &
                            + (gradu_z(i,j,k)- gradw_x(i,j,k))**2 &
                            + (gradv_x(i,j,k)- gradu_y(i,j,k))**2 &
                             )&
                       )

            CD = 2*density*sigma_w2*(gradtk_x(i,j,k)*gradtw_x(i,j,k)&
                                   + gradtk_y(i,j,k)*gradtw_y(i,j,k)&
                                   + gradtk_z(i,j,k)*gradtw_z(i,j,k)&
                                    )/tw
            CD = max(CD, 10.0**(-limiter))
            F1 = sst_F1(i,j,k)


            sigma_k     =    sigma_k1*F1  +    sigma_k2*(1. - F1)
            sigma_w     =    sigma_w1*F1  +    sigma_w1*(1. - F1)
            gama        =       gama1*F1  +       gama2*(1. - F1)
            beta        =       beta1*F1  +       beta2*(1. - F1)



            ! ____ Dissipation term ___
            D_k = bstar*density*tw*tk
            D_w = beta*density*tw**2

            ! ____ PRODUCTION term____ 
            divergence = gradu_x(i,j,k) + gradv_y(i,j,k) + gradw_z(i,j,k)
            P_k = mu_t(i,j,k)*(vort**2) -((2.0/3.0)*density*tk*divergence)
            P_k = min(P_k,limiter*D_k)
            P_w = (density*gama/mu_t(i,j,k))*P_k

            ! ____ cross diffusion term ___
            lamda = (1. - F1)*CD

            S_k = P_k - D_k           !Source term TKE
            S_w = P_w - D_w  +lamda   !source term omega

            S_k = S_k * cells(i, j, k)%volume
            S_w = S_w * cells(i, j, k)%volume

            residue(i, j, k, 6) = residue(i, j, k, 6) - S_k
            residue(i, j, k, 7) = residue(i, j, k, 7) - S_w

          end do
        end do
      end do

    end subroutine add_sst_source