compute_face_state Subroutine

        subroutine compute_face_state(qp, f_qp_left, f_qp_right, flags, lam_switch, turb_switch, dims)
          !< Subroutine to calculate state at the face, generalized for
          !< all direction : I,J, and K.
            implicit none
            type(extent), intent(in) :: dims
            !< Extent of the domain:imx,jmx,kmx
            integer, dimension(3), intent(in) :: flags
            !< Flags for direction 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
            real(wp), dimension(1-flags(1):dims%imx-1+2*flags(1), 1-flags(2):dims%jmx-1+2*flags(2), 1-flags(3):dims%kmx-1+2*flags(3), 1:dims%n_var), intent(inout) :: f_qp_left, f_qp_right
            !< primitive variable at cell faces
            integer, intent(in) :: lam_switch
            !< Limiter switch for laminar variables
            integer, intent(in) :: turb_switch
            !< Limiter switch for turbulent variables
            integer :: i, j, k, l
            !< integer used for DO loop
            integer :: ii, jj, kk
            !< Variable for ALFA family limiter
            real(wp) :: alpha
            !< Flags to determine face direction
            real(wp) :: psi1, psi2
            !< limiters
            real(wp) :: fd
            !< forward difference
            real(wp) :: bd
            !< backward difference
            real(wp) :: r
            !< ratio of differences

            DebugCall('compute_face_state')


            alpha = 2./3. !Koren limiter 
            phi = 1.0
            kappa = 1./3.
            switch_L=lam_switch

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

            do l = 1, dims%n_var
              if(l>=6)then
                switch_L=turb_switch
              end if
             do k = 1-kk, dims%kmx - 1 + kk
              do j = 1-jj, dims%jmx - 1 + jj
               do i = 1-ii, dims%imx - 1 + ii
                ! Cell based
                ! Koren limiter for now
                ! From paper: delta: forward difference 'fd'
                !             nabla: backward difference 'bd'
                fd = qp(i+ii, j+jj, k+kk, l) - qp(i, j, k, l)
                bd = qp(i, j, k, l) - qp(i-ii, j-jj, k-kk, l)

                r = fd / max(bd,1e-10)
                psi1 = max(0., min(2*r, alpha*(r-1.0) + 1.0, 2.))  !alpha limiter
!                psi1 = max(0., min(2*r,1.), min(r,2.))    ! superbee
!                psi1 = ((r*r) + r)/((r*r) + 1.0)          ! Van-Albda 
!                psi1 = (abs(r) + r)/(abs(r) + 1.0)          ! Van-Leer

                r = bd / max(fd, 1e-10)
                psi2 = max(0., min(2*r, alpha*(r-1.0) + 1.0, 2.))
!                psi2 = max(0., min(2*r,1.), min(r,2.))
!                psi2 = ((r*r) + r)/((r*r) + 1.0)          ! Van-Albda 
!                psi2 = (abs(r) + r)/(abs(r) + 1.0)          ! Van-Leer

                psi1 = (1 - (1 - psi1)*switch_L )
                psi2 = (1 - (1 - psi2)*switch_L )

                f_qp_left(i+ii, j+jj, k+kk, l) = qp(i, j, k, l) + 0.25*phi* &
                    (((1.-kappa) * psi1 * bd) + ((1.+kappa) * psi2 * fd))
                f_qp_right(i, j, k, l) = qp(i, j, k, l) - 0.25*phi* &
                    (((1.+kappa) * psi1 * bd) + ((1.-kappa) * psi2 * fd))
               end do
              end do
             end do
            end do


        end subroutine compute_face_state