A Fortran implementation of Isogeometric Analysis (IGA) for thin shell problems of various geometries was presented. Non-Uniform Rational B-Splines (NURBS) basis functions were used for both the approximations of the solution and description of the geometry. The necessary continuities can be easily met with NURBS basis functions in the Kirchhoff-Love shell theory. For thin shell problems, low basis functions may produce inaccurate results and it can be improved by h-refinement and k-refinement. However, the computational cost increases. The high performance of the Fortran implementation can handle the issue. Due to the non-interpolated property of NURBS basis functions, the penalty method was applied for imposing boundary conditions. Several typical examples showed the efficiency and accuracy of the Fortran implementation.