25 subroutine r8tcspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
26 > ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
27 > ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th,
28 > ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph,
29 > wk,nwk,ilinx,ilinth,ilinph,ier)
34 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
35 INTEGER inth,inph,inf4,inf5,ibcxmin,ibcxmax,inb1x,ibcthmin
36 INTEGER ibcthmax,inb1th,ibcphmin,ibcphmax,inb1ph,nwk,ilinx
37 INTEGER ilinth,ilinph,ier,inx,iflg,ith,ix,itest,ierx,ierth
38 INTEGER ierph,iaspl2,iabcx1,iabcx2,iabcth1,iabcth2,iawk,inwk
39 INTEGER iph,inpho,ic1,ic2,ibcphmina,ibcphmaxa,iabcph1,iabcph2
40 INTEGER iaccoef,i,iskip1,iskip2,ia1,ia2
45 real*8 x(inx),th(inth),ph(inph)
46 real*8 fspl(4,4,4,inf4,inf5,inph),wk(nwk)
47 real*8 bcxmin(inb1x,*),bcxmax(inb1x,*)
48 real*8 bcthmin(inb1th,*),bcthmax(inb1th,*)
49 real*8 bcphmin(inb1ph,*),bcphmax(inb1ph,*)
167 if(ibcphmin.ne.-1)
then
168 if((ibcphmin.eq.1).or.(ibcphmin.eq.2))
then
171 if(bcphmin(ix,ith).ne.z0) iflg=1
175 if((ibcphmax.eq.1).or.(ibcphmax.eq.2))
then
178 if(bcphmax(ix,ith).ne.z0) iflg=1
184 itest=10*max(inx,inth,inph) + 20*inx*inth
186 itest=16*inx*inth*inph
189 if(nwk.lt.itest)
then
190 write(6,9901) nwk,itest
192 9901
format(
' ?tcspline: workspace too small.'/
193 >
' user supplied nwk=',i7,
'; need at least: ',i7/
194 >
' If no explicit df/dph boundary condition is set,'/
195 >
' nwk = 20*inx*inth + 10*max(inx,inth,inph) can be used.'/
196 >
' If an explicit df/dph or d2f/dph2 boundary condition'/
197 >
' is set, nwk=16*inx*inth*inph is required.')
200 write(6,
'('' ?tcspline: at least 2 x points required.'')')
204 write(6,
'('' ?tcspline: need at least 2 theta points.'')')
208 write(6,
'('' ?tcspline: need at least 2 phi points.'')')
212 if((ibcxmin.eq.1).or.(ibcxmax.eq.1).or.(ibcxmin.eq.2).or.
213 > (ibcxmax.eq.2))
then
214 if(inb1x.lt.inth)
then
217 >
'('.lt.
' ?tcspline: 1st dim of bcxmin/max arrays inth'')')
221 if((ibcthmin.eq.1).or.(ibcthmax.eq.1).or.(ibcthmin.eq.2).or.
222 > (ibcthmax.eq.2))
then
223 if(inb1th.lt.inx)
then
226 >
'('.lt.
' ?tcspline: 1st dim of bcthmin/max arrays inx'')')
230 if((ibcphmin.eq.1).or.(ibcphmax.eq.1).or.(ibcphmin.eq.2).or.
231 > (ibcphmax.eq.2))
then
232 if(inb1ph.lt.inx)
then
235 >
'('.lt.
' ?tcspline: 1st dim of bphmin/max arrays inx'')')
239 call ibc_ck(ibcxmin,
'tcspline',
'xmin',-1,7,ier)
240 if(ibcxmin.ge.0)
call ibc_ck(ibcxmax,
'tcspline',
'xmax',0,7,ier)
242 call ibc_ck(ibcthmin,
'tcspline',
'thmin',-1,7,ier)
243 if(ibcthmin.ge.0)
call ibc_ck(ibcthmax,
'tcspline',
'thmax',0,7,ier)
245 call ibc_ck(ibcphmin,
'tcspline',
'phmin',-1,7,ier)
246 if(ibcphmax.ge.0)
call ibc_ck(ibcphmax,
'tcspline',
'phmax',0,7,ier)
250 call r8splinck(x,inx,ilinx,ztol,ierx)
254 write(6,
'('' ?tcspline: x axis not strict ascending'')')
259 call r8splinck(th,inth,ilinth,ztol,ierth)
263 write(6,
'('' ?tcspline: theta axis not strict ascending'')')
268 call r8splinck(ph,inph,ilinph,ztol,ierph)
272 write(6,
'('' ?tcspline: phi axis not strict ascending'')')
285 iabcx1=iaspl2+16*inx*inth
297 wk(iabcth1+ix-1)=0.0_r8
298 wk(iabcth2+ix-1)=0.0_r8
299 if((ibcthmin.eq.1).or.(ibcthmin.eq.2))
then
300 wk(iabcth1+ix-1)=bcthmin(ix,iph)
302 if((ibcthmin.ne.-1).and.
303 > ((ibcthmax.eq.1).or.(ibcthmax.eq.2)))
then
304 wk(iabcth2+ix-1)=bcthmax(ix,iph)
308 wk(iabcx1+ith-1)=0.0_r8
309 wk(iabcx2+ith-1)=0.0_r8
310 if((ibcxmin.eq.1).or.(ibcxmin.eq.2))
then
311 wk(iabcx1+ith-1)=bcxmin(ith,iph)
313 if((ibcxmin.ne.-1).and.
314 > ((ibcxmax.eq.1).or.(ibcxmax.eq.2)))
then
315 wk(iabcx2+ith-1)=bcxmax(ith,iph)
321 call r8tcsp23(x,inx,th,inth,fspl(1,1,1,1,1,iph),inf4,
322 > ibcxmin,wk(iabcx1),ibcxmax,wk(iabcx2),
323 > ibcthmin,wk(iabcth1),ibcthmax,wk(iabcth2),
324 > wk(iaspl2),wk(iawk),inwk,ilinx,ilinth,ier)
327 write(6,*)
' ?tcspline: error in 2d spline, exiting.'
351 wk(4*(iph-1)+1)=fspl(ic1,ic2,1,ix,ith,iph)
365 if((ibcphmin.eq.1).or.(ibcphmin.eq.2)) ibcphmina=0
366 if((ibcphmax.eq.1).or.(ibcphmax.eq.2)) ibcphmaxa=0
369 call r8v_spline(ibcphmina,ibcphmaxa,inph,ph,wk,
375 fspl(ic1,ic2,2,ix,ith,iph)=wk(4*(iph-1)+2)
376 fspl(ic1,ic2,3,ix,ith,iph)=wk(4*(iph-1)+3)*xo2
377 fspl(ic1,ic2,4,ix,ith,iph)=wk(4*(iph-1)+4)*xo6
396 iabcph2=iabcph1+inx*inth
397 iaccoef=iabcph2+inx*inth
398 iawk=iaccoef+12*inx*inth*inph
409 if(ibcphmin.eq.1)
then
411 else if(ibcphmin.eq.2)
then
418 if(ibcphmax.eq.1)
then
420 else if(ibcphmax.eq.2)
then
434 call r8tcspeval(x(ix),th(ith),ph(1),iselect1, zcur,
435 > x,inx,th,inth,ph,inph,ilinx,ilinth,ilinph,
436 > fspl,inf4,inf5,ier)
438 write(6,*)
' ?? tcspline: error in tcspeval call'
441 wk(ia1)=bcphmin(ix,ith)-zcur(1)
447 call r8tcspeval(x(ix),th(ith),ph(inph),iselect2, zcur,
448 > x,inx,th,inth,ph,inph,ilinx,ilinth,ilinph,
449 > fspl,inf4,inf5,ier)
451 write(6,*)
' ?? tcspline: error in tcspeval call'
454 wk(ia2)=bcphmax(ix,ith)-zcur(1)
461 call r8tcspcorr(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
462 > ibcxmin,ibcxmax,ibcthmin,ibcthmax,
463 > ibcphmin,wk(iabcph1),ibcphmax,wk(iabcph2),
464 > wk(iaccoef),wk(iawk),inwk)
471 subroutine r8tcspcorr(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
472 > ibcxmin,ibcxmax,ibcthmin,ibcthmax,
473 > ibcphmin,bcph1,ibcphmax,bcph2,ccorr,wk,nwk)
484 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
485 INTEGER inth,inph,inf4,inf5,nwk,inx,iawk2,inpho,ith,ix,iph
486 INTEGER inxo,icph,intho,icx,icth
489 real*8 xo2,xo6,z0,zfac
495 real*8 fspl(4,4,4,inf4,inf5,inph)
497 integer ibcxmin,ibcxmax
498 integer ibcthmin,ibcthmax
499 integer ibcphmin,ibcphmax
501 real*8 bcph1(inx,inth)
502 real*8 bcph2(inx,inth)
506 real*8 ccorr(3,4,inx,inth,inph)
517 if(nwk.lt.10*max(inx,inth,inph))
then
518 write(6,*)
' ?? programming error in tcspcorr (tcspline)'
535 if(ibcphmin.eq.1)
then
537 else if(ibcphmin.eq.2)
then
541 if(ibcphmax.eq.1)
then
542 wk(inpho+2)=bcph2(ix,ith)
543 else if(ibcphmax.eq.2)
then
544 wk(inpho+3)=bcph2(ix,ith)
547 call r8v_spline(ibcphmin,ibcphmax,inph,ph,wk,wk(iawk2))
552 ccorr(1,1,ix,ith,iph)=wk(4*(iph-1)+2)
553 ccorr(2,1,ix,ith,iph)=wk(4*(iph-1)+3)*xo2
554 ccorr(3,1,ix,ith,iph)=wk(4*(iph-1)+4)*xo6
572 wk(4*(ix-1)+1)=ccorr(icph,1,ix,ith,iph)
582 call r8v_spline(ibcxmin,ibcxmax,inx,x,wk,wk(iawk2))
585 ccorr(icph,2,ix,ith,iph)=wk(4*(ix-1)+2)
586 ccorr(icph,3,ix,ith,iph)=wk(4*(ix-1)+3)*xo2
587 ccorr(icph,4,ix,ith,iph)=wk(4*(ix-1)+4)*xo6
610 wk(4*(ith-1)+1)=ccorr(icph,icx,ix,ith,iph)
620 call r8v_spline(ibcthmin,ibcthmax,inth,th,wk,
626 if(icth.eq.3) zfac=xo2
627 if(icth.eq.4) zfac=xo6
628 fspl(icx,icth,icph+1,ix,ith,iph)=
629 > fspl(icx,icth,icph+1,ix,ith,iph)+
630 > wk(4*(ith-1)+icth)*zfac
643 subroutine r8tcsp23(x,inx,th,inth,fspl,inf4,
644 > ibcxmin,bcxmin,ibcxmax,bcxmax,
645 > ibcthmin,bcthmin,ibcthmax,bcthmax,
646 > fspl2,wk,nwk,ilinx,ilinth,ier)
655 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
656 INTEGER inth,inf4,ibcxmin,ibcxmax,ibcthmin,ibcthmax,nwk,ilinx
657 INTEGER ilinth,ier,inx,ith,ix,j,i
662 REAL*8 fspl(4,4,4,inf4,inth)
663 real*8 fspl2(4,4,inx,inth)
665 real*8 bcxmin(inth),bcxmax(inth)
666 real*8 bcthmin(inx),bcthmax(inx)
676 fspl2(1,1,ix,ith)=fspl(1,1,1,ix,ith)
682 call r8bcspline(x,inx,th,inth,fspl2,inx,
683 > ibcxmin,bcxmin,ibcxmax,bcxmax,
684 > ibcthmin,bcthmin,ibcthmax,bcthmax,
685 > wk,nwk,ilinx,ilinth,ier)
694 fspl(i,j,1,ix,ith)=fspl2(i,j,ix,ith)