compute_viscous_fluxes_kkl Subroutine

    subroutine compute_viscous_fluxes_kkl(F, qp, cells, faces, flags, dims)
      !< Compute viscous fluxes for additianal equations due to k-kL turbulence model
      implicit none
      real(wp), dimension(:, :, :, :), intent(inout) :: F 
      !< Flux array
      type(extent), intent(in) :: dims
      !< Extent of the domain: imx,jmx,kmx
      integer, dimension(3), intent(in) :: flags
      !< flags for directions switch
      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
      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+2+flags(1),-2:dims%jmx+2+flags(2),-2:dims%kmx+2+flags(3)), intent(in) :: faces
      !< Face quantities: area and unit normal
      ! local variables
      real(wp) :: tkface
      real(wp) :: rhoface
      real(wp) :: normal_comp
      real(wp) :: d_LR
      real(wp) :: mu_f
      real(wp) :: mut_f
      real(wp) :: total_mu
      real(wp) :: delx
      real(wp) :: dely
      real(wp) :: delz
      real(wp) :: deltk
      real(wp) :: deltkl
      real(wp) :: Tau_xx
      real(wp) :: Tau_yy
      real(wp) :: Tau_zz
      real(wp) :: nx
      real(wp) :: ny
      real(wp) :: nz
      real(wp) :: area
      integer :: i, j, k
      integer :: ii, jj, kk
      !--- kkl variable requirement  ---!
      real(wp) :: dtkdx, dtkdy, dtkdz
      real(wp) :: dtkldx, dtkldy, dtkldz


      ii = flags(1)
      jj = flags(2)
      kk = flags(3)


      !---------------------------------------------------------------------
      ! Calculating the turbulent viscous fluxes at the faces
      !--------------------------------------------------------------------
      do k = 1, dims%kmx - 1 + kk
       do j = 1, dims%jmx - 1 + jj
        do i = 1, dims%imx - 1 + ii

          !--- FACE Gradients ---!
          ! Gradients at face as average of gradients at cell centres
          dtkdx  = 0.5 * (gradtk_x(i-ii, j-jj, k-kk) + gradtk_x(i, j, k))
          dtkdy  = 0.5 * (gradtk_y(i-ii, j-jj, k-kk) + gradtk_y(i, j, k))
          dtkdz  = 0.5 * (gradtk_z(i-ii, j-jj, k-kk) + gradtk_z(i, j, k))
          dtkldx = 0.5 * (gradtkl_x(i-ii, j-jj, k-kk) + gradtkl_x(i, j, k))
          dtkldy = 0.5 * (gradtkl_y(i-ii, j-jj, k-kk) + gradtkl_y(i, j, k))
          dtkldz = 0.5 * (gradtkl_z(i-ii, j-jj, k-kk) + gradtkl_z(i, j, k))
          !--- For ODD-EVEN coupling error ---!
          ! distance between cell center of adjacent cell for the i,j,k face
          delx = cells(i, j, k)%centerx - cells(i-ii, j-jj, k-kk)%centerx
          dely = cells(i, j, k)%centery - cells(i-ii, j-jj, k-kk)%centery
          delz = cells(i, j, k)%centerz - cells(i-ii, j-jj, k-kk)%centerz

          d_LR = sqrt(delx*delx + dely*dely + delz*delz)

          ! difference in state across face
          deltk  = qp(i, j, k, 6) - qp(i-ii, j-jj, k-kk, 6)   !TKE
          deltkl = qp(i, j, k, 7) - qp(i-ii, j-jj, k-kk, 7)   !Kl

          !normal_comp   = ( delta(phi) - (grad(phi).dot.delR) )/magnitudeR
          !new grad(phi) =  grad(phi) + correction(normal_comp.dot.delR/magnitudeR)

          normal_comp = (deltk - (dtkdx*delx + dtkdy*dely + dtkdz*delz))/d_LR
          dtkdx       =  dtkdx + (normal_comp * delx / d_LR)
          dtkdy       =  dtkdy + (normal_comp * dely / d_LR)
          dtkdz       =  dtkdz + (normal_comp * delz / d_LR)

          normal_comp  = (deltkl - (dtkldx*delx + dtkldy*dely + dtkldz*delz))/d_LR
          dtkldx       =  dtkldx + (normal_comp * delx / d_LR)
          dtkldy       =  dtkldy + (normal_comp * dely / d_LR)
          dtkldz       =  dtkldz + (normal_comp * delz / d_LR)
          !--- end of ODD-EVEN coupling correction ---!

          mu_f     = 0.5*(  mu(i-ii, j-jj, k-kk) +  mu(i, j, k))
          mut_f    = 0.5*(mu_t(i-ii, j-jj, k-kk) + mu_t(i, j, k))

          total_mu = mu_f + mut_f
          rhoface  = 0.5 * (qp(i-ii, j-jj, k-kk, 1) + qp(i, j, k, 1))
          tkface   = 0.5 * (qp(i-ii, j-jj, k-kk, 6) + qp(i, j, k, 6))
          ! k in reynolds stress
          Tau_xx = -2.0*rhoface*tkface/3.0
          Tau_yy = Tau_xx
          Tau_zz = Tau_xx

          ! calling some element from memory and keep them handy for calculation
          nx    = faces(i,j,k)%nx
          ny    = faces(i,j,k)%ny
          nz    = faces(i,j,k)%nz
          area  = faces(i,j,k)%A

          ! adding viscous fluxes to stored convective flux
          F(i, j, k, 2) = F(i, j, k, 2) - (Tau_xx * nx * area)
          F(i, j, k, 3) = F(i, j, k, 3) - (Tau_yy * ny * area)
          F(i, j, k, 4) = F(i, j, k, 4) - (Tau_zz * nz * area)
          F(i, j, k, 5) = F(i, j, k, 5) - (area*((mu_f + sigma_k*mut_f)*(dtkdx*nx + dtkdy*ny + dtkdz*nz)))
          F(i, j, k, 6) = F(i, j, k, 6) - (area*((mu_f + sigma_k*mut_f)*(dtkdx*nx + dtkdy*ny + dtkdz*nz)))
          F(i, j, k, 7) = F(i, j, k, 7) - (area*((mu_f + sigma_phi*mut_f)*(dtkldx*nx + dtkldy*ny + dtkldz*nz)))
         
        end do
       end do
      end do

    end subroutine compute_viscous_fluxes_kkl