V3FIT
mkspl3pb.f
1  subroutine mkspl3pb(fun,x,nx,th,nth,ph,nph,fspl,nf4,nf5,
2  > ibcxmin,bcxmin,ibcxmax,bcxmax,nb1,
3  > wk,inwk,ilinx,ilinth,ilinph,ier)
4 C
5 C create a tricubic biperiodic spline with knots at grid points and
6 C function values from the callable function `fun' passed.
7 C
8 C use tpsplinb to setup the spline coefficients
9 C
10 C periodic boundary condition for theta & phi grids.
11 C boundary condition data may be specified at x(1) and x(nx) for each
12 C theta & phi:
13 C ibcxmin=0 ==> "not a knot" boundary condition, cf cubspl.for
14 C ibcxmin=1 ==> match slope, bcxmin(ith,iph) gives df/dx at
15 C x(1),th(ith),ph(iph)
16 C ibcxmin=2 ==> match 2nd derivative, given at x(1),th(ith),ph(iph)
17 C by bcxmin(ith,iph)
18 C
19 C ibcxmax,bcxmax have analogous interpretation -- at x(nx).
20 C
21  external fun ! passed real function(x,th)
22  real x(nx) ! x coordinate array
23  real th(nth) ! th coordinate array
24  real ph(nph) ! ph coordinate array
25 C
26  real fspl(4,4,4,nf4,nf5,nph) ! function data / spline coeff array
27  real wk(inwk) ! workspace for bpsplinb
28 C
29  integer ibcxmin ! boundary condition indicator
30  real bcxmin(nb1,nph) ! boundary condition data
31  integer ibcxmax ! boundary condition indicator
32  real bcxmax(nb1,nph) ! boundary condition data
33 C
34  integer ilinx ! output =1 if x(...) evenly spaced
35  integer ilinth ! output =1 if th(...) evenly spaced
36  integer ilinph ! output =1 if ph(...) evenly spaced
37 C
38  integer ier ! completion code from bpspline
39 C
40 C----------------------------
41 C
42  ier=0
43  if(nf4.lt.nx) then
44  write(6,'('.lt.' ?mkspl3pb -- array dim error, nf4 nx'')')
45  ier=1
46  endif
47  if(nf5.lt.nth) then
48  write(6,'('.lt.' ?mkspl3pb -- array dim error, nf5 nth'')')
49  ier=2
50  endif
51  if(nb1.lt.nth) then
52  write(6,'('.lt.' ?mkspl3pb -- array dim error, nb1 nth'')')
53  ier=3
54  endif
55  if(inwk.lt.5*max(nx,nth,nph)) then
56  write(6,'('' ?mkspl3pb -- array dim error, inwk too small'')')
57  ier=4
58  endif
59 
60  if(ier.ne.0) return
61 C
62  do iph=1,nph
63  do ith=1,nth
64  do ix=1,nx
65  fspl(1,1,1,ix,ith,iph)=fun(x(ix),th(ith),ph(iph))
66  enddo
67  enddo
68  enddo
69 C
70  call tpsplinb(x,nx,th,nth,ph,nph,fspl,nf4,nf5,
71  > ibcxmin,bcxmin,ibcxmax,bcxmax,nb1,
72  > wk,inwk,ilinx,ilinth,ilinph,ier)
73 C
74  return
75  end