C	WESTERN MICHIGAN UNIVERSITY
C	INTERP.FOR (FILENAME ON LIBRARY DECTAPE)
C	INTERP, 2.14.1 (CALLING NAME, SUBLST NO.)
C	INTERPOLATION AND CURVE FITTING
C	ADAPTED BY B. GRANET FROM (A) CACM 15,10 (OCT. 1972) 914-918
C	 (B) CACM 17,1 (JAN. 1974) "BIVARIATE INTERPOLATION AND
C	 SMOOTH SURFACE FITTING BASED ON LOCAL PROCEDURES.
C	REPRINTING PRIVILEGES WERE GRANTED BY PERMISSION OF THE
C	 ASSOCIATION FOR COMPUTING MACHINERY BUT NOT FOR PROFIT.
C	LIBRARY DECTAPE PROGS. USED:  USAGE.MAC
C	FORWMU PROGS. USED:  TTYPTY, ALLCOR, DEVCHR, EXISTS,
C	EXIST, GES, GETPPN, JOBNUM, PRINTS, RUNUUO
C	APLB10 PROGS. USED:  IO, GETFOR
C	INTERNAL SUBR. USED:  INTRPL, SFCFIT, CRVFIT, ITPLBV
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
	COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
	1,IRSP,IDUMMY,JDUMMY
	COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
	INTEGER ARRAY(8)
	DIMENSION A(1)
	DATA ARRAY/'UN1CF','UN2CF','BIVCF','BIVIN','UNINT',
	1'HEAD','FORM','HELP'/
	IPAGCT=-1
	IDLG=-1
	IRSP=-4
	NDEVI=4
	NDEVO=6
	NCOL=5
	FMT=0
	WRITE(IDLG,1)
C	CALL USAGE('INTERP')
C---------------1, 0, NDEVO, NDEVI ARE INPUT.  OTHER ARGS. ARE
C--------------- RETURNED.  1 MEANS OUTPUT? PRINTS, 0 MEANS INPUT? PRINTS
	CALL IO(1,NDEVO,DEVNAM,IDVO,IFLNMO,IPJ,IPG,IBNK)
33	CALL IO(0,NDEVI,DEVNAM,IDVI,IFLNMI,IPJ,IPG,IBNK)
C---------------ICODE RETURNED, =0 MEANS TERMINAL JOB, =1 MEANS BATCH
	CALL TTYPTY(ICODE)
46	WRITE(IDLG,3)
11	WRITE(IDLG,4)
	READ(IRSP,5)JTYPE
	DO 6 II=1,8
	IF(ARRAY(II).EQ.JTYPE)
	1GO TO(7,7,8,8,9,97,98,112),II
6	CONTINUE
	GO TO 22
97	WRITE(IDLG,99)
54	READ(IRSP,100)(IDENT(I),I=1,4)
	GO TO 46
98	FMT=1.0
	GO TO 46
7	WRITE(IDLG,102)
103	READ(IRSP,29,ERR=106)L,M
	N=(L-1)*M+1
	IF(II.EQ.1)MD=1
	IF(II.EQ.2)MD=2
	MAX=2*(L+N)
	CALL ALLCOR(MAX,IERR,I1,A)
	IF(IERR.EQ.0)GO TO 114
	WRITE(IDLG,115)
	GO TO 7
114	I2=I1+L
	I3=I2+L
	I4=I3+N
	CALL CRVFIT(NDEVO,MD,L,M,N,A(I1),A(I2),A(I3),A(I4))
	GO TO 33
8	IF(II.EQ.3)WRITE(IDLG,104)
	IF(II.EQ.4)WRITE(IDLG,109)
60	READ(IRSP,29,ERR=105)LX,LY,MX,MY
	IF(II.EQ.4)GO TO 107
	NU=(LX-1)*MX+1
	NV=(LY-1)*MY+1
	MAX=LX+LY+LX*LY+NU+NV+NU*NV
	CALL ALLCOR(MAX,IERR,I1,A)
	IF(IERR.EQ.0)GO TO 116
	WRITE(IDLG,115)
	GO TO 8
116	I2=I1+LX
	I3=I2+LY
	I4=I3+LX*LY
	I5=I4+NU
	I6=I5+NV
	CALL SFCFIT(NDEVO,LX,LY,MX,MY,NU,NV,A(I1),A(I2),A(I3),
	1A(I4),A(I5),A(I6))
	GO TO 33
107	WRITE(IDLG,108)
111	READ(IRSP,29,ERR=63)NINTRP
	MAX=LX+LY+LX*LY+3*NINTRP
	CALL ALLCOR(MAX,IERR,I1,A)
	IF(IERR.EQ.0)GO TO 117
	WRITE(IDLG,115)
	GO TO 8
117	I2=I1+LX
	I3=I2+LY
	I4=I3+LX*LY
	I5=I4+NINTRP
	I6=I5+NINTRP
	CALL ITPLBV(NDEVO,LX,LY,NINTRP,A(I1),A(I2),A(I3),A(I4),
	1A(I5),A(I6))
	GO TO 33
9	WRITE(IDLG,110)
55	READ(IRSP,29,ERR=71)L,NINTRP
	MAX=2*(L+NINTRP)
	CALL ALLCOR(MAX,IERR,I1,A)
	IF(IERR.EQ.0)GO TO 118
	WRITE(IDLG,115)
	GO TO 9
118	I2=I1+L
	I3=I2+L
	I4=I3+NINTRP
	CALL INTRPL(L,NINTRP,A(I1),A(I2),A(I3),A(I4))
	GO TO 33
112	WRITE(IDLG,113)
113	FORMAT(1X,'THE AVAILABLE OPTIONS ARE: '/
	11X,'UN1CF-UNIVARIATE,SINGLE VALUED FUNCTION,CURVE FITTING'/
	21X,'UN2CF-UNIVARIATE MULTIVALUED FUNCTION,CURVE FITTING'/
	31X,'BIVCF-BIVARIATE CURVE FITTING'/
	41X,'BIVIN- BIVARIATE INTERPOLATION'/
	51X,'UNINT-UNIVARIATE,SINGLE VALUED FUNCTION,INTERPOLATION'/
	61X,'HEAD-ENTER IDENTIFICATION FOR OUTPUT.'/
	71X,'FORM-ENTER USER SPECIFIED FORMAT.'/)
	GO TO 46
105	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,53)
	GO TO 60
106	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,53)
	GO TO 103
63	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,53)
	GO TO 111
71	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,53)
	GO TO 55
22	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,47)JTYPE
	GO TO 11
1	FORMAT(1X,'WMU INTERPOLATION'/)
3	FORMAT(1X,'ENTER OPTION.'/)
4	FORMAT(/' *',$)
5	FORMAT(A5)
29	FORMAT(4I)
47	FORMAT(1X,A5,3X,' IS NOT VALID. TRY AGAIN.'/)
53	FORMAT(1X,'ERROR IN INPUT,TRY AGAIN.'/)
99	FORMAT(1X,'ENTER IDENTIFICATION.'/)
100	FORMAT(16A5)
102	FORMAT(1X,'ENTER NO. OF INPUT PTS. AND SUBINTERVALS.'/)
104	FORMAT(1X,'ENTER NO. PTS. IN X AND Y COORDINATES AND NO. '/
	11X,'OF SUBINTERVALS IN X AND Y COORD. SEPARATED BY COMMAS.'/)
108	FORMAT(1X,'ENTER NO. OF PTS. OF INTERP.'/)
109	FORMAT(1X,'ENTER NO. OF PTS. IN X AND Y COORD.',
	1' SEPARATED BY COMMA.'/)
110	FORMAT(1X,'ENTER NO. INPUT PTS. AND NO. PTS. OF INTERP.'/)
115	FORMAT(1X,'NOT ABLE TO ALLOCATE CORE'/)
	END

C---------------L, N ARE INPUT. OTHER ARGS. ARE RETURNED.
      SUBROUTINE  INTRPL(L,N,X,Y,U,V)
C INTERPOLATION OF A SINGLED-VALUED FUNCTION
C THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION
C GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE
C AND  FOR A GIVEN SET OF X VALUES (ABSCISSAS), THE VALUES OF
C A SINGLE-VALUED FUNCTION Y = Y(X).
C THE  INPUT PARAMETERS ARE
C     IU = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT
C     L  = NUMBER OF INPUT DATA POINTS
C          (MUST BE 2 OR GREATER)
C     X  = ARRAY OF DIMENSION L STORING THE X VALUES
C          (ABSCISSAS) OF INPUT DATA POINTS
C          (IN ASCENDING ORDER)
C     Y  = ARRAY OF DIMENSION L STORING THE Y VALUES
C          (ORDINATES) OF INPUT DATA POINTS
C     N  = NUMBER OF POINTS AT WHICH INTERPOLATION OF THE
C          Y VALUE (ORDINATE) IS DESIRED
C          (MUST BE 1 OR GREATER)
C     U  = ARRAY OF DIMENSION N STORING THE X VALUES
C          (ABSCISSAS) OF DESIRED POINTS
C THE  OUTPUT PARAMETER IS
C     V  = ARRAY OF DIMENSION N WHERE THE INTERPOLATED Y
C          VALUES (ORDINATES) ARE TO BE DISPLAYED
C DECLARATION STATEMENTS
C---------------IDLG, IRSP ARE INPUT THRU COMMON /IOB. IDENT, 
C--------------- FMT, NDEVI, NDEVO, ICODE, IDVI ARE INPUT THRU COMMON 
C--------------- /ALL/. IFMT IS RETURNED THRU COMMON /ALL/
	COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
	1,IRSP,IDUMMY,JDUMMY
	COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
      DIMENSION    X(1),Y(1),U(1),V(1)
      EQUIVALENCE  (P0,X3),(Q0,Y3),(Q1,T3)
      REAL         M1,M2,M3,M4,M5
      EQUIVALENCE  (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5),
     1             (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3)
C PRELIMINARY PROCESSING
109	FORMAT(1X,'DATA BEING PROCESSED.'/)
270	FORMAT(1X,'ENTER PTS. OF INTERP.'/)
300	FORMAT(1X,'  I',T10,'X(I)',T23,'Y(I)'/)
320	FORMAT(1X,I3,3X,F10.3,3X,F10.3)
490	FORMAT(1X,'ENTER I,XI,YI IN THIS ORDER.'/)
530	FORMAT(1X,'ERROR IN INPUT,TRY AGAIN.'/)
810	FORMAT(F)
830	FORMAT(/)
1000	FORMAT(16A5)
	J=0
	K=0
	ISTD=1
	IF(FMT.NE.1)GO TO 101
