V3FIT
r8genxpkg.f
1  subroutine r8genxpkg(nx,x,xpkg,iper,imsg,itol,ztol,ialg,ier)
2 c
3 c from an x axis assemble a "package":
4 c
5 c MOD DMC Feb 2010: handle very small grids: nx=2, nx=3...
6 c interchanged meaning of xpkg(1,4) and xpkg(3,4);
7 c xpkg(3,4) only set if nx.ge.3;
8 c xpkg(4,4) only set if nx.ge.4.
9 c
10 c there are corresponding changes in xlookup: simplified code lookup
11 c code for the cases nx=2 and nx=3.
12 c
13 c xpkg(j,1) = x(j), j = 1 to nx ... nx.ge.2
14 c
15 c if(abs(ialg).ne.3) then...
16 c for j=1:nx-1
17 c xpkg(j,2) = h(j) = x(j+1)-x(j), j=1 to nx-1
18 c else
19 c for j=1:nx-1
20 c xpkg(j,2) = index location, with linear offset, of
21 c xpkg(1,1)+<h>*(j-1) in the original x(1:nx)
22 c (piecewise linear indexing function)
23 c endif
24 c xpkg(nx,2) = <h> = (x(nx)-x(1))/(nx-1)
25 c
26 c xpkg(j,3) = 1/h(j)
27 c xpkg(nx,3) = 1/<h>
28 c
29 c xpkg(1,4) = +/- tolerance epsilon for out-of-range warnings/errors
30 c +if message is to be written on out of range errors
31 c -if no message to be written. In either case, a
32 c warning flag is set.
33 c
34 c xpkg(2,4) = 1.0 if x is *periodic*, else 0.0
35 c
36 c only set if nx.ge.3:
37 c xpkg(3,4) = 0.0 if x is *evenly spaced* (to ztol or 5.e-7 rel) else:
38 c = 1.0 => use (1/h) newton's method like search algorithm
39 c = 2.0 => use binary search algorithm
40 c = 3.0 => use piecewise-linear indexing function...
41 c
42 c only set if nx.ge.4:
43 c xpkg(4,4) = 0.0 -- do not use past search result as start point
44 c for next search;
45 c = 1.0 -- do use past search result as start point for
46 c next search.
47 c
48 c tolerance epsilon means:
49 c if xtest is within epsilon*max(abs(x(1)),abs(x(nx))) outside the
50 c range [x(1),x(nx)] treat xtest as being at the endpoint x(1) or
51 c x(nx) whichever is closer.
52 c
53 c input arguments:
54 c
55 !============
56 ! idecl: explicitize implicit INTEGER declarations:
57  IMPLICIT NONE
58  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
59  INTEGER ialgu,iabs,ix,ixp,itest,i
60 !============
61 ! idecl: explicitize implicit REAL declarations:
62  real*8 ztolr,ztola,zh,xtest
63 !============
64  integer nx ! length of x, .ge.4
65  real*8 x(nx) ! x axis vector, strict ascending
66 c
67  integer iper ! =1 if x is periodic
68  integer imsg ! =1 for range error messages
69  integer itol ! =1 to specify tolerance, else default
70 c
71 c default tolerance is 5.0e-7
72 c
73  real*8 ztol ! range tolerance, if itol=1
74 c
75 c lookup algorithm selection for uneven grids:
76 c
77  integer ialg ! = +/- 1: <1/h(j)> "Newton" method
78 c ! = +/- 2: binary search
79 c ! = +/- 3: indexing function
80 c
81 c to use past search result as init. cond. for next search, give
82 c ialg .lt. 0; otherwise give ialg .gt. 0.
83 c
84 c output arguments:
85 c
86  real*8 xpkg(nx,4) ! xpkg, as described above
87  integer ier ! completion code, 0=OK
88 c
89 c------------------------------------------------
90 c
91  if(nx.lt.2) then
92  write(6,*) .ge.' %genxpkg: nx2 required!'
93  ier=1
94  return
95  else
96  ier=0
97  endif
98 c
99  ialgu=ialg
100  if(ialgu.eq.0) ialgu=3
101  if(iabs(ialgu).gt.3) ialgu=3
102 c
103 c get tolerance for end point range check & even spacing check
104 c
105  if(itol.eq.1) then
106  ztolr=abs(ztol)
107  else
108  ztolr=5.0e-7_r8
109  endif
110 c
111  ztola=max(abs(x(1)),abs(x(nx)))*ztolr
112 c
113 c assume even spacing for now...
114 c
115  if(nx.ge.3) then
116  xpkg(3,4)=0.0_r8
117  endif
118 c
119 c mark if x axis is a periodic coordinate
120 c
121  if(iper.eq.1) then
122  xpkg(2,4)=1.0_r8
123  else
124  xpkg(2,4)=0.0_r8
125  endif
126 c
127 c store tolerance parameter
128 c
129  xpkg(1,4)=ztola
130 c
131 c mark if messages are to be written if range lookup errors occur
132 c
133  if(imsg.eq.1) then
134  continue ! xpkg(1,4) left .gt. 0.0
135  else
136  xpkg(1,4)=-xpkg(1,4)
137  endif
138 c
139 c OK check linearity and spacing
140 c
141  ier=0
142 c
143  xpkg(nx,2)=(x(nx)-x(1))/(nx-1) ! average spacing
144 c
145  do ix=1,nx
146  xpkg(ix,1)=x(ix)
147  if((ier.eq.0).and.(ix.lt.nx)) then
148  if(x(ix+1).le.x(ix)) then
149  ier=1
150  write(6,*) ' %genxpkg: x axis not strict ascending!'
151  else
152  zh=x(ix+1)-x(ix)
153 c
154 c checking even spacing now...
155 c
156  if(nx.ge.3) then
157  if(abs(zh-xpkg(nx,2)).gt.ztola) xpkg(3,4)=1.0_r8
158  endif
159 c
160  xpkg(ix,2)=zh
161  xpkg(ix,3)=1.0_r8/zh
162 c
163  endif
164  endif
165  enddo
166 c
167  if(ier.ne.0) return
168 c
169 c OK, store inverse average spacing
170 c
171  xpkg(nx,3)=1.0_r8/xpkg(nx,2)
172 c
173 c if even spacing is detected, redefine x axis slightly, for
174 c improved regularity of behaviour
175 c
176  if(nx.ge.3) then
177 c
178  if(xpkg(3,4).eq.0.0_r8) then
179  do ix=1,nx-2
180  ixp=ix+1
181  xpkg(ixp,1)=xpkg(1,1)+ix*xpkg(nx,2)
182  enddo
183  if(nx.gt.3) then
184  if(ialgu.lt.0) then
185  xpkg(4,4)=1.0_r8 ! check init guess
186  else
187  xpkg(4,4)=0.0_r8
188  endif
189  endif
190  endif
191 c
192 c if uneven spacing is detected, must use an algorithm...
193 c
194  if(xpkg(3,4).ne.0.0_r8) then
195  xpkg(3,4)=abs(ialgu)
196  if(nx.gt.3) then
197  if(ialgu.lt.0) then
198  xpkg(4,4)=1.0_r8 ! check init guess
199  else
200  xpkg(4,4)=0.0_r8
201  endif
202  endif
203 c
204  if(abs(ialgu).eq.3) then
205 c
206 c construct a piecewise linear indexing function
207 c
208  xpkg(1,2)=1.0_r8
209  xtest=xpkg(1,1)
210  itest=1
211  do i=2,nx-1
212  xtest=xtest+xpkg(nx,2) ! x1 + (i-1)*<h>
213  10 continue
214  if((xpkg(itest,1).le.xtest).and.
215  > (xtest.le.xpkg(itest+1,1))) then
216  xpkg(i,2)=itest+(xtest-xpkg(itest,1))/
217  > (xpkg(itest+1,1)-xpkg(itest,1))
218  else
219  itest=itest+1
220  go to 10
221  endif
222  enddo
223 c
224 c (implicitly, xpkg(nx,2)=nx*1.0; but we leave <h> in xpkg(nx,2)
225 c instead...)
226 c
227  endif
228  endif
229 c
230  endif ! nx.ge.3
231 c
232  return
233  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17