V3FIT
r8akherm1.f
1  subroutine r8akherm1(x,nx,fherm,ilinx,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 1d routine -- default boundary condition (based on divided differences)
7 C
8 C input:
9 C
10 !============
11 ! idecl: explicitize implicit INTEGER declarations:
12  IMPLICIT NONE
13  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
14  INTEGER ier,ilinx
15 !============
16  integer nx ! array dimensions
17  real*8 x(nx) ! x coordinate array
18  real*8 fherm(0:1,nx) ! data/Hermite array
19 C
20 C fherm(0,i) = function value f at x(i) **on input**
21 C
22 C fherm(1,i) = derivative df/dx at x(i) **on output**
23 C
24 C addl output:
25 C ilinx=1 if x is "evenly spaced" ier=0 if no errors
26 C
27 C ** x must be strict ascending **
28 C
29 C work by calling akherm1p; no periodic boundary condition
30 C
31  call r8akherm1p(x,nx,fherm,ilinx,0,ier)
32 C
33  return
34  end
35 C----------------------------
36  subroutine r8akherm1p(x,nx,fherm,ilinx,ipx,ier)
37 C
38 C Akima subroutine with boundary condition option:
39 C
40 C =>ipx=1 for periodic boundary condition
41 C
42 C =>ipx=2 for user specified boundary condition:
43 C in which case fherm(1,1) and fherm(1,nx) must contain
44 C the derivatives fx(x(1)) and fx(x(nx)), respectively, and these
45 C will not be modified on output.
46 C
47 C =>ipx=0: default boundary conditions
48 C
49 C other arguments as with akherm1.
50 C
51 C input:
52 C
53 !============
54 ! idecl: explicitize implicit INTEGER declarations:
55  IMPLICIT NONE
56  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
57  INTEGER ier,ilinx,ierbc,ix,ixm2,ixm1,ixmm2,ixmm1,ixp2,ixp1
58  INTEGER ixpp2,ixpp1
59 !============
60 ! idecl: explicitize implicit REAL declarations:
61  real*8 ztol,cxp,cxm,cxpp,cxmm,cxtrap0,cxtrap1
62 !============
63  integer nx ! array dimensions
64  real*8 x(nx) ! x coordinate array
65  real*8 fherm(0:1,nx) ! data/Hermite array
66  integer ipx ! =1: f' periodic in x
67 C
68 C----------------------------
69 C
70  real*8 wx(2)
71 C
72 C error checks...
73 C
74  ztol=1.0e-3_r8
75  ier=0
76 C
77  call r8splinck(x,nx,ilinx,ztol,ier)
78  if(ier.ne.0) then
79  write(6,*) '?akherm1: x axis not strict ascending.'
80  ier=ier+1
81  endif
82 C
83  ierbc=0
84  call ibc_ck(ipx,'akherm1','Bdy Cond',0,2,ierbc)
85  ier=ier+ierbc
86  if(ier.gt.0) return
87 C
88 C deal with boundary region. The boundary derivatives are set
89 C and "ghost" points are provided...
90 C
91 C all paths need the 1st div. diffs at the bdy...
92 C
93  cxp=(fherm(0,2)-fherm(0,1))/(x(2)-x(1))
94  cxm=(fherm(0,nx)-fherm(0,nx-1))/(x(nx)-x(nx-1))
95 C
96  if(ipx.eq.1) then
97 C
98 C periodic BC
99 C
100 C LHS/RHS
101 C
102  if(nx.gt.2) then
103  cxpp=(fherm(0,3)-fherm(0,2))/(x(3)-x(2))
104  cxmm=(fherm(0,nx-1)-fherm(0,nx-2))/(x(nx-1)-x(nx-2))
105 C
106  call r8akherm0(cxmm,cxm,cxp,cxpp,wx,fherm(1,1))
107  fherm(1,nx)=fherm(1,1)
108  else
109  ! nx=2
110  fherm(1,1)=cxp ! =cxm
111  fherm(1,nx)=cxm ! =cxp
112  endif
113 c
114  cxtrap0=cxm
115  cxtrap1=cxp
116 C
117  else if(ipx.eq.0) then
118 C
119 C default BC -- standard numeric extrapolation
120 C
121  if(nx.gt.2) then
122  cxpp=(fherm(0,3)-fherm(0,2))/(x(3)-x(2))
123  fherm(1,1)=1.5_r8*cxp-0.5_r8*cxpp
124 C
125  cxmm=(fherm(0,nx-1)-fherm(0,nx-2))/(x(nx-1)-x(nx-2))
126  fherm(1,nx)=1.5_r8*cxm-0.5_r8*cxmm
127 C
128  else
129  ! nx=2
130  fherm(1,1)=cxp ! =cxm
131  fherm(1,nx)=cxm ! =cxp
132  endif
133 C
134 C extrapolate to slope to ghost points just past bdy...
135 C
136  cxtrap0=2.0_r8*fherm(1,1)-cxp
137  cxtrap1=2.0_r8*fherm(1,nx)-cxm
138 C
139  else
140 C
141 C BC supplied by user. Also use this for extrapolation...
142 C extrapolate to slope to ghost points just past bdy...
143 C
144  cxtrap0=2.0_r8*fherm(1,1)-cxp
145  cxtrap1=2.0_r8*fherm(1,nx)-cxm
146 C
147  endif
148 C
149 C NOTE: this loop is inactive if nx=2
150 C
151  do ix=2,nx-1
152 c
153 c x div. diffs in vicinity
154 c
155  ixm2=ix
156  ixm1=ixm2-1
157  ixmm2=ixm1
158  ixmm1=ixm1-1
159 c
160  ixp2=ix+1
161  ixp1=ix
162  ixpp2=ix+2
163  ixpp1=ixp2
164 c
165  if(ix.eq.2) then
166  cxmm=cxtrap0
167  else
168  cxmm=(fherm(0,ixmm2)-fherm(0,ixmm1))/(x(ixmm2)-x(ixmm1))
169  endif
170 c
171  if(ix.eq.nx-1) then
172  cxpp=cxtrap1
173  else
174  cxpp=(fherm(0,ixpp2)-fherm(0,ixpp1))/(x(ixpp2)-x(ixpp1))
175  endif
176 c
177  cxm=(fherm(0,ixm2)-fherm(0,ixm1))/(x(ixm2)-x(ixm1))
178  cxp=(fherm(0,ixp2)-fherm(0,ixp1))/(x(ixp2)-x(ixp1))
179 C
180  call r8akherm0(cxmm,cxm,cxp,cxpp,wx,fherm(1,ix))
181 c
182  enddo
183 C
184  return
185  end
186 C--------------------------------------
187  subroutine r8akherm0(cxmm,cxm,cxp,cxpp,wx,result)
188 c
189 c basic akima formula for 1st derivative at pt p:
190 c
191 c cxmm = numerical slope 2 zones to left
192 c cxm = numerical slope 1 zone to left
193 c cxp = numerical slope 1 zone to right
194 c cxpp = numerical slope 2 zones to right
195 c
196 c return slope at point p and weighting factors
197 c
198  IMPLICIT NONE
199  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
200  real*8 cxmm,cxm,cxp,cxpp ! nearby slopes (in)
201  real*8 wx(2) ! weights (out)
202  real*8 result ! Akima nominal slope (out)
203 c
204  wx(2)=abs(cxm-cxmm)
205  wx(1)=abs(cxp-cxpp)
206  if(wx(1)+wx(2).eq.0.0_r8) then
207  wx(1)=1.0_r8
208  wx(2)=1.0_r8
209  endif
210 c
211 c the derivative -- weighted average of neighbouring slopes
212 c
213  result=(wx(1)*cxm+wx(2)*cxp)/(wx(1)+wx(2))
214 c
215  return
216  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17