A new algorithm for constructing orthogonal curvilinear grids on a sphere for a fairly general geometric shape of the modeling region is implemented as a “compile-once - use forever” software package. It is based on the numerical solution of the inverse problem by an iterative procedure -- finding such distribution of grid points along its perimeter, so that the conformal transformation of the perimeter into a rectangle turns this distribution into uniform one. The iterative procedure itself turns out to be multilevel - i.e. an iterative loop built around another, internal iterative procedure. Thereafter, knowing this distribution, the grid nodes inside the region are obtained solving an elliptic problem. It is shown that it was possible to obtain the exact orthogonality of the perimeter at the corners of the grid, to achieve very small, previously unattainable level of orthogonality errors, as well as make it isotropic -- local distances between grid nodes about both directions are equal to each other.