V3FIT
genxpkg.f
1  subroutine genxpkg(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  integer nx ! length of x, .ge.4
56  real x(nx) ! x axis vector, strict ascending
57 c
58  integer iper ! =1 if x is periodic
59  integer imsg ! =1 for range error messages
60  integer itol ! =1 to specify tolerance, else default
61 c
62 c default tolerance is 5.0e-7
63 c
64  real ztol ! range tolerance, if itol=1
65 c
66 c lookup algorithm selection for uneven grids:
67 c
68  integer ialg ! = +/- 1: <1/h(j)> "Newton" method
69 c ! = +/- 2: binary search
70 c ! = +/- 3: indexing function
71 c
72 c to use past search result as init. cond. for next search, give
73 c ialg .lt. 0; otherwise give ialg .gt. 0.
74 c
75 c output arguments:
76 c
77  real xpkg(nx,4) ! xpkg, as described above
78  integer ier ! completion code, 0=OK
79 c
80 c------------------------------------------------
81 c
82  if(nx.lt.2) then
83  write(6,*) .ge.' %genxpkg: nx2 required!'
84  ier=1
85  return
86  else
87  ier=0
88  endif
89 c
90  ialgu=ialg
91  if(ialgu.eq.0) ialgu=3
92  if(iabs(ialgu).gt.3) ialgu=3
93 c
94 c get tolerance for end point range check & even spacing check
95 c
96  if(itol.eq.1) then
97  ztolr=abs(ztol)
98  else
99  ztolr=5.0e-7
100  endif
101 c
102  ztola=max(abs(x(1)),abs(x(nx)))*ztolr
103 c
104 c assume even spacing for now...
105 c
106  if(nx.ge.3) then
107  xpkg(3,4)=0.0
108  endif
109 c
110 c mark if x axis is a periodic coordinate
111 c
112  if(iper.eq.1) then
113  xpkg(2,4)=1.0
114  else
115  xpkg(2,4)=0.0
116  endif
117 c
118 c store tolerance parameter
119 c
120  xpkg(1,4)=ztola
121 c
122 c mark if messages are to be written if range lookup errors occur
123 c
124  if(imsg.eq.1) then
125  continue ! xpkg(1,4) left .gt. 0.0
126  else
127  xpkg(1,4)=-xpkg(1,4)
128  endif
129 c
130 c OK check linearity and spacing
131 c
132  ier=0
133 c
134  xpkg(nx,2)=(x(nx)-x(1))/(nx-1) ! average spacing
135 c
136  do ix=1,nx
137  xpkg(ix,1)=x(ix)
138  if((ier.eq.0).and.(ix.lt.nx)) then
139  if(x(ix+1).le.x(ix)) then
140  ier=1
141  write(6,*) ' %genxpkg: x axis not strict ascending!'
142  else
143  zh=x(ix+1)-x(ix)
144 c
145 c checking even spacing now...
146 c
147  if(nx.ge.3) then
148  if(abs(zh-xpkg(nx,2)).gt.ztola) xpkg(3,4)=1.0
149  endif
150 c
151  xpkg(ix,2)=zh
152  xpkg(ix,3)=1.0/zh
153 c
154  endif
155  endif
156  enddo
157 c
158  if(ier.ne.0) return
159 c
160 c OK, store inverse average spacing
161 c
162  xpkg(nx,3)=1.0/xpkg(nx,2)
163 c
164 c if even spacing is detected, redefine x axis slightly, for
165 c improved regularity of behaviour
166 c
167  if(nx.ge.3) then
168 c
169  if(xpkg(3,4).eq.0.0) then
170  do ix=1,nx-2
171  ixp=ix+1
172  xpkg(ixp,1)=xpkg(1,1)+ix*xpkg(nx,2)
173  enddo
174  if(nx.gt.3) then
175  if(ialgu.lt.0) then
176  xpkg(4,4)=1.0 ! check init guess
177  else
178  xpkg(4,4)=0.0
179  endif
180  endif
181  endif
182 c
183 c if uneven spacing is detected, must use an algorithm...
184 c
185  if(xpkg(3,4).ne.0.0) then
186  xpkg(3,4)=abs(ialgu)
187  if(nx.gt.3) then
188  if(ialgu.lt.0) then
189  xpkg(4,4)=1.0 ! check init guess
190  else
191  xpkg(4,4)=0.0
192  endif
193  endif
194 c
195  if(abs(ialgu).eq.3) then
196 c
197 c construct a piecewise linear indexing function
198 c
199  xpkg(1,2)=1.0
200  xtest=xpkg(1,1)
201  itest=1
202  do i=2,nx-1
203  xtest=xtest+xpkg(nx,2) ! x1 + (i-1)*<h>
204  10 continue
205  if((xpkg(itest,1).le.xtest).and.
206  > (xtest.le.xpkg(itest+1,1))) then
207  xpkg(i,2)=itest+(xtest-xpkg(itest,1))/
208  > (xpkg(itest+1,1)-xpkg(itest,1))
209  else
210  itest=itest+1
211  go to 10
212  endif
213  enddo
214 c
215 c (implicitly, xpkg(nx,2)=nx*1.0; but we leave <h> in xpkg(nx,2)
216 c instead...)
217 c
218  endif
219  endif
220 c
221  endif ! nx.ge.3
222 c
223  return
224  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17