C---------------IFMT, ISTD ARE RETURNED.  OTHER ARGS. ARE INPUT.
C--------------- 96=NO. OF FMT. WORDS FOR OBJ. TIME FORMAT (6 LINES).
C--------------- 4 MEANS UNRESTRICTED FORMAT.
	CALL GETFOR(IRSP,IDLG,IFMT,ISTD,96,4)
101	IF(ISTD.EQ.1)IFMT(1)='(I,2F'
	IF(ISTD.EQ.1)IFMT(2)=')'
	WRITE(IDLG,270)
105	K=K+1
104	READ(IRSP,810,ERR=710)U(K)
	IF(K.NE.N)GO TO 105
	IF(IDVI.EQ.'TTY')GO TO 108
	WRITE(IDLG,109)
	GO TO 103
108	WRITE(IDLG,490)
100	J=J+1
103	READ(NDEVI,IFMT,ERR=520)I1,XI,YI
	J=I1
	X(J)=XI
	Y(J)=YI
	IF(J.NE.L)GO TO 100
	GO TO 10
520	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,530)
	GO TO 103
710	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,530)
	GO TO 104
C---------------COME HERE FROM ST. 80+1 (NEAR END OF SUBR.)
110	WRITE(NDEVO,1000)(IDENT(I),I=1,4)
	WRITE(NDEVO,300)
	DO 310 I=1,N
310	WRITE(NDEVO,320)I,U(I),V(I)
	WRITE(NDEVO,830)
	RETURN
10    L0=L
      LM1=L0-1
      LM2=LM1-1
      LP1=L0+1
      N0=N
      IF(LM2.LT.0)         GO TO 90
      IF(N0.LE.0)          GO TO 91
      DO 11  I=2,L0
        IF(X(I-1)-X(I))    11,95,96
11      CONTINUE
      IPV=0
C MAIN DO-LOOP
      DO 80  K=1,N0
        UK=U(K)
C ROUTINE TO LOCATE THE DESIRED POINT
20      IF(LM2.EQ.0)       GO TO 27
        IF(UK.GE.X(L0))    GO TO 26
        IF(UK.LT.X(1))     GO TO 25
        IMN=2
        IMX=L0
21      I=(IMN+IMX)/2
        IF(UK.GE.X(I))     GO TO 23
22      IMX=I
        GO TO 24
23      IMN=I+1
24      IF(IMX.GT.IMN)     GO TO 21
        I=IMX
        GO TO 30
25      I=1
        GO TO 30
26      I=LP1
        GO TO 30
27      I=2
C CHECK IF I = IPV
30      IF(I.EQ.IPV)       GO TO 70
        IPV=I
C ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND
C          TO ESTIMATE THEM IF NECESSARY
40      J=I
        IF(J.EQ.1)         J=2
        IF(J.EQ.LP1)       J=L0
        X3=X(J-1)
        Y3=Y(J-1)
        X4=X(J)
        Y4=Y(J)
        A3=X4-X3
        M3=(Y4-Y3)/A3
        IF(LM2.EQ.0)       GO TO 43
        IF(J.EQ.2)         GO TO 41
        X2=X(J-2)
        Y2=Y(J-2)
        A2=X3-X2
        M2=(Y3-Y2)/A2
        IF(J.EQ.L0)        GO TO 42
41      X5=X(J+1)
        Y5=Y(J+1)
        A4=X5-X4
        M4=(Y5-Y4)/A4
        IF(J.EQ.2)         M2=M3+M3-M4
        GO TO 45
42      M4=M3+M3-M2
        GO TO 45
43      M2=M3
        M4=M3
45    IF(J.LE.3)           GO TO 46
        A1=X2-X(J-3)
        M1=(Y2-Y(J-3))/A1
        GO TO 47
46      M1=M2+M2-M3
47      IF(J.GE.LM1)       GO TO 48
        A5=X(J+2)-X5
        M5=(Y(J+2)-Y5)/A5
        GO TO 50
48      M5=M4+M4-M3
C NUMERCIAL DIFFERENTIATION
50      IF(I.EQ.LP1)       GO TO 52
        W2=ABS(M4-M3)
        W3=ABS(M2-M1)
        SW=W2+W3
        IF(SW.NE.0.0)      GO TO 51
        W2=0.5
        W3=0.5
        SW=1.0
51      T3=(W2*M2+W3*M3)/SW
        IF(I.EQ.1)         GO TO 54
52      W3=ABS(M5-M4)
        W4=ABS(M3-M2)
        SW=W3+W4
        IF(SW.NE.0.0)      GO TO 53
        W3=0.5
        W4=0.5
        SW=1.0
53      T4=(W3*M3+W4*M4)/SW
        IF(I.NE.LP1)       GO TO 60
        T3=T4
        SA=A2+A3
        T4=0.5*(M4+M5-A2*(A2-A3)*(M2-M3)/(SA*SA))
        X3=X4
        Y3=Y4
        A3=A2
        M3=M4
        GO TO 60
54      T4=T3
        SA=A3+A4
        T3=0.5*(M1+M2-A4*(A3-A4)*(M3-M4)/(SA*SA))
        X3=X3-A4
        Y3=Y3-M2*A4
        A3=A4
        M3=M2
C DETERMINATION OF THE COEFFICIENTS
60      Q2=(2.0*(M3-T3)+M3-T4)/A3
        Q3=(-M3-M3+T3+T4)/(A3*A3)
C COMPUTATION OF THE POLYNOMIAL
70      DX=UK-P0
80      V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3))
	GO TO 110
C ERROR EXIT
90    WRITE (NDEVO,2090)
      GO TO 99
91    WRITE (NDEVO,2091)
      GO TO 99
95    WRITE (IU,2095)
      GO TO 97
96    WRITE (NDEVO,2096)
97    WRITE (NDEVO,2097)  I,X(I)
99    WRITE (NDEVO,2099)  L0,N0
	RETURN
