1      SUBROUTINE TBFUNX(X,Y,DYDX,NP,XA,YA,C,IN,
2     1           MI,NG,LEXL,LEXU,MESS,NMSS,ROUT)
3C
4C**
5C** INPUTS
6C**
7C**    X    = X COORDINATE OF POINT OF INTEREST
8C**    NP   = NUMBER OF INPUT DATA POINTS
9C**    XA   = ARRAY OF X COORDINATES FOR DATA POINTS
10C**    YA   = ARRAY OF Y COORDINATES FOR DATA POINTS
11C**    C    = DUMMY ARRAY    -- NOT USED
12C**    IN   = DUMMY VARIABLE -- NOT USED
13C**    MI   = DUMMY VARIABLE -- NOT USED
14C**    NG   = DUMMY VARIABLE -- NOT USED
15C**    LEXL = INDICATOR IF BELOW LOWER LIMIT OF DATA
16C**           IF .LE. 0, Y IS YA(1) AND DYDX DETERMINED
17C**           BY 2ND-ORDER CURVE FIT
18C**           IF .EQ. 1, Y IS EXTRAPOLATED LINEARLY USING
19C**           SLOPE OF FIRST 2 DATA POINTS
20C**           IF .GT. 1, Y AND DYDX ARE DETERMINED BY
21C**           2ND-ORDER CURVE FIT
22C**    LEXU = INDICATOR IF ABOVE UPPER LIMIT OF DATA
23C**           IF .LE. 0, Y IS YA(1) AND DYDX DETERMINED
24C**           BY 2ND-ORDER CURVE FIT
25C**           IF .EQ. 1, Y IS EXTRAPOLATED LINEARLY USING
26C**           SLOPE OF LAST 2 DATA POINTS
27C**           IF .GT. 1, Y AND DYDX ARE DETERMINED BY
28C**           2ND-ORDER CURVE FIT
29C**    MESS = FIGURE NUMBER IN EXTRAPOLATION MESSAGE
30C**    NMSS = DIMENSION OF MESS
31C**    ROUT = ROUTINE NAME CALLING TBFUNX
32C**
33C** OUTPUTS
34C**
35C**    Y    = Y COORDINATE AT POINT OF INTEREST
36C**    DYDX = FIRST DERIVATIVE AT POINT OF INTEREST
37C**       *          *          *          *          *
38      COMMON /CONSNT/ PI,DEG,UNUSED,RAD
39      DIMENSION XA(NP),YA(NP),C(6),MESS(2),ROUT(2),MSSCL(9)
40      DIMENSION XX(3),YY(3)
41      EQUIVALENCE (QSSCL4,MSSCL(4)),(QSSCL6,MSSCL(6))
42      DATA MSSCL/4HTBFU, 4HNX  ,2*0,1,4*0/
43C**
44C**    SET FLAGS FOR USING SUBROUTINE QUAD...
45C**
46      DATA ADER/4HDERI/,BDER/4HYNOD/
47
48C AJT INITIALISE L
49      L=1
50
51      NP1=NP-1
52      IF(NP.LT.3)GO TO 1070
53
54C**
55C**    CHECK IF EXTRAPOLATING OR INTERPOLATING
56C**
57C AJT      IF(X.GE.XA(NP).OR.X.LE.XA(1))GO TO 1020
58C AJT FIND MAX AN MIN VALUES OF XA()
59      XAMAX=XA(1)
60      XAMIN=XA(1)
61      DO 500 I=1,NP
62       IF(XA(I).GT. XAMAX) XAMAX=XA(I)
63       IF(XA(I).LT. XAMIN) XAMIN=XA(I)
64 500  CONTINUE
65C AJT      IF(X.GE.XA(NP).OR.X.LE.XA(1))GO TO 1020
66C AJT CHECK IF X IS BETWEEN XAMAX AND XAMIN
67      IF(X.GE.XAMAX.OR.X.LE.XAMIN)GO TO 1020
68C AJT  END OF PATCH
69
70C**
71C**    HERE IF INTERPOLATING...
72C**
73C**    DETERMINE POINTS XA(I) SURROUNDING X...
74C**
75      DO 1000 I=1,NP1
76         IF(X.GE.XA(I))L=I
77 1000 CONTINUE
78      IF(L.EQ.1)L=2
79      L2=L-2
80C**
81C**    FILL VECTORS FOR SUBROUTINE QUAD...
82C**
83      DO 1010 I=1,3
84         XX(I)=XA(L2+I)
85 1010 YY(I)=YA(L2+I)
86C**
87C**    CALCULATE Y BY LINEAR INTERPOLATION...
88C**
89      DEL=XX(3)-XX(2)
90      XDEL=XX(3)
91      IF(DEL.EQ.0.)DEL=UNUSED
92      Y=YY(2)+(YY(3)-YY(2))*(X-XX(2))/DEL
93      DEL=XX(2)-XX(1)
94      XDEL=XX(2)
95      IF(DEL.EQ.0.)DEL=UNUSED
96      IF(X.LT.XA(2))Y=YY(1)+(YY(2)-YY(1))*(X-XX(1))/DEL
97C**
98C**    SET FLAG FOR CALCULATING DERIVATIVE AND
99C**    CALL QUAD TO CALCULATE DYDX...
100C**
101      YIND=ADER
102      CALL QUAD(XX(1),YY(1),X,YIND)
103      DYDX=YIND
104      GO TO 1050
105C**
106C**    HERE IF EXTRAPOLATING...
107C**
108 1020 CONTINUE
109      NP3=NP-3
110C**
111C**    DETERMINE IF ABOVE OR BELOW DATA LIMITS
112C**
113      LE=0
114      LLE=1
115      LIND=LEXL
116      IF(X.GE.XA(NP))LE=NP3
117      IF(X.GE.XA(NP))LLE=NP
118      IF(X.GE.XA(NP))LIND=LEXU
119C**
120C**    FILL VECTORS FOR SUBROUTINE QUAD...
121C**
122      DO 1030 I=1,3
123         XX(I)=XA(LE+I)
124 1030 YY(I)=YA(LE+I)
125C**
126C**    DETERMINE WHICH METHODS FOR CALCULATING Y AND DYDX
127C**
128C**    IF USING LINEAR EXTRAPOLATION, CALCULATE
129C**    Y AND DYDX...
130C**
131      DEL=XA(2)-XA(1)
132      XDEL=XA(2)
133      IF(DEL.EQ.0.)DEL=UNUSED
134      IF(LEXL.EQ.1)DYDX=(YA(2)-YA(1))/DEL
135      DEL=XA(NP)-XA(NP1)
136      XDEL=XA(NP)
137      IF(DEL.EQ.0.)DEL=UNUSED
138      IF(LEXU.EQ.1.AND.LE.EQ.NP3)DYDX=(YA(NP)-YA(NP1))/DEL
139      IF(LEXL.EQ.1)Y=(X-XA(1))*DYDX+YA(1)
140      IF(LEXU.EQ.1.AND.LE.EQ.NP3)Y=(X-XA(NP))*DYDX+YA(NP)
141      IF(LIND.EQ.1)GO TO 1050
142C**
143C**    HERE FOR RETURNING VALUES AT ENDPOINTS FOR Y...
144C**
145      IF(LIND.LE.0)Y=YA(LLE)
146      IF(LIND.LE.0)GO TO 1040
147C**
148C**    IF USING SUBROUTINE QUAD, CALCULATE Y...
149C**
150      YIND=BDER
151      CALL QUAD(XX(1),YY(1),X,YIND)
152      Y=YIND
153C**
154C**    IF USING SUBROUTINE QUAD, CALCULATE DYDX...
155C**
156 1040 YIND=ADER
157      CALL QUAD(XX(1),YY(1),X,YIND)
158      DYDX=YIND
159 1045 CONTINUE
160      IF(X.GE.XA(1).AND.X.LE.XA(NP))GO TO 1050
161      MSSCL(3)=NMSS
162      QSSCL4=Y
163      QSSCL6=X
164      MSSCL(7)=NP
165      MSSCL(8)=LEXL
166      MSSCL(9)=LEXU
167      CALL MESSGE(ROUT,MESS,XA,MJ,MJ,MJ,MSSCL)
168 1050 CONTINUE
169      IF(DEL.EQ.UNUSED)WRITE(6,1060)XDEL,ROUT(1),ROUT(2),(MESS(I),
170     1                      I=1,NMSS)
171 1060 FORMAT(37H WARNING - DUPLICATE X VALUES AT X = ,1PE14.7,
172     1 24H IN TBFUNX, CALLED FROM ,2A4,12H FOR FIGURE ,4A4)
173      RETURN
174 1070 IF(NP.LE.1)Y=YA(1)
175      IF(NP.LE.1)DYDX=0.
176      IF(NP.LE.1)RETURN
177      DEL=0.
178      DYDX=(YA(2)-YA(1))/(XA(2)-XA(1))
179      Y=YA(1)+DYDX*(X-XA(1))
180      GO TO 1045
181      END
182