V3FIT
akherm2.f
1  subroutine akherm2(x,nx,y,ny,fherm,nf2,ilinx,iliny,ier)
2 C
3 C create a data set for Hermite interpolation, based on Akima's method
4 C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
5 C
6 C input:
7 C
8  integer nx,ny,nf1 ! array dimensions
9  real x(nx) ! x coordinate array
10  real y(ny) ! y coordinate array
11  real fherm(0:3,nf2,ny) ! data/Hermite array
12 C
13 C fherm(0,i,j) = function value f at x(i),y(j) **on input**
14 C
15 C fherm(1,i,j) = derivative df/dx at x(i),y(j) **on output**
16 C fherm(2,i,j) = derivative df/dy at x(i),y(j) **on output**
17 C fherm(3,i,j) = derivative d2f/dxdy at x(i),y(j) **on output**
18 C
19 C addl output:
20 C ilinx=1 if x axis is evenly spaced
21 C iliny=1 if y axis is evenly spaced
22 C ier=0 if no error:
23 C x, y must both be strict ascending
24 C nf2.ge.nx is required.
25 C
26 C a default boundary condition is used, based on divided differences
27 C in the edge region. For more control of BC, use akherm2p...
28 C
29  call akherm2p(x,nx,y,ny,fherm,nf2,ilinx,iliny,0,0,ier)
30 C
31  return
32  end
33 C----------------------------
34  subroutine akherm2p(x,nx,y,ny,fherm,nf2,ilinx,iliny,ipx,ipy,ier)
35 C
36 C create a data set for Hermite interpolation, based on Akima's method
37 C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
38 C
39 C with independently settable boundary condition options:
40 C ipx or ipy =0 -- default, boundary conditions from divided diffs
41 C ipx or ipy =1 -- periodic boundary condition
42 C ipx or ipy =2 -- user supplied df/dx or df/dy
43 C input:
44 C
45  integer nx,ny,nf1 ! array dimensions
46  real x(nx) ! x coordinate array
47  real y(ny) ! y coordinate array
48  real fherm(0:3,nf2,ny) ! data/Hermite array
49 C
50  integer ipx ! =1 if df/dx periodic in x
51  integer ipy ! =1 if df/dy periodic in y
52 C
53 C fherm(0,1:nx,1:ny) supplied by user; this routine computes the
54 C rest of the elements, but note:
55 C
56 C if ipx=2: fherm(1,1,1:ny) & fherm(1,nx,1:ny) = INPUT df/dx BCs
57 C if ipy=2: fherm(2,1:nx,1) & fherm(2,1:nx,ny) = INPUT df/dy BCs
58 C
59 C on output, at all grid points (i,j) covering (1:nx,1:ny):
60 C fherm1(1,i,j) -- df/dx at the grid point
61 C fherm1(2,i,j) -- df/dy at the grid point
62 C fherm1(3,i,j) -- d2f/dxdy at the grid point
63 C
64 C---------------------
65 C local...
66 C
67  real wx(2),wy(2),e(2,2)
68 c
69  real, dimension(:,:), allocatable :: ftmp
70 c
71  real xx(0:nx+1)
72  real yy(0:ny+1)
73 c
74 c---------------------
75 c
76 c error checks
77 c
78  ier=0
79 c
80  call splinck(x,nx,ilinx,1.0e-3,ierx)
81  if(ierx.ne.0) ier=ier+1
82 c
83  if(ierx.ne.0) then
84  write(6,'('' ?akherm2: x axis not strict ascending'')')
85  endif
86 c
87  call splinck(y,ny,iliny,1.0e-3,iery)
88  if(iery.ne.0) ier=ier+1
89 c
90  if(iery.ne.0) then
91  write(6,'('' ?akherm2: y axis not strict ascending'')')
92  endif
93 c
94  if(nf2.lt.nx) then
95  ier=ier+1
96  write(6,*) '?akherm2: fherm array dimension too small.'
97  endif
98 C
99  ierbc=0
100  call ibc_ck(ipx,'akherm2','X Bdy Cond',0,2,ierbc)
101  ier=ier+ierbc
102 C
103  ierbc=0
104  call ibc_ck(ipy,'akherm2','Y Bdy Cond',0,2,ierbc)
105  ier=ier+ierbc
106 C
107  if(ier.ne.0) return
108 C
109 C---------------------------------------
110 C
111 C get a temporary array for f -- will extend out by 1 zone in
112 C each direction so that numerical derivative evaluation can be
113 C done without a lot of special case logic near the edges...
114 C
115  allocate(ftmp(0:nx+1,0:ny+1))
116 C
117  do iy=1,ny
118  do ix=1,nx
119  ftmp(ix,iy)=fherm(0,ix,iy)
120  enddo
121  enddo
122 C
123 C also create expanded axes grids...
124 C
125  xx(1:nx)=x
126  yy(1:ny)=y
127  xx(0)=2*x(1)-x(2)
128  xx(nx+1)=2*x(nx)-x(nx-1)
129  yy(0)=2*y(1)-y(2)
130  yy(ny+1)=2*y(ny)-y(ny-1)
131 C
132 C---------------------------------------
133 C
134 C handle boundary conditions and create rows of extrapolated points
135 C in ftmp. first do ftmp(0,1:ny), ftmp(nx+1,1:ny),
136 C then ftmp(1:nx,0), ftmp(1:nx,ny+1),
137 c then ... fill in the corners
138 c
139 c also, for ipx.le.1 fill in the bdy fherm(1,*) values;
140 c for ipy.le.1 fill in the bdy fherm(2,*) values
141 c
142 c x bc's
143 c
144  do iy=1,ny
145 c
146  cxp=(ftmp(2,iy)-ftmp(1,iy))/(xx(2)-xx(1))
147  cxm=(ftmp(nx,iy)-ftmp(nx-1,iy))/(xx(nx)-xx(nx-1))
148 c
149  if(ipx.eq.1) then
150 c
151 c periodic BC
152 c
153  if(nx.gt.2) then
154  cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
155  cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/(xx(nx-1)-xx(nx-2))
156 c
157  call akherm0(cxmm,cxm,cxp,cxpp,wx,fherm(1,1,iy))
158  fherm(1,nx,iy)=fherm(1,1,iy)
159  else
160  fherm(1,1,iy) = cxp ! =cxm, nx=2
161  fherm(1,nx,iy) = fherm(1,1,iy)
162  endif
163 c
164  cxtrap0=cxm
165  cxtrapn=cxp
166 c
167  else if(ipx.eq.0) then
168 C
169 C default BC -- standard numeric extrapolation
170 C
171  if(nx.gt.2) then
172  cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
173  fherm(1,1,iy)=1.5*cxp-0.5*cxpp
174 C
175  cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/(xx(nx-1)-xx(nx-2))
176  fherm(1,nx,iy)=1.5*cxm-0.5*cxmm
177 c
178  else
179  fherm(1,1,iy) = cxp ! =cxm, nx=2
180  fherm(1,nx,iy) = fherm(1,1,iy)
181  endif
182 C
183 C extrapolate to slope to ghost points just past bdy...
184 C
185  cxtrap0=2.0*fherm(1,1,iy)-cxp
186  cxtrapn=2.0*fherm(1,nx,iy)-cxm
187 C
188  else
189 C
190 C BC supplied by user. Also use this for extrapolation...
191 C
192  cxtrap0=2.0*fherm(1,1,iy)-cxp
193  cxtrapn=2.0*fherm(1,nx,iy)-cxm
194 C
195  endif
196 C
197  ftmp(0,iy)=ftmp(1,iy)-cxtrap0*(xx(1)-xx(0))
198  ftmp(nx+1,iy)=ftmp(nx,iy)+cxtrapn*(xx(nx+1)-xx(nx))
199 C
200  enddo
201 c
202 c y bc's
203 c
204  do ix=1,nx
205 c
206  cyp=(ftmp(ix,2)-ftmp(ix,1))/(yy(2)-yy(1))
207  cym=(ftmp(ix,ny)-ftmp(ix,ny-1))/(yy(ny)-yy(ny-1))
208 c
209  if(ipy.eq.1) then
210 c
211 c periodic BC
212 c
213  if(ny.gt.2) then
214  cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
215  cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/(yy(ny-1)-yy(ny-2))
216 c
217  call akherm0(cymm,cym,cyp,cypp,wy,fherm(2,ix,1))
218  fherm(2,ix,ny)=fherm(2,ix,1)
219 c
220  else
221  fherm(2,ix,1) = cyp ! =cym, ny=2
222  fherm(2,ix,ny)=fherm(2,ix,1)
223  endif
224 c
225  cytrap0=cym
226  cytrapn=cyp
227 c
228  else if(ipy.eq.0) then
229 C
230 C default BC -- standard numeric extrapolation
231 C
232  if(ny.gt.2) then
233  cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
234  fherm(2,ix,1)=1.5*cyp-0.5*cypp
235 C
236  cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/(yy(ny-1)-yy(ny-2))
237  fherm(2,ix,ny)=1.5*cym-0.5*cymm
238 c
239  else
240  fherm(2,ix,1) = cyp ! =cym, ny=2
241  fherm(2,ix,ny)=fherm(2,ix,1)
242  endif
243 C
244 C extrapolate to slope to ghost points just past bdy...
245 C
246  cytrap0=2.0*fherm(2,ix,1)-cyp
247  cytrapn=2.0*fherm(2,ix,ny)-cym
248 C
249  else
250 C
251 C BC supplied by user. Also use this for extrapolation...
252 C
253  cytrap0=2.0*fherm(2,ix,1)-cyp
254  cytrapn=2.0*fherm(2,ix,ny)-cym
255 C
256  endif
257 C
258  ftmp(ix,0)=ftmp(ix,1)-cytrap0*(yy(1)-yy(0))
259  ftmp(ix,ny+1)=ftmp(ix,ny)+cytrapn*(yy(ny+1)-yy(ny))
260 C
261  enddo
262 C
263 C and do something for the corners...
264 C
265  do ix=0,1
266  do iy=0,1
267  icx=ix*(nx+1)
268  icy=iy*(ny+1)
269  incx=1-2*ix
270  incy=1-2*iy
271 
272  ix1=icx+incx
273  iy1=icy+incy
274 
275  ftmp(icx,icy)=ftmp(icx,iy1)+ftmp(ix1,icy)-ftmp(ix1,iy1)
276 C
277  enddo
278  enddo
279 C
280 C----------------------------------------------------------------
281 C OK, now ready to compute all the interior coefficients and the
282 C rest of the edge coefficients as well...
283 C
284  do iy=1,ny
285  iym2=iy
286  iym1=iym2-1
287 c
288  iymm2=iy-1
289  iymm1=iymm2-1
290 c
291  iyp2=iy+1
292  iyp1=iyp2-1
293 c
294  iypp2=iy+2
295  iypp1=iypp2-1
296 c
297  do ix=1,nx
298 c
299 c x div. diffs in vicinity
300 c
301  ixm2=ix
302  ixm1=ixm2-1
303 c
304  ixmm2=ix-1
305  ixmm1=ixmm2-1
306 c
307  iflagx=0
308  cxm=(ftmp(ixm2,iy)-ftmp(ixm1,iy))/(xx(ixm2)-xx(ixm1))
309  if(ix.gt.1) then
310  cxmm=(ftmp(ixmm2,iy)-ftmp(ixmm1,iy))/
311  > (xx(ixmm2)-xx(ixmm1))
312  else
313  if(ipx.eq.1) then
314  cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/
315  > (xx(nx-1)-xx(nx-2))
316  else
317  iflagx=1
318  endif
319  endif
320 c
321  ixp2=ix+1
322  ixp1=ixp2-1
323 c
324  ixpp2=ix+2
325  ixpp1=ixpp2-1
326 c
327  cxp=(ftmp(ixp2,iy)-ftmp(ixp1,iy))/(xx(ixp2)-xx(ixp1))
328  if(ix.lt.nx) then
329  cxpp=(ftmp(ixpp2,iy)-ftmp(ixpp1,iy))/
330  > (xx(ixpp2)-xx(ixpp1))
331  else
332  if(ipx.eq.1) then
333  cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
334  else
335  cxpp=cxp+(cxm-cxmm)
336  endif
337  endif
338 c
339  if(iflagx.eq.1) then
340  cxmm=cxm+(cxp-cxpp)
341  endif
342 c
343 c Akima weightings + df/dx for interior pts
344 c
345  call akherm0(cxmm,cxm,cxp,cxpp,wx,zansr)
346  if((ix.gt.1).and.(ix.lt.nx)) fherm(1,ix,iy)=zansr
347 c
348 c y div. diffs in vicinity
349 c
350  iflagy=0
351  cym=(ftmp(ix,iym2)-ftmp(ix,iym1))/(yy(iym2)-yy(iym1))
352  if(iy.gt.1) then
353  cymm=(ftmp(ix,iymm2)-ftmp(ix,iymm1))/
354  > (yy(iymm2)-yy(iymm1))
355  else
356  if(ipy.eq.1) then
357  cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/
358  > (yy(ny-1)-yy(ny-2))
359  else
360  iflagy=1
361  endif
362  endif
363 c
364  cyp=(ftmp(ix,iyp2)-ftmp(ix,iyp1))/(yy(iyp2)-yy(iyp1))
365  if(iy.lt.ny) then
366  cypp=(ftmp(ix,iypp2)-ftmp(ix,iypp1))/
367  > (yy(iypp2)-yy(iypp1))
368  else
369  if(ipy.eq.1) then
370  cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
371  else
372  cypp=cyp+(cym-cymm)
373  endif
374  endif
375 c
376  if(iflagy.eq.1) then
377  cymm=cym+(cyp-cypp)
378  endif
379 c
380 c Akima weightings + df/dy for interior pts
381 c
382  call akherm0(cymm,cym,cyp,cypp,wy,zansr)
383  if((iy.gt.1).and.(iy.lt.ny)) fherm(2,ix,iy)=zansr
384 c
385 c cross derivatives (2nd order divided differences)
386 c
387  cxm2=(ftmp(ixm2,iym1)-ftmp(ixm1,iym1))/
388  > (xx(ixm2)-xx(ixm1))
389  e(1,1)=(cxm-cxm2)/(yy(iym2)-yy(iym1))
390 c
391  cxm2=(ftmp(ixm2,iyp2)-ftmp(ixm1,iyp2))/
392  > (xx(ixm2)-xx(ixm1))
393  e(1,2)=(cxm2-cxm)/(yy(iyp2)-yy(iyp1))
394 c
395  cxp2=(ftmp(ixp2,iym1)-ftmp(ixp1,iym1))/
396  > (xx(ixp2)-xx(ixp1))
397  e(2,1)=(cxp-cxp2)/(yy(iym2)-yy(iym1))
398 c
399  cxp2=(ftmp(ixp2,iyp2)-ftmp(ixp1,iyp2))/
400  > (xx(ixp2)-xx(ixp1))
401  e(2,2)=(cxp2-cxp)/(yy(iyp2)-yy(iyp1))
402 c
403 c the values
404 c
405  fherm(3,ix,iy)=(wx(1)*(wy(1)*e(1,1)+wy(2)*e(1,2))+
406  > wx(2)*(wy(1)*e(2,1)+wy(2)*e(2,2)))/
407  > ((wx(1)+wx(2))*(wy(1)+wy(2)))
408 c
409  enddo
410  enddo
411 C
412  deallocate (ftmp)
413  return
414  end