diff --git a/src/fitpack_core.f90 b/src/fitpack_core.f90
index 6f3c094..791df5c 100644
--- a/src/fitpack_core.f90
+++ b/src/fitpack_core.f90
@@ -12232,7 +12232,7 @@ pure subroutine fpspgr(iopt,ider,u,mu,v,mv,r,mr,r0,r1,s, &
       end subroutine fpspgr
 
 
-      subroutine fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eta,tol,maxit, &
+      pure subroutine fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eta,tol,maxit, &
                              ib1,ib3,nc,ncc,intest,nrest,nt,tt,np,tp,c,fp,sup,fpint,coord,f, &
                              ff,row,coco,cosi,a,q,bt,bp,spt,spp,h,index,nummer,wrk,lwrk,ier)
 
@@ -12271,8 +12271,6 @@ subroutine fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eta,tol,maxit, &
       ntt    = 0
       iband1 = 0
       
-      print *, 'ib3=',ib3,' size(h)=',size(h)
-
       bootstrap: if (iopt>=0) then
 
           !  if iopt=0 we begin by computing the weighted least-squares polynomial of the form
@@ -17238,7 +17236,7 @@ pure subroutine spgrid(iopt,ider,mu,u,mv,v,r,r0,r1,s, &
       end subroutine spgrid
 
 
-      subroutine sphere(iopt,m,teta,phi,r,w,s,ntest,npest, &
+      pure subroutine sphere(iopt,m,teta,phi,r,w,s,ntest,npest, &
                              eps,nt,tt,np,tp,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
 
       !  subroutine sphere determines a smooth bicubic spherical spline
diff --git a/src/fitpack_core_c.f90 b/src/fitpack_core_c.f90
index 2b7950c..c9c7cf2 100644
--- a/src/fitpack_core_c.f90
+++ b/src/fitpack_core_c.f90
@@ -301,7 +301,7 @@ subroutine pogrid_c(iopt,ider,mu,u,mv,v,z,z0,r,s,nuest,nvest,nu,tu,nv,tv, &
           
       end subroutine pogrid_c
 
-      subroutine sphere_c(iopt,m,teta,phi,r,w,s,ntest,npest,eps,nt,tt,np,tp,c, &
+      pure subroutine sphere_c(iopt,m,teta,phi,r,w,s,ntest,npest,eps,nt,tt,np,tp,c, &
                                fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier) bind(C,name='sphere_c')
           integer(FP_SIZE), intent(in), value :: iopt,m,ntest,npest,lwrk1,lwrk2,kwrk
           real(FP_REAL),    intent(in), value :: s,eps