We present a python implementation of the implicit level set algorithm to model propagation of planar hydraulic fractures. This algorithm was first introduced by Peirce & Detournay in (2008). The growth of a fluid driven planar fracture is modelled by coupling elasticity, described by a hypersingular integral equation (assuming a homogeneous medium) with fluid flow inside the fracture, described by the lubrication equation together with the mode I fracture propagation condition. The solution of the coupled system is numerically challenging due to a number of reasons. First, the lubrication equation describing the viscous fluid flow degenerate at the fracture tip (see e.g. Detournay and Peirce, 2014). Secondly, the solution of the elasticity equation is not straight forward as it involves a boundary integral equation with a hypersingular kernel. Finally, the footprint of the fracture is not known as the resolution of the fracture front requires the velocity field, which is a-priori unknown. The Implicit Level Set Algorithm (ILSA) tracks the free boundary by enforcing implicitely the hydraulic fracture opening tip assymptote at the computational domain cells in the vicinity off the tip. The regular rectangular computational grid is divided into three types of cells. The tip cells, i.e. the cells containing the fracture tip, the ribbon cells comprising of a narrow band of cells adjacent to the tip cells inside the fracture and the channel cells containing the remaining cells inside the fracture. The boundary conditions for the Eikonal equations are prescribed by inverting the tip assymptote in the ribbon cells from the knowledge of the trial frature opening in those cells. The fracture front is then located by solving the Eikonal equation using the Fast Marching Method (Sethian, 1996). Since the method is implicit, the Eikonal equation is solved iteratively in a time step until the correct position of the fracture front satisfying both the mass and momentum conservation of the fluid and the elasticity equation is located. We have verified the algorithm by testing the propagation of a radial fracture in the toughness dominated regime for which an analytical solution exist (Abé et al., 1976). In this propagation regime, the viscous loses are negligible and the pressure is uniform throughout the fracture. Our results show that the numerical solution matches with the analytical solution. We also present a number of others simulations for the case where the confining stresses is spatially varying, resulting in deviation of the fracture footprint from the initially radial geometry.