      INTEGER NT,NX,NTP1,STEPDOWN,NPAD,SLOP,SLIP,NTPS
      REAL DT,DX,BETA,VC                                                        
      LOGICAL TAPER                                                             
      COMMON CORE(2560000)                                                      
C                                                                               
      CALL SETPRS(NT,NX,STEPDOWN,DT,DX,TRICK,VC,TAPER,NPAD,SLOP,SLIP)     
      NX=NX+NPAD*2                                                              
      NTP1=NT+1                                                                 
      NTPS=NT+SLIP+1
      NTAU=(NT-1)/STEPDOWN                                                      
C                                                                               
      LP = 1                                                                    
      LALP = LP + NTPS*NX*2+SLOP                                            
      LA = LALP + NTPS*NX+SLOP                                                  
      LB = LA + NTPS*NX+SLOP                                                  
      LC = LB + NTPS*NX+SLOP                                                  
      LR1 = LC + NTPS*NX+SLOP                                                 
      LR2 = LR1 + NTPS*NX+SLOP                                                
      LEFAC = LR2 + NTPS*NX+SLOP                                              
      LEND = LEFAC + NX                                                         
      IF((LEND-LP).GT.2560000) STOP 1                                           
      CALL PARMIG(CORE(LP),CORE(LALP),CORE(LA),CORE(LB),CORE(LC),CORE(LR        
     11),CORE(LR2),CORE(LEFAC),DT,DX,BETA,NT,NX,VC,NTP1,TAPER,STEPDOWN,N        
     2TAU,NPAD,NTPS)                                                                 
      STOP                                                                      
      END                                                                       
      SUBROUTINE PARMIG(P,ALP,A,B,C,R1,R2,EFAC,DT,DX,BETA,NT,NX,VC,NTP1,        
     1TAPER,STEPDOWN,NTAU,NPAD,NTPS)
      INTEGER NT,NX,NTP1,IX,JX,IT,ITAU,STEPDOWN,NTAU,NTPS
      INTEGER ITBOT,ITTOP,MIN0,NPAD                                             
      REAL P(NX,NTPS,2),ALP(NX,NTPS),A(NX,NTPS),B(NX,NTPS),C(NX,NTPS)           
      REAL R1(NX,NTPS),R2(NX,NTPS),SC,BETA,DT,DX,VC,FACTOR,EXP,EFAC(NX)
      REAL FLOAT                                                                
      REAL ETIME,TARRAY(2)
      INTEGER ITEMP,IPREV,INOW                                                  
      INTEGER N1VEL,N2VEL,NXM1,NXD32,IXBOT,IXTOP,RITE
      LOGICAL TAPER                                                             
      NXM1=NX-1                                                                 
      NXD32=NX/32                                                               
      IXBOT=1+NPAD                                                              
      IXTOP=NX-NPAD                                                             
C                                                                               
C INTIALIZE A SPIKE INPUT SECTION                                               
C                                                                               
           DO 23015 IT=1,NTP1                                                   
      DO 23016 IX=IXBOT,IXTOP                                                   
              P(IX,IT,1)=0.0                                                    
23016    CONTINUE                                                               
23015       CONTINUE                                                            
      IT=NTP1/2                                                                 
      IX=NX/2                                                                   
      P(IX,IT,1)=1.0                                                            
      P(IX,NT-1,1)=1.0
C                                                                               
C TAPER EDGES TO REDUCE BOUNDARY REFLECTIONS                                    
C                                                                               
      IF(.NOT.TAPER)GOTO 23018                                                  
         DO 23020 IX=IXBOT,NXD32                                                
            FACTOR=FLOAT(IX-1-NXD32)/FLOAT(NXD32)                               
            EFAC(IX)=EXP(-FACTOR**2)                                            
23020       CONTINUE                                                            
            DO 23026 IT=1,NT                                                    
         DO 23022 IX=IXBOT,NXD32                                                
              JX=NX+1-IX                                                        
               P(IX,IT,1) = P(IX,IT,1)*EFAC(IX)                                 
               P(JX,IT,1) = P(JX,IT,1)*EFAC(IX)                                 
23022       CONTINUE                                                            
23026          CONTINUE                                                         
23018 CONTINUE                                                                  
C                                                                               
C  PRINT OUT TIMES
C                                                                               
      ELPSTM=ETIME(TARRAY)*1000.                                                
      WRITE(0,113)ELPSTM                                                  
 113  FORMAT(' AFTER SECINIT:  WALLMSEC=',1PE16.8/)        
C                                                                               
      SC = - (1./32.)*(DT/DX)**2                                                
         SC = SC * STEPDOWN * VC**2                                             
            DO 23046 ITAU=1,NTAU                                                
         DO 23044 IX=IXBOT,IXTOP                                                
               ALP(IX,ITAU) = SC                                                
