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