reconstruct_kmax Subroutine

    subroutine reconstruct_kmax(qp, scheme, bc)
      !< Reconstruct state at the KMAX boundary face with MUSCL scheme

      implicit none
      real(wp), dimension(-2:imx+2, -2:jmx+2, -2:kmx+2, 1:n_var), intent(in) :: qp
      type(schemetype), intent(in) :: scheme
      type(boundarytype), intent(in) :: bc
      real(wp) :: psi1, psi2, fd, bd, r
      integer :: i, j, k, l
      real(wp) :: kappa, phi
    
      phi = 1.0
      kappa = 1./3.
      switch_L=scheme%klimiter_switch

      do k = kmx-1, kmx-1
       do l = 1, n_var
         if(l>=6) switch_L=scheme%ktlimiter_switch
         if(l<6) switch_L=scheme%klimiter_switch
        do j = 1, jmx - 1
         do i = 1, imx - 1
          ! left face of kmx ghost cell
          z_qp_right(i, j, k+1, l) = 0.5 * (qp(i, j, k, l) + qp(i, j, k+1, l))

          if (ppm_flag==1) then

            fd = qp(i, j, k+1, l) - qp(i, j, k, l)
            bd = qp(i, j, k, l) - qp(i, j, k-1, l)

            r = fd / bd
            psi1 = max(0., min(2*r, (2 + r)/3., 2.))
            psi1 = (1 - (1 - psi1)*switch_L )
            r = bd / fd
            psi2 = max(0., min(2*r, (2 + r)/3., 2.))
            psi2 = (1 - (1 - psi2)*switch_L )

            ! right face of last k interior cell
            z_qp_left(i, j, k+1, l) = qp(i, j, k, l) + 0.25*phi* &
                (((1-kappa) * psi1 * bd) + ((1+kappa) * psi2 * fd))

            ! left face of last k cell
            z_qp_right(i, j, k, l) = qp(i, j, k, l) - 0.25*phi* &
                (((1+kappa) * psi1 * bd) + ((1-kappa) * psi2 * fd))
         end if
         end do
        end do
       end do
      end do
      if(bc%kmax_id==-8 .or. bc%kmax_id==-9)then
         z_qp_left(1:imx-1,1:jmx-1,kmx,1:n_var) = qp(1:imx-1,1:jmx-1,kmx,1:n_var) 
        z_qp_right(1:imx-1,1:jmx-1,kmx,1:n_var) = qp(1:imx-1,1:jmx-1,kmx,1:n_var) 
      else
        z_qp_right(1:imx-1,1:jmx-1,kmx,1:n_var) = 0.5*(qp(1:imx-1,1:jmx-1,kmx-1,1:n_var)&
                                                      +qp(1:imx-1,1:jmx-1,kmx  ,1:n_var))
      end if

    end subroutine reconstruct_kmax