23044       CONTINUE                                                            
23046          CONTINUE                                                         
      ITEMP=IXBOT-1                                                             
         DO 23050 IT=1,NTP1                                                     
      DO 23048 IX=ITEMP,1,-1                                                    
           JX=IX+1                                                              
            P(IX,IT,1) = P(JX,IT,1)                                             
23048    CONTINUE                                                               
23050       CONTINUE                                                            
         DO 23052 ITAU=1,NTAU                                                   
	 DO 23053 IX=ITEMP,1,-1
	    JX=IX+1
            ALP(IX,ITAU) = ALP(JX,ITAU)                                         
23053    CONTINUE
23052       CONTINUE                                                            
      ITEMP=IXTOP+1                                                             
         DO 23056 IT=1,NTP1                                                     
      DO 23054 IX=ITEMP,NX,1                                                    
           JX=IX-1                                                              
            P(IX,IT,1) = P(JX,IT,1)                                             
23054    CONTINUE                                                               
23056       CONTINUE                                                            
         DO 23058 ITAU=1,NTAU                                                   
	 DO 23062 IX=ITEMP,NX,1
	    JX=IX-1
            ALP(IX,ITAU) = ALP(JX,ITAU)                                         
23062    CONTINUE                                                               
23058       CONTINUE                                                            
      DO 23066 ITAU=1,NTAU                                                      
         B(1,ITAU)=1.0                                                          
         C(1,ITAU)=-1.0                                                         
      DO 23060 IX=2,NXM1                                                        
            A(IX,ITAU)=BETA+ALP(IX,ITAU)                                        
            B(IX,ITAU)=1.0-2.0*A(IX,ITAU)                                       
            C(IX,ITAU)=A(IX,ITAU)                                               
23060    CONTINUE                                                               
         A(NX,ITAU)=-1.0                                                        
         B(NX,ITAU)=1.0                                                         
23066    CONTINUE                                                               
C                                                                               
C  PRINT OUT TIMES
C                                                                               
      ELPSTM=ETIME(TARRAY)*1000.                                                
      WRITE(0,114)ELPSTM                                                  
 114  FORMAT(' BEFORE TRIFAC: WALLMSEC=',1PE16.8/)             
C                                                                               
      CALL TRIFAC(A,NTPS,B,NTPS,C,NTPS,NTAU,NX)
         DO 23070 ITAU=1,NTAU                                                   
      DO 23068 IX=1,NX                                                          
            ALP(IX,ITAU) = BETA - ALP(IX,ITAU)                                  
23068    CONTINUE                                                               
23070       CONTINUE                                                            
C                                                                               
C COPY INPUT TO PARALLEL SUBDIAGONAL                                            
C                                                                               
         DO 23074 IT=1,NTP1                                                     
      DO 23072 IX=1,NX                                                          
            P(IX,IT,2)=P(IX,IT,1)                                               
23072    CONTINUE                                                               
23074       CONTINUE                                                            
C                                                                               
C  PRINT OUT TIMES
C                                                                               
      ELPSTM=ETIME(TARRAY)*1000.                                                
      WRITE(0,115)ELPSTM                                                  
 115  FORMAT(' BEF PMIG LOOP: WALLMSEC=',1PE16.8/)          
C                                                                               
C                                                                               
C NOW FOR PARALLEL MIGRATION                                                    
C                                                                               
      INOW=1                                                                    
      IPREV=2                                                                   
      ITBOT=1+STEPDOWN                                                          
      DO 23076 IT=NT,ITBOT,-1                                                   
         ITTOP = NTP1-IT                                                        
         IF(.NOT.(STEPDOWN.GT.1))GOTO 23078                                     
            ITTOP = MIN0(ITTOP,(IT-1)/(STEPDOWN-1))                             
23078    CONTINUE                                                               
         CALL RHS(NTAU,NT,NX,NTP1,P(1,IT,INOW),P(1,IT+1,INOW),ALP,              
     &    R1,R2,ITTOP,NTPS)                                                          
         CALL TRIPLY(A,NTPS,B,NTPS,C,NTPS,R2,NTPS,ITTOP,NX)
         CALL LHS(NTAU,NX,NTP1,P(1,IT+1,IPREV),P(1,IT,IPREV),R2,ITTOP,
     &            NTPS)          
         ITEMP=IPREV                                                            
         IPREV=INOW                                                             
         INOW=ITEMP                                                             
C                                                                               
C WRITE OUT PARTIALLY MIGRATED SECTION                              
C                                                                               
CC    DO 23080 IX=IXBOT,IXTOP                                                   
CC       ICNT=RITE(OUTFD,P(IX,1,INOW),NT*4)                                     
C23080    CONTINUE                                                               
23076    CONTINUE                                                               
C                                                                               
C  PRINT OUT TIMES
C                                                                               
      ELPSTM=ETIME(TARRAY)*1000.0                                               
      WRITE(0,116)ELPSTM                                                  
 116  FORMAT(' AFT PMIG LOOP: WALLMSEC=',1PE16.8/)           
