C
C
C     **************************************************************
C
      SUBROUTINE MODEL(LU,LOU,MM,MTEXT)
C
C     APPROXIMATION OF INTERFACES AND VELOCITY DISTRIBUTION
C     IN INDIVIDUAL LAYERS (BICUBIC SPLINE APPROXIMATION OF VELOCITY)
C
      CHARACTER*80 MTEXT
      DIMENSION IPR(101),X(100),FX(100),AUX(12),VEL(6),Y(4)
      COMMON/MEDIM/V(1000),SX(350),SY(350),NX(20),NY(20),NVS(19),
     1PTOS(19)
      COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20),
     1BRD(2),III(30,20),NPNT(20),NINT
      COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT,
     1NTR
      COMMON/DEN/RHO1(19),RHO2(19),NRHO
      COMMON/AUXX/MMX(20),MMY(20),MMXY(20),MAUX
      COMMON/ABSR/QP1(19),QP2(19),QP3(19),QS1(19),QS2(19),QS3(19),
     1NQP(19),NQS(19),NABS
C
      MI=0
      IF(MM.NE.MI)MM=-1
      IF(MM.LT.0)RETURN
C
C     INPUT DATA - INTERFACES
C
      READ(LU,*)NINT,(NPNT(I),I=1,NINT)
      IF(MPRINT.GE.0)WRITE(LOU,102)NINT,(NPNT(I),I=1,NINT)
      DO 11 I=1,NINT
      NC=NPNT(I)
      READ(LU,*)(X(J),FX(J),III(J,I),J=1,NC)
      IF(MPRINT.GE.0)WRITE(LOU,103)(X(J),FX(J),III(J,I),J=1,NC)
C
C     DETERMINATION OF COEFFICIENTS OF CUBIC PARABOLAS
C     APPROXIMATING INTERFACES
C
      DO 1 J=1,NC
      X1(J,I)=X(J)
    1 A1(J,I)=FX(J)
      J1=1
      NMIN=1
    2 DO 3 J=J1,NC
      J2=J
      IF(III(J,I))4,3,6
    3 CONTINUE
    4 IF(NMIN.EQ.J2)GO TO 5
      FX(NMIN)=A1(NMIN,I)
      CALL SPLIN(X,FX,NMIN,J2)
      KEY=0
      GO TO 8
    5 IF(J2.EQ.NC)GO TO 11
      J1=J2+1
      NMIN=J2
      GO TO 2
    6 IF(NMIN.EQ.J2)GO TO 7
      FX(NMIN)=A1(NMIN,I)
      CALL SPLIN(X,FX,NMIN,J2)
      KEY=1
      GO TO 8
    7 IN=III(J2,I)
      X1(J2,I)=X1(IN,I-1)
      A1(J2,I)=A1(IN,I-1)
      B1(J2,I)=B1(IN,I-1)
      C1(J2,I)=C1(IN,I-1)
      D1(J2,I)=D1(IN,I-1)
      IF(J2.EQ.(NC-1))GO TO 11
      J1=J2+1
      NMIN=J1
      GO TO 2
    8 IF((J2-NMIN).EQ.1)GO TO 10
      J3=J2-1
      DO 9 J=NMIN,J3
      H=X(J+1)-X(J)
      D=(A1(J+1,I)-A1(J,I))/H
      D1(J,I)=(FX(J+1)-FX(J))/(3.*H)
      C1(J,I)=FX(J)
    9 B1(J,I)=D-.333333*H*(FX(J+1)+2.*FX(J))
      IF(KEY)5,5,7
   10 D1(NMIN,I)=0.
      C1(NMIN,I)=0.
      B1(NMIN,I)=(A1(J2,I)-A1(NMIN,I))/(X(J2)-X(NMIN))
      IF(KEY)5,5,7
   11 CONTINUE
