pressure_based_switching Subroutine

        subroutine pressure_based_switching(qp, f_qp_left, f_qp_right, pdif,  flags,  flow, dims)
          !< Pressure based switching. 
          !< User x,y, or z for I,J,or K face respectively
          !----------------------------------------------

            implicit none
            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):: qp
            !< Store primitive variable at cell center
            type(flowtype), intent(in) :: flow
            !< Information about fluid flow: freestream-speed, ref-viscosity,etc.
            integer, dimension(3), intent(in) :: flags
            !< flags for direction switch
            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
            real(wp), dimension(0:dims%imx,0:dims%jmx,0:dims%kmx), intent(inout) :: pdif
            !< pressure difference 
            integer :: i, j, k, i_end, j_end, k_end
            integer :: i_f, j_f, k_f  ! Flags to determine face direction
            real(wp) :: pd2

            DebugCall('pressure_based_switching')


            i_f = flags(1)
            j_f = flags(2)
            k_f = flags(3)
            i_end = dims%imx - 1 +i_f
            j_end = dims%jmx - 1 +j_f
            k_end = dims%kmx - 1 +k_f

            ! i_end and j_end denote number of faces
            ! Total number of cells including ghost_cells is
            ! (i_end+1) * j_end for xi faces and i_end*(j_end+1) for
            ! eta faces. 

            ! Loop over cells (physical)
            do k = 1, dims%kmx - 1
             do j = 1, dims%jmx - 1
              do i = 1, dims%imx - 1
                pd2 = abs(qp(i + i_f*1, j + j_f*1, k + k_f*1, 5) - &  !pressure
                          qp(i - i_f*1, j - j_f*1, k - k_f*1, 5))
                pdif(i, j, k) = 1 - (pd2/(pd2 + flow%pressure_inf))
              end do
             end do
            end do

            ! Update at ghost cells
            pdif((1-i_f):(1-i_f)*(dims%imx-1), (1-j_f):(1-j_f)*(dims%jmx-1), &
                 (1-k_f):(1-k_f)*(dims%kmx-1)) = &
                pdif(1:dims%imx-1 - i_f*(dims%imx-2), 1:dims%jmx-1 - j_f*(dims%jmx-2), &
                     1:dims%kmx-1 - k_f*(dims%kmx-2))
                    
            pdif((1-i_f*(-dims%imx+1)):(i_f)+(dims%imx-1), (1-j_f*(-dims%jmx+1)):(j_f)+(dims%jmx-1), &
                 (1-k_f*(-dims%kmx+1)):(k_f)+(dims%kmx-1)) = &
                pdif(1+i_f*(dims%imx-2):dims%imx-1 , 1+j_f*(dims%jmx-2):dims%jmx-1, &
                     1+k_f*(dims%kmx-2):dims%kmx-1)
             

            ! Loop over faces
            do k = 1, dims%kmx - (1 - k_f)            
             do j = 1, dims%jmx - (1 - j_f)
              do i = 1, dims%imx - (1 - i_f)
                f_qp_left(i, j, k, :) = qp(i - i_f*1, j - j_f*1, k - k_f*1, :) + (&
                    pdif(i - i_f*1, j - j_f*1, k - k_f*1) * ( &
                    f_qp_left(i, j, k, :) - qp(i - i_f*1, j - j_f*1, k - k_f*1, :)))

                f_qp_right(i, j, k, :) = qp(i, j, k, :) - (&
                    pdif(i, j, k) * ( &
                    qp(i, j, k, :) - f_qp_right(i, j, k, :)))
              end do
             end do
            end do

        end subroutine pressure_based_switching