      INTEGER NT,NX,NTP1,STEPDOWN,NPAD,SLOP                            
      REAL DT,DX,BETA,VC                                                        
      LOGICAL TAPER                                                             
      COMMON CORE(2560000)                                                      
C                                                                               
      CALL LINK("unit6=(timetabl,create,hc)//")                                 
      CALL TIMEF(TINIT)
      CALL SETPRS(NT,NX,STEPDOWN,DT,DX,TRICK,VC,TAPER,NPAD,SLOP)     
      NX=NX+NPAD*2                                                              
      NTP1=NT+1                                                                 
      NTAU=(NT-1)/STEPDOWN                                                      
C                                                                               
      LP = 1                                                                    
      LALP = LP + NTP1*NX*2+SLOP                                            
      LA = LALP + NTP1*NX+SLOP                                                  
      LB = LA + NTP1*NX+SLOP                                                  
      LC = LB + NTP1*NX+SLOP                                                  
      LR1 = LC + NTP1*NX+SLOP                                                 
      LR2 = LR1 + NTP1*NX+SLOP                                                
      LEFAC = LR2 + NTP1*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)                                                                 
      STOP                                                                      
      END                                                                       
      SUBROUTINE PARMIG(P,ALP,A,B,C,R1,R2,EFAC,DT,DX,BETA,NT,NX,VC,NTP1,        
     1TAPER,STEPDOWN,NTAU,NPAD)                                                 
      INTEGER NT,NX,NTP1,IX,JX,IT,ITAU,STEPDOWN,NTAU                            
      INTEGER ITBOT,ITTOP,MIN0,NPAD                                             
      REAL P(NTP1,NX,2),ALP(NTP1,NX),A(NTP1,NX),B(NTP1,NX),C(NTP1,NX)           
      REAL R1(NTP1,NX),R2(NTP1,NX),SC,BETA,DT,DX,VC,FACTOR,EXP,EFAC(NX)         
      REAL FLOAT                                                                
      INTEGER ITEMP,IPREV,INOW                                                  
      INTEGER N1VEL,N2VEL,NXM1,NXD32,IXBOT,IXTOP                                
      LOGICAL TAPER                                                             
      NXM1=NX-1                                                                 
      NXD32=NX/32                                                               
      IXBOT=1+NPAD                                                              
      IXTOP=NX-NPAD                                                             
C                                                                               
C INTIALIZE A SPIKE INPUT SECTION                                               
C                                                                               
      DO 23016 IX=IXBOT,IXTOP                                                   
           DO 23015 IT=1,NTP1                                                   
              P(IT,IX,1)=0.0                                                    
23015       CONTINUE                                                            
         P(NTP1,IX,1)=0.0                                                       
