@@ -3615,12 +3615,12 @@ pure subroutine fpclos(iopt,idim,m,u,mx,x,w,k,s,nest,tol, &
36153615 end subroutine fpclos
36163616
36173617 pure subroutine fpclos_reset_interp(idim,k,m,mx,n,nc,nest,kk,kk1,u,x,t,c,fp,per,fp0,s,fpint,nrdata,done)
3618- integer(FP_SIZE), intent(in) :: idim,k,m,mx,n,nc,nest
3618+ integer(FP_SIZE), intent(in) :: idim,k,m,mx,n,nc,nest
36193619 integer(FP_SIZE), intent(inout) :: kk,kk1
3620- real(FP_REAL), intent(in) :: u(m),x(mx),per,fp0,s
3620+ real(FP_REAL), intent(in) :: u(m),x(mx),per,fp0,s
36213621 real(FP_REAL), intent(inout) :: t(nest),c(nc),fp,fpint(nest)
36223622 integer(FP_SIZE), intent(inout) :: nrdata(nest)
3623- logical(FP_BOOL), intent(out) :: done
3623+ logical(FP_BOOL), intent(out) :: done
36243624
36253625 integer(FP_SIZE) :: i,j,j1,jj,m1
36263626
@@ -4765,11 +4765,11 @@ pure subroutine fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol, &
47654765
47664766 ! ..
47674767 ! ..scalar arguments..
4768- real(FP_REAL), intent(in) :: xb,xe,s,tol
4769- real(FP_REAL), intent(out) :: fp
4770- integer(FP_SIZE), intent(in) :: iopt,m,k,nest,maxit,k1,k2
4771- integer(FP_SIZE), intent(out ) :: n
4772- integer(FP_FLAG), intent(out) :: ier
4768+ real(FP_REAL), intent(in) :: xb,xe,s,tol
4769+ real(FP_REAL), intent(out) :: fp
4770+ integer(FP_SIZE), intent(in) :: iopt,m,k,nest,maxit,k1,k2
4771+ integer(FP_SIZE), intent(inout ) :: n
4772+ integer(FP_FLAG), intent(out) :: ier
47734773
47744774 ! ..array arguments..
47754775 real(FP_REAL), intent(in) :: x(m),y(m),w(m)
@@ -4893,7 +4893,7 @@ pure subroutine fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol, &
48934893 fp = zero
48944894
48954895 ! initialize the observation matrix a.
4896- z(1:nk1) = zero
4896+ z(1:nk1) = zero
48974897 a(1:nk1,1:k1) = zero
48984898 l = k1
48994899 coefs: do it=1,m
@@ -5198,10 +5198,10 @@ pure subroutine fpcuro(a,b,c,d,x,n)
51985198
51995199 else
52005200
5201- u = sign(sqrt(abs(q)),r)
5201+ u = sign(sqrt(abs(q)),r)
52025202 p3 = atan2(sqrt(-disc),abs(r))*e3
52035203 u2 = u+u
5204- n = 3
5204+ n = 3
52055205 x(1) = -u2*cos(p3)-b1
52065206 x(2) = u2*cos(pi3-p3)-b1
52075207 x(3) = u2*cos(pi3+p3)-b1
@@ -12232,7 +12232,7 @@ pure subroutine fpspgr(iopt,ider,u,mu,v,mv,r,mr,r0,r1,s, &
1223212232 end subroutine fpspgr
1223312233
1223412234
12235- pure subroutine fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eta,tol,maxit, &
12235+ subroutine fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eta,tol,maxit, &
1223612236 ib1,ib3,nc,ncc,intest,nrest,nt,tt,np,tp,c,fp,sup,fpint,coord,f, &
1223712237 ff,row,coco,cosi,a,q,bt,bp,spt,spp,h,index,nummer,wrk,lwrk,ier)
1223812238
@@ -12270,6 +12270,8 @@ pure subroutine fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eta,tol,maxit, &
1227012270 lwest = 0
1227112271 ntt = 0
1227212272 iband1 = 0
12273+
12274+ print *, 'ib3=',ib3,' size(h)=',size(h)
1227312275
1227412276 bootstrap: if (iopt>=0) then
1227512277
@@ -12515,7 +12517,9 @@ pure subroutine fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eta,tol,maxit, &
1251512517 jlt = j+lt
1251612518 htj = ht(j)
1251712519 if (jlt==3) then
12518- h(1:3) = [h(1)+htj,facc*htj,facs*htj]
12520+ h(1) = h(1)+htj
12521+ h(2) = facc*htj
12522+ h(3) = facs*htj
1251912523 j1 = 3
1252012524 elseif (jlt==nt4) then
1252112525 h(j1+1:j1+3) = htj*[facc,facs,one]
@@ -12811,7 +12815,8 @@ pure subroutine fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eta,tol,maxit, &
1281112815 if (j>1 .and. j<nt6) then
1281212816 h(1:npp) = row(1:npp)
1281312817 else
12814- h(1:2) = [facc,facs]
12818+ h(1) = facc
12819+ h(2) = facs
1281512820 if(j==1) jrot = 2
1281612821 endif
1281712822
@@ -14201,12 +14206,22 @@ pure subroutine insert(iopt,t,n,c,k,x,tt,nn,cc,nest,ier)
1420114206 integer(FP_SIZE) :: kk,k1,l,nk,nk1
1420214207 ! ..
1420314208 ! before starting computations a data check is made. if the input data
14204- ! are invalid control is immediately repassed to the calling program.
14209+ ! are invalid control is immediately repassed to the calling program.
1420514210 ier = FITPACK_INPUT_ERROR
14206- if (nest<=n) return
14211+ if (nest<=n) then
14212+ tt = t
14213+ cc = c
14214+ nn = n
14215+ return
14216+ endif
1420714217 k1 = k+1
1420814218 nk = n-k
14209- if (x<t(k1) .or. x>t(nk)) return
14219+ if (x<t(k1) .or. x>t(nk)) then
14220+ tt = t
14221+ cc = c
14222+ nn = n
14223+ return
14224+ endif
1421014225 ! search for knot interval t(l) <= x < t(l+1).
1421114226 nk1 = nk-1
1421214227 l = k1
@@ -14215,11 +14230,21 @@ pure subroutine insert(iopt,t,n,c,k,x,tt,nn,cc,nest,ier)
1421514230 end do
1421614231
1421714232 ! no interval found in whole range
14218- if(t(l)>=t(l+1)) return
14233+ if (t(l)>=t(l+1)) then
14234+ tt = t
14235+ cc = c
14236+ nn = n
14237+ return
14238+ endif
1421914239
14220- if(iopt/=0) then
14240+ if (iopt/=0) then
1422114241 kk = 2*k
14222- if (l<=kk .and. l>=(n-kk)) return
14242+ if (l<=kk .and. l>=(n-kk)) then
14243+ tt = t
14244+ cc = c
14245+ nn = n
14246+ return
14247+ end if
1422314248 endif
1422414249
1422514250 ier = FITPACK_OK
@@ -14232,10 +14257,10 @@ end subroutine insert
1423214257 ! that called subroutine insert with the same variables as input and output arguments
1423314258 pure subroutine insert_inplace(iopt,t,n,c,k,x,nest,ier)
1423414259
14235- integer(FP_SIZE), intent(in) :: iopt,k,nest
14236- integer(FP_FLAG), intent(out) :: ier
14237- real(FP_REAL), intent(in) :: x
14238- integer(FP_SIZE), intent(inout) :: n
14260+ integer(FP_SIZE), intent(in) :: iopt,k,nest
14261+ integer(FP_FLAG), intent(out) :: ier
14262+ real(FP_REAL), intent(in) :: x
14263+ integer(FP_SIZE), intent(inout) :: n
1423914264 ! ..array arguments..
1424014265 real(FP_REAL), intent(inout) :: t(nest),c(nest)
1424114266
@@ -17213,7 +17238,7 @@ pure subroutine spgrid(iopt,ider,mu,u,mv,v,r,r0,r1,s, &
1721317238 end subroutine spgrid
1721417239
1721517240
17216- pure subroutine sphere(iopt,m,teta,phi,r,w,s,ntest,npest, &
17241+ subroutine sphere(iopt,m,teta,phi,r,w,s,ntest,npest, &
1721717242 eps,nt,tt,np,tp,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
1721817243
1721917244 ! subroutine sphere determines a smooth bicubic spherical spline
0 commit comments