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
            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 state at faces
            real(wp), dimension(0:dims%imx,0:dims%jmx,0:dims%kmx), intent(inout) :: pdif
            !< pressure difference
            type(flowtype), intent(in) :: flow
            ! Character can be x or y or z
            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) - &
                          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(((dims%imx-1)*i_f)+1:dims%imx-1+i_f, &
                 ((dims%jmx-1)*j_f)+1:dims%jmx-1+j_f, &
                 ((dims%kmx-1)*k_f)+1:dims%kmx-1+k_f) &
                                        =   &
                pdif(i_f*(dims%imx-2)+1:dims%imx-1, &
                     j_f*(dims%jmx-2)+1:dims%jmx-1, &
                     k_f*(dims%kmx-2)+1: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