C
C     INPUT DATA - VELOCITY DISTRIBUTION IN INDIVIDUAL LAYERS
      WRITE(LOU,117)
      MX2=0
      MY2=0
      MXY2=0
      NC=NINT-1
      DO 13 I=1,NC
      MX1=MX2+1
      MY1=MY2+1
      MXY1=MXY2+1
      READ(LU,*)MX,MY
      IF(MPRINT.GE.0)WRITE(LOU,102)MX,MY
      NX(I)=MX
      NY(I)=MY
      MX2=MX1+MX-1
      MY2=MY1+MY-1
      MXY2=MXY1+MX*MY-1
      READ(LU,*)(SX(J),J=MX1,MX2)
      READ(LU,*)(SY(J),J=MY1,MY2)
      IF(MPRINT.GE.0)WRITE(LOU,104)(SX(J),J=MX1,MX2)
      IF(MPRINT.GE.0)WRITE(LOU,104)(SY(J),J=MY1,MY2)
      M1=MXY1
      DO 12 L=1,MX
      M2=M1+MY-1
      READ(LU,*)(V(J),J=M1,M2)
      IF(MPRINT.GE.0)WRITE(LOU,104)(V(J),J=M1,M2)
   12 M1=M2+1
C
C     DETERMINATION OF COEFFICIENTS OF BICUBIC PARABOLAS
C     APPROXIMATING VELOCITY DISTRIBUTION
C
      CALL BIAP(MX1,MX,MY1,MY,MXY1)
      MMX(I)=MX1
      MMY(I)=MY1
      MMXY(I)=MXY1
   13 CONTINUE
C
C     INPUT DATA - DENSITIES AND S VELOCITIES
C
      READ(LU,*)NRO,(NVS(I),I=1,NC),NABS
      IF(MPRINT.GE.0)WRITE(LOU,102)NRO,(NVS(I),I=1,NC),NABS
      IF(NRO.EQ.0)GO TO 30
      READ(LU,*)(RHO1(I),RHO2(I),I=1,NC)
      IF(MPRINT.GE.0)WRITE(LOU,104)(RHO1(I),RHO2(I),I=1,NC)
   30 IF(NABS.EQ.0)GO TO 34
      DO 35 I=1,NC
      READ(LU,*)NQP(I),NQS(I),QP1(I),QP2(I),QP3(I),QS1(I),QS2(I),
     1QS3(I)
      IF(MPRINT.GE.0)WRITE(LOU,106)
     1NQP(I),NQS(I),QP1(I),QP2(I),QP3(I),QS1(I),QS2(I),QS3(I)
   35 CONTINUE
   34 IF(NRO.NE.0)GO TO 31
      DO 32 IRHO=1,NC
      RHO1(IRHO)=1.7
      RHO2(IRHO)=.2
   32 CONTINUE
   31 CONTINUE
      READ(LU,*)(PTOS(I),I=1,NC)
      IF(MPRINT.GE.0)WRITE(LOU,104)(PTOS(I),I=1,NC)
      DO 33 IRHO=1,NC
   33 IF(PTOS(IRHO).LT.1.41421356)PTOS(IRHO)=1.732
      IF(PTOS(1).GE.100.)RHO1(1)=1.
      IF(PTOS(1).GE.100.)RHO2(1)=0.
      NRHO=0
      IF(PTOS(1).GE.100.)NRHO=1
C
C     INPUT DATA - PRINTER PLOT OF THE VELOCITY DISTRIBUTION
C                  X COORDINATES OF VERTICAL BOUNDARIES OF THE MODEL
C
      READ(LU,*)VMIN,VMAX,BMIN,BMAX,BLEFT,BRIGHT
      IF(MPRINT.GE.0)WRITE(LOU,104)VMIN,VMAX,BMIN,BMAX,BLEFT,BRIGHT
      IF(ABS(BMIN).LT..00001)BMIN=A1(1,1)
      IF(ABS(BMAX).LT..00001)BMAX=A1(1,NINT)
      IF(ABS(BLEFT-BRIGHT).GT..001)GO TO 14
      NC=NPNT(1)
      BLEFT=X1(1,1)
      BRIGHT=X1(NC,1)
