add_turbulent_time Subroutine

      subroutine add_turbulent_time(qp,delta_t,cells,Ifaces,Jfaces,Kfaces,CFL,flow,dims)
        !< Addition to local time step due to turbulence 
        implicit none
        real(wp), intent(in) :: CFL
        !< CFL number
        type(flowtype), intent(in) :: flow
        !< Information about fluid flow: freestream-speed, ref-viscosity,etc.
        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(in), target :: qp
        !< Store primitive variable at cell center
        real(wp) , dimension(1:dims%imx-1, 1:dims%jmx-1, 1:dims%kmx-1), intent(inout) :: 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
        !< Store face quantites for I faces 
        type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+3,-2:dims%kmx+2), intent(in) :: Jfaces
        !< Store face quantites for J faces 
        type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+3), intent(in) :: Kfaces
        !< Store face quantites for K faces 
        real(wp) :: lmx1, lmx2, lmx3, lmx4, lmx5, lmx6, lmxsum
        integer :: i, j, k

        DebugCall('add_viscous_time_step')

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

           ! Faces with lower index

           
           ! For left face: i.e., lower index face along xi direction
           lmx1 = mu_t(i,j,k)/(qp(i,j,k,1)*abs( &
                ((cells(i-1,j,k)%centerx - cells(i,j,k)%centerx) * Ifaces(i, j, k)%nx) + &
                ((cells(i-1,j,k)%centery - cells(i,j,k)%centery) * Ifaces(i, j, k)%ny) + &
                ((cells(i-1,j,k)%centerz - cells(i,j,k)%centerz) * Ifaces(i, j, k)%nz)))
           ! For front face, i.e., lower index face along eta direction
           lmx2 = mu_t(i,j,k)/(qp(i,j,k,1)*abs( &
                ((cells(i,j-1,k)%centerx - cells(i,j,k)%centerx) * Jfaces(i, j, k)%nx) + &
                ((cells(i,j-1,k)%centery - cells(i,j,k)%centery) * Jfaces(i, j, k)%ny) + &
                ((cells(i,j-1,k)%centerz - cells(i,j,k)%centerz) * Jfaces(i, j, k)%nz)))
           ! For bottom face, i.e., lower index face along zeta direction
           lmx3 = mu_t(i,j,k)/(qp(i,j,k,1)*abs( &
                ((cells(i,j,k-1)%centerx - cells(i,j,k)%centerx) * Kfaces(i, j, k)%nx) + &
                ((cells(i,j,k-1)%centery - cells(i,j,k)%centery) * Kfaces(i, j, k)%ny) + &
                ((cells(i,j,k-1)%centerz - cells(i,j,k)%centerz) * Kfaces(i, j, k)%nz)))

           
           ! For right face, i.e., higher index face along xi direction
           lmx4 = mu_t(i+1,j,k)/(qp(i+1,j,k,1)*abs( &
                ((cells(i,j,k)%centerx - cells(i+1,j,k)%centerx) * Ifaces(i+1, j, k)%nx) + &
                ((cells(i,j,k)%centery - cells(i+1,j,k)%centery) * Ifaces(i+1, j, k)%ny) + &
                ((cells(i,j,k)%centerz - cells(i+1,j,k)%centerz) * Ifaces(i+1, j, k)%nz)))
           ! For back face, i.e., higher index face along eta direction
           lmx5 = mu_t(i,j+1,k)/(qp(i,j+1,k,1)*abs( &
                ((cells(i,j,k)%centerx - cells(i,j+1,k)%centerx) * Jfaces(i, j+1, k)%nx) + &
                ((cells(i,j,k)%centery - cells(i,j+1,k)%centery) * Jfaces(i, j+1, k)%ny) + &
                ((cells(i,j,k)%centerz - cells(i,j+1,k)%centerz) * Jfaces(i, j+1, k)%nz)))
           ! For top face, i.e., higher index face along zeta direction
           lmx6 = mu_t(i,j,k+1)/(qp(i,j,k+1,1)*abs( &
                ((cells(i,j,k)%centerx - cells(i,j,k+1)%centerx) * Kfaces(i, j, k+1)%nx) + &
                ((cells(i,j,k)%centery - cells(i,j,k+1)%centery) * Kfaces(i, j, k+1)%ny) + &
                ((cells(i,j,k)%centerz - cells(i,j,k+1)%centerz) * Kfaces(i, j, k+1)%nz)))

               lmxsum = (Ifaces(i, j, k)%A * lmx1) + &
                        (Jfaces(i, j, k)%A * lmx2) + &
                        (Kfaces(i, j, k)%A * lmx3) + &
                        (Ifaces(i+1, j, k)%A * lmx4) + &
                        (Jfaces(i, j+1, k)%A * lmx5) + &
                        (Kfaces(i, j, k+1)%A * lmx6)

           lmxsum = flow%gm*lmxsum/flow%tPr

           lmxsum = 2./(lmxsum + (2.*CFL*cells(i,j,k)%volume/delta_t(i,j,k)))
        
           delta_t(i, j, k) = CFL*( lmxsum * cells(i, j, k)%volume)
          end do
         end do
        end do
      end subroutine add_turbulent_time