25 subroutine tcspline(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)
31 real x(inx),th(inth),ph(inph)
32 real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
33 real bcxmin(inb1x,*),bcxmax(inb1x,*)
34 real bcthmin(inb1th,*),bcthmax(inb1th,*)
35 real bcphmin(inb1ph,*),bcphmax(inb1ph,*)
153 if(ibcphmin.ne.-1)
then
154 if((ibcphmin.eq.1).or.(ibcphmin.eq.2))
then
157 if(bcphmin(ix,ith).ne.z0) iflg=1
161 if((ibcphmax.eq.1).or.(ibcphmax.eq.2))
then
164 if(bcphmax(ix,ith).ne.z0) iflg=1
170 itest=10*max(inx,inth,inph) + 20*inx*inth
172 itest=16*inx*inth*inph
175 if(nwk.lt.itest)
then
176 write(6,9901) nwk,itest
178 9901
format(
' ?tcspline: workspace too small.'/
179 >
' user supplied nwk=',i7,
'; need at least: ',i7/
180 >
' If no explicit df/dph boundary condition is set,'/
181 >
' nwk = 20*inx*inth + 10*max(inx,inth,inph) can be used.'/
182 >
' If an explicit df/dph or d2f/dph2 boundary condition'/
183 >
' is set, nwk=16*inx*inth*inph is required.')
186 write(6,
'('' ?tcspline: at least 2 x points required.'')')
190 write(6,
'('' ?tcspline: need at least 2 theta points.'')')
194 write(6,
'('' ?tcspline: need at least 2 phi points.'')')
198 if((ibcxmin.eq.1).or.(ibcxmax.eq.1).or.(ibcxmin.eq.2).or.
199 > (ibcxmax.eq.2))
then
200 if(inb1x.lt.inth)
then
203 >
'('.lt.
' ?tcspline: 1st dim of bcxmin/max arrays inth'')')
207 if((ibcthmin.eq.1).or.(ibcthmax.eq.1).or.(ibcthmin.eq.2).or.
208 > (ibcthmax.eq.2))
then
209 if(inb1th.lt.inx)
then
212 >
'('.lt.
' ?tcspline: 1st dim of bcthmin/max arrays inx'')')
216 if((ibcphmin.eq.1).or.(ibcphmax.eq.1).or.(ibcphmin.eq.2).or.
217 > (ibcphmax.eq.2))
then
218 if(inb1ph.lt.inx)
then
221 >
'('.lt.
' ?tcspline: 1st dim of bphmin/max arrays inx'')')
225 call ibc_ck(ibcxmin,
'tcspline',
'xmin',-1,7,ier)
226 if(ibcxmin.ge.0)
call ibc_ck(ibcxmax,
'tcspline',
'xmax',0,7,ier)
228 call ibc_ck(ibcthmin,
'tcspline',
'thmin',-1,7,ier)
229 if(ibcthmin.ge.0)
call ibc_ck(ibcthmax,
'tcspline',
'thmax',0,7,ier)
231 call ibc_ck(ibcphmin,
'tcspline',
'phmin',-1,7,ier)
232 if(ibcphmax.ge.0)
call ibc_ck(ibcphmax,
'tcspline',
'phmax',0,7,ier)
236 call splinck(x,inx,ilinx,ztol,ierx)
240 write(6,
'('' ?tcspline: x axis not strict ascending'')')
245 call splinck(th,inth,ilinth,ztol,ierth)
249 write(6,
'('' ?tcspline: theta axis not strict ascending'')')
254 call splinck(ph,inph,ilinph,ztol,ierph)
258 write(6,
'('' ?tcspline: phi axis not strict ascending'')')
271 iabcx1=iaspl2+16*inx*inth
285 if((ibcthmin.eq.1).or.(ibcthmin.eq.2))
then
286 wk(iabcth1+ix-1)=bcthmin(ix,iph)
288 if((ibcthmin.ne.-1).and.
289 > ((ibcthmax.eq.1).or.(ibcthmax.eq.2)))
then
290 wk(iabcth2+ix-1)=bcthmax(ix,iph)
296 if((ibcxmin.eq.1).or.(ibcxmin.eq.2))
then
297 wk(iabcx1+ith-1)=bcxmin(ith,iph)
299 if((ibcxmin.ne.-1).and.
300 > ((ibcxmax.eq.1).or.(ibcxmax.eq.2)))
then
301 wk(iabcx2+ith-1)=bcxmax(ith,iph)
307 call tcsp23(x,inx,th,inth,fspl(1,1,1,1,1,iph),inf4,
308 > ibcxmin,wk(iabcx1),ibcxmax,wk(iabcx2),
309 > ibcthmin,wk(iabcth1),ibcthmax,wk(iabcth2),
310 > wk(iaspl2),wk(iawk),inwk,ilinx,ilinth,ier)
313 write(6,*)
' ?tcspline: error in 2d spline, exiting.'
337 wk(4*(iph-1)+1)=fspl(ic1,ic2,1,ix,ith,iph)
351 if((ibcphmin.eq.1).or.(ibcphmin.eq.2)) ibcphmina=0
352 if((ibcphmax.eq.1).or.(ibcphmax.eq.2)) ibcphmaxa=0
355 call v_spline(ibcphmina,ibcphmaxa,inph,ph,wk,
361 fspl(ic1,ic2,2,ix,ith,iph)=wk(4*(iph-1)+2)
362 fspl(ic1,ic2,3,ix,ith,iph)=wk(4*(iph-1)+3)*xo2
363 fspl(ic1,ic2,4,ix,ith,iph)=wk(4*(iph-1)+4)*xo6
382 iabcph2=iabcph1+inx*inth
383 iaccoef=iabcph2+inx*inth
384 iawk=iaccoef+12*inx*inth*inph
395 if(ibcphmin.eq.1)
then
397 else if(ibcphmin.eq.2)
then
404 if(ibcphmax.eq.1)
then
406 else if(ibcphmax.eq.2)
then
420 call tcspeval(x(ix),th(ith),ph(1),iselect1, zcur,
421 > x,inx,th,inth,ph,inph,ilinx,ilinth,ilinph,
422 > fspl,inf4,inf5,ier)
424 write(6,*)
' ?? tcspline: error in tcspeval call'
427 wk(ia1)=bcphmin(ix,ith)-zcur(1)
433 call tcspeval(x(ix),th(ith),ph(inph),iselect2, zcur,
434 > x,inx,th,inth,ph,inph,ilinx,ilinth,ilinph,
435 > fspl,inf4,inf5,ier)
437 write(6,*)
' ?? tcspline: error in tcspeval call'
440 wk(ia2)=bcphmax(ix,ith)-zcur(1)
447 call tcspcorr(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
448 > ibcxmin,ibcxmax,ibcthmin,ibcthmax,
449 > ibcphmin,wk(iabcph1),ibcphmax,wk(iabcph2),
450 > wk(iaccoef),wk(iawk),inwk)
457 subroutine tcspcorr(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
458 > ibcxmin,ibcxmax,ibcthmin,ibcthmax,
459 > ibcphmin,bcph1,ibcphmax,bcph2,ccorr,wk,nwk)
471 real fspl(4,4,4,inf4,inf5,inph)
473 integer ibcxmin,ibcxmax
474 integer ibcthmin,ibcthmax
475 integer ibcphmin,ibcphmax
482 real ccorr(3,4,inx,inth,inph)
493 if(nwk.lt.10*max(inx,inth,inph))
then
494 write(6,*)
' ?? programming error in tcspcorr (tcspline)'
511 if(ibcphmin.eq.1)
then
513 else if(ibcphmin.eq.2)
then
517 if(ibcphmax.eq.1)
then
518 wk(inpho+2)=bcph2(ix,ith)
519 else if(ibcphmax.eq.2)
then
520 wk(inpho+3)=bcph2(ix,ith)
523 call v_spline(ibcphmin,ibcphmax,inph,ph,wk,wk(iawk2))
528 ccorr(1,1,ix,ith,iph)=wk(4*(iph-1)+2)
529 ccorr(2,1,ix,ith,iph)=wk(4*(iph-1)+3)*xo2
530 ccorr(3,1,ix,ith,iph)=wk(4*(iph-1)+4)*xo6
548 wk(4*(ix-1)+1)=ccorr(icph,1,ix,ith,iph)
558 call v_spline(ibcxmin,ibcxmax,inx,x,wk,wk(iawk2))
561 ccorr(icph,2,ix,ith,iph)=wk(4*(ix-1)+2)
562 ccorr(icph,3,ix,ith,iph)=wk(4*(ix-1)+3)*xo2
563 ccorr(icph,4,ix,ith,iph)=wk(4*(ix-1)+4)*xo6
586 wk(4*(ith-1)+1)=ccorr(icph,icx,ix,ith,iph)
596 call v_spline(ibcthmin,ibcthmax,inth,th,wk,
602 if(icth.eq.3) zfac=xo2
603 if(icth.eq.4) zfac=xo6
604 fspl(icx,icth,icph+1,ix,ith,iph)=
605 > fspl(icx,icth,icph+1,ix,ith,iph)+
606 > wk(4*(ith-1)+icth)*zfac
619 subroutine tcsp23(x,inx,th,inth,fspl,inf4,
620 > ibcxmin,bcxmin,ibcxmax,bcxmax,
621 > ibcthmin,bcthmin,ibcthmax,bcthmax,
622 > fspl2,wk,nwk,ilinx,ilinth,ier)
631 real fspl(4,4,4,inf4,inth)
632 real fspl2(4,4,inx,inth)
634 real bcxmin(inth),bcxmax(inth)
635 real bcthmin(inx),bcthmax(inx)
645 fspl2(1,1,ix,ith)=fspl(1,1,1,ix,ith)
651 call bcspline(x,inx,th,inth,fspl2,inx,
652 > ibcxmin,bcxmin,ibcxmax,bcxmax,
653 > ibcthmin,bcthmin,ibcthmax,bcthmax,
654 > wk,nwk,ilinx,ilinth,ier)
663 fspl(i,j,1,ix,ith)=fspl2(i,j,ix,ith)