C
C     NUMERICAL FORM OF INTERFACES
C
   14 AUX1=(BRIGHT-BLEFT)/25.
      IF(MPRINT.GT.0)WRITE(LOU,101)MTEXT
      IF(MPRINT.GE.2)WRITE(LOU,105)NINT
      IF(MPRINT.GE.2)WRITE(LOU,107)
      DO 19 I=1,NINT
      NC=NPNT(I)-1
      IF(MPRINT.GE.2)WRITE(LOU,112)
      IF(MPRINT.GE.2)WRITE(LOU,108)I
      DO 15 J=1,NC
      IF(MPRINT.GE.2)
     1WRITE(LOU,109)D1(J,I),C1(J,I),B1(J,I),A1(J,I),X1(J,I),X1(J+1,I),
     2III(J,I)
   15 CONTINUE
      IF(MPRINT.GE.2)WRITE(LOU,110)
      AUX2=BLEFT
      AUX3=AUX2
      AUX(1)=AUX2
      L=1
      J=1
   16 AUX3=AUX3+AUX1
      IF(AUX3.GT.BRIGHT)GO TO 17
      L=L+1
      AUX(L)=AUX3
      IF(L.LT.12)GO TO 16
   17 IF(MPRINT.GE.2)WRITE(LOU,111)(AUX(M),M=1,L)
      K=1
   18 IF(AUX(K).GT.X1(J+1,I))J=J+1
      IF(AUX(K).GT.X1(J+1,I))GO TO 18
      AUX4=AUX(K)-X1(J,I)
      AUX(K)=((D1(J,I)*AUX4+C1(J,I))*AUX4+B1(J,I))*AUX4+A1(J,I)
      K=K+1
      IF(K.LE.L)GO TO 18
      IF(MPRINT.GE.2)WRITE(LOU,111)(AUX(M),M=1,L)
      IF(MPRINT.GE.2)WRITE(LOU,112)
      IF((AUX3+AUX1).GT.BRIGHT)GO TO 19
      L=0
      GO TO 16
   19 CONTINUE
C
C     NUMERICAL FORM OF VELOCITY DISTRIBUTION
C
      LAY=1
      ITYPE=1
      Y(3)=1.
      VMM=VMAX-VMIN
      VMMM=VMM/10.
      IF(MPRINT.GE.1)WRITE(LOU,114)VMIN,VMAX,VMMM
      DY=(BMAX-BMIN)/50.
      DX=(BRIGHT-BLEFT)/5.
      Y(2)=BMIN-DY
      K1=1
      AUX(1)=BLEFT
      DO 20 I=2,6
   20 AUX(I)=AUX(I-1)+DX
      IF(MPRINT.GE.1)WRITE(LOU,113)(AUX(I),I=1,6)
      DX=(BRIGHT-BLEFT)/100.
      MAUX=0
      NDER=0
      DO 29 K=1,51
      Y(2)=Y(2)+DY
      Y(1)=BLEFT-DX
      DO 28 L=1,101
      Y(1)=Y(1)+DX
      IF(LAY.GE.NINT)GOTO 24
   21 LAY1=LAY+1
      NC=NPNT(LAY1)-1
      DO 22 I=1,NC
      J=I
      IF(Y(1).LT.X1(I+1,LAY1))GO TO 23
   22 CONTINUE
   23 A2=A1(J,LAY1)
      B2=B1(J,LAY1)
      C2=C1(J,LAY1)
      D2=D1(J,LAY1)
      X2=X1(J,LAY1)
      AUX1=Y(1)-X2
      ZINT=((D2*AUX1+C2)*AUX1+B2)*AUX1+A2
      IF(Y(2).GE.ZINT)LAY=LAY+1
      IF(LAY.GE.NINT)GO TO 27
      IF(Y(2).GT.ZINT)GO TO 21
      IF(LAY.LE.0)GO TO 27
   24 NC=NPNT(LAY)-1
      DO 25 I=1,NC
      J=I
      IF(Y(1).LT.X1(I+1,LAY))GO TO 26
   25 CONTINUE
   26 A2=A1(J,LAY)
      B2=B1(J,LAY)
      C2=C1(J,LAY)
      D2=D1(J,LAY)
      X2=X1(J,LAY)
      AUX1=Y(1)-X2
      ZINT=((D2*AUX1+C2)*AUX1+B2)*AUX1+A2
      IF(Y(2).LT.ZINT)LAY=LAY-1
      IF(LAY.LE.0)GO TO 27
      IF(Y(2).LT.ZINT)GO TO 24
   27 IF(LAY.LE.0.OR.LAY.GE.NINT)IPR(L)=11
      IF(LAY.LE.0.OR.LAY.GE.NINT)GO TO 28
      CALL VELOC(Y,VEL)
      AUX1=10.*(VEL(1)-VMIN)/VMM
      IPR(L)=IFIX(AUX1)
      IF(AUX1.LT.0..OR.AUX1.GT.10.)IPR(L)=11
      ITYPE=1
   28 CONTINUE
