V3FIT
r8dnherm2.f
1  subroutine r8dnherm2(x,nx,y,ny,fherm,nf2,ilinx,iliny,ier)
2 C
3 C create a data set for Hermite interpolation, based on simple
4 C numerical differentiation using the given grid. The grid should
5 C be "approximately" evenly spaced.
6 C
7 C 2d routine
8 C
9 C input:
10 C
11 C nf2.gt.nx expected!
12 C
13 !============
14 ! idecl: explicitize implicit INTEGER declarations:
15  IMPLICIT NONE
16  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
17  INTEGER iliny,ier,ilinx,ierx,iery,iy,iyp,iym,ix,ixp,ixm
18 !============
19 ! idecl: explicitize implicit REAL declarations:
20  real*8 zd
21 !============
22  integer nx,ny,nf2 ! array dimensions
23  real*8 x(nx) ! x coordinate array
24  real*8 y(ny) ! y coordinate array
25  real*8 fherm(0:3,nf2,ny) ! data/Hermite array
26 C
27 C fherm(0,i,j) = function value f at x(i),y(j) **on input**
28 C
29 C fherm(1,i,j) = derivative df/dx at x(i),y(j) **on output**
30 C fherm(2,i,j) = derivative df/dy at x(i),y(j) **on output**
31 C fherm(3,i,j) = derivative d2f/dxdy at x(i),y(j) **on output**
32 C
33 C addl output:
34 C ilinx=1 if x axis is evenly spaced
35 C iliny=1 if y axis is evenly spaced
36 C ier=0 if no error:
37 C x, y must both be strict ascending
38 C nf2.ge.nx is required.
39 C
40 C----------------------------
41 c
42  ier=0
43 c
44  call r8splinck(x,nx,ilinx,1.0e-3_r8,ierx)
45  if(ierx.ne.0) ier=2
46 c
47  if(ier.eq.2) then
48  write(6,'('' ?dnherm2: x axis not strict ascending'')')
49  endif
50 c
51  call r8splinck(y,ny,iliny,1.0e-3_r8,iery)
52  if(iery.ne.0) ier=3
53 c
54  if(ier.eq.3) then
55  write(6,'('' ?dnherm2: y axis not strict ascending'')')
56  endif
57 c
58  if(nf2.lt.nx) then
59  ier=4
60  write(6,*) '?dnherm2: fherm array dimension too small.'
61  endif
62 C
63  if(ier.ne.0) return
64 C
65  do iy=1,ny
66 c
67  iyp=min(ny,iy+1)
68  iym=max(1,iy-1)
69 c
70  do ix=1,nx
71 c
72  ixp=min(nx,ix+1)
73  ixm=max(1,ix-1)
74 c
75 c x diffs in vicinity
76 c
77  zd=(fherm(0,ixp,iy)-fherm(0,ixm,iy))/(x(ixp)-x(ixm))
78 c
79  fherm(1,ix,iy)=zd
80 c
81 c y diffs in vicinity
82 c
83  zd=(fherm(0,ix,iyp)-fherm(0,ix,iym))/(y(iyp)-y(iym))
84 c
85  fherm(2,ix,iy)=zd
86 c
87 c xy cross derivs (except at corners, iedge=2)
88 c
89  fherm(3,ix,iy)=(fherm(0,ixp,iyp)-fherm(0,ixm,iyp)
90  > -fherm(0,ixp,iym)+fherm(0,ixm,iym))/
91  > ((x(ixp)-x(ixm))*(y(iyp)-y(iym)))
92 c
93  enddo
94  enddo
95 C
96  return
97  end