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