C
C     PRINTER PLOT OF THE VELOCITY DISTRIBUTION
C
      IF(K1.EQ.K.AND.MPRINT.GE.1)WRITE(LOU,115)Y(2),IPR
      IF(K1.NE.K.AND.MPRINT.GE.1)WRITE(LOU,116)IPR
      IF(K1.EQ.K)K1=K1+10
   29 CONTINUE
C
  101 FORMAT(1H1/A)
  102 FORMAT(16I5)
  103 FORMAT(3(2F10.5,I5))
  104 FORMAT(8F10.5)
  105 FORMAT(////1X,'M O D E L   O F   M E D I U M'/2X,'NUMBER OF INTERF
     1ACES - ',I3)
  106 FORMAT(2I5,6F10.5)
  107 FORMAT(///////1X,'I N T E R F A C E S'///1X,'INTERFACES ARE APPROX
     1IMATED BY CUBIC PARABOLAS  Z=D*(X-X1)**3+C*(X-X1)**2+B*(X-X1)+A  B
     2ETWEEN  X1  AND  X2'/1X,'COEFFICIENTS OF PARABOLAS ARE DETERMINED
     3BY CUBIC SPLINE INTERPOLATION'///)
  108 FORMAT(/1X,'COEFFICIENTS OF PARABOLAS APPROXIMATING INTERFACE',I3/
     1/15X,'D',14X,'C',14X,'B',14X,'A',2X,'IN INTERVAL',5X,'FROM X1',
     27X,'TO X2',5X,'INDEX')
  109 FORMAT(1X,4E15.5,15X,F10.5,F12.5,I10)
  110 FORMAT(////1X,'NUMERICAL FORM OF INTERFACE'/1X,'UPPER ROW - X-COOR
     1DINATES OF POINTS OF INTERFACE, LOWER ROW - CORRESPONDING Z-COORDI
     2NATES OF POINTS OF INTERFACE'//)
  111 FORMAT(1X,13F9.4)
  112 FORMAT(/)
  113 FORMAT(5X,5(F7.3,13X),F7.3)
  114 FORMAT(////1X,'VELOCITY DISTRIBUTION'/1X,'ISOLINES CONSTRUCTED FRO
     1M',F10.5,'  TO ',F10.5,'  WITH INCREMENT',F10.5//)
  115 FORMAT(1X,F6.2,2X,101I1)
  116 FORMAT(9X,101I1)
  117 FORMAT(1X,'BICUBIC SPLINE INTERPOLATION')
C
      RETURN
      END
C
C     *****************************************************************
C
      SUBROUTINE BIAP(MX1,MX,MY1,MY,MXY1)
C
      DIMENSION X(100),FX(100)
      COMMON/MEDIM/V(1000),SX(350),SY(350),NX(20),NY(20),NVS(19),
     1PTOS(19)
      COMMON/VCOEF/A02(1000),A20(1000),A22(1000)
C
C
C     THIS ROUTINE PERFORMS THE DETERMINATION OF COEFFICIENTS
C     OF THE BICUBIC SPLINE INTERPOLATION
C
      DO 1 J=1,MX
      L=MX1+J-1
    1 X(J)=SX(L)
      DO 3 I=1,MY
      DO 2 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    2 FX(J)=V(K)
      CALL SPLIN(X,FX,1,MX)
      DO 3 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    3 A20(K)=FX(J)
C
      DO 4 I=1,MY
      L=MY1+I-1
    4 X(I)=SY(L)
      DO 6 J=1,MX
      DO 5 I=1,MY
      K=MXY1+(J-1)*MY+I-1
    5 FX(I)=V(K)
      CALL SPLIN(X,FX,1,MY)
      DO 6 I=1,MY
      K=MXY1+(J-1)*MY+I-1
    6 A02(K)=FX(I)
C
      DO 7 J=1,MX
      L=MX1+J-1
    7 X(J)=SX(L)
      DO 9 I=1,MY
      DO 8 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    8 FX(J)=A02(K)
      CALL SPLIN(X,FX,1,MX)
      DO 9 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    9 A22(K)=FX(J)
C
      RETURN
      END
C
C     *************************************************************
C
      SUBROUTINE VELOC(Y,VEL)
C
C     DETERMINATION OF VELOCITY AND ITS DERIVATIVES
C     FOR BICUBIC POLYNOMIAL APPROXIMATION
C
      DIMENSION Y(4),VEL(6)
      COMMON/MEDIM/V(1000),SX(350),SY(350),NX(20),NY(20),NVS(19),
     1PTOS(19)
      COMMON/VCOEF/A02(1000),A20(1000),A22(1000)
      COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20),
     1BRD(2),III(30,20),NPNT(20),NINT
      COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1
      COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT,
     1NTR
      COMMON/AUXX/MMX(20),MMY(20),MMXY(20),MAUX
      COMMON/ABSR/QP1(19),QP2(19),QP3(19),QS1(19),QS2(19),QS3(19),
     1NQP(19),NQS(19),NABS
C
      MX=NX(LAY)
      MY=NY(LAY)
      DO 1 I=2,MX
      K=MMX(LAY)+I-1
      IF(Y(1).LT.SX(K))GO TO 2
    1 CONTINUE
    2 I1=K
      DO 3 I=2,MY
      K=MMY(LAY)+I-1
      IF(Y(2).LT.SY(K))GO TO 4
    3 CONTINUE
    4 J1=K
      IF(MAUX.EQ.0)GO TO 5
      IF(I1.EQ.I2.AND.J1.EQ.J2)GO TO 6
    5 I=I1-MMX(LAY)
      J=J1-MMY(LAY)
      I2=I1
      J2=J1
      XX=SX(I1-1)
      YY=SY(J1-1)
      MM=MMXY(LAY)-1
      K=MM+(I-1)*MY+J
      B20=A20(K)
      B02=A02(K)
      B22=A22(K)
      B00=V(K)
      K=MM+I*MY+J
      C20=A20(K)
      C02=A02(K)
      C22=A22(K)
      C00=V(K)
      K=MM+(I-1)*MY+J+1
      D20=A20(K)
      D02=A02(K)
      D22=A22(K)
      D00=V(K)
      K=MM+I*MY+J+1
      E20=A20(K)
      E02=A02(K)
      E22=A22(K)
      E00=V(K)
      HX=SX(I1)-XX
      HY=SY(J1)-YY
      XA=3.*HX
      YA=3.*HY
      XB=HX/3.
      YB=HY/3.
      D32=(E22-D22)/XA
      D30=(E20-D20)/XA
      B30=(C20-B20)/XA
      B32=(C22-B22)/XA
      D12=(E02-D02)/HX-XB*(E22+2.*D22)
      D10=(E00-D00)/HX-XB*(E20+2.*D20)
      B10=(C00-B00)/HX-XB*(C20+2.*B20)
      B12=(C02-B02)/HX-XB*(C22+2.*B22)
      B03=(D02-B02)/YA
      B13=(D12-B12)/YA
      B23=(D22-B22)/YA
      B33=(D32-B32)/YA
      B01=(D00-B00)/HY-YB*(D02+2.*B02)
      B11=(D10-B10)/HY-YB*(D12+2.*B12)
      B21=(D20-B20)/HY-YB*(D22+2.*B22)
      B31=(D30-B30)/HY-YB*(D32+2.*B32)
      MAUX=1
    6 AX=Y(1)-XX
      AZ=Y(2)-YY
      IF(AX.GE.0..AND.AX.LE.HX.AND.AZ.GE.0..AND.AZ.LE.HY)GO TO 11
      IF(AZ.LT.0..AND.LAY.NE.1)IND=21
      IF(AZ.GT.HY.AND.LAY.NE.(NINT-1))IND=21
C     DETERMINATION OF VELOCITY AT A POINT SITUATED
C     OUTSIDE THE VELOCITY NETWORK
      DO 13 I=2,6
   13 VEL(I)=0.
      RETURN
   11 AUX1=((B33*AZ+B32)*AZ+B31)*AZ+B30
      AUX2=((B23*AZ+B22)*AZ+B21)*AZ+B20
      AUX3=((B13*AZ+B12)*AZ+B11)*AZ+B10
      AUX4=((B03*AZ+B02)*AZ+B01)*AZ+B00
      VEL(1)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4
      VEL(2)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3
      IF(NDER.EQ.0)GO TO 7
      VEL(4)=6.*AUX1*AX+2.*AUX2
    7 AUX1=(3.*B33*AZ+2.*B32)*AZ+B31
      AUX2=(3.*B23*AZ+2.*B22)*AZ+B21
      AUX3=(3.*B13*AZ+2.*B12)*AZ+B11
      AUX4=(3.*B03*AZ+2.*B02)*AZ+B01
      VEL(3)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4
      IF(NDER.EQ.0)GO TO 8
      VEL(5)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3
      AUX1=3.*B33*AZ+B32
      AUX2=3.*B23*AZ+B22
      AUX3=3.*B13*AZ+B12
      AUX4=3.*B03*AZ+B02
      VEL(6)=2.*(((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4)
    8 IF(ITYPE.GT.0.AND.NVS(LAY).EQ.0)GO TO 10
      IF(ITYPE.LT.0.AND.NVS(LAY).EQ.1)GO TO 10
      IF(ITYPE.GT.0.AND.NVS(LAY).EQ.1)GO TO 9
      IF(PTOS(LAY).GE.100.)GO TO 12
      VEL(1)=VEL(1)/PTOS(LAY)
      VEL(2)=VEL(2)/PTOS(LAY)
      VEL(3)=VEL(3)/PTOS(LAY)
      IF(NDER.EQ.0)GO TO 10
      VEL(4)=VEL(4)/PTOS(LAY)
      VEL(5)=VEL(5)/PTOS(LAY)
      VEL(6)=VEL(6)/PTOS(LAY)
      GO TO 10
   12 VEL(1)=0.
      VEL(2)=0.
      VEL(3)=0.
      IF(NDER.EQ.0)GO TO 10
      VEL(4)=0.
      VEL(5)=0.
      VEL(6)=0.
      GO TO 10
    9 VEL(1)=PTOS(LAY)*VEL(1)
      VEL(2)=PTOS(LAY)*VEL(2)
      VEL(3)=PTOS(LAY)*VEL(3)
      IF(NDER.EQ.0)GO TO 10
      VEL(4)=PTOS(LAY)*VEL(4)
      VEL(5)=PTOS(LAY)*VEL(5)
      VEL(6)=PTOS(LAY)*VEL(6)
   10 IF(Y(3).LT.100.)RETURN
      V1=VEL(1)
      AY(6,N)=VEL(1)
      AY(7,N)=VEL(2)
      AY(8,N)=VEL(3)
      IF(NDER.EQ.0)RETURN
      AY(9,N)=VEL(4)
      AY(10,N)=VEL(5)
      AY(11,N)=VEL(6)
      IF(NABS.EQ.0)RETURN
      IF(ITYPE.LT.0)GO TO 14
      Q1=QP1(LAY)
      Q2=QP2(LAY)
      Q3=QP3(LAY)
      NQ=NQP(LAY)
      GO TO 15
   14 Q1=QS1(LAY)
      Q2=QS2(LAY)
      Q3=QS3(LAY)
      NQ=NQS(LAY)
   15 IF(NQ.EQ.0)AY(12,N)=Q1+Q2*V1+Q3*V1*V1
      IF(NQ.EQ.1)AY(12,N)=Q1+Q2*V1+Q3/V1
      RETURN
      END
C
C     *************************************************************
C
      SUBROUTINE SPLIN(X,FX,NMIN,NMAX)
C
      DIMENSION A(100),B(100),H(100),F(100),X(100),FX(100)
C
C     CUBIC SPLINE INTERPOLATION WITH ZERO CURVATURES AT
C     THE END POINTS
C
      IF((NMAX-NMIN).EQ.1)GO TO 4
      NMIN1=NMIN+1
      NMAX1=NMAX-1
      H(NMIN1)=X(NMIN1)-X(NMIN)
      D2=(FX(NMIN1)-FX(NMIN))/H(NMIN1)
      DO 1 I=NMIN1,NMAX1
      H(I+1)=X(I+1)-X(I)
      D1=D2
      D2=(FX(I+1)-FX(I))/H(I+1)
      B(I)=H(I)+H(I+1)
      FX(I)=3.*(D2-D1)/B(I)
      A(I)=H(I)/B(I)
    1 B(I)=H(I+1)/B(I)
    4 FX(NMIN)=0.
      FX(NMAX)=0.
      IF((NMAX-NMIN).EQ.1)RETURN
      H(NMIN)=0.
      F(NMIN)=0.
      DO 2 I=NMIN1,NMAX1
      XPOM=2.+A(I)*H(I-1)
      H(I)=-B(I)/XPOM
    2 F(I)=(FX(I)-A(I)*F(I-1))/XPOM
      DO 3 I=NMIN,NMAX1
      J=NMAX1-(I-NMIN)
    3 FX(J)=H(J)*FX(J+1)+F(J)
      RETURN
      END
C