This paper deals with numerical solution of singularly perturbed parabolic partial differential equations with large negative shift on the spatial variable and integral boundary condition on the right side of the domain. For small values of perturbation parameter ɛ, solution of the problem exhibits an interior layer besides a parabolic boundary layers. The numerical methods combine the Crank–Nicolson difference scheme for the temporal direction on the uniform mesh and exponentially fitted operator finite difference method in spatial discretization. The integral boundary condition is treated using numerical integration techniques. The rate of convergence of the scheme is optimized by applying Richardson extrapolation technique. The stability and uniform convergence analysis has been carried out. To validate the applicability of the scheme, numerical examples are presented and the method is more accurate and shown to be second order uniformly convergent in both direction, and it also improves the results of the methods (Elango et al., 2021; Gobena and Duressa, 2021) existing in the literature.