C                                                                               
C                                                                               
C ALL DONE                                                                      
C                                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE RHS(NTAU,NT,NX,NTP1,P1,P3,ALP,R1,R2,ITTOP,NTPS)
      INTEGER NTAU,NT,NX,NTP1,IX,IT,ITAU,ITTOP                                  
      REAL P1(NX,NTPS),P3(NX,NTPS),ALP(NX,NTPS)
      REAL R1(NX,NTPS),R2(NX,NTPS)                                              
      INTEGER IXTOP,MIN0,MOD,ITEMP,NXM1                                         
      NXM1=NX-1                                                                 
      DO 23096 IT=1,ITTOP                                                       
      DO 23082 IX=1,NX                                                          
                R1(IX,IT) = P1(IX,IT)+P3(IX,IT)                                 
23082    CONTINUE                                                               
            R2(1,IT)=0.0                                                        
      IXTOP=MIN0((NX-1),(2-1)+MOD(1+(NX-1)-(2),5))                              
      DO 23088 IX=2,IXTOP                                                       
                R2(IX,IT) = R1(IX,IT) +                                         
     1          ALP(IX,IT)*(R1(IX-1,IT)-2.0*R1(IX,IT)+R1(IX+1,IT))              
23088    CONTINUE                                                               
      ITEMP=IXTOP+1                                                             
      DO 23092 IX=ITEMP,NXM1,5                                                  
                R2(IX,IT) = R1(IX,IT) +                                         
     1          ALP(IX,IT)*(R1(IX-1,IT)-2.0*R1(IX,IT)+R1(IX+1,IT))              
                R2(IX+1,IT) = R1(IX+1,IT) +                                     
     1          ALP(IX+1,IT)*(R1(IX,IT)-2.0*R1(IX+1,IT)+R1(IX+2,IT))            
                R2(IX+2,IT) = R1(IX+2,IT) +                                     
     1          ALP(IX+2,IT)*(R1(IX+1,IT)-2.0*R1(IX+2,IT)+R1(IX+3,IT))          
                R2(IX+3,IT) = R1(IX+3,IT) +                                     
     1          ALP(IX+3,IT)*(R1(IX+2,IT)-2.0*R1(IX+3,IT)+R1(IX+4,IT))          
                R2(IX+4,IT) = R1(IX+4,IT) +                                     
     1          ALP(IX+4,IT)*(R1(IX+3,IT)-2.0*R1(IX+4,IT)+R1(IX+5,IT))          
23092    CONTINUE                                                               
            R2(NX,IT)=0.0                                                       
23096    CONTINUE                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE LHS(NTAU,NX,NTP1,P2,P4,R2,ITTOP,NTPS)
      INTEGER NTAU,NX,NTP1,IX,IT,ITTOP,NTPS
      REAL P4(NX,NTPS),P2(NX,NTPS)
      REAL R2(NX,NTPS)                                                          
         DO 23100 IT=1,ITTOP                                                    
      DO 23098 IX=1,NX                                                          
              P4(IX,IT) = R2(IX,IT) - P2(IX,IT)                                 
23098    CONTINUE                                                               
23100       CONTINUE                                                            
      RETURN                                                                    
      END                                                                       
      SUBROUTINE TRIFAC(A,I,B,J,C,K,M,N)
      INTEGER I,J,K,M,N
      REAL A(N,I),B(N,J),C(N,K)                                                 
      INTEGER II,IM,ITEMP                                                       
      IF(.NOT.(N.LE.0))GOTO 23102                                               
         RETURN                                                                 
23102 CONTINUE                                                                  
      DO 23110 II=1,M                                                           
             B(1,II)=1.0/ B(1,II)                                               
             C(1,II)= C(1,II)* B(1,II)                                          
      ITEMP=N-1                                                                 
      DO 23106 IM=2,ITEMP,1                                                     
                  B(IM,II)=1.0/( B(IM,II)- A(IM,II)* C(IM-1,II))                
                  A(IM,II)= A(IM,II)* B(IM-1,II)                                
                  C(IM,II)= C(IM,II)* B(IM,II)                                  
23106    CONTINUE                                                               
      IF(.NOT.(N.EQ.1))GOTO 23110                                               
             B(N,II)=1.0/( B(N,II)- A(N,II)* C(N-1,II))                         
             A(N,II)= A(N,II)* B(N-1,II)                                       
23110 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       
      INTEGER FUNCTION RITE(IDUM,BUF,NBYTE)
      INTEGER IDUM,NBYTE
      BYTE BUF(NBYTE)
      INTEGER I,RC,PUTC
      DO 10 I=1,NBYTE
      RC=PUTC(BUF(I))
 10   CONTINUE
      RITE=RC
      RETURN
      END
