V3FIT
mkbicub.f
1  subroutine mkbicub(x,nx,y,ny,f,nf2,
2  > ibcxmin,bcxmin,ibcxmax,bcxmax,
3  > ibcymin,bcymin,ibcymax,bcymax,
4  > ilinx,iliny,ier)
5 C
6 C setup bicubic spline, dynamic allocation of workspace
7 C fortran-90 fixed form source
8 C
9 C ** see mkbicubw.for **
10 C arguments are identical, except the workspace argument is
11 C omitted (created here)
12 C
13 C --NOTE-- dmc 22 Feb 2004 -- rewrite for direct calculation of
14 C coefficients, to avoid large transient use of memory. mkbicubw is
15 C no longer called (but reference to mkbicubw comments is still correct).
16 C
17 C abbreviated description of arguments, mkbicubw has details:
18 C
19  implicit NONE
20 C
21 C input:
22  integer nx ! length of x vector
23  integer ny ! length of y vector
24  real x(nx) ! x vector, strict ascending
25  real y(ny) ! y vector, strict ascending
26 C
27  integer nf2 ! 2nd dimension of f, nf2.ge.nx
28 C input/output:
29  real f(4,nf2,ny) ! data & spline coefficients
30 C
31 C input bdy condition data:
32  integer ibcxmin ! bc flag for x=xmin
33  real bcxmin(ny) ! bc data vs. y at x=xmin
34  integer ibcxmax ! bc flag for x=xmax
35  real bcxmax(ny) ! bc data vs. y at x=xmax
36 C
37  integer ibcymin ! bc flag for y=ymin
38  real bcymin(nx) ! bc data vs. x at y=ymin
39  integer ibcymax ! bc flag for y=ymax
40  real bcymax(nx) ! bc data vs. x at y=ymax
41 C
42 C output linear grid flags and completion code (ier=0 is normal):
43 C
44  integer ilinx ! =1: x grid is "nearly" equally spaced
45  integer iliny ! =1: y grid is "nearly" equally spaced
46 C
47  integer ier ! =0 on exit if there is no error.
48 C
49 C-----------------------
50  integer ierx,iery
51 C
52  real, dimension(:,:), allocatable :: fwk
53  real :: zbcmin,zbcmax
54  integer ix,iy,ibcmin,ibcmax
55 c
56  real, dimension(:,:,:), allocatable :: fcorr
57  integer iflg2
58  real zdiff(2),hy
59 C
60 C-----------------------
61 c
62 c see if 2nd pass is needed due to inhomogeneous d/dy bdy cond.
63 c
64  iflg2=0
65  if(ibcymin.ne.-1) then
66  if((ibcymin.eq.1).or.(ibcymin.eq.2)) then
67  do ix=1,nx
68  if (bcymin(ix).ne.0.0) iflg2=1
69  enddo
70  endif
71  if((ibcymax.eq.1).or.(ibcymax.eq.2)) then
72  do ix=1,nx
73  if (bcymax(ix).ne.0.0) iflg2=1
74  enddo
75  endif
76  endif
77 c
78 c check boundary condition specifications
79 c
80  ier=0
81 c
82  call ibc_ck(ibcxmin,'bcspline','xmin',-1,7,ier)
83  if(ibcxmin.ge.0) call ibc_ck(ibcxmax,'bcspline','xmax',0,7,ier)
84  call ibc_ck(ibcymin,'bcspline','ymin',-1,7,ier)
85  if(ibcymin.ge.0) call ibc_ck(ibcymax,'bcspline','ymax',0,7,ier)
86 c
87 c check ilinx & x vector
88 c
89  call splinck(x,nx,ilinx,1.0e-3,ierx)
90  if(ierx.ne.0) ier=2
91 c
92  if(ier.eq.2) then
93  write(6,'('' ?bcspline: x axis not strict ascending'')')
94  endif
95 c
96 c check iliny & y vector
97 c
98  call splinck(y,ny,iliny,1.0e-3,iery)
99  if(iery.ne.0) ier=3
100 c
101  if(ier.eq.3) then
102  write(6,'('' ?bcspline: y axis not strict ascending'')')
103  endif
104 c
105  if(ier.ne.0) return
106 c
107 c------------------------------------
108  allocate(fwk(2,max(nx,ny)))
109 c
110 c evaluate fxx (spline in x direction)
111 c
112  zbcmin=0
113  zbcmax=0
114  do iy=1,ny
115  fwk(1,1:nx) = f(1,1:nx,iy)
116  if((ibcxmin.eq.1).or.(ibcxmin.eq.2)) zbcmin=bcxmin(iy)
117  if((ibcxmax.eq.1).or.(ibcxmax.eq.2)) zbcmax=bcxmax(iy)
118  call mkspline(x,nx,fwk,
119  > ibcxmin,zbcmin,ibcxmax,zbcmax,ilinx,ier)
120  if(ier.ne.0) return
121  f(2,1:nx,iy)=fwk(2,1:nx)
122  enddo
123 c
124 c evaluate fyy (spline in y direction)
125 c use homogeneous boundary condition; correction done later if necessary
126 c
127  zbcmin=0
128  zbcmax=0
129  ibcmin=ibcymin
130  ibcmax=ibcymax
131  do ix=1,nx
132  fwk(1,1:ny) = f(1,ix,1:ny)
133  if(iflg2.eq.1) then
134  if((ibcymin.eq.1).or.(ibcymin.eq.2)) ibcmin=0
135  if((ibcymax.eq.1).or.(ibcymax.eq.2)) ibcmax=0
136  endif
137  call mkspline(y,ny,fwk,
138  > ibcmin,zbcmin,ibcmax,zbcmax,iliny,ier)
139  if(ier.ne.0) return
140  f(3,ix,1:ny)=fwk(2,1:ny)
141  enddo
142 c
143 c evaluate fxxyy (spline fxx in y direction; BC simplified; avg
144 c d2(d2f/dx2)/dy2 and d2(df2/dy2)/dx2
145 c
146  zbcmin=0
147  zbcmax=0
148  ibcmin=ibcymin
149  ibcmax=ibcymax
150  do ix=1,nx
151  fwk(1,1:ny) = f(2,ix,1:ny)
152  if(iflg2.eq.1) then
153  if((ibcymin.eq.1).or.(ibcymin.eq.2)) ibcmin=0
154  if((ibcymax.eq.1).or.(ibcymax.eq.2)) ibcmax=0
155  endif
156  call mkspline(y,ny,fwk,
157  > ibcmin,zbcmin,ibcmax,zbcmax,iliny,ier)
158  if(ier.ne.0) return
159  f(4,ix,1:ny)= fwk(2,1:ny)
160  enddo
161 c
162  if(iflg2.eq.1) then
163  allocate(fcorr(2,nx,ny))
164 c
165 c correct for inhomogeneous y boundary condition
166 c
167  do ix=1,nx
168  ! the desired inhomogenous BC is the difference btw the
169  ! requested derivative (1st or 2nd) and the current value
170 
171  zdiff(1)=0.0
172  if(ibcymin.eq.1) then
173  hy=y(2)-y(1)
174  zdiff(1)=(f(1,ix,2)-f(1,ix,1))/hy +
175  > hy*(-2*f(3,ix,1)-f(3,ix,2))/6
176  zdiff(1)=bcymin(ix)-zdiff(1)
177  else if(ibcymin.eq.2) then
178  zdiff(1)=bcymin(ix)-f(3,ix,1)
179  endif
180 
181  zdiff(2)=0.0
182  if(ibcymax.eq.1) then
183  hy=y(ny)-y(ny-1)
184  zdiff(2)=(f(1,ix,ny)-f(1,ix,ny-1))/hy +
185  > hy*(2*f(3,ix,ny)+f(3,ix,ny-1))/6
186  zdiff(2)=bcymax(ix)-zdiff(2)
187  else if(ibcymax.eq.2) then
188  zdiff(2)=bcymax(ix)-f(3,ix,ny)
189  endif
190 c
191  fwk(1,1:ny)=0.0 ! values are zero; only BC is not
192  call mkspline(y,ny,fwk,ibcymin,zdiff(1),ibcymax,zdiff(2),
193  > iliny,ier)
194  if(ier.ne.0) return
195  fcorr(1,ix,1:ny)=fwk(2,1:ny) ! fyy-correction
196  enddo
197 c
198  zbcmin=0
199  zbcmax=0
200  do iy=1,ny
201  fwk(1,1:nx)=fcorr(1,1:nx,iy)
202  call mkspline(x,nx,fwk,ibcxmin,zbcmin,ibcxmax,zbcmax,
203  > ilinx,ier)
204  if(ier.ne.0) return
205  fcorr(2,1:nx,iy)=fwk(2,1:nx) ! fxxyy-correction
206  enddo
207 c
208  f(3:4,1:nx,1:ny)=f(3:4,1:nx,1:ny)+fcorr(1:2,1:nx,1:ny)
209 c
210  deallocate(fcorr)
211  endif
212 c
213 c correction spline -- f=fxx=zero; fyy & fxxyy are affected
214 c
215  deallocate(fwk)
216 c------------------------------------
217 C
218 C thats all
219 C
220  return
221  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17