C FORMAT SATATEMENTS
2090  FORMAT(1X/22H  ***   L = 1 OR LESS./)
2091  FORMAT(1X/22H  ***   N = 0 OR LESS./)
2095  FORMAT(1X/27H  ***   IDENTICAL X VALUES./)
2096  FORMAT(1X/33H  ***   X VALUES OUT OF SEQUENCE./)
2097  FORMAT(6H   I =,I7,10X,6HX(I) =,E12.3)
2099  FORMAT(6H   L =,I7,10X,3H( =,I7/
     1       36H ERROR DETECTED IN ROUTINE    INTRPL)
      END
C---------------X,Y,U,V ARE RETURNED.  OTHER ARGS. ARE INPUT.
      SUBROUTINE  CRVFIT(IU,MD,L,M,N,X,Y,U,V)
C SMOOTH CURVE FITTING
C THIS SUBROUTINE FITS A SMOOTH CURVE TO A GIVEN SET OF IN-
C PUT DATA POINTS IN AN X-Y PLANE.  IT INTERPOLATES POINTS
C IN EACH INTERVAL BETWEEN A PAIR OF DATA POINTS AND GENER-
C ATES A SET OF OUTPUT POINTS CONSISTING OF THE INPUT DATA
C POINTS AND THE INTERPOLATED POINTS.  IT CAN PROCESS EITHER
C A SINGLE-VALUED FUNCTION OR A MULTIPLE-VALUED FUNCTION.
C THE INPUT PARAMETERS ARE
C     IU = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT (FOR ERROR 
C	   MESSAGES
C     MD = MODE OF THE CURVE (MUST BE 1 OR 2)
C        = 1 FOR A SINGLE-VALUED FUNCTION
C        = 2 FOR A MULTIPLE-VALUED FUNCTION
C     L  = NUMBER OF INPUT DATA POINTS
C          (MUST BE 2 OR GREATER)
C     X  = ARRAY OF DIMENSION L STORING THE ABSCISSAS OF
C          INPUT DATA POINTS (IN ASCENDING OR DESCENDING
C          ORDER FOR MD = 1)
C     Y  = ARRAY OF DIMENSION L STORING THE ORDINATES OF
C          INPUT DATA POINTS
C     M  = NUMBER OF SUBINTERVALS BETWEEN EACH PAIR OF
C          INPUT DATA POINTS (MUST BE 2 OR GREATER)
C     N  = NUMBER OF OUTPUT POINTS
C        = (L-1)*M+1
C THE OUTPUT PARAMETERS ARE
C     U  = ARRAY OF DIMENSION N WHERE THE ABSCISSAS OF
C          OUTPUT POINTS ARE TO BE DISPLAYED
C     V  = ARRAY OF DIMENSION N WHERE THE ORDINATES OF
C          OUTPUT POINTS ARE TO BE DISPLAYED
C DECLARATION STATEMENTS
C---------------IDLG, IRSP ARE INPUT THRU COMMON /IOB/.  FMT, NDEVI,
C--------------- NDEVO, IDENT, ICODE, IDVI ARE INPUT THRU COMMON
C--------------- /ALL/.  IFMT IS RETURNED THRU COMMON /ALL/.
	COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
	1,IRSP,IDUMMY,JDUMMY
	COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
      DIMENSION    X(1),Y(1),U(1),V(1)
      EQUIVALENCE  (M1,B1),(M2,B2),(M3,B3),(M4,B4),
     1             (X2,P0),(Y2,Q0),(T2,Q1)
      REAL         M1,M2,M3,M4
      EQUIVALENCE  (W2,Q2),(W3,Q3),(A1,P2),(B1,P3),
     1             (A2,DZ),(SW,R,Z)
C PRELIMINARY PROCESSING
109	FORMAT(1X,'DATA BEING PROCESSED.'/)
300	FORMAT(1X,'  I',T10,'X(I)',T23,'Y(I)'/)
320	FORMAT(1X,I3,3X,F10.3,3X,F10.3)
490	FORMAT(1X,'ENTER I,XI,YI IN THIS ORDER.'/)
530	FORMAT(1X,'ERROR IN INPUT ,TRY AGAIN.'/)
830	FORMAT(/)
1000	FORMAT(16A5)
	I=0
	K=0
	ISTD=1
	IF(FMT.NE.1.0)GO TO 105
	CALL GETFOR(IRSP,IDLG,IFMT,ISTD,96,4)
105	IF(ISTD.EQ.1)IFMT(1)='(I,2F'
	IF(ISTD.EQ.1)IFMT(2)=')'
	IF(IDVI.EQ.'TTY')GO TO 108
	WRITE(IDLG,109)
	GO TO 100
108	WRITE(IDLG,490)
100	I=I+1
103	READ(NDEVI,IFMT,ERR=520)I1,XI,YI
	I=I1
	X(I)=XI
	Y(I)=YI
	IF(I.NE.L)GO TO 100
	GO TO 10
520	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,530)
	GO TO 103
110	WRITE(NDEVO,1000)(IDENT(I),I=1,4)
	WRITE(NDEVO,300)
	DO 310 I=1,N
310	WRITE(NDEVO,320)I,U(I),V(I)
	WRITE(NDEVO,830)
	RETURN
10    MD0=MD
      MDM1=MD0-1
      L0=L
      LM1=L0-1
      M0=M
      MM1=M0-1
      N0=N
      IF(MD0.LE.0)        GO TO 90
      IF(MD0.GE.3)        GO TO 90
      IF(LM1.LE.0)        GO TO 91
      IF(MM1.LE.0)        GO TO 92
      IF(N0.NE.LM1*M0+1)  GO TO 93
      GO TO (11,16), MD0
11    I=2
      IF(X(1)-X(2))       12,95,14
12    DO 13  I=3,L0
        IF(X(I-1)-X(I))   13,95,96
13      CONTINUE
      GO TO 18
14    DO 15  I=3,L0
        IF(X(I-1)-X(I))   96,95,15
15      CONTINUE
      GO TO 18
16    DO 17  I=2,L0
        IF(X(I-1).NE.X(I))  GO TO 17
        IF(Y(I- 1).EQ.Y(I))  GO TO 97
17      CONTINUE
18    K=N0+M0
      I=L0+1
      DO 19  J=1,L0
        K=K-M0
        I=I-1
        U(K)=X(I)
19      V(K)=Y(I)
      RM=M0
      RM=1.0/RM
C MAIN DO-LOOP
20    K5=M0+1
      DO 80  I=1,L0
C ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND
C          TO ESTIMATE THEM IF NECESSARY
        IF(I.GT.1)        GO TO 40
30      X3=U(1)
        Y3=V(1)
        X4=U(M0+1)
        Y4=V(M0+1)
        A3=X4-X3
        B3=Y4-Y3
        IF(MDM1.EQ.0)     M3=B3/A3
        IF(L0.NE.2)       GO TO 41
        A4=A3
        B4=B3
31      GO TO (33,32), MD0
32      A2=A3+A3-A4
        A1=A2+A2-A3
33      B2=B3+B3-B4
        B1=B2+B2-B3
        GO TO (51,56), MD0
40      X2=X3
        Y2=Y3
        X3=X4
        Y3=Y4
        X4=X5
        Y4=Y5
        A1=A2
        B1=B2
        A2=A3
        B2=B3
        A3=A4
        B3=B4
        IF(I.GE.LM1)      GO TO 42
41      K5=K5+M0
        X5=U(K5)
        Y5=V(K5)
        A4=X5-X4
        B4=Y5-Y4
        IF(MDM1.EQ.0)     M4=B4/A4
        GO TO 43
42      IF(MDM1.NE.0)     A4=A3+A3-A2
        B4=B3+B3-B2
43      IF(I.EQ.1)        GO TO 31
        GO TO (50,55), MD0
C NUMERICAL DIFFERENTIATION
50      T2=T3
51      W2=ABS(M4-M3)
        W3=ABS(M2-M1)
        SW=W2+W3
        IF(SW.NE.0.0)     GO TO 52
        W2=0.5
        W3=0.5
        SW=1.0
52      T3=(W2*M2+W3*M3)/SW
        IF(I-1) 80,80,60
55      COS2=COS3
        SIN2=SIN3
56      W2=ABS(A3*B4-A4*B3)
        W3=ABS(A1*B2-A2*B1)
        IF(W2+W3.NE.0.0)  GO TO 57
        W2=SQRT(A3*A3+B3*B3)
        W3=SQRT(A2*A2+B2*B2)
57      COS3=W2*A2+W3*A3
        SIN3=W2*B2+W3*B3
        R=COS3*COS3+SIN3*SIN3
        IF(R.EQ.0.0)      GO TO 58
        R=SQRT(R)
        COS3=COS3/R
        SIN3=SIN3/R
58      IF(I-1) 80,80,65
C DETERMINATION OF THE COEFFICIENTS
60      Q2=(2.0*(M2-T2)+M2-T3)/A2
        Q3=(-M2-M2+T2+T3)/(A2*A2)
        GO TO 70
65      R=SQRT(A2*A2+B2*B2)
        P1=R*COS2
        P2=3.0*A2-R*(COS2+COS2+COS3)
        P3=A2-P1-P2
        Q1=R*SIN2
        Q2=3.0*B2-R*(SIN2+SIN2+SIN3)
        Q3=B2-Q1-Q2
        GO TO 75
C COMPUTATION OF THE POLYNOMIALS
70      DZ=A2*RM
        Z=0.0
        DO 71  J=1,MM1
          K=K+1
          Z=Z+DZ
          U(K)=P0+Z
71        V(K)=Q0+Z*(Q1+Z*(Q2+Z*Q3))
        GO TO 79
75      Z=0.0
        DO 76  J=1,MM1
          K=K+1
          Z=Z+RM
          U(K)=P0+Z*(P1+Z*(P2+Z*P3))
76        V(K)=Q0+Z*(Q1+Z*(Q2+Z*Q3))
79      K=K+1
80      CONTINUE
	GO TO 110
C ERROR EXIT
90    WRITE (IU,2090)
      GO TO 99
91    WRITE (IU,2091)
      GO TO 99
92    WRITE (IU,2092)
      GO TO 99
93    WRITE (IU,2093)
      GO TO 99
95    WRITE (IU,2095)
      GO TO 98
96    WRITE (IU,2096)
      GO TO 98
97    WRITE (IU,2097)
98    WRITE (IU,2098)  I,X(I),Y(I)
99    WRITE (IU,2099)  MD0,L0,M0,N0
	RETURN
C FORMAT STATEMENTS
2090  FORMAT(1X/31H  ***   MD OUT OF PROPER RANGE./)
2091  FORMAT(1X/22H  ***   L = 1 OR LESS./)
2092  FORMAT(1X/22H  ***   M = 1 OR LESS./)
2093  FORMAT(1X/25H  ***   IMPROPER N VALUE./)
2095  FORMAT(1X/27H  ***   IDENTICAL X VALUES./)
2096  FORMAT(1X/33H  ***   X VALUES OUT OF SEQUENCE./)
2097  FORMAT(1X/33H  ***   IDENTICAL X AND Y VALUES./)
2098  FORMAT(7H   I  =,I4,10X,6HX(I) =,E12.3,
     1                    10X,6HY(I) =,E12.3)
2099  FORMAT(7H   MD =,I4,8X,3HL =,I5,8X,
     1       3HM =,I5,8X,3HN =,I5/
     2       36H ERROR DETECTED IN ROUTINE    CRVFIT)
      END
C---------------IU,LX, LY ARE INPUT.  OTHER ARGS. ARE RETURNED.
      SUBROUTINE ITPLBV(IU,LX,LY,N,X,Y,Z,U,V,W)
C BIVARIATE INTERPOLATION
C THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION
C GIVEN AT INPUT GRID POINTS IN AN X-Y PLANE AND FOR A GIVEN
C SET OF POINTS IN THE PLANE, THE VALUES OF A SINGLE-VALUED
C BIVARIATE FUNCTION Z = Z(X,Y).
C THE METHOD IS BASED ON A PIECE-WISE FUNCTION COMPOSED OF
C A SET OF BICUBIC POLYNOMIALS IN X AND Y.  EACH POLYNOMIAL
C IS APPLICABLE TO A RECTANGLE OF THE INPUT GRID IN THE X-Y
C PLANE.  EACH POLYNOMIAL IS DETERMINED LOCALLY.
C THE INPUT PARAMETERS ARE
C IU  = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT
C LX  = NUMBER OF INPUT GRID POINTS IN THE X COORDINATE
C       (MUST BE 2 OR GREATER)
C LY  = NUMBER OF INPUT GRID POINTS IN THE Y COORDINATE
C       (MUST BE 2 OR GREATER)
C X   = ARRAY OF DIMENSION LX STORING THE X COORDINATES
C       OF INPUT GRID POINTS (IN ASCENDING ORDER)
C Y   = ARRAY OF DIMENSION LY STORING THE Y COORDINATES
C       OF INPUT GRID POINTS (IN ASENDING ORDER)
C Z   = DOUBLY-DIMENSIONED ARRAY OF DIMENSION (LX,LY)
C     STORING THE VALUES OF THE FUNCTION (Z VALUES)
C     AT INPUT GRID POINTS
C N   = NUMBER OF POINTS AT WHICH INTERPOLATION OF THE
C       Z VALUE IS DESIRED (MUST BE 1 OR GREATER)
C U   = ARRAY OF DIMENSION N STORING THE X COORDINATES
C       OF DESIRED POINTS
C V   = ARRAY OF DIMENSION N STORING THE Y COORDINATES
C       OF DESIRED POINTS
C THE OUTPUT PARAMETER IS
C W   = ARRAY OF DIMENSION N WHERE THE INTERPOLATED Z
C       VALUES AT DESIRED POINTS ARE TO BE DISPLAYED
C       SOME VARIABLES INTERNALLY USED ARE
C ZA  = DIVIDED DIFFERENCE OF Z WITH RESPECT TO X
C ZB  = DIVIDED DIFFERENCE OF Z WITH RESPECT TO Y
C ZAB = SECOND ORDER DIVIDED DIFFERENCE OF Z WITH
C       RESPECT TO X AND Y
C ZX  = PARTIAL DERIVATIVE OF Z WITH RESPECT TO X
C ZY  = PARTIAL DERIVATIVE OF Z WITH RESPECT TO Y
C ZXY = SECOND ORDER PARTIAL DERIVATIVE OF Z WITH
C       RESPECT TO X AND Y
C DECLARATION STATEMENTS
C---------------IDLG, IRSP ARE INPUT THRU COMMON /IOB/.  IDENT,
C--------------- FMT, NDEVI, NDEVO, ICODE, IDVI ARE INPUT THRU COMMON 
C--------------- /ALL/.  IFMT IS RETURNED THRU COMMON /ALL/.
	COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
	1,IRSP,IDUMMY,JDUMMY
	COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
      DIMENSION X(1), Y(1), Z(1), U(1), V(1), W(1)
      DIMENSION ZA(5,2), ZB(2,5), ZAB(3,3), ZX(4,4), ZY(4,4),
     * ZXY(4,4)
      EQUIVALENCE (Z3A1,ZA(1)), (Z3A2,ZA(2)), (Z3A3,ZA(3)),
     * (Z3A4,ZA(4)), (Z3A5,ZA(5)), (Z4A1,ZA(6)), (Z4A2,ZA(7)),
     * (Z4A3,ZA(8)), (Z4A4,ZA(9)), (Z4A5,ZA(10)), (Z3B1,ZB(1)),
     * (Z3B2,ZB(3)), (Z3B3,ZB(5)), (Z3B4,ZB(7)), (Z3B5,ZB(9)),
     * (Z4B1,ZB(2)), (Z4B2,ZB(4)), (Z4B3,ZB(6)), (Z4B4,ZB(8)),
     * (Z4B5,ZB(10)), (ZA2B2,ZAB(1)), (ZA3B2,ZAB(2)),
     * (ZA4B2,ZAB(3)), (ZA2B3,ZAB(4)), (ZA3B3,ZAB(5)),
     * (ZA4B3,ZAB(6)), (ZA2B4,ZAB(7)), (ZA3B4,ZAB(8)),
     * (ZA4B4,ZAB(9)), (ZX33,ZX(6)), (ZX43,ZX(7)),
     * (ZX34,ZX(10)), (ZX44,ZX(11)), (ZY33,ZY(6)),
     * (ZY43,ZY(7)), (ZY34,ZY(10)), (ZY44,ZY(11)),
     * (ZXY33,ZXY(6)), (ZXY43,ZXY(7)), (ZXY34,ZXY(10)),
     * (ZXY44,ZXY(11)), (P00,Z33), (P01,ZY33), (P10,ZX33),
     * (P11,ZXY33)
      EQUIVALENCE (LX0,ZX(1)), (LXM1,ZX(4)), (LXM2,ZX(13)),
     * (LXP1,ZX(16)), (LY0,ZY(1)), (LYM1,ZY(4)), (LYM2,ZY(13)),
     * (LYP1,ZY(16)), (IX,ZXY(1)), (IY,ZXY(4)), (IXPV,ZXY(13)),
     * (IYPV,ZXY(16)), (IMN,JX), (IMX,JY), (JXM2,JX1),
     * (JYM2,JY1), (UK,DX), (VK,DY), (A1,A5,B1,B5,ZX(2),A,Q0),
     * (A2,ZX(5),B,Q1), (A4,ZX(8),C,Q2), (B2,ZY(2),D,Q3),
     * (B4,ZY(14),E), (X2,ZX(3),A3SQ), (X4,ZX(9)), (X5,ZX(12)),
     * (Y2,ZX(14)), (Y4,ZY(3),B3SQ), (Y5,ZX(15),P02),
     * (Z23,ZY(5),P03), (Z24,ZY(8),P12), (Z32,ZY(9),P13),
     * (Z34,ZY(12),P20), (Z35,ZY(15),P21), (Z42,ZXY(2),P22),
     * (Z43,ZXY(5),P23), (Z44,ZXY(3),P30), (Z45,ZXY(8),P31),
     * (Z53,ZXY(9),P32), (Z54,ZXY(12),P33), (W2,WY2,W4),
     * (W3,WY3,W1,W5), (WX2,ZXY(14)), (WX3,ZXY(15))
C PRELIMINARY PROCESSING
C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES
917	FORMAT(1X,'ENTER I,J,XI,YJ,ZIJ SEPARATED BY COMMAS.'/)
916	FORMAT(/)
900	FORMAT(1X,'ENTER PTS. OF INTERP. X,Y.'/)
901	FORMAT(2F)
904	FORMAT(1X,'ERROR IN INPUT ,TRY AGAIN.'/)
906	FORMAT(1X,'DATA BEING PROCESSED.'/)
912	FORMAT(16A5)
913	FORMAT(1X,'  I',T10,'X(I)',T23,'Y(J)',T36,'Z(I,J)'/)
915	FORMAT(1X,I3,3X,F10.3,3X,F10.3,3X,F10.3)
	DO 918 I=1,4
	DO 918 J=1,4
	ZX(I,J)=0
	ZY(I,J)=0
	ZXY(I,J)=0
918	CONTINUE
	DO 919 I=1,5
	DO 919  J=1,2
919	ZA(I,J)=0
	DO 920 I=1,2
	DO 920 J=1,5
920	ZB(I,J)=0
	DO 921 I=1,3
	DO 921 J=1,3
921	ZAB(I,J)=0
 	K=0
	NOPTS=LX*LY
	ISTD=1
	IF(FMT.NE.1)GO TO 905
	CALL GETFOR(IRSP,IDLG,IFMT,ISTD,96,4)
905	IF(ISTD.EQ.1)IFMT(1)='(2I,3'
	IF(ISTD.EQ.1)IFMT(2)='F)'
	WRITE(IDLG,900)
903	K=K+1
	READ(IRSP,901,ERR=902)U(K),V(K)
	IF(K.NE.N)GO TO 903
	K=0
	IF(IDVI.EQ.'TTY')GO TO 907
	WRITE(IDLG,906)
	GO TO 908
907	WRITE(IDLG,917)
908	K=K+1
909	READ(NDEVI,IFMT,ERR=910)I1,J1,XI,YJ,ZIJ
	I=I1
	J=J1
	X(I)=XI
	Y(J)=YJ
	IJ=(I-1)*LY+J
	Z(IJ)=ZIJ
	IF(K.NE.NOPTS)GO TO 908
	GO TO 1
902	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,904)
	GO TO 905
