V3FIT
mkintrp2d.f
1  subroutine mkintrp2d(x,nx,y,ny,jspline,
2  > f,icoeff,ixdim,iydim,
3  > ibcxmin,bcxmin,ibcxmax,bcxmax,
4  > ibcymin,bcymin,ibcymax,bcymax,
5  > ier)
6 C
7 C setup bicubic spline, or bicubic hermite, or Hybrid linear/zonal with
8 C 1d cubic Hermite or Spline interpolation in 1 dimension
9 C
10  implicit NONE
11 C
12 C input:
13  integer nx ! length of x vector
14  integer ny ! length of y vector
15  real x(nx) ! x vector, strict ascending
16  real y(ny) ! y vector, strict ascending
17 C
18  integer :: jspline(2) ! interpolation method control
19 C (1) -- 1st dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
20 C (2) -- 2nd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
21 C
22 C Standard interpolation-- all jspline values match
23 C e.g. jspline(1)=jspline(2)=2 for bicubic spline
24 C Hybrid interpolation occurs when multiple jspline values are specified.
25 C
26 C RESTRICTION: jspline = (/ 1, 2 /) or (/ 2, 1 /) not allowed. I.e.
27 C Hermite and Spline interpolation cannot be mixed in a hybrid interpolant.
28 C This restriction exists because of technical issues in the
29 C implementation (it could be removed in principle but the work to do
30 C this has not been scheduled).
31 C
32 C coefficient buffer dimensioning:
33  integer :: icoeff ! #coefficients per data point
34  integer :: ixdim ! nx; nx-1 if jspline(1)==-1
35  integer :: iydim ! ny; ny-1 if jspline(2)==-1
36 C input/output:
37  real f(icoeff,ixdim,iydim) ! data & spline coefficients
38 C
39 C input bdy condition data:
40 C meaning of arguments as described in "mkbicubw.for"
41 C *ignored* for dimensions with piecewise linear or zonal interpolation;
42 C For Hermite interpolation only the values {-1,0,1} are accepted.
43 C
44  integer ibcxmin ! bc flag for x=xmin
45  real bcxmin(iydim) ! bc data vs. y at x=xmin
46  integer ibcxmax ! bc flag for x=xmax
47  real bcxmax(iydim) ! bc data vs. y at x=xmax
48 C
49  integer ibcymin ! bc flag for y=ymin
50  real bcymin(ixdim) ! bc data vs. x at y=ymin
51  integer ibcymax ! bc flag for y=ymax
52  real bcymax(ixdim) ! bc data vs. x at y=ymax
53 C
54  integer ier ! =0 on exit if there is no error.
55 c ier -- completion code, 0 for normal
56 C
57 C-----------------------
58  integer :: kspline
59  integer :: ii,imul,imin,imax,ickx,icky
60  integer :: idum1,idum2
61  integer :: ipx,ipy
62  real :: ztol = 0.0001
63  real, dimension(:,:), allocatable :: wk2
64 C-----------------------
65 c
66  ier = 0
67 c
68  imin=3
69  imax=-2
70  imul=1
71  do ii=1,2
72  imin=min(imin,jspline(ii))
73  imax=max(imax,jspline(ii))
74  if(jspline(ii).gt.0) imul=imul*2
75  enddo
76 
77  if((imin.lt.-1).or.(imax.gt.2)) then
78  ier = 1
79  write(6,*)
80  > ' ?mkintrp2d: spline type control out of range -1 to 2: ',
81  > jspline
82  endif
83  if(ier.ne.0) return
84 
85  if(imin.eq.imax) then
86  kspline=imin ! same interp type on all dimensions
87  else
88  kspline=-99 ! hybrid
89  if(imin.gt.0) then
90  ier = 1
91  write(6,*)
92  > ' ?mkintrp2d: spline/Hermite hybrid not supported (',
93  > jspline,')'
94  else if(imax.lt.1) then
95  ier = 1
96  write(6,*)
97  > ' ?mkintrp2d: zonal/linear hybrid not supported (',
98  > jspline,')'
99  endif
100  endif
101  if(ier.ne.0) return
102 c
103  if(imul.ne.icoeff) then
104  write(6,*)
105  > ' ?coeff dimension inconsistency for spline type codes ',
106  > jspline
107  write(6,*) ' in mkintrp2d: expected: ',imul,' got: ',icoeff
108  ier=1
109  return
110  endif
111 c
112 c check dimensioning consistency
113 c
114  if(jspline(1).eq.-1) then
115  ickx=nx-1
116  else
117  ickx=nx
118  endif
119 c
120  if(jspline(2).eq.-1) then
121  icky=ny-1
122  else
123  icky=ny
124  endif
125 c
126  if((ickx.ne.ixdim).or.(icky.ne.iydim)) then
127  write(6,*)
128  > ' ?mkintrp2d: dimensioning inconsistent with '//
129  > 'interpolation controls: ',jspline
130  write(6,*) ' expected: ',ickx,icky,'; got: ',ixdim,iydim
131  ier=1
132  return
133  endif
134 c
135  if(jspline(1).le.0) then
136  call splinck(x,nx,idum1,ztol,ier)
137  if(ier.ne.0) then
138  write(6,*) ' ?mkintrp2d: x axis not strict ascending.'
139  return
140  endif
141  endif
142 c
143  if(jspline(2).le.0) then
144  call splinck(y,ny,idum1,ztol,ier)
145  if(ier.ne.0) then
146  write(6,*) ' ?mkintrp2d: y axis not strict ascending.'
147  return
148  endif
149  endif
150 c
151 c if no work to be done: exit now
152  if(imul.eq.1) return
153 c
154 c check Hermite BCs if necessary
155 c
156  if(jspline(1).eq.1) then
157  if((min(ibcxmin,ibcxmax).lt.-1).or.
158  > (max(ibcxmin,ibcxmax).gt.1)) then
159  write(6,*) ' ?mkintrp2d: Bdy Cond code out of range for'
160  write(6,*) ' Hermite interpolation; (-1:1) allowed, '//
161  > 'found: ',ibcxmin,ibcxmax
162  ier=1
163  return
164  endif
165  ipx=0
166  if(ibcxmin.eq.-1) then
167  ipx=1
168  else if((ibcxmin.eq.1).or.(ibcxmax.eq.1)) then
169  ipx=2
170  endif
171  endif
172 c
173  if(jspline(2).eq.1) then
174  if((min(ibcymin,ibcymax).lt.-1).or.
175  > (max(ibcymin,ibcymax).gt.1)) then
176  write(6,*) ' ?mkintrp2d: Bdy Cond code out of range for'
177  write(6,*) ' Hermite interpolation; (-1:1) allowed, '//
178  > 'found: ',ibcymin,ibcymax
179  ier=1
180  return
181  endif
182  ipy=0
183  if(ibcymin.eq.-1) then
184  ipy=1
185  else if((ibcymin.eq.1).or.(ibcymax.eq.1)) then
186  ipy=2
187  endif
188  endif
189 
190  if(kspline.eq.1) then
191  ! bicubic Hermite
192 
193  ! stuff the BCs inside the function data at the right locations:
194  call util_bcherm2(f, ixdim, iydim,
195  > ibcxmin, ibcxmax, ibcymin, ibcymax,
196  > bcxmin, bcxmax, bcymin, bcymax,
197  > x, y)
198 
199  call r8akherm2p(x,ixdim,y,iydim,
200  > f,ixdim,idum1,idum2,ipx,ipy,ier)
201 
202  else if(kspline.eq.2) then
203  ! bicubic Spline
204 
205  call mkbicub(x,nx, y,ny, f,ixdim,
206  > ibcxmin,bcxmin,ibcxmax,bcxmax,
207  > ibcymin,bcymin,ibcymax,bcymax,
208  > idum1,idum2,ier)
209 
210  else
211  ! Hybrid
212 
213  if(jspline(1).gt.0) then
214  ! cubic in x direction
215  do ii=1,iydim
216  if(jspline(1).eq.1) then
217  call util_bcherm1(f(1,1,ii), ixdim,
218  > ibcxmin, ibcxmax, bcxmin(ii), bcxmax(ii), x)
219  call akherm1p(x,ixdim,f(1,1,ii),idum1,ipx,ier)
220 
221  else if(jspline(1).eq.2) then
222  call mkspline(x,ixdim,f(1,1,ii),
223  > ibcxmin,bcxmin(ii), ibcxmax,bcxmax(ii),
224  > idum1,ier)
225  endif
226  if(ier.ne.0) exit
227  enddo
228 
229  else
230  ! cubic in y direction
231  allocate(wk2(2,ny))
232 
233  do ii=1,ixdim
234  wk2(1,1:iydim) = f(1,ii,1:iydim)
235  wk2(2,1:iydim) = 0.0
236  if(jspline(2).eq.1) then
237  call util_bcherm1(wk2, iydim,
238  > ibcymin, ibcymax, bcymin(ii), bcymax(ii), y)
239  call akherm1p(y,iydim,wk2,idum2,ipy,ier)
240 
241  else if(jspline(2).eq.2) then
242  call mkspline(y,iydim,wk2,
243  > ibcymin,bcymin(ii), ibcymax,bcymax(ii),
244  > idum2,ier)
245 
246  endif
247  if(ier.ne.0) exit
248  f(1:2,ii,1:iydim) = wk2(1:2,1:iydim)
249  enddo
250 
251  deallocate(wk2)
252  endif
253 
254  endif
255 
256  return
257  end