C  * * * * * * * * * * * * * * * * * * * * * * *                        
C     SUBROUTINE FOR LAGRANGEAN INTERPOLATION                           
C  * * * * * * * * * * * * * * * * * * * * * * *                        
C                                                                       
C     PURPOSE                                                           
C       TO INTERPOLATE THE GIVEN, EQUI-SPACED DATA BY LAGRANGE'S        
C       FOUR-POINTS METHOD                                              
C                                                                       
C     USAGE                                                             
C       CALL LAIN(DX,NN,A,ND1,AA,ND2,NDIV)                              
C                                                                       
C     DESCRIPTION OF ARGUMENTS                                          
C       DX      - INTERVAL OF ORIGINAL/INTERPOLATED DATA                
C                 AT CALL/RETURN                                        
C       NN      - TOTAL NUMBER OF ORIGINAL/INTERPOLATED DATA            
C                 AT CALL/RETURN  NN(AT CALL).GE.4                      
C       A(ND1)  - ORIGINAL EQUI-SPACED DATA                             
C       ND1     - DIMENSION OF A IN CALLING PROGRAM                     
C       AA(ND2) - INTERPOLATED EQUI-SPACED DATA                         
C       ND2     - DIMENSION OF AA IN CALLING PROGRAM                    
C                 ND2.GE.NN(AT CALL)*NDIV                               
C       NDIV    - NUMBER OF SUBDIVISION  NDIV.LE.20                     
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       NONE                                                            
C                                                                       
      SUBROUTINE LAIN(DX,NN,A,ND1,AA,ND2,NDIV)                          
C                                                                       
      DIMENSION A(ND1),AA(ND2),XD(4),D(4),C(80,4),AD(4)                 
      DATA      XD/-1.,0.,1.,2./,D/-6.,2.,-2.,6./                       
C                                                                       
      DIV=REAL(NDIV)                                                    
      I4=NDIV*3+1                                                       
      DO 130 I=1,I4                                                     
      X=REAL(I-NDIV-1)/DIV                                              
      DO 120 J=1,4                                                      
      P=1.                                                              
      DO 110 K=1,4                                                      
      IF(K.EQ.J) GO TO 110                                              
      P=P*(X-XD(K))                                                     
  110 CONTINUE                                                          
      C(I,J)=P/D(J)                                                     
  120 CONTINUE                                                          
  130 CONTINUE                                                          
      N=0                                                               
      IS=1                                                              
      IE=NDIV*2                                                         
      DO 170 M=2,NN-2                                                   
      DO 140 J=1,4                                                      
      AD(J)=A(M+J-2)                                                    
  140 CONTINUE                                                          
      DO 160 I=IS,IE                                                    
      AA1=0.                                                            
      DO 150 J=1,4                                                      
      AA1=AA1+AD(J)*C(I,J)                                              
  150 CONTINUE                                                          
      N=N+1                                                             
      AA(N)=AA1                                                         
  160 CONTINUE                                                          
      IF(M.EQ.2) IS=NDIV+1                                              
      IF(M.EQ.NN-3) IE=I4                                               
  170 CONTINUE                                                          
      NN=N                                                              
      DX=DX/DIV                                                         
      RETURN                                                            
      END                                                               
