Microbubble enhanced HIFU is of great interest to tissue ablation for solid tumor treatments, where contrast agents/microbubbles are injected into the targeted region to promote heating and reduce pre-focal tissue damage. A compressible Euler-Lagrange model has been developed to accurately characterize the acoustic and thermal fields during this process. This employs a compressible Navier-Stokes solver for the ultrasound acoustic field and a discrete singularities model for bubble dynamics. To address the demanding computational cost in practical medical applications, a hybrid MPI-OpenMP parallelization scheme is developed to take advantage of both scalability of MPI and load balancing of OpenMP. Firstly the Eulerian computational domain is divided into multiple subdomains and the bubbles are subdivided in groups based on which subdomain they fall into. Then, in each subdomain containing bubbles, multiple OpenMP threads are activated to speed up the computations of the bubble dynamics. OpenMP threads are more heavily distributed to subdomains where the bubbles are clustered. Thus MPI load imbalance due to uneven bubble distribution is mitigated by OpenMP speedup locally for those subdomains hosting more bubbles. The is used to conduct simulations and physical studies of bubble-enhanced HIFU containing a large number of microbubbles. The phenomenon of acoustic shadowing caused by the bubble cloud is then analyzed and discussed. Efficiency tests on two different machines with 48 processors are conducted and indicate 2 to 3 times of speedup with the same hardware by introducing an OpenMP parallelization in combination to the MPI parallelization.