23016    CONTINUE                                                               
      IT=NTP1/2                                                                 
      IX=NX/2                                                                   
      P(IT,IX,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 23022 IX=IXBOT,NXD32                                                
              JX=NX+1-IX                                                        
            DO 23024 IT=1,NT                                                    
               P(IT,IX,1) = P(IT,IX,1)*EFAC(IX)                                 
23024          CONTINUE                                                         
            DO 23026 IT=1,NT                                                    
               P(IT,JX,1) = P(IT,JX,1)*EFAC(IX)                                 
23026          CONTINUE                                                         
23022       CONTINUE                                                            
23018 CONTINUE                                                                  
C                                                                               
C  PRINT OUT TIMES
C                                                                               
      CALL SECOND(CPUTM)                                                        
      ELPSTM=TIMEF(ELPSTM)                                                      
      WRITE(6,113)CPUTM,ELPSTM                                                  
 113  FORMAT(' AFTER SECINIT: CPUSEC=',1PE16.8,'   WALLMSEC=',1PE16.8/)        
C                                                                               
      SC = - (1./32.)*(DT/DX)**2                                                
      IF(.NOT.(VC .EQ. 0.0))GOTO 23028                                          
C        VELFD = AUXIN('VEL')                                                   
         N1VEL = NT                                                             
         N2VEL = NX                                                             
C        CALL AUXPAR('N1','D',N1VEL,'VEL')                                      
C        CALL AUXPAR('N2','D',N2VEL,'VEL')                                      
         IF(.NOT.(N1VEL .LT. NT))GOTO 23030                                     
C           CALL EREXIT('TRACE AND VELOCITY PROFILE LENGTH MISMATCH')           
23030    CONTINUE                                                               
         IF(.NOT.(N2VEL .NE. 1 .AND. N2VEL .LT. NX))GOTO 23032                  
C           CALL EREXIT('TOO FEW VELOCITY PROFILES')                            
23032    CONTINUE                                                               
         DO 23034 IX=IXBOT,IXTOP                                                
C           ICNT=LSEEK(VELFD,(IX-1)*N1VEL*4,0)                                  
C           ICNT=REED(VELFD,ALP(1,IX),NT*4)                                     
C           IF(.NOT.(ICNT .EQ. 0))GOTO 23036                                    
               JX=IX-1                                                          
               DO 23038 ITAU=1,NTAU                                             
                  ALP(ITAU,IX) = ALP(ITAU,JX)                                   
23038             CONTINUE                                                      
               GOTO 23037                                                       
C           ELSE                                                                
23036          CONTINUE                                                         
               DO 23040 IT=1,NT                                                 
                  ALP(IT,IX)=SC*ALP(IT,IX)**2                                   
23040             CONTINUE                                                      
               IF(.NOT.(STEPDOWN .GT. 1))GOTO 23042                             
                  CALL VPACK(ALP(1,IX),NT,ALP(1,IX),NTAU,STEPDOWN)              
23042          CONTINUE                                                         
23037       CONTINUE                                                            
23034       CONTINUE                                                            
         GOTO 23029                                                             
C     ELSE                                                                      
23028    CONTINUE                                                               
         SC = SC * STEPDOWN * VC**2                                             
         DO 23044 IX=IXBOT,IXTOP                                                
            DO 23046 ITAU=1,NTAU                                                
               ALP(ITAU,IX) = SC                                                
23046          CONTINUE                                                         
23044       CONTINUE                                                            
23029 CONTINUE                                                                  
      ITEMP=IXBOT-1                                                             
      DO 23048 IX=ITEMP,1,-1                                                    
           JX=IX+1                                                              
         DO 23050 IT=1,NTP1                                                     
            P(IT,IX,1) = P(IT,JX,1)                                             
23050       CONTINUE                                                            
         DO 23052 ITAU=1,NTAU                                                   
            ALP(ITAU,IX) = ALP(ITAU,JX)                                         
23052       CONTINUE                                                            
23048    CONTINUE                                                               
      ITEMP=IXTOP+1                                                             
      DO 23054 IX=ITEMP,NX,1                                                    
           JX=IX-1                                                              
         DO 23056 IT=1,NTP1                                                     
            P(IT,IX,1) = P(IT,JX,1)                                             
23056       CONTINUE                                                            
         DO 23058 ITAU=1,NTAU                                                   
            ALP(ITAU,IX) = ALP(ITAU,JX)                                         
23058       CONTINUE                                                            
23054    CONTINUE                                                               
      DO 23060 IX=2,NXM1                                                        
         DO 23062 ITAU=1,NTAU                                                   
            A(ITAU,IX)=BETA+ALP(ITAU,IX)                                        
            B(ITAU,IX)=1.0-2.0*A(ITAU,IX)                                       
            C(ITAU,IX)=A(ITAU,IX)                                               
23062       CONTINUE                                                            
23060    CONTINUE                                                               
      DO 23064 ITAU=1,NTAU                                                      
         A(ITAU,NX)=-1.0                                                        
         B(ITAU,NX)=1.0                                                         
23064    CONTINUE                                                               
      DO 23066 ITAU=1,NTAU                                                      
         B(ITAU,1)=1.0                                                          
         C(ITAU,1)=-1.0                                                         
23066    CONTINUE                                                               
C                                                                               
C  PRINT OUT TIMES
C                                                                               
      CALL SECOND(CPUTM)                                                        
      ELPSTM=TIMEF(ELPSTM)                                                      
      WRITE(6,114)CPUTM,ELPSTM                                                  
 114  FORMAT(' BEFORE TRIFAC: CPUSEC=',1PE16.8,'   WALLMSEC=',1PE16.8/)             
C                                                                               
      CALL TRIFAC(A,NTP1,B,NTP1,C,NTP1,NTAU,NX)                                 
      DO 23068 IX=1,NX                                                          
         DO 23070 ITAU=1,NTAU                                                   
            ALP(ITAU,IX) = BETA - ALP(ITAU,IX)                                  
23070       CONTINUE                                                            
23068    CONTINUE                                                               
C                                                                               
C COPY INPUT TO PARALLEL SUBDIAGONAL                                            
C                                                                               
      DO 23072 IX=1,NX                                                          
         DO 23074 IT=1,NTP1                                                     
            P(IT,IX,2)=P(IT,IX,1)                                               
23074       CONTINUE                                                            
23072    CONTINUE                                                               
C                                                                               
C  PRINT OUT TIMES
C                                                                               
      CALL SECOND(CPUTM)                                                        
      ELPSTM=TIMEF(ELPSTM)                                                      
      WRITE(6,115)CPUTM,ELPSTM                                                  
 115  FORMAT(' BEF PMIG LOOP: CPUSEC=',1PE16.8,'   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(IT,1,INOW),P(IT+1,1,INOW),ALP,              
     &    R1,R2,ITTOP)                                                          
         CALL TRIPLY(A,NTP1,B,NTP1,C,NTP1,R2,NTP1,ITTOP,NX)                     
         CALL LHS(NTAU,NX,NTP1,P(IT+1,1,IPREV),P(IT,1,IPREV),R2,ITTOP)          
         ITEMP=IPREV                                                            
         IPREV=INOW                                                             
         INOW=ITEMP                                                             
23076    CONTINUE                                                               
C                                                                               
C  PRINT OUT TIMES
C                                                                               
      CALL SECOND(CPUTM)                                                        
      ELPSTM=TIMEF(ELPSTM)                                                      
      WRITE(6,116)CPUTM,ELPSTM                                                  
 116  FORMAT(' AFT PMIG LOOP: CPUSEC=',1PE16.8,'   WALLMSEC=',1PE16.8/)           
C                                                                               
C                                                                               
C WRITE OUT MIGRATED SECTION                                                    
C                                                                               
      DO 23080 IX=IXBOT,IXTOP                                                   
C        ICNT=RITE(OUTFD,P(1,IX,INOW),NT*4)                                     
23080    CONTINUE                                                               
C                                                                               
C ALL DONE                                                                      
C                                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE RHS(NTAU,NT,NX,NTP1,P1,P3,ALP,R1,R2,ITTOP)                     
      INTEGER NTAU,NT,NX,NTP1,IX,IT,ITAU,ITTOP                                  
      REAL P1(NTP1,NX),P3(NTP1,NX),ALP(NTP1,NX)                                 
      REAL R1(NTP1,NX),R2(NTP1,NX)                                              
      INTEGER IXTOP,MIN0,MOD,ITEMP,NXM1                                         
      NXM1=NX-1                                                                 
      DO 23082 IX=1,NX                                                          
         DO 23084 IT=1,ITTOP                                                    
                R1(IT,IX) = P1(IT,IX)+P3(IT,IX)                                 
23084       CONTINUE                                                            
23082    CONTINUE                                                               
      DO 23086 IT=1,ITTOP                                                       
            R2(IT,1)=0.0                                                        
23086    CONTINUE                                                               
      IXTOP=MIN0((NX-1),(2-1)+MOD(1+(NX-1)-(2),5))                              
      DO 23088 IX=2,IXTOP                                                       
         DO 23090 IT=1,ITTOP                                                    
                R2(IT,IX) = R1(IT,IX) +                                         
     1          ALP(IT,IX)*(R1(IT,IX-1)-2.0*R1(IT,IX)+R1(IT,IX+1))              
23090       CONTINUE                                                            
23088    CONTINUE                                                               
      ITEMP=IXTOP+1                                                             
      DO 23092 IX=ITEMP,NXM1,5                                                  
         DO 23094 IT=1,ITTOP                                                    
                R2(IT,IX) = R1(IT,IX) +                                         
     1          ALP(IT,IX)*(R1(IT,IX-1)-2.0*R1(IT,IX)+R1(IT,IX+1))              
                R2(IT,IX+1) = R1(IT,IX+1) +                                     
     1          ALP(IT,IX+1)*(R1(IT,IX)-2.0*R1(IT,IX+1)+R1(IT,IX+2))            
                R2(IT,IX+2) = R1(IT,IX+2) +                                     
     1          ALP(IT,IX+2)*(R1(IT,IX+1)-2.0*R1(IT,IX+2)+R1(IT,IX+3))          
                R2(IT,IX+3) = R1(IT,IX+3) +                                     
     1          ALP(IT,IX+3)*(R1(IT,IX+2)-2.0*R1(IT,IX+3)+R1(IT,IX+4))          
                R2(IT,IX+4) = R1(IT,IX+4) +                                     
     1          ALP(IT,IX+4)*(R1(IT,IX+3)-2.0*R1(IT,IX+4)+R1(IT,IX+5))          
23094       CONTINUE                                                            
23092    CONTINUE                                                               
      DO 23096 IT=1,ITTOP                                                       
            R2(IT,NX)=0.0                                                       
23096    CONTINUE                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE LHS(NTAU,NX,NTP1,P2,P4,R2,ITTOP)                               
      INTEGER NTAU,NX,NTP1,IX,IT,ITTOP                                          
      REAL P4(NTP1,NX),P2(NTP1,NX)                                              
      REAL R2(NTP1,NX)                                                          
      DO 23098 IX=1,NX                                                          
         DO 23100 IT=1,ITTOP                                                    
              P4(IT,IX) = R2(IT,IX) - P2(IT,IX)                                 
23100       CONTINUE                                                            
23098    CONTINUE                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE TRIFAC(A,I,B,J,C,K,M,N)                                        
      INTEGER I,J,K,M,N                                                         
      REAL A(I,N),B(J,N),C(K,N)                                                 
      INTEGER II,IM,ITEMP                                                       
      IF(.NOT.(N.LE.0))GOTO 23102                                               
         RETURN                                                                 
23102 CONTINUE                                                                  
      DO 23104 II=1,M                                                           
             B(II,1)=1.0/ B(II,1)                                               
             C(II,1)= C(II,1)* B(II,1)                                          
23104    CONTINUE                                                               
      ITEMP=N-1                                                                 
      DO 23106 IM=2,ITEMP,1                                                     
         DO 23108 II=1,M                                                        
                  B(II,IM)=1.0/( B(II,IM)- A(II,IM)* C(II,IM-1))                
                  A(II,IM)= A(II,IM)* B(II,IM-1)                                
                  C(II,IM)= C(II,IM)* B(II,IM)                                  
23108       CONTINUE                                                            
23106    CONTINUE                                                               
      IF(.NOT.(N.EQ.1))GOTO 23110                                               
         RETURN                                                                 
23110 CONTINUE                                                                  
      DO 23112 II=1,M                                                           
             B(II,N)=1.0/( B(II,N)- A(II,N)* C(II,N-1))                         
             A(II,N)= A(II,N)* B(II,N-1)                                        
23112    CONTINUE                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE VPACK(VIN,NT,VOUT,NTAU,STEPDOWN)                               
      INTEGER NT,NTAU,STEPDOWN                                                  
      REAL VIN(STEPDOWN,NTAU),VOUT(NTAU)                                        
      INTEGER ITAU,IT                                                           
CDIR$ NEXTSCALAR                                                                
      DO 23134 IT=2,STEPDOWN                                                    
CDIR$ IVDEP                                                                     
         DO 23136 ITAU=1,NTAU                                                   
            VIN(1,ITAU) = VIN(1,ITAU)+VIN(IT,ITAU)                              
23136       CONTINUE                                                            
23134    CONTINUE                                                               
      DO 23138 ITAU=1,NTAU                                                      
         VOUT(ITAU) = VIN(1,ITAU)                                               
23138    CONTINUE                                                               
      RETURN                                                                    
      END                                                                       
