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