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