910	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,904)
	GO TO 909
911	WRITE(NDEVO,912)(IDENT(I),I=1,4)
	WRITE(NDEVO,913)
	DO 914 I=1,N
914	WRITE(NDEVO,915)I,U(I),V(I),W(I)
	WRITE(NDEVO,916)
	RETURN
1      IU0 = IU
      LX0 = LX
      LXM1 = LX0 - 1
      LXM2 = LXM1 - 1
      LXP1 = LX0 + 1
      LY0 = LY
      LYM1 = LY0 - 1
      LYM2 = LYM1 - 1
      LYP1 = LY0 + 1
      N0 = N
C ERROR CHECK
      IF (LXM2.LT.0) GO TO 710
      IF (LYM2.LT.0) GO TO 720
      IF (N0.LT.1) GO TO 730
      DO 10 IX=2,LX0
        IF (X(IX-1)-X(IX)) 10, 740, 750
   10 CONTINUE
      DO 20 IY=2,LY0
        IF (Y(IY-1)-Y(IY)) 20, 770, 780
   20 CONTINUE
C INITIAL SETTING OF PREVIOUS VALUES OF IX AND IY
      IXPV = 0
      IYPV = 0
C MAIN DO-LOOP
      DO 700 K=1,N0
      UK = U(K)
      VK = V(K)
