V3FIT
r8dnherm3.f
1  subroutine r8dnherm3(x,nx,y,ny,z,nz,fherm,nf2,nf3,
2  > ilinx,iliny,ilinz,ier)
3 C
4 C create a data set for Hermite interpolation, based on simple
5 C numerical differentiation using the given grid. The grid should
6 C be "approximately" evenly spaced.
7 C
8 C 3d routine
9 C
10 C input:
11 C
12 C nf2.get.nx and nf3.ge.ny expected!
13 C
14 !============
15 ! idecl: explicitize implicit INTEGER declarations:
16  IMPLICIT NONE
17  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
18  INTEGER iliny,ilinz,ier,ilinx,ierx,iery,ierz,iz,izp,izm,iy
19  INTEGER iyp,iym,ix,ixp,ixm
20 !============
21 ! idecl: explicitize implicit REAL declarations:
22  real*8 zd
23 !============
24  integer nx,ny,nz,nf2,nf3 ! array dimensions
25  real*8 x(nx) ! x coordinate array
26  real*8 y(ny) ! y coordinate array
27  real*8 z(nz) ! z coordinate array
28  real*8 fherm(0:7,nf2,nf3,nz) ! data/Hermite array
29 C
30 C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input**
31 C
32 C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output**
33 C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output**
34 C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output**
35 C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output**
36 C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output**
37 C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output**
38 C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output**
39 C
40 C additional outputs:
41 C
42 C ilinx = 1 if x is evenly spaced (approximately)
43 C iliny = 1 if y is evenly spaced (approximately)
44 C ilinz = 1 if z is evenly spaced (approximately)
45 C
46 C ier = 0 if there are no errors
47 C
48 C note possible errors: x y or z NOT strict ascending
49 C nf2.lt.nx .or. nf3.lt.ny
50 C
51 C----------------------------
52 C
53 C
54 C error check
55 C
56  ier=0
57 c
58  call r8splinck(x,nx,ilinx,1.0e-3_r8,ierx)
59  if(ierx.ne.0) ier=2
60 c
61  if(ier.eq.2) then
62  write(6,'('' ?dnherm3: x axis not strict ascending'')')
63  endif
64 c
65  call r8splinck(y,ny,iliny,1.0e-3_r8,iery)
66  if(iery.ne.0) ier=3
67 c
68  if(ier.eq.3) then
69  write(6,'('' ?dnherm3: y axis not strict ascending'')')
70  endif
71 c
72  call r8splinck(z,nz,ilinz,1.0e-3_r8,ierz)
73  if(ierz.ne.0) ier=4
74 c
75  if(ier.eq.4) then
76  write(6,'('' ?dnherm3: z axis not strict ascending'')')
77  endif
78 c
79  if(nf2.lt.nx) then
80  ier=5
81  write(6,*) '?dnherm3: fherm (x) array dimension too small.'
82  endif
83 c
84  if(nf3.lt.ny) then
85  ier=6
86  write(6,*) '?dnherm3: fherm (y) array dimension too small.'
87  endif
88 C
89  if(ier.ne.0) return
90 C
91 C error check OK
92 C
93  do iz=1,nz
94 c
95  izp=min(nz,iz+1)
96  izm=max(1,iz-1)
97 c
98  do iy=1,ny
99 c
100  iyp=min(ny,iy+1)
101  iym=max(1,iy-1)
102 c
103  do ix=1,nx
104 c
105  ixp=min(nx,ix+1)
106  ixm=max(1,ix-1)
107 c
108 c x diffs in vicinity
109 c
110  zd=(fherm(0,ixp,iy,iz)-fherm(0,ixm,iy,iz))/
111  > (x(ixp)-x(ixm))
112 c
113  fherm(1,ix,iy,iz)=zd
114 c
115 c y diffs in vicinity
116 c
117  zd=(fherm(0,ix,iyp,iz)-fherm(0,ix,iym,iz))/
118  > (y(iyp)-y(iym))
119 c
120  fherm(2,ix,iy,iz)=zd
121 c
122 c z diffs in vicinity
123 c
124  zd=(fherm(0,ix,iy,izp)-fherm(0,ix,iy,izm))/
125  > (z(izp)-z(izm))
126 c
127  fherm(3,ix,iy,iz)=zd
128 c
129 c xy cross derivs
130 c
131  fherm(4,ix,iy,iz)=
132  > (fherm(0,ixp,iyp,iz)-fherm(0,ixm,iyp,iz)
133  > -fherm(0,ixp,iym,iz)+fherm(0,ixm,iym,iz))/
134  > ((x(ixp)-x(ixm))*(y(iyp)-y(iym)))
135 c
136 c xz cross derivs
137 c
138  fherm(5,ix,iy,iz)=
139  > (fherm(0,ixp,iy,izp)-fherm(0,ixm,iy,izp)
140  > -fherm(0,ixp,iy,izm)+fherm(0,ixm,iy,izm))/
141  > ((x(ixp)-x(ixm))*(z(izp)-z(izm)))
142 c
143 c yz cross derivs
144 c
145  fherm(6,ix,iy,iz)=
146  > (fherm(0,ix,iyp,izp)-fherm(0,ix,iym,izp)
147  > -fherm(0,ix,iyp,izm)+fherm(0,ix,iym,izm))/
148  > ((y(iyp)-y(iym))*(z(izp)-z(izm)))
149 c
150 c xyz cross deriv
151 c
152  fherm(7,ix,iy,iz)=
153  > ((fherm(0,ixp,iyp,izp)-fherm(0,ixp,iym,izp)
154  > -fherm(0,ixp,iyp,izm)+fherm(0,ixp,iym,izm))-
155  > (fherm(0,ixm,iyp,izp)-fherm(0,ixm,iym,izp)
156  > -fherm(0,ixm,iyp,izm)+fherm(0,ixm,iym,izm)))/
157  > ((x(ixp)-x(ixm))*(y(iyp)-y(iym))*(z(izp)-z(izm)))
158 c
159  enddo
160  enddo
161  enddo
162 C
163  return
164  end