23 subroutine r8bcspline(x,inx,th,inth,fspl,inf3,
24 > ibcxmin,bcxmin,ibcxmax,bcxmax,
25 > ibcthmin,bcthmin,ibcthmax,bcthmax,
26 > wk,nwk,ilinx,ilinth,ier)
31 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
32 INTEGER inth,inf3,ibcxmin,ibcxmax,ibcthmin,ibcthmax,nwk,ilinx
33 INTEGER ilinth,ier,inx,iflg2,ix,itest,ierx,ierth,inxo,ith
34 INTEGER intho,ic,ibcthmina,ibcthmaxa,iasc,iinc,iawk,jx,jth,ii
35 INTEGER iadr,ia5w,iaspl
38 real*8 xo2,xo6,zhxn,zhth,zdiff1,zdiff2
40 real*8 x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
41 real*8 bcxmin(*),bcxmax(*)
42 real*8 bcthmin(*),bcthmax(*)
125 if(ibcthmin.ne.-1)
then
126 if((ibcthmin.eq.1).or.(ibcthmin.eq.2))
then
128 if (bcthmin(ix).ne.0.0_r8) iflg2=1
131 if((ibcthmax.eq.1).or.(ibcthmax.eq.2))
then
133 if (bcthmax(ix).ne.0.0_r8) iflg2=1
139 itest=5*max(inx,inth)
141 itest=itest +4*inx*inth
144 if(nwk.lt.itest)
then
145 write(6,9901) nwk,itest
146 9901
format(
' ?bcspline: workspace too small:'/
147 >
' user supplied: nwk=',i6,
'; need at least: ',i6/
148 >
' nwk=4*nx*ny +5*max(nx,ny) will work for any user'/
149 >
' choice of bdy conditions.')
153 write(6,
'('' ?bcspline: at least 2 x points required.'')')
157 write(6,
'('' ?bcspline: need at least 2 theta points.'')')
161 call ibc_ck(ibcxmin,
'bcspline',
'xmin',-1,7,ier)
162 if(ibcxmin.ge.0)
call ibc_ck(ibcxmax,
'bcspline',
'xmax',0,7,ier)
163 call ibc_ck(ibcthmin,
'bcspline',
'thmin',-1,7,ier)
164 if(ibcthmin.ge.0)
call ibc_ck(ibcthmax,
'bcspline',
'thmax',0,7,ier)
168 call r8splinck(x,inx,ilinx,1.0e-3_r8,ierx)
172 write(6,
'('' ?bcspline: x axis not strict ascending'')')
177 call r8splinck(th,inth,ilinth,1.0e-3_r8,ierth)
181 write(6,
'('' ?bcspline: th axis not strict ascending'')')
199 wk(4*(ix-1)+1)=fspl(1,1,ix,ith)
202 if(ibcxmin.eq.1)
then
204 else if(ibcxmin.eq.2)
then
208 if(ibcxmax.eq.1)
then
209 wk(inxo+2)=bcxmax(ith)
210 else if(ibcxmax.eq.2)
then
211 wk(inxo+3)=bcxmax(ith)
216 call r8v_spline(ibcxmin,ibcxmax,inx,x,wk,wk(4*inx+1))
221 fspl(2,1,ix,ith)=wk(4*(ix-1)+2)
222 fspl(3,1,ix,ith)=wk(4*(ix-1)+3)*xo2
223 fspl(4,1,ix,ith)=wk(4*(ix-1)+4)*xo6
242 wk(4*(ith-1)+1)=fspl(ic,1,ix,ith)
256 if((ibcthmin.eq.1).or.(ibcthmin.eq.2)) ibcthmina=0
257 if((ibcthmax.eq.1).or.(ibcthmax.eq.2)) ibcthmaxa=0
260 call r8v_spline(ibcthmina,ibcthmaxa,inth,th,wk,wk(4*inth+1))
265 fspl(ic,2,ix,ith)=wk(4*(ith-1)+2)
266 fspl(ic,3,ix,ith)=wk(4*(ith-1)+3)*xo2
267 fspl(ic,4,ix,ith)=wk(4*(ith-1)+4)*xo6
286 zhth=th(inth)-th(inth-1)
293 if(ibcthmin.eq.1) iselect1(3)=1
294 if(ibcthmin.eq.2) iselect1(5)=1
295 if(ibcthmax.eq.1) iselect2(3)=1
296 if(ibcthmax.eq.2) iselect2(5)=1
304 if(ibcthmin.eq.1)
then
306 zcur(1)=fspl(1,2,ix,1)
308 zcur(1)=fspl(1,2,jx,1)+zhxn*(fspl(2,2,jx,1)+zhxn*
309 > (fspl(3,2,jx,1)+zhxn*fspl(4,2,jx,1)))
311 zdiff1=bcthmin(ix)-zcur(1)
312 else if(ibcthmin.eq.2)
then
314 zcur(1)=2.0_r8*fspl(1,3,ix,1)
317 > (fspl(1,3,jx,1)+zhxn*(fspl(2,3,jx,1)+zhxn*
318 > (fspl(3,3,jx,1)+zhxn*fspl(4,3,jx,1))))
320 zdiff1=bcthmin(ix)-zcur(1)
327 if(ibcthmax.eq.1)
then
331 > fspl(1,2,ix,jth)+zhth*(2.0_r8*fspl(1,3,ix,jth)+
332 > zhth*3.0_r8*fspl(1,4,ix,jth))
334 call r8bcspeval(x(inx),th(inth),iselect2, zcur,
335 > x,inx,th,inth,ilinx,ilinth,fspl,inf3,ier)
338 zdiff2=bcthmax(ix)-zcur(1)
339 else if(ibcthmax.eq.2)
then
342 zcur(1)=2.0_r8*fspl(1,3,ix,jth)+
343 > 6.0_r8*zhth*fspl(1,4,ix,jth)
345 call r8bcspeval(x(inx),th(inth),iselect2, zcur(1),
346 > x,inx,th,inth,ilinx,ilinth,fspl,inf3,ier)
349 zdiff2=bcthmax(ix)-zcur(1)
358 iadr=iasc+(ix-1)*iinc
360 wk(iadr+4*(ith-1))=0.0_r8
365 wk(iadr+intho+1)=0.0_r8
366 wk(iadr+intho+2)=0.0_r8
368 if(ibcthmin.eq.1)
then
370 else if(ibcthmin.eq.2)
then
374 if(ibcthmax.eq.1)
then
375 wk(iadr+intho+1)=zdiff2
376 else if(ibcthmax.eq.2)
then
377 wk(iadr+intho+2)=zdiff2
380 call r8v_spline(ibcthmin,ibcthmax,inth,th,wk(iadr),wk(iawk))
386 iadr=iasc+(ix-1)*iinc
388 wk(iadr+4*(ith-1)+2)=wk(iadr+4*(ith-1)+2)*xo2
389 wk(iadr+4*(ith-1)+3)=wk(iadr+4*(ith-1)+3)*xo6
391 fspl(1,2,ix,ith)=fspl(1,2,ix,ith)+wk(iadr+4*(ith-1)+1)
392 fspl(1,3,ix,ith)=fspl(1,3,ix,ith)+wk(iadr+4*(ith-1)+2)
393 fspl(1,4,ix,ith)=fspl(1,4,ix,ith)+wk(iadr+4*(ith-1)+3)
405 iaspl=iasc+iinc*(ix-1)
406 wk(iawk+4*(ix-1))=wk(iaspl+4*(ith-1)+(ic-1))
413 wk(iawk+inxo+1)=0.0_r8
414 wk(iawk+inxo+2)=0.0_r8
419 call r8v_spline(ibcxmin,ibcxmax,inx,x,wk(iawk),wk(ia5w))
422 fspl(2,ic,ix,ith)=fspl(2,ic,ix,ith)+
423 > wk(iawk+4*(ix-1)+1)
424 fspl(3,ic,ix,ith)=fspl(3,ic,ix,ith)+
425 > wk(iawk+4*(ix-1)+2)*xo2
426 fspl(4,ic,ix,ith)=fspl(4,ic,ix,ith)+
427 > wk(iawk+4*(ix-1)+3)*xo6