C ROUTINE TO LOCATE THE DESIRED POINT
C TO FIND OUT THE IX VALUE FOR WHICH
C (U(K).GE.X(IX-1).AND.(U(K).LT.X(IX))
        IF (LXM2.EQ.0) GO TO 80
        IF (UK.GE.X(LX0)) GO TO 70
        IF (UK.LT.X(1)) GO TO 60
        IMN = 2
        IMX = LX0
   30   IX = (IMN+IMX)/2
        IF (UK.GE.X(IX)) GO TO 40
        IMX = IX
        GO TO 50
   40 IMN = IX + 1
   50   IF (IMX.GT.IMN) GO TO 30
        IX = IMX
        GO TO 90
   60   IX = 1
        GO TO 90
   70   IX = LXP1
        GO TO 90
   80   IX = 2
C TO FIND OUT THE IY VALUE FOR WHICH
C (V(K).GE.Y(IY-1)).AND.(V(K).LT.Y(IY))
   90   IF (LYM2.EQ.0) GO TO 150
        IF (VK.GE.Y(LY0)) GO TO 140
        IF (VK.LT.Y(1)) GO TO 130
        IMN = 2
        IMX = LY0
  100   IY = (IMN+IMX)/2
        IF (VK.GE.Y(IY)) GO TO 110
        IMX = IY
        GO TO 120
  110   IMN = IY + 1
  120   IF (IMX.GT.IMN) GO TO 100
        IY = IMX
        GO TO 160
  130   IY = 1
        GO TO 160
  140   IY = LYP1
        GO TO 160
  150   IY = 2
C TO CHECK IF THE DESIRED POINT IS IN THE SAME RECTANGLE
C AS THE PREVIOUS POINT.  IF YES, SKIP TO THE COMPUTATION
C OF THE POLYNOMIAL
  160   IF (IX.EQ.IXPV .AND. IY.EQ.IYPV) GO TO 690
        IXPV = IX
        IYPV = IY
C ROUTINES TO PICK UP NECESSARY X, Y, AND Z VALUES, TO
C COMPUTE THE ZA, ZB, AND ZAB VALUES, AND TO ESTIMATE THEM
C WHEN NECESSARY
        JX = IX
        IF (JX.EQ.1) JX = 2
        IF (JX.EQ.LXP1) JX = LX0
        JY = IY
        IF (JY.EQ.1) JY = 2
        IF (JY.EQ.LYP1) JY = LY0
        JXM2 = JX - 2
        JXML = JX - LX0
        JYM2 = JY - 2
        JYML = JY - LY0
C IN THE CORE AREA, I.E., IN THE RECTANGLE THAT CONTAINS
C THE DESIRED POINT
        X3 = X(JX-1)
        X4 = X(JX)
        A3 = 1.0/(X4-X3)
        Y3 = Y(JY-1)
        Y4 = Y(JY)
        B3 = 1.0/(Y4-Y3)
	JJ1=(JX-2)*LY+JY-1
	JJ2=JJ1+LY
	JJ3=JJ1+1
	JJ4=JJ3+LY
	Z33=Z(JJ1)
	Z43=Z(JJ2)
	Z34=Z(JJ3)
	Z44=Z(JJ4)
        Z3A3 = (Z43-Z33)*A3
        Z4A3 = (Z44-Z34)*A3
        Z3B3 = (Z34-Z33)*B3
        Z4B3 = (Z44-Z43)*B3
        ZA3B3 = (Z4B3-Z3B3)*A3
C IN THE X DIRECTION
        IF (LXM2.EQ.0) GO TO 230
        IF (JXM2.EQ.0) GO TO 170
        X2 = X(JX-2)
        A2 = 1.0/(X3-X2)
	JZ23=(JX-3)*LY+JY-1
	JZ24=JZ23+1
	Z23=Z(JZ23)
	Z24=Z(JZ24)
        Z3A2 = (Z33-Z23)*A2
        Z4A2 = (Z34-Z24)*A2
        IF (JXML.EQ.0) GO TO 180
  170   X5 = X(JX+1)
        A4 = 1.0/(X5-X4)
	JZ53=JX*LY+JY-1
	JZ54=JZ53+1
	Z53=Z(JZ53)
	Z54=Z(JZ54)
        Z3A4 = (Z53-Z43)*A4
        Z4A4 = (Z54-Z44)*A4
        IF (JXM2.NE.0) GO TO 190
        Z3A2 = Z3A3 + Z3A3 - Z3A4
        Z4A2 = Z4A3 + Z4A3 - Z4A4
        GO TO 190
  180   Z3A4 = Z3A3 + Z3A3 - Z3A2
        Z4A4 = Z4A3 + Z4A3 - Z4A2
  190   ZA2B3 = (Z4A2-Z3A2)*B3
        ZA4B3 = (Z4A4-Z3A4)*B3
        IF (JX.LE.3) GO TO 200
        A1 = 1.0/(X2-X(JX-3))
	J3A1=(JX-4)*LY+JY-1
	J4A1=J3A1+1
	Z3A1=(Z23-Z(J3A1))*A1
	Z4A1=(Z24-Z(J4A1))*A1
        GO TO 210
  200   Z3A1 = Z3A2 + Z3A2 - Z3A3
        Z4A1 = Z4A2 + Z4A2 - Z4A3
  210   IF (JX.GE.LXM1) GO TO 220
        A5 = 1.0/(X(JX+2)-X5)
	J3A5=(JX+1)*LY+JY-1
	J4A5=J3A5+1
	Z3A5=(Z(J3A5)-Z53)*A5
	Z4A5=(Z(J4A5)-Z54)*A5
        GO TO 240
  220   Z3A5 = Z3A4 + Z3A4 - Z3A3
        Z4A5 = Z4A4 + Z4A4 - Z4A3
        GO TO 240
  230   Z3A2 = Z3A3
        Z4A2 = Z4A3
        GO TO 180
C IN THE Y DIRECTION
  240   IF (LYM2.EQ.0) GO TO 310
        IF (JYM2.EQ.0) GO TO 250
        Y2 = Y(JY-2)
        B2 = 1.0/(Y3-Y2)
	JZ32=(JX-2)*LY+JY-2
	JZ42=(JX-1)*LY+JY-2
	Z32=Z(JZ32)
	Z42=Z(JZ42)
        Z3B2 = (Z33-Z32)*B2
        Z4B2 = (Z43-Z42)*B2
        IF (JYML.EQ.0) GO TO 260
  250   Y5 = Y(JY+1)
        B4 = 1.0/(Y5-Y4)
	JZ35=(JX-2)*LY+JY+1
	JZ45=(JX-1)*LY+JY+1
	Z35=Z(JZ35)
	Z45=Z(JZ45)
	Z3B4=(Z35-Z34)*B4
        Z4B4 = (Z45-Z44)*B4
        IF (JYM2.NE.0) GO TO 270
        Z3B2 = Z3B3 + Z3B3 - Z3B4
        Z4B2 = Z4B3 + Z4B3 - Z4B4
        GO TO 270
  260   Z3B4 = Z3B3 + Z3B3 - Z3B2
        Z4B4 = Z4B3 + Z4B3 - Z4B2
  270   ZA3B2 = (Z4B2-Z3B2)*A3
        ZA3B4 = (Z4B4-Z3B4)*A3
        IF (JY.LE.3) GO TO 280
        B1 = 1.0/(Y2-Y(JY-3))
	J3B1=(JX-2)*LY+JY-3
	J4B1=(JX-1)*LY+JY-3
	Z3B1=(Z32-Z(J3B1))*B1
	Z4B1=(Z42-Z(J4B1))*B1
        GO TO 290
  280   Z3B1 = Z3B2 + Z3B2 - Z3B3
        Z4B1 = Z4B2 + Z4B2 - Z4B3
  290   IF (JY.GE.LYM1) GO TO 300
        B5 = 1.0/(Y(JY+2)-Y5)
	J3B5=(JX-2)*LY+JY+2
	J4B5=(JX-1)*LY+JY+2
	Z3B5=(Z(J3B5)-Z35)*B5
	Z4B5=(Z(J4B5)-Z45)*B5
        GO TO 320
  300   Z3B5 = Z3B4 + Z3B4 - Z3B3
        Z4B5 = Z4B4 + Z4B4 - Z4B3
        GO TO 320
  310   Z3B2 = Z3B3
        Z4B2 = Z4B3
        GO TO 260
C IN THE DIAGONAL DIRECTIONS
  320   IF (LXM2.EQ.0) GO TO 400
        IF (LYM2.EQ.0) GO TO 410
        IF (JXML.EQ.0) GO TO 350
        IF (JYM2.EQ.0) GO TO 330
	J4B2=JZ53-1
	ZA4B2=((Z53-Z(J4B2))*B2-Z4B2)*A4
        IF (JYML.EQ.0) GO TO 340
  330	J4B4=J4B2+3
	ZA4B4=((Z(J4B4)-Z54)*B4-Z4B4)*A4
        IF (JYM2.NE.0) GO TO 380
        ZA4B2 = ZA4B3 + ZA4B3 - ZA4B4
        GO TO 380
  340   ZA4B4 = ZA4B3 + ZA4B3 - ZA4B2
        GO TO 380
  350   IF (JYM2.EQ.0) GO TO 360
	J2B2=(JX-3)*LY+JY-2
	ZA2B2=(Z3B2-(Z23-Z(J2B2))*B2)*A2
        IF (JYML.EQ.0) GO TO 370
  360	J2B4=J2B2+3
	ZA2B4=(Z3B4-(Z(J2B4)-Z24)*B4)*A2
        IF (JYM2.NE.0) GO TO 390
        ZA2B2 = ZA2B3 + ZA2B3 - ZA2B4
        GO TO 390
  370   ZA2B4 = ZA2B3 + ZA2B3 - ZA2B2
        GO TO 390
  380   IF (JXM2.NE.0) GO TO 350
        ZA2B2 = ZA3B2 + ZA3B2 - ZA4B2
	ZA2B4=ZA3B4+ZA3B4-ZA4B4
        GO TO 420
  390   IF (JXML.NE.0) GO TO 420
        ZA4B2 = ZA3B2 + ZA3B2 - ZA2B2
        ZA4B4 = ZA3B4 + ZA3B4 - ZA2B4
        GO TO 420
  400   ZA2B2 = ZA3B2
        ZA4B2 = ZA3B2
        ZA2B4 = ZA3B4
        ZA4B4 = ZA3B4
        GO TO 420
  410   ZA2B2 = ZA2B3
        ZA2B4 = ZA2B3
        ZA4B2 = ZA4B3
        ZA4B4 = ZA4B3
C NUMERICAL DIFFERENTIATION   ---   TO DETERMINE PARTIAL
C DERIVATIVES ZX, ZY, AND ZXY AS WEIGHTED MEANS OF DIVIDED
C DIFFERENCES ZA, ZB, AND ZAB, RESPECTIVELY
  420   DO 480 JY=2,3
          DO 470 JX=2,3
            W2 = ABS(ZA(JX+2,JY-1)-ZA(JX+1,JY-1))
            W3 = ABS(ZA(JX,JY-1)-ZA(JX-1,JY-1))
            SW = W2 + W3
            IF (SW.EQ.0.0) GO TO 430
            WX2 = W2/SW
            WX3 = W3/SW
            GO TO 440
  430       WX2 = 0.5
            WX3 = 0.5
  440       ZX(JX,JY) = WX2*ZA(JX,JY-1) + WX3*ZA(JX+1,JY-1)
            W2 = ABS(ZB(JX-1,JY+2)-ZB(JX-1,JY+1))
            W3 = ABS(ZB(JX-1,JY)-ZB(JX-1,JY-1))
            SW = W2 + W3
            IF (SW.EQ.0.0) GO TO 450
            WY2 = W2/SW
            WY3 = W3/SW
            GO TO 460
  450       WY2 = 0.5
            WY3 = 0.5
  460       ZY(JX,JY) = WY2*ZB(JX-1,JY) + WY3*ZB(JX-1,JY+1)
            ZXY(JX,JY) =
     *      WY2*(WX2*ZAB(JX-1,JY-1)+WX3*ZAB(JX,JY-1)) +
     *      WY3*(WX2*ZAB(JX-1,JY)+WX3*ZAB(JX,JY))
  470     CONTINUE
  480   CONTINUE
        IF (IX.EQ.LXP1) GO TO 530
        IF (IX.NE.1) GO TO 590
        W2 = A4*(3.0*A3+A4)
        W1 = 2.0*A3*(A3-A4)+W2
        DO 500 JY = 2,3
          ZX(1,JY)=(W1*ZA(1,JY-1)+W2*ZA(2,JY-1))/(W1+W2)
          ZY(1,JY) = ZY(2,JY) + ZY(2,JY) - ZY(3,JY)
          ZXY(1,JY) = ZXY(2,JY) + ZXY(2,JY) - ZXY(3,JY)
          DO 490 JX1=2,3
            JX = 5 - JX1
            ZX(JX,JY) = ZX(JX-1,JY)
            ZY(JX,JY) = ZY(JX-1,JY)
            ZXY(JX,JY) = ZXY(JX-1,JY)
  490     CONTINUE
  500   CONTINUE
        X3 = X3 - 1.0/A4
        Z33 = Z33 - Z3A2/A4
        DO 510 JY=1,5
          ZB(2,JY) = ZB(1,JY)
  510   CONTINUE
        DO 520 JY=2,4
          ZB(1,JY) = ZB(1,JY) - ZAB(1,JY-1)/A4
  520   CONTINUE
        A3 = A4
        JX = 1
        GO TO 570
  530   W4 = A2*(3.0*A3+A2)
        W5 = 2.0*A3*(A3-A2) + W4
        DO 550 JY=2,3
          ZX(4,JY) = (W4*ZA(4,JY-1)+W5*ZA(5,JY-1))/(W4+W5)
          ZY(4,JY) = ZY(3,JY) + ZY(3,JY) - ZY(2,JY)
          ZXY(4,JY) = ZXY(3,JY) + ZXY(3,JY) - ZXY(2,JY)
          DO 540 JX=2,3
            ZX(JX,JY) = ZX(JX+1,JY)
            ZY(JX,JY) = ZY(JX+1,JY)
            ZXY(JX,JY) = ZXY(JX+1,JY)
  540     CONTINUE
  550   CONTINUE
        X3 = X4
        Z33 = Z43
        DO 560 JY=1,5
          ZB(1,JY) = ZB(2,JY)
  560   CONTINUE
        A3 = A2
        JX = 3
  570   ZA(3,1) = ZA(JX+1,1)
        DO 580 JY=1,3
          ZAB(2,JY) = ZAB(JX,JY)
  580   CONTINUE
C WHEN (V(K).LT.Y(1)).OR.(V(K).GT.Y(LY))
  590   IF (IY.EQ.LYP1) GO TO 630
        IF (IY.NE.1) GO TO 680
        W2 = B4*(3.0*B3+B4)
        W1 = 2.0*B3*(B3-B4) + W2
        DO 620 JX=2,3
          IF (JX.EQ.3 .AND. IX.EQ.LXP1) GO TO 600
          IF (JX.EQ.2 .AND. IX.EQ.1) GO TO 600
          ZY(JX,1) = (W1*ZB(JX-1,1)+W2*ZB(JX-1,2))/(W1+W2)
          ZX(JX,1) = ZX(JX,2) + ZX(JX,2) - ZX(JX,3)
          ZXY(JX,1) = ZXY(JX,2) + ZXY(JX,2) - ZXY(JX,3)
  600     DO 610 JY1=2,3
            JY = 5 - JY1
            ZY(JX,JY) = ZY(JX,JY-1)
            ZX(JX,JY) = ZX(JX,JY-1)
            ZXY(JX,JY) = ZXY(JX,JY-1)
  610     CONTINUE
  620   CONTINUE
        Y3 = Y3 - 1.0/B4
        Z33 = Z33 - Z3B2/B4
        Z3A3 = Z3A3 - ZA3B2/B4
        Z3B3 = Z3B2
        ZA3B3 = ZA3B2
        B3 = B4
        GO TO 670
  630   W4 = B2*(3.0*B3+B2)
        W5 = 2.0*B3*(B3-B2) + W4
        DO 660 JX=2,3
          IF (JX.EQ.3 .AND. IX.EQ.LXP1) GO TO 640
          IF (JX.EQ.2 .AND. IX.EQ.1) GO TO 640
          ZY(JX,4) = (W4*ZB(JX-1,4)+W5*ZB(JX-1,5))/(W4+W5)
          ZX(JX,4) = ZX(JX,3) + ZX(JX,3) - ZX(JX,2)
          ZXY(JX,4) = ZXY(JX,3) + ZXY(JX,3) - ZXY(JX,2)
  640     DO 650 JY=2,3
            ZY(JX,JY) = ZY(JX,JY+1)
            ZX(JX,JY) = ZX(JX,JY+1)
            ZXY(JX,JY) = ZXY(JX,JY+1)
  650     CONTINUE
  660   CONTINUE
        Y3 = Y4
        Z33 = Z33 + Z3B3/B3
        Z3A3 = Z3A3 + ZA3B3/B3
        Z3B3 = Z3B4
        ZA3B3 = ZA3B4
        B3 = B2
  670   IF (IX.NE.1 .AND. IX.NE.LXP1) GO TO 680
        JX = IX/LXP1 + 2
        JX1 = 5 - JX
        JY = IY/LYP1 + 2
        JY1 = 5 - JY
        ZX(JX,JY) = ZX(JX1,JY) + ZX(JX,JY1) - ZX(JX1,JY1)
        ZY(JX,JY) = ZY(JX1,JY) + ZY(JX,JY1) - ZY(JX1,JY1)
        ZXY(JX,JY) = ZXY(JX1,JY) + ZXY(JX,JY1) - ZXY(JX1,JY1)
C DETERMINATION OF THE COEFFICIENTS OF THE POLYNOMIAL
  680   ZX3B3 = (ZX34-ZX33)*B3
        ZX4B3 = (ZX44-ZX43)*B3
        ZY3A3 = (ZY43-ZY33)*A3
        ZY4A3 = (ZY44-ZY34)*A3
        A = ZA3B3 - ZX3B3 - ZY3A3 + ZXY33
        B = ZX4B3 - ZX3B3 - ZXY43 + ZXY33
        C = ZY4A3 - ZY3A3 - ZXY34 + ZXY33
        D = ZXY44 - ZXY43 - ZXY34 + ZXY33
        E = A + A - B - C
        A3SQ = A3*A3
        B3SQ = B3*B3
        P02=(2.0*(Z3B3-ZY33)+Z3B3-ZY34)*B3
        P03=(-2.0*Z3B3+ZY34+ZY33)*B3SQ
        P12=(2.0*(ZX3B3-ZXY33)+ZX3B3-ZXY34)*B3
        P13=(-2.0*ZX3B3+ZXY34+ZXY33)*B3SQ
        P20=(2.0*(Z3A3-ZX33)+Z3A3-ZX43)*A3
        P21=(2.0*(ZY3A3-ZXY33)+ZY3A3-ZXY43)*A3
        P22=(3.0*(A+E)+D)*A3*B3
        P23=(-3.0*E-B-D)*A3*B3SQ
        P30 = (-2.0*Z3A3+ZX43+ZX33)*A3SQ
        P31 = (-2.0*ZY3A3+ZXY43+ZXY33)*A3SQ
        P32 = (-3.0*E-C-D)*B3*A3SQ
        P33 = (D+E+E)*A3SQ*B3SQ
C COMPUTATION OF THE POLYNOMIAL
  690   DY = VK - Y3
        Q0 = P00 + DY*(P01+DY*(P02+DY*P03))
        Q1 = P10 + DY*(P11+DY*(P12+DY*P13))
        Q2 = P20 + DY*(P21+DY*(P22+DY*P23))
        Q3 = P30 + DY*(P31+DY*(P32+DY*P33))
        DX = UK - X3
        W(K) = Q0 + DX*(Q1+DX*(Q2+DX*Q3))
  700 CONTINUE
C NORMAL EXIT
	GO TO 911
C ERROR EXIT
  710 WRITE (IU0,99999)
      GO TO 800
  720 WRITE (IU0,99998)
     GO TO 800
  730 WRITE (IU0,99997)
      GO TO 800
  740 WRITE (IU0,99996)
      GO TO 760
  750 WRITE (IU0,99995)
  760 WRITE (IU0,99994) IX, X(IX)
      GO TO 800
  770 WRITE (IU0,99993)
      GO TO 790
  780 WRITE (IU0,99992)
  790 WRITE (IU0,99991) IY, Y(IY)
  800 WRITE (IU0,99990) LX0, LY0, N0
	RETURN
C FORMAT STATEMENTS
99999 FORMAT(1X/23H  ***  LX = 1 OR LESS./)
99998 FORMAT(1X/23H  ***  LY = 1 OR LESS./)
99997 FORMAT(1X/22H  ***  N = 0 OR LESS./)
99996 FORMAT(1X/27H  ***  IDENTICAL X VALUES./)
99995 FORMAT(1X/33H  ***  X VALUES OUT OF SEQUENCE./)
99994 FORMAT(7H   IX =, I6, 10X, 7HX(IX) =, E12.3)
99993 FORMAT(1X/27H  ***  IDENTICAL Y VALUES./)
99992 FORMAT(1X/33H  ***  Y VALUES OUT OF SEQUENCE./)
99991 FORMAT(7H   IY =, I6, 10X, 7HY(IY) =, E12.3)
99990 FORMAT(7H   LX =, I6, 10X, 4HLY =, I6, 10X, 3HN =, I7/
     *36H ERROR DETECTED IN ROUTINE    ITPLBV)
      END
C---------------IU, LX, LY, MX, MY, NU, NV ARE INPUT.
C--------------- OTHER ARGS. ARE RETURNED.
	SUBROUTINE SFCFIT(IU,LX,LY,MX,MY,NU,NV,X,Y,Z,U,V,W)
C SMOOTH SURFACE FITTING
C THIS SUBROUTINE FITS A SMOOTH SURFACE OF A SINGLE-VALUED
C BIVARIATE FUNTION Z = Z(X,Y) TO A SET OF INPUT DATA
C POINTS GIVEN AT INPUT GRID POINTS IN AN X-Y PLANE. IT
C GENERATES A SET OF OUTPUT GRID POINTS BY EQUALLY DIVIDING
C THE X AND Y COORDINATES IN EACH INTERVAL BETWEEN A PAIR
C OF INPUT GRID POINTS, INTERPOLATES THE Z VALUE FOR THE
C X AND Y VALUES OF EACH OUTPUT GRID POINT, AND GENERATES
C A SET OF OUTPUT POINTS CONSISTING OF INPUT DATA POINTS
C AND THE INTERPOLATED POINTS.
C THE METHOD IS BASED ON A PIECE-WISE FUNTION COMPOSED OF
C A SET OF BICUBIC POLYNOMIALS IN X AND Y. EACH POLYNOMIAL
C IS APPLICABLE TO A RECTANGLE OF THE INPUT GRID IN THE X-Y
C PLANE.  EACH POLYNOMIAL IS DETERMINED LOCALLY.
C THE INPUT PARAMETERS ARE
C IU  = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT
C LX  = NUMBER OF INPUT GRID POINTS IN THE X COORDINATE
C       (MUST BE 2 OR GREATER)
C LY   = NUMBER OF INPUT GRID POINTS IN THE Y COORDINATE
C       (MUST BE 2 OR GREATER)
C X   = ARRAY OF DIMENSION LX STORING THE X COORDINATES
C       OF INPUT GRID POINTS (IN ASCENDING OR DESCENDING
C       ORDER)
C Y   = ARRAY OF DIMENSION LY STORING THE Y COORDINATES
C       OF INPUT GRID POINTS (IN ASCENDING OR DESCENDING
C       ORDER)
C Z   = DOUBLY-DIMENSIONED ARRAY OF DIMENSION (LX,LY)
C       STORING THE VALUES OF THE FUNCTION AT INPUT
C       GRID POINTS
C MX  = NUMBER OF SUBINTERVALS BETWEEN EACH PAIR OF
C       INPUT GRID POINTS IN THE X COORDINATE
C       (MUST BE 2 OR GREATER)
C MY  = NUMBER OF SUBINTERVALS BETWEEN EACH PAIR OF
C       INPUT GRID POINTS IN THE Y COORDINATE
C       (MUST BE 2 OR GREATER)
C NU  NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE
C     = (LX-1)*MX+1
C NV  = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE
C     = (LY-1)*MY+1
C THE OUTPUT PARAMETERS ARE
C U   = ARRAY OF DIMENSION NU WHERE THE X COORDINATES OF
C       OUTPUT POINTS ARE TO BE DISPLAYED
C V   = ARRAY OF DIMENSION NV WHERE THE Y COORINATES OF
C       OUTPUT POINTS ARE TO BE DISPLAYED
C      WHERE THE Z COORIDANETS OF OUTPUT POINTS ARE TO
C     BE DISPLAYED
C SOME VARIABLES INTERNALLY USED ARE
C ZA  = DIVIDED DIFFERENCE OF Z WITH RESPECT TO X
C ZB  = DIVIDED DIFFERENCE OF Z WITH RESPECT TO Y
C ZAB = SECOND ORDER DIVIDED DIFFERENCE OF Z WITH
C       RESPECT TO X AND Y
C ZX  = PARTIAL DERIVATIVE OF Z WITH RESPECT TO X
C ZY  = PARTIAL DERIVATIVE OF Z WITH RESPECT TO Y
C ZXY = SECOND ORDER PARTIAL DERIVATIVE OF Z WITH
C       RESPECT TO X AND Y
C DECLARATION STATEMENTS
C---------------IDLG, IRSP ARE INPUT THRU COMMON /IOB/. IDENT, FMT, NDEVI,
C--------------- NDEVO, ICODE, IDVI ARE INPUT THRU COMMON /ALL/. IFMT
C--------------- IS RETURNED THRU COMMON /ALL/.
	COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
	1,IRSP,IDUMMY,JDUMMY
	COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
      DIMENSION X(1),Y(1),Z(1),U(1),V(1),W(1)
      DIMENSION ZA(4,2), ZB(5), ZAB(2,3), ZX(2), ZY(2), ZXY(2)
      EQUIVALENCE (Z3A2,ZA(1)), (Z3A3,ZA(2)), (Z3A4,ZA(3)),
     * (Z3A5,ZA(4)), (Z4A2,ZA(5)), (Z4A3,ZA(6)), (Z4A4,ZA(7)),
     * (Z4A5,ZA(8)), (Z4B1,ZB(1)), (Z4B2,ZB(2)), (Z4B3,ZB(3)),
     * (Z4B4,ZB(4)), (Z4B5,ZB(5)), (ZA3B2,ZAB(1)),
     * (ZA4B2,ZAB(2)), (ZA3B3,ZAB(3)), (ZA4B3,ZAB(4)),
     * (ZA3B4,ZAB(5)), (ZA4B4,ZAB(6)), (ZX43,ZX(1)),
     * (ZX44,ZX(2)), (ZY43,ZY(1)), (ZY44,ZY(2)),
     * (ZXY43,ZXY(1)), (ZXY44,ZXY(2)), (P00,Z33), (P01,ZY33),
     * (P10,ZX33), (P11,ZXY33)
      EQUIVALENCE (IXM1,JX), (IXML,JY), (DU,DV,DX,DY),
     * (FMX,RMX,FMY,RMY,SW,E), (W2,WY2,A,Q0), (W3,WY3,B,Q1),
     * (WX2,C,Q2), (WX3,D,Q3), (Z3A2,P02), (Z4A2,P03),
     * (Z4B1,P12), (Z4B2,P13), (Z4B4,P20), (Z4B5,P21),
     * (ZA3B2,P22), (ZA3B4,P23)
C PRELIMINARY PROCESSING
C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES
919	FORMAT(1X,T14,1HI,5(F6.1,4X))
920	FORMAT(1X,60(1H-)/)
922	FORMAT(1X,I3,4X,F5.1,1HI,5(F7.2,3X))
924	FORMAT(//1X,T14,1HI,48(1H-)/1X,T14,1HI,T15,' J ='
	1,/1X,T14,1HI,T15,I3,T25,I3,T35,I3,T45,I3,
	2T55,I3)
926	FORMAT(1X,55(1H-)/1X,'I',5X,'X(I)',T14,1HI,T15,'Y(J) =')
906	FORMAT(1X,'DATA BEING PROCESSED.'/)
916	FORMAT(/)
912	FORMAT(16A5)
913	FORMAT(1X,T14,1HI,T32,'W(I,J)'/1X,T14,1HI,48(1H-)/1X,
	1T14,1HI,T15,' J= ',/1X,T14,1HI,T15,'1',T25,'2',T35,'3',
	2T45,'4',T55,'5',/1X,55(1H-),/1X,'I',5X,'X(I)',T14,
	31HI,T15,'Y(J) =')

917	FORMAT(1X,'ENTER I,J,XI,YJ,ZIJ SEPARATED BY COMMAS.'/)
904	FORMAT(1X,'ERROR IN INPUT,TRY AGAIN.'/)
	NCOL=5
	DO 928 I=1,2
	ZX(I)=0
	ZY(I)=0
	ZXY(I)=0
	DO 928 J=1,3
	ZAB(I,J)=0
928	CONTINUE
	DO 929 I=1,5
929	ZB(I)=0
	DO 930 I=1,4
	DO 930 J=1,2
930	ZA(I,J)=0
	K=0
	NOPTS=LX*LY
	ISTD=1
	IF(FMT.NE.1)GO TO 905
	CALL GETFOR(IRSP,IDLG,IFMT,ISTD,96,4)
905	IF(ISTD.EQ.1)IFMT(1)='(2I,3'
	IF(ISTD.EQ.1)IFMT(2)='F)'
	IF(IDVI.EQ.'TTY')GO TO 907
	WRITE(IDLG,906)
	GO TO 908
907	WRITE(IDLG,917)
908	K=K+1
909	READ(NDEVI,IFMT,ERR=910)I1,J1,XI,YJ,ZIJ
	I=I1
	J=J1
	X(I)=XI
	Y(J)=YJ
	IJ=(I-1)*LY+J
	Z(IJ)=ZIJ
	IF(K.NE.NOPTS)GO TO 908
	GO TO 1
911	WRITE(NDEVO,912)(IDENT(I),I=1,4)
	WRITE(NDEVO,913)
	NFIRST=1
	NLAST=NCOL
	IF(NLAST.GT.NV)NLAST=NV
925	IF(NFIRST.LE.NCOL)GO TO 918
	WRITE(NDEVO,924)(I,I=NFIRST,NLAST)
	WRITE(NDEVO,926)
918	WRITE(NDEVO,919)(V(I),I=NFIRST,NLAST)
	WRITE(NDEVO,920)
	DO 921 I=1,NU
	IJ=(I-1)*NV
	WRITE(NDEVO,922)I,U(I),(W(IJ+J),J=NFIRST,NLAST)
921	CONTINUE
	IF(NLAST.GE.NV)GO TO 923
	NFIRST=NFIRST+NCOL
	NLAST=NLAST+NCOL
	IF(NLAST.GT.NV)NLAST=NV
	GO TO 925
923	WRITE(NDEVO,916)
	RETURN
910	IF(ICODE.EQ.-1)CALL EXIT
	WRITE(IDLG,904)
	GO TO 909
1      IU0 = IU
      LX0 = LX
      LXM1 = LX0 - 1
      LXM2 = LXM1 - 1
      LY0 = LY
      LYM1 = LY0 - 1
      LYM2 = LYM1 - 1
      MX0 = MX
      MXP1 = MX0 + 1
      MXM1 = MX0 - 1
      MY0 = MY
      MYP1 = MY0 + 1
      MYM1 = MY0 - 1
      NU0 = NU
      NV0 = NV
C ERROR CHECK
      IF (LXM2.LT.0) GO TO 400
      IF (LYM2.LT.0) GO TO 410
      IF (MXM1.LE.0) GO TO 420
      IF (MYM1.LE.0) GO TO 430
      IF (NU0.NE.LXM1*MX0+1) GO TO 440
      IF (NV0.NE.LYM1*MY0+1) GO TO 450
      IX = 2
      IF (X(1)-X(2)) 10, 460, 30
   10 DO 20 IX=3,LX0
        IF (X(IX-1)-X(IX)) 20, 460, 470
   20 CONTINUE
      GO TO 50
   30 DO 40 IX=3,LX0
        IF (X(IX-1)-X(IX)) 470, 460, 40
   40 CONTINUE
   50 IY = 2
      IF (Y(1)-Y(2)) 60, 490, 80
   60 DO 70 IY=3,LY0
        IF (Y(IY-1)-Y(IY)) 70, 490, 500
   70 CONTINUE
      GO TO 100
   80 DO 90 IY=3,LY0
        IF (Y(IY-1)-Y(IY)) 500, 490, 90
   90 CONTINUE
C COMPUTATION OF THE U ARRAY
  100 FMX = MX0
      RMX = 1.0/FMX
      KU = 1
      X4 = X(1)
      U(1) = X4
      DO 120 IX=2,LX0
        X3 = X4
        X4 = X(IX)
        DU = (X4-X3)*RMX
        DO 110 JX=1,MXM1
          KU = KU + 1
          U(KU) = U(KU-1) + DU
  110   CONTINUE
        KU = KU + 1
      U(KU) = X4
  120 CONTINUE
C COMPUTATION OF THE V ARRAY
      FMY = MY0
      RMY = 1.0/FMY
      KV = 1
      Y4 = Y(1)
      V(1) = Y4
      DO 140 IY=2,LY0
        Y3 = Y4
        Y4 = Y(IY)
        DV = (Y4-Y3)*RMY
        DO 130 JY=1,MYM1
          KV = KV + 1
          V(KV) = V(KV-1) + DV
  130   CONTINUE
        KV = KV + 1
        V(KV) = Y4
  140 CONTINUE
C MAIN DO-LOOPS
      JYMX = MY0
      KV0 = 0
      DO 390 IY=2,LY0
        IYM2 = IY - 2
        IYM3 = IYM2 - 1
        IYML = IY - LY0
        IYML1 = IYML + 1
        IX6 = 0
        IF (IYML.EQ.0) JYMX = MYP1
        JXMX = MX0
        KU0 = 0
        DO 380 IX=1,LX0
          IXM1 = IX - 1
          IXML = IX - LX0
          IF (IXML.EQ.0) JXMX = MXP1
C ROUTINES TO PICK UP NECESSARY X,Y, AND Z VALUES, TO
C COMPUTE THE ZA, ZB, AND ZAB VLAUES, AND TO ESTIMATE THEM
C WHEN NECESSARY
C PRELIMINARY WHEN IX.EQ.1
          IF (IXM1.NE.0) GO TO 150
          Y3 = Y(IY-1)
          Y4 = Y(IY)
          B3 = 1.0/(Y4-Y3)
          B3SQ = B3*B3
      IF (IYM2.GT.0) B2 = 1.0/(Y3-Y(IY-2))
          IF (IYM3.GT.0) B1 = 1.0/(Y(IY-2)-Y(IY-3))
          IF (IYML.LT.0) B4 = 1.0/(Y(IY+1)-Y4)
          IF (IYML1.LT.0) B5 = 1.0/(Y(IY+2)-Y(IY+1))
          GO TO 180
C TO SAVE THE OLD VALUES
  150     Z3A2 = Z3A3
          Z4A2 = Z4A3
          X3 = X4
          Z33 = Z43
          Z3B3 = Z4B3
          A3 = A4
          A3SQ = A3*A3
          Z3A3 = Z3A4
          Z4A3 = Z4A4
          ZA3B2 = ZA4B2
          ZA3B3 = ZA4B3
          ZA3B4 = ZA4B4
  160     X4 = X5
          Z43 = Z53
          Z4B1 = Z5B1
          Z4B2 = Z5B2
          Z4B3 = Z5B3
          Z4B4 = Z5B4
          Z4B5 = Z5B5
          A4 = A5
          Z3A4 = Z3A5
          Z4A4 = Z4A5
          ZA4B2 = ZA5B2
          ZA4B3 = ZA5B3
          ZA4B4 = ZA5B4
  170     X5 = X6
          Z53 = Z63
          Z54 = Z64
          Z5B1 = Z6B1
          Z5B2 = Z6B2
          Z5B3 = Z6B3
          Z5B4 = Z6B4
          Z5B5 = Z6B5
C TO COMPUTE THE ZA,ZB, AND ZAB VALUES AND
C TO ESTIMATE THE ZB VALUES
C WHEN (IY.LE.3).OR.(IY.GE.LY-1)
  180     IX6 = IX6 + 1
          IF (IX6.GT.LX0) GO TO 260
          X6 = X(IX6)
	JZ63=(IX6-1)*LY+IY-1
	Z63=Z(JZ63)
	JZ64=JZ63+1
	Z64=Z(JZ64)
          Z6B3 = (Z64-Z63)*B3
          IF (LYM2.EQ.0) GO TO 200
          IF (IYM2.EQ.0) GO TO 190
	JZ62=JZ64-2
	Z62=Z(JZ62)
          Z6B2 = (Z63-Z62)*B2
          IF (IYML.NE.0) GO TO 190
          Z6B4 = Z6B3 + Z6B3 - Z6B2
          GO TO 210
190	JZ65=JZ64+1
	Z65=Z(JZ65)
          Z6B4 = (Z65-Z64)*B4
          IF (IYM2.NE.0) GO TO 210
          Z6B2 = Z6B3 + Z6B3 - Z6B4
          GO TO 210
  200     Z6B2 = Z6B3
          Z6B4 = Z6B3
  210     IF (IYM3.LE.0) GO TO 220
	J6B1=(IX6-1)*LY+IY-3
	Z6B1=(Z62-Z(J6B1))*B1
          GO TO 230
  220     Z6B1 = Z6B2 + Z6B2 - Z6B3
  230     IF (IYML1.GE.0) GO TO 240
	J6B5=(IX6-1)*LY+IY+2
	Z6B5=(Z(J6B5)-Z65)*B5
          GO TO 250
  240     Z6B5 = Z6B4 + Z6B4 - Z6B3
  250     IF (IX6.EQ.1) GO TO 170
          A5 = 1.0/(X6-X5)
          Z3A5 = (Z63-Z53)*A5
          Z4A5 = (Z64-Z54)*A5
          ZA5B2 = (Z6B2-Z5B2)*A5
          ZA5B3 = (Z6B3-Z5B3)*A5
          ZA5B4 = (Z6B4-Z5B4)*A5
          IF (IX6.EQ.2) GO TO 160
          GO TO 280
C TO ESTIMATE THE ZA AND ZAB VALUES
C WHEN (IX.GE.LX-1).AND.(LX.GT.2)
  260     IF (LXM2.EQ.0) GO TO 270
          Z3A5 = Z3A4 + Z3A4 - Z3A3
          Z4A5 = Z4A4 + Z4A4 -Z4A3
          IF (IXML.EQ.0) GO TO 290
          ZA5B2 = ZA4B2 + ZA4B2 - ZA3B2
          ZA5B3 = ZA4B3 + ZA4B3 - ZA3B3
          ZA5B4 = ZA4B4 + ZA4B4 - ZA3B4
          GO TO 290
C TO ESTIMATE THE ZA AND ZAB VALUES
C WHEN (IX.GE.LX-1).AND.(LX.EQ.2)
  270     Z3A5 = Z3A4
          Z4A5 = Z4A4
          IF (IXML.EQ.0) GO TO 290
          ZA5B2 = ZA4B2
          ZA5B3 = ZA4B3
          ZA5B4 = ZA4B4
C TO ESTIMATE THE ZA AND ZAB VALUES
C WHEN IX.EQ.1
  280     IF (IXM1.NE.0) GO TO 290
          Z3A3 = Z3A4 + Z3A4 - Z3A5
          Z3A2 = Z3A3 + Z3A3 - Z3A4
          Z4A3 = Z4A4 + Z4A4 - Z4A5
          Z4A2 = Z4A3 + Z4A3 - Z4A4
          ZA3B2 = ZA4B2 + ZA4B2 - ZA5B2
          ZA3B3 = ZA4B3 + ZA4B3 - ZA5B3
          ZA3B4 = ZA4B4 + ZA4B4 - ZA5B4
          GO TO 300
C NUMERICAL DIFFERENTIATION   ---   TO DETERMINE PARTIAL
C DERIVATIVES ZX, ZY, AND ZXY AS WEIGHTED MEANS OF DIVIDED
C DIFFERENCES ZA, ZB, AND ZAB, RESPECTIVELY
C TO SAVE THE OLD VALUES WHEN IX.NE.1
  290     ZX33 = ZX43
          ZX34 = ZX44
          ZY33 = ZY43
          ZY34 = ZY44
          ZXY33 = ZXY43
          ZXY34 = ZXY44
C NEW COMPUTATION
  300     DO 350 JY=1,2
            W2 = ABS(ZA(4,JY)-ZA(3,JY))
            W3 = ABS(ZA(2,JY)-ZA(1,JY))
            SW = W2 + W3
            IF (SW.EQ.0.0) GO TO 310
            WX2 = W2/SW
            WX3 = W3/SW
            GO TO 320
  310       WX2 = 0.5
            WX3 = 0.5
  320       ZX(JY) = WX2*ZA(2,JY) + WX3*ZA(3,JY)
            W2 = ABS(ZB(JY+3)-ZB(JY+2))
            W3 = ABS(ZB(JY+1)-ZB(JY))
            SW = W2 + W3
            IF (SW.EQ.0.0) GO TO 330
            WY2 = W2/SW
            WY3 = W3/SW
            GO TO 340
  330       WY2 = 0.5
            WY3 = 0.5
  340       ZY(JY) = WY2*ZB(JY+1) + WY3*ZB(JY+2)
            ZXY(JY) = WY2*(WX2*ZAB(1,JY)+WX3*ZAB(2,JY)) +
     *      WY3*(WX2*ZAB(1,JY+1)+WX3*ZAB(2,JY+1))
  350     CONTINUE
          IF (IXM1.EQ.0) GO TO 380
C DETERMINATION FO THE COEFFICIENTS OF THE POLYNOMIAL
          ZX3B3 = (ZX34-ZX33)*B3
          ZX4B3 = (ZX44-ZX43)*B3
          ZY3A3 = (ZY43-ZY33)*A3
          ZY4A3 = (ZY44-ZY34)*A3
          A = ZA3B3 - ZX3B3 - ZY3A3 + ZXY33
          B = ZX4B3 - ZX3B3 - ZXY43 + ZXY33
          C = ZY4A3 - ZY3A3 - ZXY34 + ZXY33
          D = ZXY44 - ZXY43 - ZXY34 + ZXY33
          E = A + A - B - C
          P02 = (2.0*(Z3B3-ZY33)+Z3B3-ZY34)*B3
          P03 = (-2.0*Z3B3+ZY34+ZY33)*B3SQ
          P12 = (2.0*(ZX3B3-ZXY33)+ZX3B3-ZXY34)*B3
          P13 = (-2.0*ZX3B3+ZXY34+ZXY33)*B3SQ
          P20 = (2.0*(Z3A3-ZX33)+Z3A3-ZX43)*A3
          P21 = (2.0*(ZY3A3-ZXY33)+ZY3A3-ZXY43)*A3
          P22 = (3.0*(A+E)+D)*A3*B3
          P23 = (-3.0*E-B-D)*A3*B3SQ
	P30=(-2.0*Z3A3+ZX43+ZX33)*A3SQ
          P31 = (-2.0*ZY3A3+ZXY43+ZXY33)*A3SQ
          P32 = (-3.0*E-C-D)*B3*A3SQ
          P33 = (D+E+E)*A3SQ*B3SQ
C COMPUTATION OF THE POLYNOMIAL
          DO 370 JY=1,JYMX
            KV = KV0 + JY
            DY = V(KV) - Y3
            Q0 = P00 + DY*(P01+DY*(P02+DY*P03))
            Q1 = P10 + DY*(P11+DY*(P12+DY*P13))
            Q2 = P20 + DY*(P21+DY*(P22+DY*P23))
            Q3 = P30 + DY*(P31+DY*(P32+DY*P33))
            DO 360 JX=1,JXMX
              KU = KU0 + JX
              DX = U(KU) - X3
	JKUKV=(KU-1)*NV+KV
              W(JKUKV) = Q0 + DX*(Q1+DX*(Q2+DX*Q3))
  360       CONTINUE
  370     CONTINUE
          KU0 = KU0 + MX0
  380     CONTINUE
          KV0 = KV0 + MY0
  390 CONTINUE
C NORMAL EXIT
	GO TO 911
C ERROR EXIT
  400 WRITE (IU0,99999)
      GO TO 520
  410 WRITE (IU0,99998)
      GO TO 520
  420 WRITE (IU0,99997)
      GO TO 520
  430 WRITE (IU0,99996)
      GO TO 520
  440 WRITE (IU0,99995)
      GO TO 520
  450 WRITE (IU0,99994)
      GO TO 520
  460 WRITE (IU0,99993)
      GO TO 480
  470 WRITE (IU0,99992)
  480 WRITE (IU0,99991) IX, X(IX)
      GO TO 520
  490 WRITE (IU0,99990)
      GO TO 510
  500 WRITE (IU0,99989)
  510 WRITE (IU0,99988) IY, Y(IY)
  520 WRITE (IU0,99987) LX0, MX0, NU0, LY0, MY0, NV0
      RETURN
C FORMAT STATEMENTS
99999 FORMAT(1X/23H  ***   LX = 1 OR LESS./)
99998 FORMAT(1X/23H  ***   LY = 1 OR LESS./)
99997 FORMAT(1X/22H  ***   N = 0 OR LESS./)
99996 FORMAT(1X/27H  ***   IDENTICAL X VALUES./)
99995 FORMAT(1X/33H  ***   X VALUES OUT OF SEQUENCE./)
99994 FORMAT(7H   IX =, I6, 10X, 7HX(IX) =, E12.3)
99993 FORMAT(1X/27H  ***   IDENTICAL Y VALUES./)
99992 FORMAT(1X/33H  ***   Y VALUES OUT OF SEQUENCE./)
99991 FORMAT(7H   IY =, I6, 10X, 7HY(IY) =, E12.3)
99990 FORMAT(1X/27H  ***   IDENTICAL Y VALUES./)
99989 FORMAT(1X/33H  ***   Y VALUES OUT OF SEQUENCE./)
99988 FORMAT(7H   IY =, I6, 10X, 7HY(IY) =, E12.3)
99987 FORMAT(7H   LX =, I6, 10X 4HMX =, I6, 10X, 4HNU =, I6/
     * 7H   LY =, I6, 10X, 4HMY =, I6, 10X, 4HNV =, I6/6H ERROR,
     * 30H DETECTED IN ROUTINE    SFCFIT)
      END

