V3FIT
r8mkbicubw.f
1  subroutine r8mkbicubw(x,nx,y,ny,f,nf2,
2  > ibcxmin,bcxmin,ibcxmax,bcxmax,
3  > ibcymin,bcymin,ibcymax,bcymax,
4  > wk,nwk,
5  > ilinx,iliny,ier)
6 C
7 C setup a bicubic spline; store coefficients in compact form
8 C (as per suggestion of L. Zakharov, PPPL, Feb. 1999)
9 C 4 coeffs per grid point: f,fxx,fyy,fxxyy
10 C
11 C
12 C input:
13 !============
14 ! idecl: explicitize implicit INTEGER declarations:
15  IMPLICIT NONE
16  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
17  INTEGER itest,iadfp,isiz1,iadfw,inwk
18 !============
19 ! idecl: explicitize implicit REAL declarations:
20  real*8 hxlast,hylast
21 !============
22  integer nx ! length of x vector
23  integer ny ! length of y vector
24  real*8 x(nx) ! x vector, strict ascending
25  real*8 y(ny) ! y vector, strict ascending
26 C
27  integer nf2 ! 2nd dimension of f, nf2.ge.nx
28 C input/output:
29 C
30  real*8 f(4,nf2,ny) ! data & spline coefficients
31 C
32 C on input: f(1,i,j) = f(x(i),y(j))
33 C on output: f(1,i,j) unchanged
34 C f(2,i,j) = d2f/dx2(x(i),y(j))
35 C f(3,i,j) = d2f/dy2(x(i),y(j))
36 C f(4,i,j) = d4f/dx2dy2(x(i),y(j))
37 C
38 C and the interpolation formula for (x,y) in (x(i),x(i+1))^(y(j),y(j+1))
39 C is:
40 C hx = x(i+1)-x(i) hy = y(j+1)-y(j)
41 C dxp= (x-x(i))/hx dxm= 1-dxp dxp,dxm in (0,1)
42 C dyp= (x-x(i))/hx dym= 1-dyp dyp,dym in (0,1)
43 C dx3p = dxp**3-dxp dx3m = dxm**3-dxm dxp3,dxm3 in (0,1)
44 C
45 C finterp = dxm*(dym*f(1,i,j)+dyp*f(1,i,j+1))
46 C +dxp*(dym*f(1,i+1,j)+dyp*f(1,i+1,j+1))
47 C +1/6*hx**2*
48 C dx3m*(dym*f(2,i,j)+dyp*f(2,i,j+1))
49 C +dx3p*(dym*f(2,i+1,j)+dyp*f(2,i+1,j+1))
50 C +1/6*hy**2*
51 C dxm*(dy3m*f(3,i,j)+dy3p*f(3,i,j+1))
52 C +dxp*(dy3m*f(3,i+1,j)+dy3p*f(3,i+1,j+1))
53 C +1/36*hx**2*hy**2*
54 C dx3m*(dym*f(4,i,j)+dyp*f(4,i,j+1))
55 C +dx3p*(dym*f(4,i+1,j)+dyp*f(4,i+1,j+1))
56 C
57 C where the f(2:4,*,*) are cleverly evaluated to assure
58 C (a) finterp is continuous and twice differentiable across all
59 C grid cell boundaries, and
60 C (b) all boundary conditions are satisfied.
61 C
62 C the boundary conditions specification options are:
63 C inputs:
64 C
65  integer ibcxmin ! bc flag for x=xmin
66  REAL*8 bcxmin(ny) ! bc data vs. y at x=xmin
67  integer ibcxmax ! bc flag for x=xmax
68  REAL*8 bcxmax(ny) ! bc data vs. y at x=xmax
69 C
70  integer ibcymin ! bc flag for y=ymin
71  real*8 bcymin(nx) ! bc data vs. x at y=ymin
72  integer ibcymax ! bc flag for y=ymax
73  REAL*8 bcymax(nx) ! bc data vs. x at y=ymax
74 C
75 C with interpretation:
76 c ibcxmin -- indicator for boundary condition at x(1):
77 c bcxmin(...) -- boundary condition data
78 c =-1 -- periodic boundary condition
79 c =0 -- use "not a knot"
80 c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
81 c =2 -- match 2nd derivative, specified at x(1),th(ith) by bcxmin(ith)
82 c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
83 c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
84 c =5 -- match 1st derivative to 1st divided difference
85 c =6 -- match 2nd derivative to 2nd divided difference
86 c =7 -- match 3rd derivative to 3rd divided difference
87 c (for more detailed definition of BCs 5-7, see the
88 c comments of subroutine mkspline)
89 c NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
90 c
91 c ibcxmax -- indicator for boundary condition at x(nx):
92 c bcxmax(...) -- boundary condition data
93 c (interpretation as with ibcxmin, bcxmin)
94 c NOTE: if ibcxmin=-1, ibcxmax is ignored! ...and the BC is periodic.
95 c
96 C and analogous interpretation for ibcymin,bcymin,ibcymax,bcymax
97 C (df/dy or d2f/dy2 boundary conditions at y=ymin and y=ymax).
98 C
99 C *** workspace ***
100 C This FORTRAN-77 routine requires a workspace of
101 C nwk = (AT LEAST) 21*nx*ny
102 C *** for dynamic allocation of this workspace use subroutine mkbicub
103 C which has the same arguments as mkbicubw, except that the workspace
104 C does not have to be provided -- mkbicub is FORTRAN-90.
105 C
106 C input:
107  integer nwk ! size of workspace
108 C modified on output:
109  REAL*8 wk(nwk) ! workspace
110 C
111 C and output arguments
112 C
113  integer ilinx ! =1: x grid is "nearly" equally spaced
114  integer iliny ! =1: y grid is "nearly" equally spaced
115 C
116 C ilinx and iliny are set to zero if corresponding grids are not equally
117 C spaced
118 C
119 C and completion code
120 C
121  integer ier ! =0 on exit if there is no error.
122 C
123 C if there is an error, ier is set and a message is output on unit 6.
124 C these are considered programming errors in the calling routine.
125 C
126 C possible errors:
127 C x(...) not strict ascending
128 C y(...) not strict ascending
129 C nx .lt. 4
130 C ny .lt. 4
131 C invalid boundary condition flag
132 C workspace too small
133 C
134 C-----------------------------------------------------
135 C local arrays
136 C
137  integer iselect(10)
138  integer zvalues(10)
139 C
140  data iselect/-1,0,0,0,0,0,0,0,0,0/
141 C
142 C-----------------------------------------------------
143 C
144  itest=21*nx*ny
145  if(nwk.lt.itest) then
146  write(6,9901) nwk,itest
147  9901 format(' ?mkbicubw: workspace too small:'/
148  > ' user supplied: nwk=',i7,'; need at least: ',i7/
149  > ' nwk = at least 21*nx*ny is required.')
150  ier=1
151  return
152  endif
153 C
154  iadfp=1
155  isiz1=16*nx*ny
156  iadfw=iadfp+isiz1
157  inwk = nwk-isiz1
158 C
159  call r8mkbicop(f,nf2,wk(iadfp),nx,ny)
160 C
161 C evaluate 4x4 continuous bicubic spline
162 C
163  call r8bcspline(x,nx,y,ny,wk(iadfp),nx,
164  > ibcxmin,bcxmin,ibcxmax,bcxmax,
165  > ibcymin,bcymin,ibcymax,bcymax,
166  > wk(iadfw),inwk,
167  > ilinx,iliny,ier)
168 C
169  if(ier.ne.0) return
170 C
171 C convert to 4-coefficient (2nd/4th partials form)
172 C
173  hxlast=x(nx)-x(nx-1)
174  hylast=y(ny)-y(ny-1)
175  call r8mkbicon(f,nf2,wk(iadfp),nx,ny,hxlast,hylast)
176 C
177  return
178  end
179 C
180 C----------------------------------------------------------------
181 C mkbicop -- copy spline function input data
182 C
183  subroutine r8mkbicop(fin,nf2,fwk,nx,ny)
184 C
185 !============
186 ! idecl: explicitize implicit INTEGER declarations:
187  IMPLICIT NONE
188  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
189  INTEGER nx,ny,nf2,iy,ix
190 !============
191  real*8 fin(4,nf2,ny)
192  real*8 fwk(4,4,nx,ny)
193 C
194  do iy=1,ny
195  do ix=1,nx
196  fwk(1,1,ix,iy)=fin(1,ix,iy)
197  enddo
198  enddo
199 C
200  return
201  end
202 C----------------------------------------------------------------
203 C mkbicon -- create compact spline representation from 4x4
204 C (bcspline) representation
205 C
206  subroutine r8mkbicon(fin,nf2,fwk,nx,ny,hxlast,hylast)
207 C
208 !============
209 ! idecl: explicitize implicit INTEGER declarations:
210  IMPLICIT NONE
211  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
212  INTEGER nx,ny,nf2,iy,ix,iflag,ixuse,iyuse,j
213 !============
214 ! idecl: explicitize implicit REAL declarations:
215  real*8 hxlast,hylast,dxuse,dyuse
216 !============
217  real*8 fin(4,nf2,ny)
218  real*8 fwk(4,4,nx,ny)
219 C
220 C-----------------------------------------------------
221 C local arrays
222 C
223  integer iselect(10)
224  REAL*8 zvalues(10)
225 C
226  data iselect/-1,0,0,0,0,0,0,0,0,0/
227 C
228 C-----------------------------------------------------
229 C
230  do iy=1,ny
231  do ix=1,nx
232 C
233 C copy derivatives from result. Special treatment needed for end zones
234 C
235  iflag=0
236  dxuse=0.0_r8
237  dyuse=0.0_r8
238  ixuse=ix
239  iyuse=iy
240  if(ix.eq.nx) then
241  iflag=1
242  dxuse=hxlast
243  ixuse=ix-1
244  endif
245  if(iy.eq.ny) then
246  iflag=1
247  dyuse=hylast
248  iyuse=iy-1
249  endif
250 C
251  if(iflag.eq.1) then
252  call r8bcspevfn(iselect,1,1,zvalues,
253  > ixuse,iyuse,dxuse,dyuse,
254  > fwk,nx,ny)
255  do j=2,4
256  fin(j,ix,iy)=zvalues(j)
257  enddo
258  else
259  fin(2,ix,iy)=2.0_r8*fwk(3,1,ix,iy)
260  fin(3,ix,iy)=2.0_r8*fwk(1,3,ix,iy)
261  fin(4,ix,iy)=4.0_r8*fwk(3,3,ix,iy)
262  endif
263 C
264  enddo
265  enddo
266 C
267  return
268  end