      SUBROUTINE BKGPTE ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ,NUMPTS,ARRAY1,ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN, CNSTNT,IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC, 
     * IPARAM )                                                       
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200),           
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1),IPARAM(5)                              
      INTEGER FACFR4                                                    
C                                                                       
C     BALTZ, KAUFFMANN, GLENDENNING AND PRUESS  2+ TELP                 
C                                                                       
C     COMPUTES THE IMAGINARY POTENTIAL AS                               
C                                                                       
C        VI, RI, AI, IMAPOWER   WOODS-SAXON                             
C      + VSI, RSI, ASI          SURFACE TERM                            
C      + BKGP TELP FOR A 2+ STATE.                                      
C                                                                       
C     THE TELP PART IS CONTROLLED BY                                    
C                                                                       
C     IPARAM5 - 1 = PROJECTILE EXCITATION                               
C               2 = TARGET EXCITATION                                   
C               IF NOT DEFINED, THEN  JA @= JB  OR  JBIGA @= JBIGB      
C          IS USED TO DETERMIN WHICH NUCLEUS IS EXCITED.                
C     PARAM20 - CONTAINS THE G2(ETA,XI) FUDGE FACTOR.  IF NOT DEFINED,  
C               THE VALUE 1 IS USED.                                    
C     BETACOULOMB                                                       
C     WE STORE SOME CONSTANTS IN THE BEGINNING OF THE WORK AREA         
C     ALLOC(KLOC(IWRK)+I):                                               
C     I   STORED                                                        
C     0   R2FAC   R2 = R2FAC*RC                                         
C     1   BETACOULOMB                                                   
C     2   ETA                                                           
C     3   WPP = (3/8) WP                                                
C     4   AC  = A PARAMETER FOR L = L(CRITICAL)                         
C     5   BC                                                            
C     6   CC                                                            
C     7   N1 = BOUNDARY AT RC.                                          
C     8   AK (K-VALUE)                                                  
C     9   RC                                                            
C                                                                       
C     3/13/79 - MADE FROM LTSTELP - S.P.                                
C     7/30/79 - IPARAM ADDED TO ARGUMENT LIST - S.P.                    
C     1/16/80 - STANDARD MORTRAN - RPG                                  
C                                                                       
      CHARACTER*4 ID                                                    
      CHARACTER*4 WORKNM(2,2)                                           
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                       
      CHARACTER*8 WORDS(2)                                              
      DATA WORDS/'PROJCTIL', '  TARGET'/                                
C                                                                       
C     DEFINITIONS OF THE COEFFICENTS IN THE 1/R EXPANSION               
C                                                                       
      ATERM(DL) = ETA*AK**2/DL**2 *(ETA*(3.0D0*DL**2+ETASQ)/
     * (DL**2+ETASQ)**2  - DATAN(DL/ETA)/DL )                                         
      BTERM(DL) = 4.0D0*ETA*AK*DL**2/(DL**2+ETASQ)**2                       
      CTERM(DL) = 2.0D0*DL**4/(DL**2+ETASQ)**2                              
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
C                                                                       
      IF ( IREQUE .EQ. 4 ) GO TO 700                                    
C                                                                       
      UNDEF = FINTRN(1)                                                 
      NOTDEF = INTRNL(3)                                                
C                                                                       
      VI = TEMPVS(2)                                                    
      AI = TEMPVS(6)                                                    
      RI = FLOATK(47)                                                    
      POWER = FLOATK(103)                                                
      VSI = FLOATK(71)                                                   
      RSI = FLOATK(69)                                                   
      ASI = FLOATK(68)                                                   
      G2 = PARAM(20)                                                    
      IF ( G2 .EQ. UNDEF ) G2 = 1.0D0                                       
      ICHANW = INTRNL(5)                                                
      AK = FKANDM(ICHANW)                                               
      RC = FLOATK(45)                                                    
      ZP = INTGER(24)                                                   
      ZT = INTGER(25)                                                   
      ETA = FKANDM(ICHANW+16)                                           
      ALPHA = 1.0D0/CNSTNT(8)                                               
      PI = CNSTNT(1)                                                    
      HBARC = CNSTNT(6)                                                 
      E = FLOATK(12)                                                     
      AM = FLOATK(24)                                                    
C                                                                       
      IF ( IREQUE .GT. 1 ) GO TO 200                                    
C                                                                       
C     IS THIS PROJECTILE OR TARGET EXCITATION.                          
C     WE USE THIS INFORMATION TO CONSTRUCT R2FAC WHICH MAKES THE        
C     SINGLE NUCLEON RADIUS OUT OF THE SUM (RC).                        
C                                                                       
      IPORT = IPARAM(5)                                                 
      IF ( IPORT .NE. NOTDEF ) GO TO 50                                 
      IF ( JBLOCK(2) .NE. JBLOCK(3) ) IPORT = 1                         
      IF ( JBLOCK(4) .NE. JBLOCK(5) ) IPORT = 2                         
 50   IF ( IPORT .NE. 1 .AND. IPORT .NE. 2 ) GO TO 930                  
C      AMP = FLOATK(40)                                                  
C      AMT = FLOATK(41)                                                  
      R0MASS = FINTRN(10)                                               
      R2FAC = FLOATK(39+IPORT)**(1.D0/3.D0) / R0MASS                     
C                                                                       
      IB = NAMLOC('BETACOUL')                                           
      IF ( IB .EQ. 0 ) GO TO 940                                        
C                                                                       
      IF ( VI .EQ. 0.0D0 ) GO TO 65                                         
      IF ( RI .EQ. UNDEF .OR. AI .EQ. UNDEF ) GO TO 950                 
  65  IF ( VSI .EQ. 0.0D0 ) GO TO 70                                        
      IF ( RSI .EQ. UNDEF .OR. ASI .EQ. UNDEF ) GO TO 945               
C                                                                       
C     GENERATE WORK ARRAY TO SAVE THE W.S.                              
C                                                                       
  70  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
C                                                                       
      IWRK = NALLOC( NUMPTS+10, WORKNM(1,1) )                           
      LWRK = KLOC(IWRK)                                                  
      ALLOC(LWRK) = R2FAC                                               
      ALLOC(LWRK+1) = ALLOC(KLOC(IB))                                    
      GO TO 300                                                         
C                                                                       
C     COMPUTE THE CONSTANTS                                             
C     THIS PART IS DONE FOR IREQUE = 2 (PRINTING) AND 3 (NEW POTENTIAL) 
C                                                                       
 200  IWRK = MYINTS(1)                                                  
      LWRK = KLOC(IWRK)                                                  
      IF ( VI.NE.0.D0.AND.(AI.LE.0.0D0.OR.RI.LE.0.0D0 ))  GO TO 960    
      IF ( VSI.NE.0.0D0.AND.(ASI.LE.0.0D0.OR.RSI.LE.0.0D0))  GO TO 960 
      RC2 = RC*ALLOC(LWRK)                                              
      BETA = ALLOC(LWRK+1)                                              
      RD = ZP*ZT*ALPHA*HBARC/E                                          
      WP = 3.0D0*(ZP*ZT*ALPHA*BETA)**2*AM*G2*RC2**4/(50.0D0*PI*AK)  
      WPP = 3.0D0*WP/8.0D0                                                   
      ALLOC(LWRK+2) = ETA                                               
      ALLOC(LWRK+3) = WPP                                               
      ALLOC(LWRK+8) = AK                                                
      ETASQ = ETA**2                                                    
C                                                                       
C     ESTIMATE CRITICAL L (IN VERY SIMPLE FASHION) AND GET CONSTANT     
C     (IN L) PARAT OF THE POTENTIAL PARAMETERS.                         
C                                                                       
      DLCSQ = AK*RC*(AK*RC - 2.0D0*ETA)                                     
      IF ( DLCSQ .LT. 0.25D0 ) DLCSQ = 0.25D0                                 
      DLC = DSQRT(DLCSQ)                                                
      AC = ATERM(DLC)                                                   
      BC = BTERM(DLC)                                                   
      CC = CTERM(DLC)                                                   
      ALLOC(LWRK+4) = AC                                                
      ALLOC(LWRK+5) = BC                                                
      ALLOC(LWRK+6) = CC                                                
      FATRC = ( (CC/RC + BC)/RC + AC )/RC**3                            
C                                                                       
 300  GO TO ( 400, 500, 600, 700 ), IREQUE                              
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = IWRK                                                  
      MYINTS(2) = II                                                    
C                                                                       
C     HERE WE INDICATE THAT THIS IS A L-DEPENDENT POTENTIAL             
C                                                                       
      IRETUR = +1                                                       
C                                                                       
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  IF ( VI .EQ. 0 ) GO TO 510                                        
      R0 = RI/R0MASS                                                    
      WRITE (IUNIT, 503) VI, RI, AI, R0, POWER                          
 503  FORMAT (' W.S. WELL:', T20, G14.3, T37, F7.4,T49, F7.4,T67,F7.4, 
     * T80, F7.3 )                                                   
 510  IF ( VSI .EQ.0.0D0 ) GO TO 520                                       
      R0 = RSI/R0MASS                                                   
      WRITE (IUNIT, 513) VSI, RSI, ASI, R0                              
 513  FORMAT ( ' SURFACE ABSORPTION:', T20, G14.3, T37, F7.4,  T49, F7.4
     *, T67, F7.4, T80, F7.3 )                                          
C                                                                       
 520  WPRC5 = WP/RC**5                                                  
      WRITE (IUNIT, 523) WORDS(IPORT), BETA, G2, RC, RD, WPRC5          
 523  FORMAT ( 1X, A8, ' 2+ BKGP TELP WITH BETACOULOMB =', F7.4,  5X, 'G
     *2 =', G12.5, 5X, 'RC =', F7.4,  5X, 'RD =', F7.4, 5X, 'WP/RC**5 ='
     *, G12.5, ' MEV' )                                                 
      FFC = FATRC*RC**5                                                 
      AAC = AC*RC**2                                                    
      BBC = BC*RC                                                       
      WRITE (IUNIT, 527) DLC, FFC, AAC, BBC, CC                         
 527  FORMAT ( ' SIMPLE L-HAT CRITICAL =', F8.1,  5X, 'F(LCRIT, RC)*RC**
     *5 =', G13.5,  5X, 'TERMS =', 23G14.5 )                            
      GO TO 900                                                         
C                                                                       
C     SET UP THE TWO PIECES                                             
C                                                                       
 600  LL = LWRK+10                                                      
      CALL WOODSX ( NUMPTS, RSTART, STEPSZ, ALLOC(LL), N1, N2,1, VI,RI, 
     * AI, POWER )                                                    
      DO 619 I = 1, NUMPTS                                              
      ARRAY1(I) = ALLOC(LL-1+I)                                         
 619  CONTINUE                                                          
C                                                                       
      IF ( VSI .EQ.0.0D0 ) GO TO 650                                       
      CALL WOODSX ( NUMPTS, RSTART, STEPSZ, ALLOC(LL), N1, N2,  3, VSI, 
     * RSI, ASI, 1.D0 )                                                  
      DO 629 I = N1, N2                                                 
      ARRAY1(I) = ARRAY1(I) + ALLOC(LL-1+I)                             
 629  CONTINUE                                                          
C                                                                       
C     NOW THE BKGP TELP FOR L = LCRIT                                   
C                                                                       
 650  N1 = (RC-RSTART)/STEPSZ + 1.5D0                                     
      N1 = MIN0(N1, NUMPTS)                                             
      ALLOC(LWRK+7) = N1                                                
      ALLOC(LWRK+9) = RC                                                
      DO 679 I = 1, N1                                                  
      ARRAY1(I) = ARRAY1(I) - WPP*FATRC                                 
 679  CONTINUE                                                          
C                                                                       
      N1 = N1+1                                                         
      IF ( N1 .GT. NUMPTS ) RETURN                                      
      R = RSTART + (N1-1)*STEPSZ                                        
      DO 689 I = N1, NUMPTS                                             
      RINV = 1/R                                                        
      ARRAY1(I) = ARRAY1(I) - WPP*((CC*RINV +BC)*RINV + AC )*RINV**3                                                                  
      R = R+STEPSZ                                                      
 689  CONTINUE                                                          
      RETURN                                                            
C                                                                       
C     ADD IN THE L-DEPENDANT PART                                       
C                                                                       
C     NOTE THAT WE ARE REQUIRED TO ADD IN TO ARRAY1 THE QUANTITY        
C                                                                       
C        - H**2/(12*E) * DELTA V(L)                                     
C                                                                       
C     WHERE  DELTA V(L)  HAS ITS TRUE SIGN (NEGATIVE FOR MORE ATTRACTION
C     -H**2/(12*E) IS TO BE FOUND IN ARRAY2(1)                          
C                                                                       
 700  LWRK = KLOC(MYINTS(1))                                             
      LL = LWRK+10                                                      
      ETA = ALLOC(LWRK+2)                                               
      WPP = ALLOC(LWRK+3) * ARRAY2(1)                                   
      AC = ALLOC(LWRK+4)                                                
      BC = ALLOC(LWRK+5)                                                
      CC = ALLOC(LWRK+6)                                                
      N1 = ALLOC(LWRK+7)                                                
      AK = ALLOC(LWRK+8)                                                
      RC = ALLOC(LWRK+9)                                                
      ETASQ = ETA**2                                                    
      DL = L+.50D0                                                      
      A = ATERM(DL)                                                     
      B = BTERM(DL)                                                     
      C = CTERM(DL)                                                     
      FATRC = ( ( (C-CC)/RC + (B-BC) )/RC + (A-AC) ) / RC**3            
      DO 779 I = 1, N1                                                  
      ARRAY1(I) = ARRAY1(I) - WPP*FATRC                                 
 779  CONTINUE                                                          
      R = RSTART + N1*STEPSZ                                            
      N1 = N1+1                                                         
      IF ( N1 .GT. NUMPTS ) RETURN                                      
      DO 789 I = N1, NUMPTS                                             
      RINV = 1/R                                                        
      ARRAY1(I) = ARRAY1(I)  - WPP * ( ( (C-CC)*RINV + (B-BC) )*RINV +
     * (A-AC) ) * RINV**3                                                 
      R = R + STEPSZ                                                    
 789  CONTINUE                                                          
      RETURN                                                            
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 930  WRITE (IUNIT, 933) IPORT                                          
 933  FORMAT ( '0**** IPARAM5 =', I3, ' MUST BE 1 OR 2 TO INDICATE',  ' 
     *PROJECTILE OR TARGET EXCITATION FOR BKGPTELP.' )                  
      GO TO 955                                                         
 940  WRITE (IUNIT, 943)                                                
 943  FORMAT ( '0**** BETACOULOMB MUST BE DEFINED FOR BKGPTELP.' )      
      GO TO 955                                                         
 945  WRITE (IUNIT, 947)                                                
 947  FORMAT ( '0**** RSI (OR RSI0) AND ASI MUST BE DEFINED.' )         
      GO TO 955                                                         
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ('0**** RI (OR RI0) AND AI MUST BE DEFINED.' )             
 955  IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 960  IRETUR = -1                                                       
      RETURN                                                            
C                                                                       
      END                                                               

      SUBROUTINE DEFORM ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ,NUMPTS,ARRAY1,ARRAY2, FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,  CNSTNT, IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC,
     * IPARAM )                                                       
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1),IPARAM(5)                              
      INTEGER FACFR4                                                    
C                                                                       
C     DEFORMED - DEFORMED WOODS-SAXON OPTICAL POTENTIAL                 
C                                                                       
C     THE POTENTIAL IS DEFINED BY  V, R, A, AND POWER                   
C     THE DEFORMATION IS DEFINED BY THE ARRAY BETAP OR BETAT, WHICH     
C     CONTAIN BETAS FOR  LXMIN, LXMIN+1, LXMIN+2, ...                   
C     IPARAM4 = FIRST LX (LXMIN) FOR WHICH THE BETA ARRAYS ARE DEFINED. 
C               THE DEFAULT IS 2. IF BOTH BETAP AND BETAT ARE           
C               DEFINED, BOTH MUST START FROM THE SAME LXMIN.           
C     IF EXCITATION IS OCCURING (I.E. NOT OPTICAL POTENTIAL)            
C     THEN JBLOCK(1) IS USED TO DETERMIN WHICH NUCLEUS IS BEING         
C     EXCITED.  JJ IS MOD(JBLOCK(1), 100) - 3 AND HAS THE MEANINGS      
C          JJ = 1 FOR PROJECTILE DEFORMATION                            
C               2 FOR TARGET                                            
C               3 FOR MUTUAL DEFORMATION                                
C     FURTHERMORE, FOR MUTUAL EXCITATION THE INDIVIDUAL MULTIPOLES      
C     ARE FOUND AS                                                      
C          LXP = MOD( JBLOCK(1)/100, 100 )   AND                        
C          LXT = MOD( JBLOCK(1)/10000, 100 ) ,                          
C     IN THIS CASE LX IS IGNORED                                        
C     NCOSINE = NUMBER OF GAUSS PTS ON (0,1) (10=DEFAULT)               
C     LX IS USED AS THE DESIRED PROJECTION                              
C     LX = -1 MEANS THAT THE DEFORMED OPTICAL POT IS TO BE FOUND        
C                                                                       
C     7/27/81 - BASED ON FIXEDWOO - S.P.                                
C     9/22/81 - SET LXOUT=-1 TO 0 - S.P.                                
C     10/16/81 - ALLOW BOTH NUCLEI TO BE DEFORMED - S.P.                
C     5/14/82 - CONSIDER LX=NOTDEF TO BE LX=-1 - S.P.                   
C                                                                       
      DIMENSION BETAS(2,20), RBETAS(2,20), IBETA(2), LDEFMX(2), 
     * RDEFS(2), LXS(2)                                                         
      CHARACTER*8 DEFNAM(3)                                             
      DATA DEFNAM/'PROJECTI', 'TARGET', 'SUMMED'/                       
C                                                                       
      COMMON /CNSTNT/ PI, RT4PI, PIINV, RADIAN, DEGREE, HBARC, AMUMEV, 
     * AFINE, BIGEST, SMLEST, BIGNUM, SMLNUM, BIGLOG                      
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
      NOTDEF = INTRNL(3)                                                
      PI = CNSTNT(1)                                                    
      BIGLOG = CNSTNT(13)                                               
C                                                                       
      V = TEMPVS(IPOTYP)                                                
      VV = 0.0D0                                                            
      IF ( V .EQ.0.0D0 ) GO TO 100                                         
      A = TEMPVS(4+IPOTYP)                                              
      R = FLOATK(43)                                                     
      IF ( IPOTYP .EQ. 2 ) R = FLOATK(47)                                
      POWER = FLOATK(101+IPOTYP)                                         
      IF (R.EQ.UNDEF.OR.A .EQ.UNDEF.OR.V.EQ.UNDEF) GO TO 950
      IF ( R .LE.0.0D0 .OR. A .LE.0.0D0 ) GO TO 960                           
      LXOUT = INTGER(5)                                                 
      IF ( LXOUT .EQ. -1 .OR. LXOUT .EQ. NOTDEF ) LXOUT = 0             
      LXS(1) = 0                                                        
      LXS(2) = 0                                                        
      LXMIN = IPARAM(4)                                                 
      IF (LXMIN .EQ. NOTDEF ) LXMIN = 2                                 
      IBETA(1) = MYINTS(1)                                              
      IBETA(2) = MYINTS(2)                                              
      IF ( IREQUE .NE. 1 ) GO TO 40                                     
C                                                                       
C     FIND BETAP, BETAT                                                 
C                                                                       
      IBETA(1) = NAMLOC( 'BETAP ' )                                     
      IBETA(2) = NAMLOC( 'BETAT ' )                                     
      IF ( IBETA(1)+IBETA(2) .EQ. 0 ) GO TO 970                         
C                                                                       
C     GET THE BETAS AND COMPUTE R*BETA                                  
C                                                                       
  40  R0MASS = FINTRN(10)                                               
      DO 49 I = 1, 2                                                    
      LDEFMX(I) = 0                                                     
      IF ( IBETA(I) .EQ. 0 ) GO TO 49                                   
      LDEFMX(I) = LXMIN + LENGTH(IBETA(I))-1                            
      LXMAX = LDEFMX(I)                                                 
      RDEF = R*FLOATK(39+I)**(1.D0/3.D0) / R0MASS                        
      RDEFS(I) = RDEF                                                   
      LBETA = KLOC(IBETA(I))                                             
      DO 44 LX = LXMIN, LXMAX                                           
      BETA = ALLOC(LBETA)                                               
      BETAS(I,LX) = BETA                                                
      RBETAS(I,LX) = RDEF*BETA                                          
      LBETA = LBETA+1                                                   
 44   CONTINUE                                                          
 49   CONTINUE                                                          
      NCOSIN = INTGER(11)                                               
      IF ( NCOSIN .EQ . NOTDEF ) NCOSIN = 10                            
      IF ( NCOSIN .GT. 100 ) NCOSIN = 10                                
      VV = V                                                            
      IPORT = 0                                                         
      IF ( INTGER(5) .EQ. -1 .OR. INTGER(5) .EQ. NOTDEF ) GO TO 100     
C                                                                       
C     DETERMIN WHAT PROJECTION IS NEEDED FOR EXCITATION                 
C                                                                       
      IF ( IREQUE .NE. 3 ) GO TO 100                                    
      IPORT = MOD( JBLOCK(1), 100 ) - 2                                 
      IF ( IPORT .GE. 3 ) GO TO 80                                      
      IF ( IPORT .LT. 1 .OR. IBETA(IPORT) .EQ. 0 ) GO TO 965            
      LXS(IPORT) = LXOUT                                                
      VV = -V/RDEFS(IPORT)                                              
      GO TO 100                                                         
  80  LXS(1) = MOD( JBLOCK(1)/100, 100 )                                
      LXS(2) = MOD( JBLOCK(1)/10000, 100 )                              
      VV = -V / ( RDEFS(1)*RDEFS(2) )                                   
C                                                                       
 100  GO TO ( 200, 500, 600 ), IREQUE                                   
 200  MYINTS(1) = IBETA(1)                                              
      MYINTS(2) = IBETA(2)                                              
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  IF ( V .EQ. 0.0D0 ) GO TO 550                                         
      WRITE (IUNIT, 503) V, R, A, POWER                                 
 503  FORMAT ( ' THE POTENTIAL IS A DEFORMED WOODS-SAXON:' /  ' V =', 
     * F7.2, '   R =', F7.4, '   A = ', F7.4,  5X, 'POWER =', F8.4 )       
      DO 519 I = 1, 2                                                   
      IF ( IBETA(I) .EQ. 0 ) GO TO 519                                  
      LXMAX = LDEFMX(I)                                                 
      WRITE (IUNIT, 513) DEFNAM(I), RDEFS(I),  ( LX, BETAS(I,LX), RBETAS
     *(I,LX), LX = LXMIN, LXMAX )                                       
  513 FORMAT (' THE ', A8, ' RADIUS OF', F8.4,  ' FM. IS BEING DEFORMED
     *.' /  ' LX   BETA   DEFORMATION LENGTH' /  ( I3, F8.4, F15.4 ) )  
  519 CONTINUE                                                          
      WRITE ( IUNIT, 523 ) NCOSIN                                       
 523  FORMAT ( I5, '-POINT GAUSS INTEGRATION IS BEING USED.' )          
      GO TO 900                                                         
 550  WRITE ( IUNIT, 553 )                                              
 553  FORMAT ( ' THE POTENTIAL IS IDENTICALLY ZERO.' )                  
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
 600  CALL DEFWOO ( NUMPTS, RSTART, STEPSZ, ARRAY1, N1, N2,  1, VV, R, A
     *, POWER, IPORT, LXMIN, LDEFMX,  RBETAS, LXS, NCOSIN )             
C                                                                       
      IF ( VV .EQ.0.0D0 ) GO TO 900                                        
      VV = VV/V                                                         
      IF ( IPORT .GT. 0 ) GO TO 610                                     
      WRITE ( IUNIT, 603 ) VV, LDEFMX                                   
 603  FORMAT ( ' DEFORMED OPTICAL POTENTIAL WEIGHTED BY', G15.5,  ' WAS 
     *COMPUTED WITH LXP AND LXT MAXIMUM =', 2I3 )                       
      GO TO 680                                                         
 610  IF ( IPORT .GE. 3 ) GO TO 620                                     
      WRITE ( IUNIT, 613 ) LXOUT, VV, DEFNAM(IPORT)                     
 613  FORMAT ( ' THE LAMBDA =', I2,  ' PROJECTION OF THE DEFORMED POTENT
     *IAL WEIGHTED BY',  G15.5, ' WAS COMPUTED WITH A DEFORMED ', A8,  '
     * RADIUS.' )                                                       
      GO TO 680                                                         
 620  WRITE ( IUNIT, 623 ) VV, LXS                                      
 623  FORMAT ( ' THE MUTUAL DEFORMED POTENTIAL WEIGHTED BY',  G15.5, ' W
     *AS COMPUTED FOR LXP, LXT =', 2I3 )                                
C                                                                       
C     INDICATE THAT NO DERIVATIVE IS NEEDED IN THE FORM FACTOR          
C                                                                       
  680 IF ( IPORT .GT. 0 ) ISWTCH(30) = 1                                
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ( '0**** R, V, AND A MUST BE DEFINED.' )                   
C                                                                       
 955  IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 960  WRITE ( IUNIT, 963 ) R, A                                         
 963  FORMAT ( '0*** R OR A IS INVALID:', 2G20.10 )                     
      GO TO 955                                                         
C                                                                       
 965  WRITE ( IUNIT, 968 ) JBLOCK(1), IPORT, IREQUE, LXS, IBETA         
 968  FORMAT ( '0*** WRONG NUCLEUS DEFORMED:', I15, 8I8 )               
      GO TO 955                                                         
C                                                                       
 970  WRITE ( IUNIT, 973 )                                              
 973  FORMAT ( '0*** BETAP OR BETAT MUST BE DEFINED' )                  
      GO TO 955                                                         
      END                                                               

      SUBROUTINE FIXEDW ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ,NUMPTS,ARRAY1,ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,  CNSTNT, IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC)
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
C                                                                       
C     FIXEDWOODS - WOODS SAXON FIXED AT A POINT                         
C                                                                       
C     THE POTENTIAL IS DEFINED BY  R, A, AND POWER AND HAS DEPTH        
C     V AT RADIUS PARAM1                                                
C                                                                       
C     4/6/78 - BASED ON SHAPE - S.P.                                    
C     7/7/78 - CDC VERSION                                              
C
      INTEGER WORKNM(2,2)                                                
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                        
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
C                                                                       
      RFIX = PARAM(IPOTYP)                                              
      VFIX = TEMPVS(IPOTYP)                                             
      A = TEMPVS(4+IPOTYP)                                              
      R = FLOATK(43)                                                     
      IF ( IPOTYP .EQ. 2 ) R = FLOATK(47)                                
      POWER = FLOATK(101+IPOTYP)                                         
      IF ( R.EQ. UNDEF .OR. A .EQ. UNDEF .OR. VFIX .EQ. UNDEF.OR. 
     * RFIX .EQ. UNDEF ) GO TO 950                                          
      IF ( R.LE.0.0D0 .OR. A.LE.0.0D0.OR.RFIX.LE.0.0D0 ) GO TO 960          
      V = VFIX * (1.0D0 + DEXP((RFIX-R)/A) )**POWER                        
C                                                                       
      GO TO ( 100, 500, 600 ), IREQUE                                   
C                                                                       
C     GENERATE RADIUS ARRAY FOR INTRPC                                  
C                                                                       
 100  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = I                                                     
      MYINTS(2) = II                                                    
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503) V, R, A, POWER, VFIX, RFIX                     
 503  FORMAT ( ' THE POTENTIAL IS A WOODS-SAXON:', 3G15.5 ,  5X, 'POWER 
     *=', G14.5 /  ' WITH V =', G15.5, 5X, 'AT R =', G15.5 )            
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
 600  CALL WOODSX ( NUMPTS, RSTART, STEPSZ, ARRAY1, N1, N2,  1, V, R, A,
     * POWER )                                                          
C                                                                       
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ( '0**** R, V, A, AND PARAM1 MUST BE DEFINED.' )           
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 960  IRETUR = -1                                                       
      RETURN                                                            
C                                                                       
      END                                                               

      SUBROUTINE GAUSSI ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ, NUMPTS, ARRAY1, ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,  CNSTNT, IWAVBL,
     * FKANDM,ID,ALLOC,ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC)       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
C                                                                       
C     GAUSSIAN - GENERALIZED GAUSSIAN POTENTIAL                         
C                                                                       
C     THE POTENTIAL IS DEFINED BY  V, R, A, AND REALPOWER,  OR BY       
C     W, RI, AI, AND IMAGPOWER.  THE FORM IS                            
C                                                                       
C     POT(X) = - V EXP( ( R**(2P) - X**(2P) )/A**(2P) )                 
C                                                                       
C     WHERE P IS ____POWER.  THUS POT(X=R) IS V.  NOTE THAT V AND R     
C     ARE NOT INDEPENDANT PARAMETERS SINCE THE SAME POTENTIAL OBTAINS   
C     FOR  V EXP((R/A)**(2P))  A CONSTANT.  THE PARAMETER R IS PROVIDED 
C     TO AVOID OVERFLOW BY ALLOWING THE DEFINITION OF THE DEPTH         
C     OF THE POTENTIAL AT A PHYSICALLY MEANINGFUL POINT.                
C                                                                       
C     5/30/78 - MADED FROM FIXEDWOOD - S.P.                             
C     7/7/78 - CDC VERSION                                              
C                                                                       
      INTEGER WORKNM(2,2)                                                
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                        
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
C                                                                       
      V = TEMPVS(IPOTYP)                                                
      A = TEMPVS(4+IPOTYP)                                              
      R = FLOATK(43)                                                     
      IF ( IPOTYP .EQ. 2 ) R = FLOATK(47)                                
      POWER = FLOATK(101+IPOTYP)                                         
      IF ( R .EQ. UNDEF .OR. A .EQ. UNDEF .OR. V .EQ. UNDEF .OR. POWER 
     *.EQ. UNDEF ) GO TO 950                                            
      IF ( R.LE.0.0D0 .OR. A.LE.0.0D0 .OR. POWER.LE.0.0D0 ) GO TO 960         
C                                                                       
      GO TO ( 100, 500, 600 ), IREQUE                                   
C                                                                       
C     GENERATE RADIUS ARRAY FOR INTRPC                                  
C                                                                       
 100  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = I                                                     
      MYINTS(2) = II                                                    
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503) V, R, A, POWER                                 
 503  FORMAT ( ' THE POTENTIAL IS A GENERALIZED GAUSSIAN:' /  ' V =', G1
     *3.5, ' AT R =', F6.2,' FM;   A =', F8.3,' FM,   POWER =', F8.4) 
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
 600  TWOP = POWER+POWER                                                
C                                                                       
C     IN THE FOLLOWING  46 IS LN(1E+20)                                 
C     AND 23 IS LN(1E+10)                                               
C                                                                       
      NMAX = (R**TWOP + A**TWOP*(46+DLOG(DABS(V))) )**(1/TWOP) / STEPSZ
      NMAX = MIN0( NUMPTS, NMAX )                                       
      NMIN = 1                                                          
      TEMP = R**TWOP - A**TWOP*(23.0D0-DLOG(DABS(V)))                       
      IF ( TEMP .GT. 0.0D0 ) NMIN = TEMP**(1/TWOP) / STEPSZ                 
      NMIN = MAX0( NMIN, 1 )                                            
C                                                                       
      IF ( NMIN .EQ. 1 ) GO TO 650                                      
      DO 639 I = 1, NMIN                                                
      ARRAY1(I) = -1.D+10                                               
 639  CONTINUE                                                          
C                                                                       
 650  RVAL = (NMIN-1)*STEPSZ/A                                          
      DO 679 I = NMIN, NMAX                                             
      ARRAY1(I) = - V * DEXP( (R/A)**TWOP - RVAL**TWOP )                
      RVAL = RVAL + (STEPSZ/A)                                          
 679  CONTINUE                                                          
C                                                                       
      IF ( NMAX .GE. NUMPTS ) RETURN                                    
      NMAX = NMAX + 1                                                   
      DO 699 I = NMAX, NUMPTS                                           
      ARRAY1(I) = 0.0D0                                                     
 699  CONTINUE                                                          
C                                                                       
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ( '0**** R, V, AND A MUST BE DEFINED.' )                   
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 960  IRETUR = -1                                                       
      RETURN                                                            
C                                                                       
      END    
      
      SUBROUTINE JDEPEN ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ, NUMPTS, ARRAY1, ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,CNSTNT,IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC)                                                                
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
C                                                                       
C     JDEPEN - J-DEPENDANT DEPTH FOR POTENTIAL                          
C                                                                       
C     THE POTENTIAL IS A WOODS SAXON WITH J DEPENDANT DEPTH:            
C                                                                       
C       V = (1+PARAM1*J+PARAM2*J**2) * WOODS SAXON( V R A )             
C                                                                       
C     6/13/78 - BASED ON FIXEDWOOD - S.P.                               
C     7/11/78 - CDC VERSION AND IRETUR CORRECTION                       
C     11/16/81 - JDEPEN BASED ON PARITWOOD                              
C                                                                       
      COMMON /CNSTNT/ PI, RT4PI, PIINV, RADIAN, DEGREE, HBARC, AMUMEV, A
     *FINE, BIGEST, SMLEST, BIGNUM, SMLNUM, BIGLOG                      
C                                                                       
      INTEGER WORKNM(2,2)                                                
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                        
      LOGICAL TBIT                                                      
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
      BIGLOG = CNSTNT(13)                                               
C                                                                       
      PARAM1 = PARAM(1)                                                 
      PARAM2 = PARAM(2)                                                 
      IF ( PARAM2 .EQ. UNDEF ) PARAM2 = 0.0D0                               
      V = TEMPVS(IPOTYP)                                                
      A = TEMPVS(4+IPOTYP)                                              
      R = FLOATK(43)                                                     
      IF ( IPOTYP .EQ. 2 ) R = FLOATK(47)                                
      POWER = FLOATK(101+IPOTYP)                                         
      IF ( R .EQ. UNDEF .OR. A .EQ. UNDEF .OR. V .EQ. UNDEF .OR. PARAM1
     * .EQ. UNDEF ) GO TO 950                                           
      IF ( R .LE.0.0D0 .OR. A .LE.0.0D0 ) GO TO 960                           
C                                                                       
      GO TO ( 100, 500, 600, 700 ), IREQUE                              
C                                                                       
C     GENERATE WORK ARRAY TO SAVE THE W.S.                              
C                                                                       
 100  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
C                                                                       
      I = NALLOC( NUMPTS, WORKNM(1,1) )                                 
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = I                                                     
      MYINTS(2) = II                                                    
C                                                                       
C     HERE WE INDICATE THAT THIS IS A L-DEPENDENT POTENTIAL             
C                                                                       
      IRETUR = +1                                                       
C                                                                       
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503) V, R, A, POWER, PARAM1, PARAM2                 
 503  FORMAT ( ' THE POTENTIAL IS A J-DEPENDANT WOODS-SAXON:', 3G15.5 , 
     * 5X, 'POWER =', G14.5 /  ' WITH THE J-DEPENDANT DEPTH FACTOR',  ' 
     * 1 + J *', G13.5, ' + J**2 *', G13.5 )                            
      GO TO 900                                                         
C                                                                       
C     SET UP THE TWO PIECES                                             
C                                                                       
 600  LL = KLOC(MYINTS(1))                                               
      CALL WOODSX ( NUMPTS, RSTART, STEPSZ, ALLOC(LL), N1, N2,  1, V, R,
     * A, POWER )                                                       
      DO 619 I = 1, NUMPTS                                              
      ARRAY1(I) = ALLOC(LL-1+I)                                         
 619  CONTINUE                                                          
C                                                                       
      RETURN                                                            
C                                                                       
C     ADD IN THE J-DEPENDANT PART                                       
C                                                                       
 700  LL = KLOC(MYINTS(1))                                               
      DJ = 0.5D0*JBLOCK(1)                                                 
      FAC = ARRAY2(1)*DJ*(PARAM1 + DJ*PARAM2)                           
      N1 = RSTART/STEPSZ + 0.5D0                                           
      LL = LL - 1 + N1                                                  
      DO 719 I = 1, NUMPTS                                              
      ARRAY1(I) = ARRAY1(I) + FAC*ALLOC(LL+I)                           
 719  CONTINUE                                                          
CCC      WRITE ( IUNIT, 723 )  DJ, FAC, RSTART, N1, NUMPTS              
CCC  723 FORMAT ( ' J-DEPENDANT DEPTH:  J, FAC =', 2G15.8,              
CCC     1  5X, 'RSTART, N1 =', G15.8, 2I6 )                             
C                                                                       
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ( '0**** R, V, A, AND PARAM1 MUST BE DEFINED.' )           
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 960  IRETUR = -1                                                       
      RETURN                                                            
C                                                                       
      END                                                               

      SUBROUTINE JDEPE1 ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ, NUMPTS, ARRAY1, ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,CNSTNT,IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC)
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
C                                                                       
C     JDEPENWS - J-DEPENDANT W.S. DEPTH FOR POTENTIAL                   
C                                                                       
C     THE POTENTIAL IS A WOODS SAXON WITH J-DEPENDANT W.S. DEPTH:       
C                                                                       
C       V = F(J) * WOODS SAXON( V R A )                                 
C       F(J) = 1 / ( 1 + EXP( (J-PARM100)/PARM101 ) )                   
C                                                                       
C     6/13/78 - BASED ON FIXEDWOOD - S.P.                               
C     7/11/78 - CDC VERSION AND IRETUR CORRECTION                       
C     11/16/81 - JDEPEN BASED ON PARITWOOD                              
C     5/11/82 - JDEPENWS MADE FROM JDEPEN - S.P.                        
C                                                                       
      COMMON /CNSTNT/ PI, RT4PI, PIINV, RADIAN, DEGREE, HBARC, AMUMEV, A
     *FINE, BIGEST, SMLEST, BIGNUM, SMLNUM, BIGLOG                      
C                                                                       
      INTEGER WORKNM(2,2)                                                
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                        
      LOGICAL TBIT                                                      
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
      BIGLOG = CNSTNT(13)                                               
C                                                                       
      PARM10 = PARAM(10)                                                
      PARM11 = PARAM(11)                                                
      V = TEMPVS(IPOTYP)                                                
      A = TEMPVS(4+IPOTYP)                                              
      R = FLOATK(43)                                                     
      IF ( IPOTYP .EQ. 2 ) R = FLOATK(47)                                
      POWER = FLOATK(101+IPOTYP)                                         
      IF ( R .EQ. UNDEF .OR. A .EQ. UNDEF .OR. V .EQ. UNDEF .OR. PARM10
     * .EQ. UNDEF .OR. PARM11 .EQ. UNDEF ) GO TO 950                    
      IF ( R .LE. 0.0D0 .OR. A .LE. 0.0D0 ) GO TO 960                           
C                                                                       
      GO TO ( 100, 500, 600, 700 ), IREQUE                              
C                                                                       
C     GENERATE WORK ARRAY TO SAVE THE W.S.                              
C                                                                       
 100  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
C                                                                       
      I = NALLOC( NUMPTS, WORKNM(1,1) )                                 
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = I                                                     
      MYINTS(2) = II                                                    
C                                                                       
C     HERE WE INDICATE THAT THIS IS A L-DEPENDENT POTENTIAL             
C                                                                       
      IRETUR = +1                                                       
C                                                                       
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503) V, R, A, POWER, PARM10, PARM11                 
 503  FORMAT ( ' THE POTENTIAL IS A J-DEPENDANT WOODS-SAXON:', 3G15.5 , 
     * 5X, 'POWER =', G14.5 /  ' WITH THE J-DEPENDANT DEPTH FACTOR',  ' 
     * 1 / ( 1 + EXP( ( J - ', G12.5, ') / ', G12.5, ') )' )            
      GO TO 900                                                         
C                                                                       
C     SET UP THE TWO PIECES                                             
C                                                                       
 600  LL = KLOC(MYINTS(1))                                               
      CALL WOODSX ( NUMPTS, RSTART, STEPSZ, ALLOC(LL), N1, N2, 1, V, R,
     * A, POWER )                                                       
      DO 619 I = 1, NUMPTS                                              
      ARRAY1(I) = ALLOC(LL-1+I)                                         
 619  CONTINUE                                                          
C                                                                       
      RETURN                                                            
C                                                                       
C     ADD IN THE J-DEPENDANT PART                                       
C                                                                       
 700  LL = KLOC(MYINTS(1))                                               
      DJ = 0.5D0*JBLOCK(1)                                                 
      FAC = 1.0D0                                                           
      X = (DJ-PARM10)/PARM11                                            
      IF ( X .LT. -30.0D0 ) GO TO 710                                       
      FAC = 0.0D0                                                           
      IF ( X .GT. 30.0D0 ) GO TO 710                                        
      FAC = 1.0D0/(1.0D0+DEXP(X))                                               
 710  FAC = ARRAY2(1) * ( FAC - 1.0D0 )                                     
      N1 = RSTART/STEPSZ + 0.5D0                                           
      LL = LL - 1 + N1                                                  
      DO 719 I = 1, NUMPTS                                              
      ARRAY1(I) = ARRAY1(I) + FAC*ALLOC(LL+I)                           
 719  CONTINUE                                                          
CCC      WRITE ( IUNIT, 723 )  DJ, FAC, RSTART, N1, NUMPTS              
CCC  723 FORMAT ( ' J-DEPENDANT DEPTH:  J, FAC =', 2G15.8,              
CCC     1  5X, 'RSTART, N1 =', G15.8, 2I6 )                             
C                                                                       
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ( '0**** R, V, A, PARAM10 & PARAM11 MUST BE DEFINED.' )    
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 960  IRETUR = -1                                                       
      RETURN                                                            
C                                                                       
      END                                                               

      SUBROUTINE LAGRAN ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J, RSTART, STEPSZ,NUMPTS, ARRAY1, ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,CNSTNT,IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC)                                                                
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(2,10),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
C                                                                       
C     LAGRANGE - CONSTRUCT POTENTIAL BY LAGRANGE INTERPOLATION          
C                                                                       
C     THIS LINKULE ASSUMES THAT THE PAIRS OF NUMBERS:                   
C       (PARAM1, PARAM2), (PARAM3, PARAM4), ..., (PARAM19, PARAM20)     
C     ARE TO BE INTERPERTED AS                                          
C       (R1, V1), (R2, V2), ..., (R10, V10)                             
C     AND DEFINE A POTENTIAL V(R).  FEWER THAN 10 POINTS MAY            
C     BE USED BY LEAVING THE REMAINING V'S OR R'S UNDEFINED.            
C     THE R'S DO NOT HAVE TO BE ORDERED.                                
C     THE POTENTIAL IS DEFINED AS (NOTE SIGN CHANGE)                    
C                                                                       
C       V(R) =  V1    R < R1 ,                                          
C               V'(R)   R1 < R < RLAST                                  
C               -A*EXP(-B*R) R > RLAST                                  
C            WHERE A AND B ARE FIT TO THE LAST POINTS.                  
C                                                                       
C     WHERE V'(R) IS THE RESULT OF LAGRANGE INTERPOLATION OR ORDER IPARA
C                                                                       
C     IF THE POTENTIAL DOES NOT CHANGE SIGN, THEN THE LAGRANGE INTERPOLA
C     IS APPLIED TO LOG!V!.  IF IT DOES CHANGE SIGN, THEN LET N BE      
C     THE LARGEST INDEX SUCH THAT  SIGN(V(N-1)) ^= SIGN(V(N)).          
C     LAGRANGE INTERPOLATION IS APPLIED TO V FOR  R < R(N)  AND TO      
C     LOG!V! FOR R > R(N).                                              
C                                                                       
C     4/10/78 - BASED ON SPLINE - S.P.                                  
C     2/28/79 - ALLOW CHANGES OF SIGN - S.P.                            
C     3/14/79 - CDC VERSION - S.P.                                      
C                                                                       
      INTEGER WORKNM(2,2)                                                
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                        
      DIMENSION R(20), V(20), A(30)                                     
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
C                                                                       
C     INTERPOLATION ORDER                                               
C                                                                       
      IORDER = INTGER(38)                                               
CIB      IF ( IORDER .EQ. INTRNL(3) )  GO TO 960                        
CDC      IF ( IORDER .EQ. INTRNL(2) )  GO TO 960                        
C                                                                       
C     GET THE DEFINED POINTS AND VALUES AND ORDER THEM                  
C                                                                       
      NV = 1                                                            
      IV = 1                                                            
  40  R(1) = PARAM(1,IV)                                                
      V(1) = PARAM(2,IV)                                                
      IF ( R(1) .NE. UNDEF .AND. V(1) .NE. UNDEF ) GO TO 50             
      IV = IV+1                                                         
      IF ( IV .LT. 10 ) GO TO 40                                        
      GO TO 950                                                         
C                                                                       
 50   IV = IV + 1                                                       
      IF ( IV .GT. 10 ) GO TO 90                                        
      IF (PARAM(1,IV).EQ.UNDEF.OR.PARAM(2,IV).EQ.UNDEF )  GO TO 50                                                                 
      DO 59 I = 1, NV                                                   
      IF ( PARAM(1,IV) .LT. R(I) ) GO TO 70                             
 59   CONTINUE                                                          
C                                                                       
C     ADD TO END                                                        
C                                                                       
      NV = NV+1                                                         
      R(NV) = PARAM(1,IV)                                               
      V(NV) = PARAM(2,IV)                                               
      GO TO 50                                                          
C                                                                       
C     INSERT POINT; FIRST MOVE REST RIGHT ONE                           
C                                                                       
 70   DO 74 II = I, NV                                                  
      R(NV+I-II+1) = R(NV+I-II)                                         
      V(NV+I-II+1) = V(NV+I-II)                                         
 74   CONTINUE                                                          
      NV = NV+1                                                         
      R(I) = PARAM(1,IV)                                                
      V(I) = PARAM(2,IV)                                                
      GO TO 50                                                          
C                                                                       
C     ALL DEFINITIONS FOUND AND ORDERED                                 
C                                                                       
 90   IF ( NV .LE. 2 ) GO TO 950                                        
C                                                                       
      GO TO ( 100, 500, 600 ), IREQUE                                   
C                                                                       
C     GENERATE RADIUS ARRAY FOR INTRPC                                  
C                                                                       
 100  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
      I = NALLOC( NUMPTS, WORKNM(1,1), KLOC )                            
C                                                                       
      RR = RSTART                                                       
      DO 159 II = 1, NUMPTS                                             
      ALLOC(KLOC(I)-1+II) = RR                                           
      RR = RR + STEPSZ                                                  
 159  CONTINUE                                                          
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = I                                                     
      MYINTS(2) = II                                                    
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503) IORDER, ( R(I), V(I), I = 1, NV )              
 503  FORMAT ( ' THE POTENTIAL IS DEFINED BY ',  'LAGRANGE INTERPOLATION
     * OF ORDER',I2,' BASED ON THE POINTS:' / (5(F10.3,G13.5))) 
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
 600  LL = KLOC(MYINTS(1)) - 1                                           
C                                                                       
C     SETUP OUTPUT X'S, DO NOT EXTRAPOLATE; RATHER USE                  
C     F(1) FOR R < R1  AND  EXPONENTIAL DECAY FOR R > R2 .              
C                                                                       
 270  RR = RSTART                                                       
      NSTRT = 1                                                         
      DO 299 N = 1, NUMPTS                                              
      IF ( RR .GE. R(1) ) GO TO 280                                     
      ARRAY1(N) = -V(1)                                                 
      NSTRT = N+1                                                       
      GO TO 290                                                         
 280  IF ( RR .LE. R(NV) ) GO TO 285                                    
      ARRAY1(N) = 0                                                     
      GO TO 290                                                         
 285  NLAST = N                                                         
 290  RR = RR + STEPSZ                                                  
 299  CONTINUE                                                          
      NOUT = NLAST+1-NSTRT                                              
C                                                                       
C     DOES THE POTENTIAL CHANGE SIGN                                    
C                                                                       
      LNORDR = IORDER                                                   
      VSIGN = -DSIGN(1.0D0, V(NV) )                                     
      LNSTRT = 1                                                        
      DO 659 II = 2, NV                                                 
      I = NV+2-II                                                       
      IF ( VSIGN*V(I-1) .GE.0.0D0 ) GO TO 660                              
 659  CONTINUE                                                          
      GO TO 700                                                         
C                                                                       
C     YES, USE STRAIGHT INTERPOLATION UP TO THE LAST CHANGE.            
C     FIND THIS TRANSITION POINT IN THE OUTPUT ARRAYS.                 
C                                                                       
 660  LNSTRT = I                                                        
      NSTRT2 = ( R(LNSTRT)-RSTART )/STEPSZ + 2                          
      CALL AITKEN ( MIN0(IORDER,LNSTRT-1), 0.D0,0.D0, LNSTRT, R, V,
     * NSTRT2-NSTRT, ALLOC(LL+NSTRT), ARRAY1(NSTRT), NFAIL, WORST, A)    
      DO 679 I = NSTRT, NSTRT2                                          
      ARRAY1(I) = -ARRAY1(I)                                            
 679  CONTINUE                                                          
C                                                                       
      LNORDR = MIN0( IORDER, NV-LNSTRT )                                
      IF ( LNORDR .LE. 0 ) RETURN                                       
      NSTRT = NSTRT2                                                    
      NOUT = NLAST+1-NSTRT2                                             
C                                                                       
C     INTERPOLATE THE LOGS                                              
C                                                                       
 700  DO 719 I = LNSTRT, NV                                             
      V(I) = DLOG( DMAX1( DABS(V(I)), 1.D-15 ) )                        
 719  CONTINUE                                                          
C                                                                       
      CALL AITKEN ( LNORDR,0.D0,0.D0,NV-LNSTRT+1,R(LNSTRT),V(LNSTRT), 
     * NOUT, ALLOC(LL+NSTRT), ARRAY1(NSTRT),  NFAIL, WORST, A )     
C                                                                       
      DO 749 I = NSTRT, NLAST                                           
      X = ARRAY1(I)                                                     
      X = DMAX1( X, -69.D0 )                                            
      X = DMIN1( X, 69.D0 )                                             
      ARRAY1(I) = VSIGN*DEXP(X)                                         
 749  CONTINUE                                                          
C                                                                       
C     FILL IN THE END WITH AN EXPONENTIAL DECAY                         
C                                                                       
      TERM = ARRAY1(NLAST)                                              
      NLAST = NLAST+1                                                   
      IF ( NLAST .GT. NUMPTS ) RETURN                                   
      STEP = TERM/ARRAY1(NLAST-2)                                       
      TERM = TERM*STEP                                                  
      IF ( 1.0D0-STEP .LT. 1.D-5 ) RETURN                                   
      DO 769 I = NLAST, NUMPTS                                          
      ARRAY1(I) = TERM                                                  
      IF ( TERM .GT. -1.D-30 ) RETURN                                   
      TERM = TERM*STEP                                                  
 769  CONTINUE                                                          
C                                                                       
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ( '0**** POTENTIAL MUST BE DEFINED ON AT LEAST',  ' TWO POI
     *NTS FOR LAGRANGE.' )                                              
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 960  WRITE (IUNIT, 963)                                                
 963  FORMAT ( ' IPARAM1 MUST BE DEFINED AS ORDER OF INTERPOLATION',  ' 
     *FOR LAGRANGE.' )                                                  
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
      END  
      
      SUBROUTINE LTSTEL ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J, RSTART, STEPSZ, NUMPTS,ARRAY1,ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,CNSTNT, IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC,
     * IPARAM )                                                       
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1),IPARAM(5)                              
      INTEGER FACFR4                                                    
      CHARACTER*4 ID                                                    
C                                                                       
C     LOVE, TERASAWA AND SATCHLER  2+ TELP                              
C                                                                       
C     COMPUTES THE IMAGINARY POTENTIAL AS                               
C                                                                       
C        VI, RI, AI, IMAPOWER   WOODS-SAXON                             
C      + VSI, RSI, ASI          SURFACE TERM                            
C      + LTS TELP FOR A 2+ STATE.                                       
C                                                                       
C     THE TELP PART IS CONTROLLED BY                                    
C                                                                       
C     IPARAM5 - 1 = PROJECTILE EXCITATION                               
C               2 = TARGET EXCITATION                                   
C               IF NOT DEFINED, THEN  JA @= JB  OR  JBIGA @= JBIGB      
C          IS USED TO DETERMIN WHICH NUCLEUS IS EXCITED.                
C     PARAM20 - CONTAINS THE G2(ETA,XI) FUDGE FACTOR.  IF NOT DEFINED,  
C               THE VALUE 1 IS USED.                                    
C     BETACOULOMB                                                       
C                                                                       
C     WE STORE SOME CONSTANTS IN THE BEGINNING OF THE WORK AREA         
C     ALLOC(KLOC(IWRK)+I):                                               
C     I   STORED                                                        
C     0   R2FAC   R2 = R2FAC*RC                                         
C     1   BETACOULOMB                                                   
C                                                                       
C     3/11/79 - MADE FROM PARITWOOD - S.P.                              
C     8/1/79 - IPARAM IS A SEPARATE ARRAY - S.P.                        
C     1/18/80 - STANDARD MORTRAN - RPG                                  
C                                                                       
      CHARACTER*4 WORKNM(2,2)                                           
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                       
CCC      LOGICAL  TBIT                                                  
      CHARACTER*8 WORDS(2)                                              
      DATA WORDS/'PROJCTIL', '  TARGET'/                                
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
      NOTDEF = INTRNL(3)                                                
C                                                                       
      VI = TEMPVS(2)                                                    
      AI = TEMPVS(6)                                                    
      RI = FLOATK(47)                                                    
      POWER = FLOATK(103)                                                
      VSI = FLOATK(71)                                                   
      RSI = FLOATK(69)                                                   
      ASI = FLOATK(68)                                                   
      G2 = PARAM(20)                                                    
      IF ( G2 .EQ. UNDEF ) G2 = 1                                       
      ICHANW = INTRNL(5)                                                
      AK = FKANDM(ICHANW)                                               
      RC = FLOATK(45)                                                    
      ZP = INTGER(24)                                                   
      ZT = INTGER(25)                                                   
      ETA = FKANDM(ICHANW+16)                                           
      ALPHA = 1/CNSTNT(8)                                               
      PI = CNSTNT(1)                                                    
      HBARC = CNSTNT(6)                                                 
      E = FLOATK(12)                                                     
      AM = FLOATK(24)                                                    
C                                                                       
      IF ( IREQUE .GT. 1 ) GO TO 200                                    
C                                                                       
C     IS THIS PROJECTILE OR TARGET EXCITATION.                          
C     WE USE THIS INFORMATION TO CONSTRUCT R2FAC WHICH MAKES THE        
C     SINGLE NUCLEON RADIUS OUT OF THE SUM (RC).                        
C                                                                       
      IPORT = IPARAM(5)                                                 
      IF ( IPORT .NE. NOTDEF ) GO TO 50                                 
      IF ( JBLOCK(2) .NE. JBLOCK(3) ) IPORT = 1                         
      IF ( JBLOCK(4) .NE. JBLOCK(5) ) IPORT = 2                         
 50   IF ( IPORT .NE. 1 .AND. IPORT .NE. 2 ) GO TO 930                  
C      AMP = FLOATK(40)                                                  
C      AMT = FLOATK(41)                                                  
      R0MASS = FINTRN(10)                                               
      R2FAC = FLOATK(39+IPORT)**(1.D0/3.D0) / R0MASS                     
C                                                                       
      IB = NAMLOC('BETACOUL')                                           
      IF ( IB .EQ. 0 ) GO TO 940                                        
C                                                                       
      IF ( VI .EQ. 0.0D0 ) GO TO 65                                         
      IF ( RI .EQ. UNDEF .OR. AI .EQ. UNDEF ) GO TO 950                 
  65  IF ( VSI .EQ. 0.0D0 ) GO TO 70                                        
      IF ( RSI .EQ. UNDEF .OR. ASI .EQ. UNDEF ) GO TO 945               
C                                                                       
C     GENERATE WORK ARRAY TO SAVE THE W.S.                              
C                                                                       
  70  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
C                                                                       
      IWRK = NALLOC( NUMPTS+10, WORKNM(1,1) )                           
      LWRK = KLOC(IWRK)                                                  
      ALLOC(LWRK) = R2FAC                                               
      ALLOC(LWRK+1) = ALLOC(KLOC(IB))                                    
      GO TO 300                                                         
C                                                                       
C     COMPUTE THE CONSTANTS                                             
C     THIS PART IS DONE FOR IREQUE = 2 (PRINTING) AND 3 (NEW POTENTIAL) 
C                                                                       
 200  IWRK = MYINTS(1)                                                  
      LWRK = KLOC(IWRK)                                                  
      IF ( VI.NE.0.0D0.AND.(AI.LE.0.0D0.OR.RI.LE.0.0D0))  GO TO 960    
      IF ( VSI.NE.0.0D0 .AND. (ASI.LE.0.0D0.OR.RSI.LE.0.0D0)) GO TO 960 
      RC2 = RC*ALLOC(LWRK)                                              
      BETA = ALLOC(LWRK+1)                                              
      RD = ZP*ZT*ALPHA*HBARC/E                                          
      WP = 3.0D0*(ZP*ZT*ALPHA*BETA)**2 *AM*G2*RC2**4/(50.0D0*PI*AK)  
C                                                                       
 300  GO TO ( 400, 500, 600, 700 ), IREQUE                              
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = IWRK                                                  
      MYINTS(2) = II                                                    
CCCC                                                                    
CCCC     HERE WE INDICATE THAT THIS IS A L-DEPENDENT POTENTIAL          
CCCC                                                                    
CCC      IRETUR = +1                                                    
      IRETUR = 0                                                        
C                                                                       
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  IF ( VI .EQ. 0.0D0 ) GO TO 510                                        
      R0 = RI/R0MASS                                                    
      WRITE (IUNIT, 503) VI, RI, AI, R0, POWER                          
 503  FORMAT ( ' W.S. WELL:', T20, G14.3, T37, F7.4,  T49, F7.4, T67, F7
     *.4, T80, F7.3 )                                                   
 510  IF ( VSI .EQ. 0.0D0 ) GO TO 520                                       
      R0 = RSI/R0MASS                                                   
      WRITE (IUNIT, 513) VSI, RSI, ASI, R0                              
 513  FORMAT ( ' SURFACE ABSORPTION:', T20, G14.3, T37, F7.4,  T49, F7.4
     *, T67, F7.4, T80, F7.3 )                                          
C                                                                       
 520  WPRC5 = WP/RC**5                                                  
      WRITE (IUNIT, 523) WORDS(IPORT), BETA, G2, RC, RD, WPRC5          
 523  FORMAT ( 1X, A8, ' 2+ LTS TELP WITH BETACOULOMB =', F7.4,  5X, 'G2
     * =', G12.5, 5X, 'RC =', F7.4,  5X, 'RD =', F7.4, 5X, 'WP/RC**5 =',
     * G12.5, ' MEV' )                                                  
      GO TO 900                                                         
C                                                                       
C     SET UP THE TWO PIECES                                             
C                                                                       
 600  LL = LWRK+10                                                      
      CALL WOODSX ( NUMPTS, RSTART, STEPSZ, ALLOC(LL), N1, N2,1,VI,RI, 
     * AI, POWER )                                                    
      DO 619 I = 1, NUMPTS                                              
      ARRAY1(I) = ALLOC(LL-1+I)                                         
 619  CONTINUE                                                          
C                                                                       
      IF ( VSI .EQ. 0 ) GO TO 650                                       
      CALL WOODSX ( NUMPTS, RSTART, STEPSZ, ALLOC(LL), N1, N2,  3, VSI, 
     *RSI, ASI, 1.D0 )                                                  
      DO 629 I = N1, N2                                                 
      ARRAY1(I) = ARRAY1(I) + ALLOC(LL-1+I)                             
 629  CONTINUE                                                          
C                                                                       
C     NOW THE LTS TELP                                                  
C      FIRST THE "COULOMB BRAKING TERM" WHICH IS ARBITRARILY TRUNCATED  
C                                                                       
 650  X = 1.0D0/DSQRT(.10D0)                                                
      N1 = (RD-RSTART)/(.90D0*STEPSZ) + 1.5D0                             
      N1 = MIN0(N1, NUMPTS)                                             
      DO 654 I = 1, N1                                                  
      ALLOC(LL-1+I) = X                                                 
 654  CONTINUE                                                          
      R = N1*STEPSZ + RSTART                                            
      DO 659 I = N1, NUMPTS                                             
      ALLOC(LL-1+I) = 1.0D0/DSQRT(1-RD/R)                                   
      R = R+STEPSZ                                                      
 659  CONTINUE                                                          
C                                                                       
C     NOW THE LTS TELP TERM WHICH HAS TWO REGIONS OF DEFINITION         
C                                                                       
      N1 = (RC-RSTART)/STEPSZ + 1.5D0                                     
      N1 = MIN0(N1, NUMPTS)                                             
      R = RSTART                                                        
      DO 679 I = 1, N1                                                  
      ARRAY1(I) = ARRAY1(I)-(2.0D0*WP/(3.0D0*RC**9))*ALLOC(LL-1+I)*R**4  
      R = R+STEPSZ                                                      
 679  CONTINUE                                                          
C                                                                       
      N1 = N1+1                                                         
      IF ( N1 .GT. NUMPTS ) RETURN                                      
      DO 689 I = N1, NUMPTS                                             
      ARRAY1(I) = ARRAY1(I)-WP*ALLOC(LL-1+I)*(1.0D0-(RC/R)**2*(2.0D0/
     * 7.D0  + (RC/R)**2 * (1.D0/21.D0) )) / R**5                     
      R = R+STEPSZ                                                      
 689  CONTINUE                                                          
      RETURN                                                            
C                                                                       
C     ADD IN THE L-DEPENDANT PART                                       
C                                                                       
 700  LWRK = KLOC(MYINTS(1))                                             
      LL = LWRK+10                                                      
CCC      FAC = ARRAY2(1)*PCOEF                                          
CCC      IF ( TBIT(L,31) )  FAC = -FAC                                  
CCC      DO 719  I = 1, NUMPTS                                          
CCC         ARRAY1(I) = ARRAY1(I) + FAC*ALLOC(LL-1+I)                   
CCC 719  CONTINUE                                                       
      RETURN                                                            
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 930  WRITE (IUNIT, 933) IPORT                                          
 933  FORMAT ( '0**** IPARAM5 =', I3, ' MUST BE 1 OR 2 TO INDICATE',  ' 
     *PROJECTILE OR TARGET EXCITATION FOR LTSTELP.' )                   
      GO TO 955                                                         
 940  WRITE (IUNIT, 943)                                                
 943  FORMAT ( '0**** BETACOULOMB MUST BE DEFINED FOR LTSTELP.' )       
      GO TO 955                                                         
 945  WRITE (IUNIT, 947)                                                
 947  FORMAT ( '0**** RSI (OR RSI0) AND ASI MUST BE DEFINED.' )         
      GO TO 955                                                         
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ('0**** RI (OR RI0) AND AI MUST BE DEFINED.' )             
 955  IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 960  IRETUR = -1                                                       
      RETURN                                                            
C                                                                       
      END 
      
      SUBROUTINE OHTA ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT,  
     * NUMOUT,L,J,RSTART, STEPSZ, NUMPTS, WAVR, WAVI, FLOATK, TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,  CNSTNT, IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC )    
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), WAVR(NUMPTS), WAVI(NUMPTS), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
C                                                                       
C     LINKULE FOR INCOMING WAVES WITH FIXED LN DERIV (TIMES K)          
C                                                                       
C     2/15/78 - FIRST VERSION                                           
C     6/1/78 - CHANGES FOR NEW WAVLJ STRUCTURE - S.P.                   
C     7/7/78 - CDC VERSION                                              
C     1/16/80 - STANDARD MORTRAN, USE SMLNUM - RPG                      
C                                                                       
      LOGICAL PBUGSW                                                    
      EQUIVALENCE ( IBUGSW, PBUGSW )                                    
C                                                                       
      INTEGER WORKNM(2)                                                  
      DATA WORKNM/'*000', 'WORK'/                                        
      DIMENSION UDERIV(2, 4), FDERIV(2, 6), BC(5),  FTERMS(2, 6, 2), 
     * FSUMS(2, 2)                                                          
      DIMENSION XS(22), BS(22), CS(22), DS(22)                          
CDC      COMMON /OHTABL/ XS, BS, CS, DS                                 
CDC      ALLOCLEVEL XS, BS, CS, DS                                      
C                                                                       
      DATA TWELTH / .0833333333333333333333 0D0 /                       
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      IBUGSW = IWAVBL(114)                                              
CDC      IBUGSW = IWAVBL(76)                                            
      IBUGSW = IWAVBL(114)                                              
      BIGNUM = CNSTNT(11)                                               
      SMLNUM = CNSTNT(12)                                               
C                                                                       
      GO TO ( 100, 500, 600 ), IREQUE                                   
C                                                                       
C     LOCATE THE BOUNDRY RADIUS (PARAM1) AND DERIVATIVE (PARAM2)        
C                                                                       
 100  IF ( PARAM(1).NE.FINTRN(1).AND.PARAM(2).NE.FINTRN(1)) GO TO 120                                                             
      WRITE ( IUNIT, 113 )                                              
 113  FORMAT ( '0**** BOTH PARAM1 AND PARAM2 MUST BE DEFINED',  ' AS THE
     * BOUNDARY RADIUS AND DERIVATIVE.' )                               
      IRETUR = -1                                                       
      GO TO 900                                                         
 120  WORKNM(1) = ID                                                    
CCCC      IWRK = NALLOC( 100,WORKNM)                                    
      UROLD = -1.D-30                                                   
      UIOLD = -1.D-30                                                   
C                                                                       
C     SAVE ADDRESS OF WORK ARRAY                                        
C                                                                       
      MYINTS(1) = IWRK                                                  
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  DERIV2 = 0.0D0                                                        
      IF ( PARAM(3) .NE. FINTRN(1) ) DERIV2 = PARAM(3)                  
      WRITE (IUNIT, 503) PARAM(1), PARAM(2), DERIV2, FKANDM(ID)         
 503  FORMAT ( ' THE WAVEFUNCTION WILL HAVE AN INCOMING',  ' BOUNDARY CO
     *NDITION AT R =', F6.2,  ' FM' / '  WITH A LOG DERIVATIVE OF -I * (
     *',  G13.6, ' +', G13.6, 'I ) * K',  '  (K =', F8.3, ' FM**-1) ',  
     *' FOR ALL VALUES OF L.' )                                         
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
C     GET SOME ARGUMENTS INTO THEIR CONVENTIONAL NAMES                  
C                                                                       
 600  NFIRST = IPOTYP                                                   
      H = STEPSZ                                                        
      ISTRT = NFIRST-2                                                  
      NSTEP = NUMPTS-1                                                  
      NWP = ID                                                          
CCCC      IWRK = MYINTS(1)                                              
CCCC      LWRK = LOC(IWRK)                                              
      STPSZ = H/FKANDM(NWP)                                             
      RB = PARAM(1)                                                     
      DERIV1 = PARAM(2)                                                 
      DERIV2 = 0.0D0                                                        
      IF ( PARAM(3) .NE. FINTRN(1) ) DERIV2 = PARAM(3)                  
      DERIV1 = FKANDM(NWP) * DERIV1                                     
      DERIV2 = FKANDM(NWP) * DERIV2                                     
C                                                                       
C     NOW GET STARTING VALUES OF THE POTENTIAL AND SAVE STARTING        
C     VALUES OF THE WAVEFUNCTION.                                       
C                                                                       
      N1 = RB/STPSZ + 0.2D0                                                
      R1 = N1*STPSZ                                                     
      X = RB - R1                                                       
      IF ( N1 .GT. ISTRT ) GO TO 250                                    
C                                                                       
C     THE BOUNDRY IS INSIDE THE FIRST POINT SO IT DOESNT COUNT          
C                                                                       
      GO TO 350                                                         
C                                                                       
C     WE WILL START AT N1                                               
C                                                                       
 250  ABACK = WAVR(N1+4)                                                
      BBACK = WAVI(N1+4)                                                
C                                                                       
C     IT MAY BE POSSIBLE TO JUST CORRECT THE PREVIOUS DERIVATIVES       
C     AND THUS SAVE THE CUBIC SPLINE TIME                               
C                                                                       
      IF ( BBACK .NE. UIOLD ) GO TO 255                                 
      DEL = L*(L+1) - LOLD*(LOLD+1)                                     
      IF (DABS(UROLD+DEL/(12.D0*N1*N1)-ABACK).GT.1.D-12)GO TO 255
C     YES, JUST ADJUST IT                                               
C                                                                       
      DEL = -DEL/RB                                                     
      DO 252 I = 1, 4                                                   
      DEL = -I*DEL/RB                                                   
      UDERIV(1,I) = UDERIV(1,I) + DEL                                   
 252  CONTINUE                                                          
      GO TO 270                                                         
C                                                                       
C     GET THE DERIVATIVES OF THE POTENTIAL FROM CUBIC SPLINES           
C                                                                       
 255  I1 = MAX0( N1-10, ISTRT+1 )                                       
      DO 254 I = 1, 21                                                  
      XS(I) = I*STPSZ                                                   
 254  CONTINUE                                                          
      I = N1-I1+1                                                       
      FACTOR = 12.0D0/STPSZ**2                                              
      IF ( PBUGSW ) WRITE ( IUNIT, 263 ) ISTRT, N1, R1, X               
 263  FORMAT ( ' BARRIER ISTRT, N1, R1, X =', 2I6, 2G15.5 )             
      IF (PBUGSW) NUMOUT = 1                                            
C                                                                       
C     FIRST THE REAL POTENTIAL DERIVATIVES                              
C                                                                       
      CALL SPLNCB ( 21, XS, WAVR(I1+4), BS, CS, DS )                    
      A = WAVR(N1+4)                                                    
      B = BS(I)                                                         
      C = CS(I)                                                         
      D = DS(I)                                                         
      UDERIV(1,1) = FACTOR*( A-1.0D0 + X*(B + X*(C + X*D )))                
      UDERIV(1,2) = FACTOR*( B + X*(2.0D0*C + 3.0D0*D*X))                       
      UDERIV(1,3) = FACTOR* 2.0D0*( C + 3.0D0*D*X )                             
      UDERIV(1,4) = FACTOR * 6.0D0*D                                        
C                                                                       
C     NOW THE IMAG POTENTIAL DERIVATIVES                                
C                                                                       
      CALL SPLNCB ( 21, XS, WAVI(I1+4), BS, CS, DS )                    
      A = WAVI(N1+4)                                                    
      B = BS(I)                                                         
      C = CS(I)                                                         
      D = DS(I)                                                         
      UDERIV(2,1) = FACTOR*( A + X*(B + X*(C + X*D )))                  
      UDERIV(2,2) = FACTOR*( B + X*(2.0D0*C + 3.0D0*D*X))                       
      UDERIV(2,3) = FACTOR* 2.0D0*( C + 3.0D0*D*X )                             
      UDERIV(2,4) = FACTOR * 6.0D0*D                                        
C                                                                       
C     SAVE STUFF TO CHECK NEXT TIME FOR CHANGES                         
C                                                                       
 270  UROLD = ABACK                                                     
      UIOLD = BBACK                                                     
      LOLD = L                                                          
      IF (PBUGSW) WRITE (IUNIT, 273) UDERIV                             
 273  FORMAT ( ' POT DERIV', ( T12, 2( 2G11.3, 3X ) ) )                 
C                                                                       
C     NOW GENERATE THE DERIVATIVES OF THE WAVEFUNCTION                  
C                                                                       
      FDERIV(1,1) = 1.0D0                                                   
      FDERIV(2,1) = 0.0D0                                                   
      FDERIV(1,2) = +DERIV2 * FDERIV(1,1)                               
      FDERIV(2,2) = -DERIV1 * FDERIV(1,1)                               
C                                                                       
C     WE GENERATE THE HIGHER DERIVATIVES BY RECURSION ON THE            
C     SCHROEDINGER EQUATION.                                            
C                                                                       
      BC(1) = 1.0D0                                                         
      DO 349 N = 2, 5                                                   
      BC(N) = 1.0D0                                                         
      NM = N-1                                                          
      FDERIV(1,N+1) = 0.0D0                                                 
      FDERIV(2,N+1) = 0.0D0                                                 
      DO 329 I = 1, NM                                                  
      FDERIV(1,N+1) = FDERIV(1,N+1) - BC(N-I) * ( FDERIV(1,I)*
     * UDERIV(1,N-I) - FDERIV(2,I)*UDERIV(2,N-I) )                                
      FDERIV(2,N+1) = FDERIV(2,N+1) - BC(N-I) * (  FDERIV(2,I)*UDERIV(1,
     *N-I) + FDERIV(1,I)*UDERIV(2,N-I) )                                
      IF ( I .NE. NM ) BC(N-I) = BC(N-I) + BC(N-I-1)                    
 329  CONTINUE                                                          
 349  CONTINUE                                                          
C                                                                       
C     NOW GENERATE THE TERMS AND SUM THEM TO GET F(R1) AND F(R2)        
C                                                                       
      X = -X                                                            
      DO 389 II = 1, 2                                                  
      DO 379 I = 1, 2                                                   
      FSUMS(I,II) = 0.0D0                                                   
      FAC = 1.0D0                                                           
      DO 379 N = 1, 6                                                   
      FTERMS(I,N,II) = FAC*FDERIV(I,N)                                  
      FSUMS(I,II) = FSUMS(I,II) + FTERMS(I,N,II)                        
      IF ( DABS(FAC) .LT. SMLNUM ) FAC = 0.0D0                              
      FAC = X*FAC/N                                                     
 379  CONTINUE                                                          
      X = X + STPSZ                                                     
 389  CONTINUE                                                          
C                                                                       
      IF (PBUGSW) WRITE (IUNIT, 393) FTERMS, FSUMS                      
 393  FORMAT ( ' FTERMS', 3( T9, 2( 2G11.3, 3X ) / )  ' FTERMS', 3( T9, 
     *2( 2G11.3, 3X ) / )  ' FSUMS', ( T9, 2( 2G11.3, 3X ) ) )          
C                                                                       
C     NOW START THE INTEGRATION FROM N1                                 
C                                                                       
      BACKR = 1.D-15 * FSUMS(1,1)                                       
      BACKI = 1.D-15 * FSUMS(2,1)                                       
      THISR = 1.D-15 * FSUMS(1,2)                                       
      THISI = 1.D-15 * FSUMS(2,2)                                       
C                                                                       
C     FILL IN THE INTERIOR WAVEFUNCTION WITH ZEROS                      
C                                                                       
      DO 399 I = ISTRT, N1                                              
      WAVR(I+1) = 0.0D0                                                     
      WAVI(I+1) = 0.0D0                                                     
 399  CONTINUE                                                          
C                                                                       
      ISTRT = N1                                                        
      NFIRST = ISTRT + 2                                                
C                                                                       
      WAVR(ISTRT+1) = BACKR                                             
      WAVI(ISTRT+1) = BACKI                                             
      WAVR(ISTRT+2)=THISR                                               
      WAVI(ISTRT+2)=THISI                                               
C                                                                       
C     COMPUTE FIRST TWO XSI'S                                           
C                                                                       
      II = ISTRT+1                                                      
      DO 309 I = ISTRT, II                                              
      WAVR(I+3) =  TWELTH*( WAVR(I+1)*WAVR(I+4) - WAVI(I+1)*WAVI(I+4) ) 
      WAVI(I+3) =  TWELTH*( WAVR(I+1)*WAVI(I+4) + WAVI(I+1)*WAVR(I+4))  
 309  CONTINUE                                                          
C                                                                       
C     WILL DO INTEGRATION HERE                                          
C                                                                       
 350  CONTINUE                                                          
C                                                                       
C     MODIFIED NUMEROV (A LA RAYNAL) IS USED IN THE FOLLOWING LOOP      
C                                                                       
C     WE DEFINE                                                         
C       W(I) = 1 - (H**2/12)*(V(I)/E - 1)                               
C       XSI(I) = W(I)*U(I)                                              
C     WHERE                                                             
C       V(I) = V(I*STEPSZ) = V(I*H/K)                                   
C     CONTAINS THE COMPLETE NUCLEAR, COULOMB AND CENTRIFIGAL POTENTIALS.
C     MOST OF W(I) IS COMPUTED AND STORED IN MAKPOT, HOWEVER            
C     THE L-DEPENDENT PIECES ARE ADDED IN JUST ABOVE.                   
C                                                                       
C     THE DESIRED WAVEFUNCTION IS                                       
C       U(I) = U(I*STEPSZ)                                              
C     EXCEPT THAT AT THIS STAGE THE NORMALIZATION IS ARBITRARY:         
C       U(R) ---> CONS*R**(L+1) .                                       
C                                                                       
C     AT STEP I IN THE LOOP, THE WAVR AND WAVI ARRAYS CONTAIN:          
C                                                                       
C     WAV_(...)   CONTENTS                                              
C          I-1    U(I-2)     NOT USED                                   
C          I      U(I-1)                                                
C          I+1    XSI(I-2),   CHANGED TO U(I) DURING STEP               
C          I+2    XSI(I-1)                                              
C          I+3    W(I-1), NOT USED, CHANGED TO XSI(I) DURING STEP       
C          I+4    W(I)                                                  
C                                                                       
C     THE SMALLEST POSSIBLE NFIRST IS 2 WHICH RESULTS IN ISTRT=0.       
C     THE U AND XSI FOR NFIRST-2, NFIRST-1 HAVE ALREADY BEEN STORED.    
C                                                                       
C     ON THE FIRST PARTIAL WAVE, WE MUST FROM TIME TO TIME CHECK        
C     FOR IMMENENT OVERFLOW AND RENORMALIZE DOWN.  AFTER THE FIRST      
C     PARTIAL WAVE WE KNOW WHERE TO START AND NEED NOT CHECK.           
C                                                                       
      N1 = NFIRST                                                       
      N2 = NSTEP                                                        
      IF ( L .GE. IWAVBL(114+NWP) ) GO TO 420                           
C                                                                       
CCCCC     GET INTERVAL INWHICH WE MAY REACH BIGNUM                      
CCCCC                                                                   
CCCC 410  N2 = N1*(BIGNUM/THISR)**(1./(DL+1))                           
CCCC      N2 = MAX0( N2, N1+5 )                                         
 410  N2 = N1 + 20                                                      
      N2 = MIN0(N2, NSTEP)                                              
C                                                                       
 420  D = ( 12.0D0 / ( WAVR(I+4)**2 + WAVI(I+4)**2 ) )                     
C                                                                       
      DO 449 I = N1, N2                                                 
C                                                                       
      WAVR(I+3) = ((-10.0D0)*WAVR(I+2)) + (WAVR(I)-WAVR(I+1))              
      WAVI(I+3) = ((-10.0D0)*WAVI(I+2)) + (WAVI(I)-WAVI(I+1))              
C                                                                       
      WAVR(I+1) = ( WAVR(I+4)*WAVR(I+3) + WAVI(I+4)*WAVI(I+3) ) * D     
      WAVI(I+1) = ( WAVR(I+4)*WAVI(I+3) - WAVI(I+4)*WAVR(I+3) ) * D     
C                                                                       
      D = ( 12.0D0 / ( WAVR(I+5)**2 + WAVI(I+5)**2 ) )                     
C                                                                       
 449  CONTINUE                                                          
      IF ( N2 .EQ. NSTEP ) GO TO 800                                    
C                                                                       
C     CHECK FOR NEARNESS OF OVERFLOW                                    
C                                                                       
      THISR = DABS(WAVR(N2+1)) + DABS(WAVI(N2+1))                       
      IF ( PBUGSW ) WRITE (IUNIT, 463) NWP, L, N1, N2, THISR            
 463  FORMAT ( ' CHECKING SIZE: NWP, L, N1, N2, MAG =',  4I5, G13.3 )   
      IF ( PBUGSW ) NUMOUT = 1                                          
C                                                                       
      IF ( THISR .LT. BIGNUM ) GO TO 490                                
C                                                                       
C     MUST SCALE DOWN.  WE WILL MULTIPLY EVERYTHING BY SMLNUM.          
C     FIRST WE FIND WHERE THIS WILL RESULT IN \U\ = SMLNUM.             
C                                                                       
      THISR = THISR*SMLNUM                                              
      N1 = ISTRT                                                        
      DO 479 ISTRT = N1, N2                                             
      IF ( DABS(WAVR(ISTRT+1)) .GE.1.0D0 ) GO TO 480                       
      WAVR(ISTRT+1) = 0.0D0                                                 
      WAVI(ISTRT+1) = 0.0D0                                                
 479  CONTINUE                                                          
C                                                                       
 480  IF ( PBUGSW ) WRITE (IUNIT, 483) ISTRT                            
 483  FORMAT ( 5X, 'RESCALING: NEW ISTRT =', I6 )                       
C                                                                       
      N1 = N2+2                                                         
      DO 489 I = ISTRT, N1                                              
      WAVR(I+1) = SMLNUM*WAVR(I+1)                                      
      WAVI(I+1) = SMLNUM*WAVI(I+1)                                      
 489  CONTINUE                                                          
C                                                                       
 490  N1 = N2+1                                                         
      GO TO 410                                                         
C                                                                       
C     END OF THE INTEGRATION                                            
C                                                                       
 800  IF ( NUMOUT .GT. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
      END 
      
      SUBROUTINE PARITW  ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR,IUNIT,
     * NUMOUT,L,J, RSTART,STEPSZ,NUMPTS, ARRAY1,ARRAY2,FLOATK,TEMPVS, 
     * PARAM,INTGER,JBLOCK,ISWTCH, INTRNL, FINTRN,CNSTNT,IWAVBL, 
     * FKANDM,ID,ALLOC,ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC) 
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
C                                                                       
C     PARITWOODSAXON - PARITY DEPENDANT WOODS SAXON                     
C                                                                       
C     THE POTENTIAL IS A WOODS SAXON WITH PARITY DEPENDANT DEPTH:       
C                                                                       
C       V = (1+PARAM1*(-1)**L) * WOODS SAXON( V R A )                   
C                                                                       
C     6/13/78 - BASED ON FIXEDWOOD - S.P.                               
C     7/11/78 - CDC VERSION AND IRETUR CORRECTION                       
C                                                                       
      INTEGER WORKNM(2,2)                                                
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                        
      LOGICAL TBIT                                                      
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
C                                                                       
      PCOEF = PARAM(IPOTYP)                                             
      V = TEMPVS(IPOTYP)                                                
      A = TEMPVS(4+IPOTYP)                                              
      R = FLOATK(43)                                                     
      IF ( IPOTYP .EQ. 2 ) R = FLOATK(47)                                
      POWER = FLOATK(101+IPOTYP)                                         
      IF ( R.LE. 0.0D0 .OR. A.LE.0.0D0 ) GO TO 960                           
C                                                                       
      GO TO ( 100, 500, 600, 700 ), IREQUE                              
C                                                                       
C     GENERATE WORK ARRAY TO SAVE THE W.S.                              
C                                                                       
 100  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
C                                                                       
      I = NALLOC( NUMPTS, WORKNM(1,1) )                                 
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = I                                                     
      MYINTS(2) = II                                                    
C                                                                       
C     HERE WE INDICATE THAT THIS IS A L-DEPENDENT POTENTIAL             
C                                                                       
      IRETUR = +1                                                       
C                                                                       
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503) V, R, A, POWER, PCOEF                          
 503  FORMAT ( ' THE POTENTIAL IS A PARITY WOODS-SAXON:', 3G15.5 ,  5X, 
     *'POWER =', G14.5 /  ' WITH THE PARITY-DEPENDANT DEPTH FACTOR',  ' 
     * 1 + (-PARAM1)**L:  PARAM1 =', G15.5 )                            
      GO TO 900                                                         
C                                                                       
C     SET UP THE TWO PIECES                                             
C                                                                       
 600  LL = KLOC(MYINTS(1))                                               
      CALL WOODSX ( NUMPTS, RSTART, STEPSZ, ALLOC(LL), N1, N2, 1, V, R,
     * A, POWER )                                                       
      DO 619 I = 1, NUMPTS                                              
      ARRAY1(I) = ALLOC(LL-1+I)                                         
 619  CONTINUE                                                          
C                                                                       
      RETURN                                                            
C                                                                       
C     ADD IN THE L-DEPENDANT PART                                       
C                                                                       
 700  LL = KLOC(MYINTS(1))                                               
      FAC = ARRAY2(1)*PCOEF                                             
      IF ( MOD(L,2).NE.0 ) FAC = -FAC                                   
      DO 719 I = 1, NUMPTS                                              
      ARRAY1(I) = ARRAY1(I) + FAC*ALLOC(LL-1+I)                         
 719  CONTINUE                                                          
      RETURN                                                            
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
C-JSW IRETUR = -1                                                       
C-JSW GO TO 900                                                         
C                                                                       
 960  IRETUR = -1                                                       
      RETURN                                                            
C                                                                       
      END
      
      SUBROUTINE RAWITS ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ, NUMPTS, WAVER, WAVEI,FLOATK,TEMPVS, 
     * PARAM, INTGER, JBLOCK, ISWTCH, INTRNL, FINTRN,  CNSTNT, IWAVBL,
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC, LENGTH, NALLOC, NAMLOC)
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), WAVER(NUMPTS), WAVEI(NUMPTS), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
      CHARACTER*4 ID                                                    
C                                                                       
C     LINKULE FOR RAWITSCHER TYPE INCOMING BOUNDARY CONDITION           
C                                                                       
C     1/22/80 - STANDARD MORTRAN, CVA%M'S - RPG                         
C                                                                       
      LOGICAL PBUGSW                                                    
      EQUIVALENCE ( IBUGSW, PBUGSW )                                    
C                                                                       
      CHARACTER*4 WORKNM(2)                                             
      DATA WORKNM/'*000', 'WORK'/                                       
      DIMENSION UDERIV(2, 4), FDERIV(2, 6), BC(5),  FSUMS(2, 2)         
      COMPLEX CSUM, CK, CF(2)                                           
      EQUIVALENCE ( CF(1), FSUMS(1) )                                   
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      IBUGSW = IWAVBL(114)                                              
C                                                                       
      GO TO ( 100, 500, 600 ), IREQUE                                   
C                                                                       
C     LOCATE THE BOUNDRY RADIUS (PARAM1)                                
C                                                                       
 100  IF ( PARAM(1) .NE. UNDEF  )                                       
     1    GO TO 120                                                     
      WRITE ( IUNIT, 113 )                                              
 113  FORMAT ( '0**** PARAM1 MUST BE DEFINED',                          
     1  ' AS THE BOUNDARY RADIUS.' )                                    
      IRETUR = -1                                                       
      GO TO 900                                                         
 120  WORKNM(1) = ID                                                    
      IWRK = NALLOC( 100,WORKNM)                                        
      ALLOC(KLOC(IWRK)) = STEPSZ                                         
C                                                                       
C     SAVE ADDRESS OF WORK ARRAY                                        
C                                                                       
      MYINTS(1) = IWRK                                                  
      IF ( NUMOUT .NE. 0 )  GO TO 900                                   
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503)  PARAM(1)                                      
 503  FORMAT ( ' THE WAVEFUNCTION WILL HAVE A',                         
     1  ' RAWITSCHER-TYPE BOUNDARY CONDITION AT R =', F6.2 )            
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
C     GET SOME ARGUMENTS INTO THEIR CONVENTIONAL NAMES                  
C                                                                       
 600  NFIRST = IPOTYP                                                   
      H = STEPSZ                                                        
      ISTRT = NFIRST-2                                                  
      NSTEP = NUMPTS-1                                                  
      IWRK = MYINTS(1)                                                  
      LWRK = KLOC(IWRK)                                                  
      STPSZ = ALLOC(LWRK)                                               
      RB = PARAM(1)                                                     
C                                                                       
C     NOW GET STARTING VALUES OF THE POTENTIAL AND SAVE STARTING        
C     VALUES OF THE WAVEFUNCTION.                                       
C                                                                       
      N1 = RB/STPSZ                                                     
      IF ( N1 .GE. ISTRT )  GO TO 250                                   
C                                                                       
C     THE BOUNDRY IS INSIDE THE FIRST POINT SO IT DOESNT COUNT          
C                                                                       
      THISR = 1.D-30                                                    
      THISI = 0.0D0                                                         
      BACKR = (RSTART/(RSTART+H))**(L+1) * THISR                        
      BACKI = 0.0D0                                                         
      ABACK = WAVER(ISTRT)                                              
      BBACK = WAVEI(ISTRT)                                              
      ATHIS = WAVER(ISTRT+1)                                            
      BTHIS = WAVEI(ISTRT+1)                                            
      GO TO 400                                                         
C                                                                       
C     WE WILL START AT N1                                               
C                                                                       
 250  ABACK = WAVER(N1)                                                 
      BBACK = WAVEI(N1)                                                 
      ATHIS = WAVER(N1+1)                                               
      BTHIS = WAVEI(N1+1)                                               
C                                                                       
C     GET CUBIC SPLINES FOR THE POTENTIAL                               
C                                                                       
      L0 = LWRK+2                                                       
      L1 = L0+21                                                        
      L2 = L1+21                                                        
      L3 = L2+21                                                        
      I1 = MAX0( N1-10, ISTRT )                                         
      DO 254  I = 1, 21                                                 
         ALLOC(L0+I-1) = I*STPSZ                                        
 254  CONTINUE                                                          
      R1 = N1*STPSZ                                                     
      X = RB - R1                                                       
      I = N1-I1                                                         
      FACTOR = 12/STPSZ**2                                              
      IF ( PBUGSW )  WRITE ( IUNIT, 263 )  ISTRT, N1, R1, X             
 263  FORMAT ( ' BARRIER ISTRT, N1, R1, X =', 2I6, 2G15.5 )             
      IF (PBUGSW)  NUMOUT = 1                                           
C                                                                       
C     FIRST THE REAL POTENTIAL DERIVATIVES                              
C                                                                       
      CALL SPLNCB ( 21, ALLOC(L0), WAVER(I1), ALLOC(L1),                
     1  ALLOC(L2), ALLOC(L3) )                                          
      A1 = WAVER(N1)                                                    
      B1 = ALLOC(L1+I)                                                  
      C1 = ALLOC(L2+I)                                                  
      D1 = ALLOC(L3+I)                                                  
C                                                                       
C     NOW THE IMAG POTENTIAL DERIVATIVES                                
C                                                                       
      CALL SPLNCB ( 21, ALLOC(L0), WAVEI(I1), ALLOC(L1),                
     1  ALLOC(L2), ALLOC(L3) )                                          
      A2 = WAVEI(N1)                                                    
      B2 = ALLOC(L1+I)                                                  
      C2 = ALLOC(L2+I)                                                  
      D2 = ALLOC(L3+I)                                                  
C                                                                       
C     INITIALIZE COMPUTATION OF EFFECTIVE K                             
C                                                                       
      CALL CKSET ( R1, STPSZ, A1, B1, C1, D1, A2, B2, C2, D2 )          
C                                                                       
C     DO RAWITSCHER INTEGRALS BY SIMPLE SIMPSON'S RULE                  
C                                                                       
      R = R1                                                            
      DO 359  I = 1, 2                                                  
         DEL = (R-RB)/4.0D0                                                 
         CSUM = CK(RB) + 4.0D0*CK(RB+DEL) + 2.0D0*CK(RB+2*DEL)                  
     1       + 4.0D0*CK(RB+3.0D0*DEL) + CK(R)                                   
         CSUM = CSUM * DEL/3.0D0                                            
CNV         CSUM = CDEXP( (0.D0,-1.D0)*CSUM )                           
         CSUM = CEXP( (0.0,-1.0)*CSUM )                                   
      CF(I) = CSUM / CSQRT(CK(R))                                       
         R = R + STPSZ                                                  
 359  CONTINUE                                                          
C                                                                       
      IF (PBUGSW)  WRITE (IUNIT, 393)  FSUMS                            
 393  FORMAT (                                                          
     1  ' FSUMS', ( T9, 2( 2G11.3, 3X ) ) )                             
C                                                                       
C     NOW START THE INTEGRATION FROM N1                                 
C                                                                       
      BACKR = 1.D-15 * FSUMS(1,1)                                       
      BACKI = 1.D-15 * FSUMS(2,1)                                       
      THISR = 1.D-15 * FSUMS(1,2)                                       
      THISI = 1.D-15 * FSUMS(2,2)                                       
C                                                                       
C     FILL IN THE INTERIOR WAVEFUNCTION WITH ZEROS                      
C                                                                       
      DO 399  I = ISTRT, N1                                             
         WAVER(I) = 0.0D0                                                   
         WAVEI(I) = 0.0D0                                                   
 399  CONTINUE                                                          
C                                                                       
      ISTRT = N1                                                        
      NFIRST = ISTRT + 2                                                
C                                                                       
 400  WAVER(ISTRT) = BACKR                                              
      WAVEI(ISTRT) = BACKI                                              
      WAVER(ISTRT+1)=THISR                                              
      WAVEI(ISTRT+1)=THISI                                              
C                                                                       
      DO 489 ISTEP = NFIRST,NSTEP                                       
C                                                                       
C     ALLOC(LWAV_+ISTEP) = WAVEFUNCTION(R = STEPSIZE*ISTEP)             
C      FOR  ISTEP = NFIRST-2,  WAVEFUNCTION IS 0                        
C                 = NFIRST-1,  WAVEFUNCTION IS RNORM*STEP_              
C                                                                       
C  INTEGRATION PARAMETER IS RHO                                         
C                                                                       
C   FOX-GOODWIN METHOD OF INTEGRATION USED HERE.                        
C                                                                       
C     HERE                                                              
C       "THIS"  MEANS  RHO = (ISTEP-1)*H                                
C       "BACK"  MEANS  RHO = (ISTEP-2)*H                                
C                                                                       
         CR = BBACK*BACKI - ABACK*BACKR + (10.0D0*BTHIS)*THISI              
     1      - (10.0D0*ATHIS-12.0D0)*THISR                                       
         CI = - ( BBACK*BACKR + ABACK*BACKI + (10.0D0*BTHIS)*THISR          
     1      + (10.0D0*ATHIS-12.0D0)*THISI )                                     
C                                                                       
         ABACK = ATHIS                                                  
         BBACK = BTHIS                                                  
         ATHIS = WAVER(ISTEP)                                           
         BTHIS = WAVEI(ISTEP)                                           
C                                                                       
C     NOW                                                               
C       "THIS"  MEANS  RHO = ISTEP*H                                    
C       "BACK"  MEANS  RHO = (ISTEP-1)*H                                
C                                                                       
         D0 = ATHIS**2 + BTHIS**2                                       
         BACKR = THISR                                                  
         BACKI = THISI                                                  
         THISR = (ATHIS*CR + BTHIS*CI)/D0                               
         WAVER(ISTEP) = THISR                                           
         THISI = (ATHIS*CI - BTHIS*CR)/D0                               
         WAVEI(ISTEP) = THISI                                           
C                                                                       
C     IF THE WAVEFUNCTION IS GETTING TOO BIG, NORMALIZE IT DOWN         
C                                                                       
         IF ( DABS(THISR) .LT. 1.D+30 )  GO TO 489                      
C                                                                       
C     WILL DIVIDE EVERYTHING BY 1.D+30.  FIND WHERE THIS WILL RESULT    
C     IN  \U\ = 1.D-30.                                                 
C                                                                       
         BACKR = 1.D-30 * BACKR                                         
         BACKI = 1.D-30 * BACKI                                         
         THISR = 1.D-30 * THISR                                         
         THISI = 1.D-30 * THISI                                         
         DO 439  I = ISTRT, ISTEP                                       
            IF ( DABS(WAVER(I)) .GE. 1.0D0 )  GO TO 440                     
            WAVER(I) = 0.0D0                                                
            WAVEI(I) = 0.0D0                                                
 439     CONTINUE                                                       
 440     IF ( PBUGSW )  WRITE ( IUNIT, 443 ) NWP, L, ISTRT, I, ISTEP    
 443     FORMAT ( ' RESCALING: NWP, L, OLD ISTRT, NEW ISTRT, ISTEP =',  
     1      5I6 )                                                       
      IF ( PBUGSW )  NUMOUT=1                                           
         ISTRT = I                                                      
         DO 449  I = ISTRT, ISTEP                                       
            WAVER(I) = 1.D-30 * WAVER(I)                                
            WAVEI(I) = 1.D-30 * WAVEI(I)                                
 449     CONTINUE                                                       
C                                                                       
 489  CONTINUE                                                          
C                                                                       
      IF ( NUMOUT .GT. 0 )  GO TO 900                                   
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
      END 
      
      FUNCTION CK( R )                                                  
C                                                                       
C     COMPUTES THE EFFECTIVE COMPLEX K                                  
C                                                                       
      RETURN                                                            
      END  
      
      SUBROUTINE REID (ALIAS, SAVINT, IPOTYP, IREQUE, IRETUR,  IUNIT, 
     * NUMOUT,L,JP,RSTART,STEPSZ, NUMPTS,ARRAY1, ARRAY2,FLT,TEMPVS, 
     * PARAM, INTGER, JBLOCK,  ISWTCH, INTRNL, FINTRN, CNSTNT,IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4,KLOC, LENGTH, NALLOC, NAMLOC, 
     * IPARAM )                                                             
C                                                                       
C     LINKULE FOR REID SOFT-CORE DEUTERON WAVE FUNCTION                 
C                                                                       
C     RELEVANT ARGUMENTS ARE:                                           
C                                                                       
C     ALIAS = LINKULE NAME, USED IN ERROR MESSAGES (INPUT).             
C     SAVINT = A TWO-ELEMENT REAL*4 ARRAY TO STORE THE INTEGRALS        
C        OF WF**2 AND WF*R**(L+1).                                      
C     IPOTYP = 6 (BOTH WF AND POTENTIAL) OR 1 (POTENTIAL ONLY) (INPUT). 
C     IREQUE = 1 FOR INITIALIZATION, 2 FOR PRINTING, 3 FOR CALCULATION  
C        (INPUT).                                                       
C     IRETUR = ERROR RETURN CODE:  <0 = ERROR (OUTPUT ).                
C     IUNIT = UNIT NO. FOR PRINTED OUTPUT (INPUT).                      
C     NUMOUT = NO. OF LINES TO BE PRINTED (OUTPUT).                     
C     L = ORBITAL ANGULAR MOMENTUM, MUST = 0 OR 2 (INPUT).              
C     JP = 2*PROJECTILE "TOTAL" ANGULAR MOMENTUM; SET TO                
C        1 OR 3 FOR L = 0 OR 2 RESPECTIVELY (OUTPUT).                   
C     RSTART = STARTING R VALUE (FERMIS) (INPUT).                       
C     STEPSZ = GRID SPACING (FERMIS) (INPUT).                           
C     NUMPTS = NO. OF GRID POINTS TO BE CALCULATED (INPUT).             
C     ARRAY1 (LENGTH NUMPTS) = PRIMARY OUTPUT ARRAY, SET TO             
C        ARRAY2(1)*THE POTENTIAL IF IPOTYP=1, OR TO THE                 
C        WAVE FUNCTION IF IPOTYP=6.                                     
C     ARRAY2: IF IPOTYP=1, ARRAY2(1) IS THE FACTOR BY WHICH THE LINKULE 
C        MUST MULTIPLY THE POTENTIAL.                                   
C        IF IPOTYP=2, ARRAY2 (LENGTH NUMPTS) IS SET TO                  
C        (-1/V) * THE POTENTIAL (OUTPUT).                               
C     FLT(62) = V, SET TO 1.0 (OUTPUT).                                 
C     FLT(43) = R, SET TO 1.0 IF UNDEFINED.                             
C     FLT(1) = A, SET TO 0.4 IF UNDEFINED.                              
C     FLT(12) = E, SET TO -2.224644 IF UNDEFINED.                       
C     FLT(24) = AM = REDUCED MASS IN MEV/C**2 (INPUT).                  
C     FLT(53) = SPAM = SPECTROSCOPIC AMPLITUDE.                         
C        SET TO SQRT(1-PD) OR SQRT(PD) IF L = 0 OR 2, RESPECTIVELY,     
C        IF INITIALLY UNDEFINED, WHERE PD = 0.064696.                   
C     FLT(54) = SPAMP = PROJECTILE SPAM.  SET TO SPAM IF NEXT = 1.      
C     FLT(55) = SPAMT = TARGET SPAM.  SET TO SPAM IF NEXT = 2.          
C     PARAM(1) = R SCALE COMPRESSION FACTOR.                            
C     INTGER(12) = NODES, SET TO 0 (OUTPUT).                            
C     INTGER(17) = IPRINT = PRINT CONTROL.  DUMPS IF THE 1 DIGIT        
C        IS >= 6.                                                       
C     JBLOCK(1) = 2*TOTAL ANGULAR MOMENTUM, SET TO 2 (OUTPUT).          
C     JBLOCK(8) = 2*PROJECTILE SPIN, SET TO 1 (OUTPUT).                 
C     JBLOCK(9) = 2*TARGET SPIN, SET TO 1 (OUTPUT).                     
C     ISWTCH(11) = NEXT = 1 IF THIS IS THE PROJECTILE, 2 IF TARGET.     
C     INTRNL(3) = NOTDEF = CONSTANT FOR UNDEFINED INTEGERS.             
C     FINTRN(1) = UNDEF = CONSTANT FOR UNDEFINED REAL NUMBERS (INPUT).  
C     CNSTNT(1) = PI (INPUT).                                           
C     CNSTNT(2) = SQRT(4*PI) (INPUT).                                   
C     CNSTNT(6) = HBARC = HBAR*C (MEV FM) (INPUT).                      
C     CNSTNT(13) = BIGLOG = NATURAL LOG OF A BIG NUMBER WHOSE SQUARE    
C        DOESN'S QUITE OVERFLOW (INPUT).                                
C                                                                       
C     REFERENCE:  R. V. REID, ANN. OF PHYSICS 50 (1968) 411-448.        
C           THIS LINKULE MAY BE USED IN EITHER OF TWO WAYS:             
C                                                                       
C           TO CALCULATE BOTH THE WAVEFUNCTION AND THE POTENTIAL, THE   
C     MINIMUM INPUT IS                                                  
C        "CHANNEL:  N + P = D"                                          
C        (OR AN APPROPRIATE "REACTION" AND "PROJECTILE" OR "TARGET")    
C        "R0 = 1   A = 0.5"      (TO AVOID COMPLAINTS)                  
C        "WAVEFUNCTION = REID"                                          
C        "L = 0;"  (OR "L = 2;")                                        
C     THEN THE LINKULE RETURNS THE WAVEFUNCTION IN ARRAY1, AND          
C     (-1/V)*THE POTENTIAL IN ARRAY2.                                   
C           TO CALCULATE ONLY THE EFFECTIVE POTENTIAL, AND LET BOUND    
C     CALCULATE THE WAVE FUNCTION, USE "REALPOTENTIAL" INSTEAD OF       
C     "WAVEFUNCTION" IN THE INPUT ABOVE, AND ADD "NODES = 0".           
C     THEN THE LINKULE RETURNS ARRAY2(1)*THE POTENTIAL IN ARRAY1.       
C     (NOTE:  THE PRESENT VERSION OF  BOUND  (10/79) WILL NOT CALCULATE 
C     THE WAVE FUCTION FROM THIS POTENTIAL BECAUSE IT CHOKES ON THE     
C     CORE.)                                                            
C           NORMALIZATION:  THE INTEGRAL OF THE SQUARE OF THE WAVE      
C     FUNCTION RETURNED IS 1.0.  SPECTROSCOPIC INFORMATION IS ENTERED   
C     USING THE KEYWORD "SPAM", WHICH IS NOT ACTUALLY USED UNTIL MUCH   
C     LATER IN THE DWBA CALCULATION.  IF "SPAM" IS NOT DEFINED, THIS    
C     LINKULE SUPPLIES DEFAULT VALUES FROM THE ACTUAL WAVE FUNCTIONS:   
C     SPAM(L=0) = SQRT(1-PD), SPAM(L=2) = SQRT(PD),                     
C     WHERE PD = 0.064696, N = 0.8776, RHO = 0.02622,                   
C     D0 = 1.251*D00 = -125.2 MEV FM**(3/2), D2 = 0.4839 FM**2.         
C                                                                       
C           NOTE:  IN DWBA CALCULATIONS, THE DEFAULT SPECTROSCOPIC      
C     AMPLITUDE IS ALWAYS USED UNLESS THE KEYWORD "SPAM" IS USED IN     
C     THE DESCRIPTION OF THIS CHANNEL.  SETTING "SPAMP" OR "SPAMT" WILL 
C     NOT WORK.  IN STAND-ALONE CALCULATIONS, YOU MUST "UNDEFINE SPAM"  
C     TO FORCE USE OF THE DEFAULT IN THE SECOND AND SUBSEQUENT          
C     CALCULATIONS.  F "SPAM" IS USED, REMEMBER THAT A PROLATE          
C     DEUTERON REQUIRES A NEGATIVE D-STATE AMPLITUDE.                   
C                                                                       
C           R COMPRESSION:  FOR DWBA TESTS, IT IS CONVENIENT TO CON-    
C     STRUCT A FAKE DEUTERON WHICH IS SMALLER THAN A REAL DEUTERON,     
C     WHILE PRESERVING ITS LEADING LOW-MOMENTUM TERMS.  THIS IS         
C     ACCOMPLISHED BY ENTERING THE COMPRESSION FACTOR USING THE         
C     "PARAM1" KEYWORD.  FOR EXAMPLE, THE FOLLOWING INPUT (IN ADDITION  
C     TO THE INPUT ABOVE) WILL CREATE A "DEUTERON" ONLY HALF AS BIG AS  
C     A REAL DEUTERON, BUT WITH THE SAME D(L) AND PROBABILITY INTEGRAL: 
C        "PARAM1 = 2                                                    
C        "STEPSIZE = 0.05"  (OR "A = 0.25" IF "STEPSPER" IS DEFINED)    
C        "BOUNDASY = 10   LOOKSTEP = 500"                               
C     THIS SHOULD GIVE DWBA RESULTS SIMILAR TO A NORMAL CALCULATION     
C     WITH STEPSIZE=0.1, BOUNDASY=20, AND LOOKSTEP=250, EXCEPT THE      
C     FINITE-RANGE EFFECTS SHOULD BE REDUCED BY A FACTOR OF 4.          
C     FOR COMPRESSED CALCULATIONS, ALWAYS USE "WAVEFUCNTION",           
C     NOT 'REALPOTENTIAL" (EVEN IF THE LATTER WORKS).                   
C     METHOD:                                                           
C           THE WAVE FUNCTIONS (BOTH S AND D) ARE CALCULATED BY         
C     INTERPOLATING BETWEEN THE VALUES TABULATED IN THE APPENDIX        
C     OF REID'S PAPER.  CUBIC INTERPOLATION IS USED, AS REID SUGGESTS.  
C     THE CENTRAL, TENSOR, AND SPIN-ORBIT POTENTIALS ARE CALCULATED     
C     USING REID'S EQUATION 30; THEY ARE SUMS OF YUKAWA TERMS.          
C     THE "EFFECTIVE POTENTIAL" FOR THE S (D) WAVE IS JUST THE          
C     DIAGONAL (CENTRAL, SPIN-ORBIT, AND TENSOR) POTENTIAL, PLUS        
C     THE OFF-DIAGONAL TENSOR POTENTIAL MULTIPLIED BY THE D/S (S/D)     
C     WAVE FUNCTION RATIO.                                              
C                                                                       
C     7/10/79 - NEW ROUTINE - RPG                                       
C     9/27/79 - REID CREATED FROM MCGEE - RPG                           
C     11/19/79 - ADD R COMPRESSION FACTOR - RPG                         
C     12/12/79 - PRINT COMPRESSION FACTOR. - RPG                        
C     1/18/80 - CNI%'LINKUL' = 'REID' - RPG                             
C     8/8/80 - FIX UNDEFINED PRINTOUT BUG - S.P.                        
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
      CHARACTER*8 ALIAS                                                 
      CHARACTER*1 SD(2)                                                 
      DATA SD/'S', 'D'/                                                 
      REAL*4 SAVINT(2)                                                  
      DIMENSION ARRAY1(1), ARRAY2(1),FLT(200),TEMPVS(6), PARAM(20), 
     *INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74),FINTRN(37), 
     *CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1),KLOC(1), 
     *LENGTH(1),IPARAM(5)                                           
CCCC      EXTERNAL NALLOC, NAMLOC                                       
C                                                                       
      LOGICAL L2SW, WFSW, PBUGSW                                        
      DIMENSION XX(34), UI(34,2), UPI(34,2), YNRMS(2), SPAMS(2)         
C                                                                       
C     THESE WAVE FUNCTION VALUES ARE FROM THE APPENDIX OF REID'S        
C     PAPER.  NOTE THAT X = 0.7*R, AND UP = DU/DX.                      
C                                                                       
      DATA XX /  .0100D0, .041250D0, .07250D0, .1350D0, .19750D0,  .2600
     *D0, .32250D0, .3850D0, .44750D0, .5100D0,  .57250D0, .6350D0, .697
     *50D0, .7600D0, .8850D0,  1.0100D0, 1.1350D0, 1.2600D0, 1.3850D0, 1
     *.5100D0,  1.7600D0, 2.0100D0, 2.5100D0, 3.0100D0, 3.5100D0,  4.010
     *0D0, 4.5100D0, 5.0100D0, 5.5100D0, 6.0100D0,  7.0100D0, 8.0100D0, 
     *9.0100D0, 10.0100D0 /                                             
      DATA UI /  .000000D0, .333730D-4, .239010D-3, .276210D-2, .127370D
     *-1,  .360620D-1, .753590D-1, .128470D0, .189930D0, .253490D0,  .31
     *3900D0, .367700D0, .413170D0, .449920D0, .499530D0,  .524060D0, .5
     *31660D0, .528640D0, .519260D0, .506210D0,  .475050D0, .442000D0, .
     *378640D0, .322490D0, .273990D0,  .232510D0, .197190D0, .167180D0, 
     *.141720D0, .120120D0,  .862900D-1, .619830D-1, .445230D-1, .319820
     *D-1,  .000000D0, .108500D-4, .840730D-4, .103690D-2, .496420D-2,  
     *.144460D-1, .307950D-1, .531570D-1, .789950D-1, .105250D0,  .12933
     *0D0, .149580D0, .165290D0, .176450D0, .187100D0,  .186540D0, .1794
     *60D0, .169100D0, .157420D0, .145530D0,  .123140D0, .103730D0, .738
     *590D-1, .532930D-1, .390770D-1,  .291150D-1, .220160D-1, .168710D-
     *1, .130790D-1, .102430D-1,  .644120D-2, .415750D-2, .273630D-2, .1
     *82770D-2 /                                                        
      DATA UPI /  .297510D-3, .261270D-2, .123350D-1, .829510D-1, .25326
     *0D0,  .500520D0, .750720D0, .933490D0, .101620D1, .100340D1,  .920
     *420D0, .796740D0, .657440D0, .519850D0, .284940D0,  .118650D0, .11
     *3970D-1, -.540900D-1, -.925230D-1, -.114180D0,  -.130970D0, -.1319
     *30D0, -.120040D0, -.104530D0, -.897090D-1,  -.765100D-1, -.650540D
     *-1, -.552290D-1, -.468500D-1, -.397260D-1, -.285470D-1, -.205080D-
     *1, -.147310D-1, -.105820D-1,  .742020D-4, .893210D-3, .447550D-2, 
     *.318920D-1, .101210D0,  .205960D0, .314770D0, .393710D0, .424660D0
     *, .408420D0,  .357720D0, .288480D0, .214280D0, .144270D0, .332940D
     *-1,  -.358990D-1, -.731510D-1, -.900910D-1, -.952990D-1, -.941580D
     *-1, -.839730D-1, -.712820D-1, -.493270D-1, -.339400D-1, -.236120D-
     *1, -.166910D-1, -.120020D-1, -.877680D-2, -.652010D-2, -.491360D-2
     *, -.289560D-2, -.177580D-2, -.112270D-2, -.726440D-3 /            
      DATA YNRMS / 0.877580D0, 0.023012780D0 /                          
C     DEFAULT SPEC. AMPLITUDES (NOTE SIGNS):                            
      DATA SPAMS / 0.9671110D0, -0.254354 /                             
      DATA D00 / -100.05650D0 /, ROOT8 / 2.8284271250D0 /               
C                                                                       
C     SKIP TO SETUP, PRINTING, OR CALCULATION CODE.                     
C                                                                       
      IRETUR = 0                                                        
      NUMOUT = 0                                                        
      IL = L/2 + 1                                                      
      UNDEF = FINTRN(1)                                                 
      IF ( IREQUE - 2 ) 100, 200, 300                                   
C                                                                       
C     SETUP:  CHECK FOR ERROR IN L (ONLY LIKELY ERROR).                 
C                                                                       
 100  IF ( L .EQ. 0 .OR. L .EQ. 2 ) GO TO 120                           
      WRITE ( IUNIT, 103 ) ALIAS, L                                     
 103  FORMAT ( '0****** ERROR IN LINKULE ', A8 /  ' ******  L =', I5, ' 
     *  L MUST = 0 OR 2' )                                              
      NUMOUT = 2                                                        
      REWIND IUNIT                                                      
      IRETUR = -1                                                       
      RETURN                                                            
C                                                                       
C     SET JP, NODES, V, R, A, E, J, JSP, JST, SPAM                      
C                                                                       
 120  JP = L + 1                                                        
      INTGER(12) = 0                                                    
C     THE POTENTIAL DEPTH IS ALWAYS SET TO 1.0                          
      FLT(62) = 1.0D0                                                      
      IF ( FLT(43) .EQ. UNDEF ) FLT(43) = 1.0D0                            
      IF ( FLT(1) .EQ. UNDEF ) FLT(1) = 0.4D0                             
      IF ( FLT(12) .EQ. UNDEF ) FLT(12) = -2.224644D0                     
      JBLOCK(1) = 2                                                     
      JBLOCK(8) = 1                                                     
      JBLOCK(9) = 1                                                     
C     THE DEFAULT SPEC. AMPLITUDES GIVE PD = 0.064696                   
      IF ( FLT(53) .EQ. UNDEF ) FLT(53) = SPAMS(IL)                     
C     "SPAMP" OR "SPAMT", AS APPROPRIATE, IS ALWAYS SET TO "SPAM".      
      I = ISWTCH(11)                                                    
      IF ( I .EQ. 1 .OR. I .EQ. 2 ) FLT(53+I) = FLT(53)                 
      RETURN                                                            
C                                                                       
C     IREQUE = 2:  PRINT STUFF ABOUT THE LINKULE.                       
C                                                                       
 200  IF ( IPOTYP .NE. 6 ) GO TO 250                                    
C                                                                       
C     PRINTOUT FOR WAVE FUNCTION CALCULATION.                           
C                                                                       
      WRITE ( IUNIT, 203 ) SD(IL), SAVINT(1)                            
 203  FORMAT ( 10X, 'REID SOFT-CORE DEUTERON WAVE FUNCTION' /  10X, A1, 
     *'-STATE PROBABILITY =', F7.4 )                                    
      NUMOUT = 1                                                        
      IF ( PARAM(1) .EQ. UNDEF ) GO TO 220                              
      WRITE ( IUNIT, 213 ) PARAM(1)                                     
 213  FORMAT ( 10X, 'R SCALE COMPRESSION FACTOR =', F7.4 )              
      NUMOUT = 2                                                        
 220  IF ( IL .EQ. 2 ) GO TO 230                                        
      T = SAVINT(2)/D00                                                 
      WRITE ( IUNIT, 223 ) SAVINT(2), T                                 
 223  FORMAT ( 10X, 'D0 =', F8.2, ' MEV FM**(3/2) =', F7.4, '*D00' )    
      GO TO 270                                                         
 230  WRITE ( IUNIT, 233 ) SAVINT(2)                                    
 233  FORMAT ( 10X, 'D2*D0 =', F9.4, ' MEV FM**(7/2)' )                 
      GO TO 270                                                         
C                                                                       
C     PRINTOUT FOR POTENTIAL CALCULATION                                
C                                                                       
 250  WRITE ( IUNIT, 253 ) FLT(62)                                      
 253  FORMAT ( 10X, 'EFFECTIVE POTENTIAL DERIVED FROM REID SOFT-CORE',  
     *' DEUTERON WAVEFUNCTION.' /  10X, 'MULTIPLIER =', F9.6 )          
C                                                                       
C     PRINTOUT FOR BOTH WF AND POTENTIAL CALCULATIONS.                  
C                                                                       
 270  WRITE (IUNIT, 273 ) FLT(53), YNRMS(IL)                            
 273  FORMAT ( 10X, 'SPECTROSCOPIC AMPLITUDE ("SPAM") =', F8.5 /  10X, '
     *ASYMPTOTIC NORMALIZATION =', F8.5, 'FM**(-1/2)' /  10X, 'REFERENCE
     *:  R. V. REID, ANNALS OF PHYSICS 50 (1968)',  ' 411-448' )        
      NUMOUT = NUMOUT + 5                                               
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
C     SET UP FOR THE MAIN CALCULATION LOOP.                             
C                                                                       
 300  L2SW = L .NE. 0                                                   
      AM = FLT(24)                                                      
      E = FLT(12)                                                       
      PBUGSW = MOD(INTGER(17), 10) .GE. 6                               
      RT4PI = CNSTNT(2)                                                 
      HBARC = CNSTNT(6)                                                 
      ALSQ = -2.0D0*AM*E/HBARC**2                                          
      ALPHA = DSQRT( ALSQ )                                             
C                                                                       
C     THE "PARAM1" KEYWORD GIVES THE RADIAL COMPRESSION FACTOR          
C     (2 MEANS FACTOR-OF-2 COMPRESSION, I.E. SMALL DEUTERON)            
C                                                                       
      SCALX = 1.0D0                                                       
      IF ( PARAM(1) .NE. UNDEF ) SCALX = PARAM(1)                       
      SCALW = DSQRT( SCALX )                                            
      SCALV = SCALW*SCALX                                               
      IF ( L2SW ) SCALV = SCALV*SCALX**2                                
      WFSW = IPOTYP .EQ. 6                                              
      ASSIGN 332 TO JUMPL                                               
      IF ( L2SW ) ASSIGN 334 TO JUMPL                                   
      ASSIGN 370 TO JUMPW                                               
      IF ( WFSW ) ASSIGN 360 TO JUMPW                                   
      ASSIGN 350 TO JUMPP                                               
      IF ( PBUGSW ) ASSIGN 340 TO JUMPP                                 
      IF ( PBUGSW ) WRITE (IUNIT, 313 )                                 
 313  FORMAT ( 6X, 'R', 7X, 'WF', 11X, 'V', 8X, 'V*PHI', 8X,'DINT')
      VMULT = ARRAY2(1)                                                 
      IF ( WFSW ) VMULT = -1.0D0/FLT(62)                                   
      DMULT = FLT(53)*RT4PI*STEPSZ/(VMULT*SPAMS(IL))                    
      IF ( L2SW ) DMULT = DMULT/15.0D0                                     
      VMULT = VMULT*SCALV                                               
      WINT = 0.0D0                                                          
      DINT = 0.0D0                                                          
      RT = RSTART                                                       
      IF ( RT .EQ. 0.0D0 ) RT = 1.D-10                                      
      X = 0.7D0*SCALX*RT                                                  
      XSTEP = 0.7D0*SCALX*STEPSZ                                          
      Y = DEXP( -X )                                                    
      YM = DEXP( -XSTEP )                                               
      IEND = (XX(34)-X)/XSTEP                                           
      IEND = MIN0( IEND, NUMPTS )                                       
      IEND = MAX0( IEND, 0 )                                            
      IF ( IEND .EQ. 0 ) GO TO 400                                      
      WF = 0.0D0                                                            
      POT = 0.0D0                                                           
      VPHI = 0.0D0                                                          
      XXP = XX(1)                                                       
      ASSIGN 338 TO JUMPZ                                               
      II = 1                                                            
C                                                                       
C     LOOP THROUGH RADII, UP TO THE HIGHEST TABULATED VALUE.            
C                                                                       
      DO 399 I = 1, IEND                                                
C                                                                       
C        IF NECESSARY, FIND THE CORRECT INTERVAL AND CALCULATE THE      
C        COEFFICIENTS.                                                  
C                                                                       
      IF ( X .LE. XXP ) GO TO 325                                       
 320  II = II + 1                                                       
      XXP = XX(II)                                                      
      IF ( X .GT. XXP ) GO TO 320                                       
      IM = II - 1                                                       
      XXM = XX(IM)                                                      
      H = XXP - XXM                                                     
      A1 = UI(IM,1)                                                     
      A2 = H*UPI(IM,1)                                                  
      UDIF = UI(II,1) - UI(IM,1)                                        
      A3 = 3.0D0*UDIF - H*( 2.0D0*UPI(IM,1) + UPI(II,1) )                     
      A4 = -2.0D0*UDIF + H*( UPI(IM,1) + UPI(II,1) )                       
      B1 = UI(IM,2)                                                     
      B2 = H*UPI(IM,2)                                                  
      UDIF = UI(II,2) - UI(IM,2)                                        
      B3 = 3.0D0*UDIF - H*(2.0D0*UPI(IM,2) + UPI(II,2) )                     
      B4 = -2.0D0*UDIF + H*( UPI(IM,2) + UPI(II,2) )                       
      ASSIGN 330 TO JUMPZ                                               
C                                                                       
 325  GO TO JUMPZ, ( 330, 338 )                                         
C                                                                       
C        COEFFICIENTS ARE CALCULATED.  CALCULATE THE DESIRED WAVE       
C        FUNCTION WF = EITHER U OR W, AND CALCULATE THE EFFECTIVE       
C        POTENTIAL POT FROM THE CENTRAL, TENSOR, AND SPIN-ORBIT         
C        PARTS (VC, VT, VLS) AND THE WAVE FUNCTION RATIO.               
C                                                                       
 330  P = (X-XXM)/H                                                     
      U = A1 + P*( A2 + P*( A3 + P*A4 ) )                               
      W = B1 + P*( B2 + P*( B3 + P*B4 ) )                               
      Y2 = Y**2                                                         
      Y4 = Y2**2                                                        
      Y6 = Y4*Y2                                                        
      VC = (-10.463D0*Y +105.468D0*Y2-3187.8D0*Y4+9924.3D0*Y6 )/X         
      VT = (-10.463D0*(Y+ 3.0D0*((Y-4.0D0*Y4)+(Y-Y4)/X)/X)+351.77D0*Y4 
     * - 1673.5D0*Y6 )/X                                                   
      GO TO JUMPL, ( 332, 334 )                                         
C                                                                       
C        S-WAVE EFFECTIVE POTENTIAL AND WAVEFUNCTION                    
C                                                                       
 332  WF = U                                                            
      POT = VC + ROOT8*VT*W/U                                           
      GO TO 336                                                         
C                                                                       
C        D-WAVE EFFECTIVE POTENTIAL AND WAVEFUNCTION                    
C                                                                       
 334  WF = W                                                            
      VLS = ( 708.91D0*Y4 - 2713.1D0*Y6 )/X                                 
      POT = VC - 2.0D0*VT - 3.0D0*VLS + ROOT8*VT*U/W                          
C                                                                       
 336  POT = POT*VMULT                                                   
      WF = WF*SCALW                                                     
C                                                                       
C        ACCUMULATE THE INTEGRALS.  WINT IS THE TOTAL                   
C        PROBABILITY INTEGRAL, AND DINT IS THE                          
C        INTEGRAL USED TO CALCULATE D0 OR D2.                           
C                                                                       
      WINT = WINT + STEPSZ*WF*WF                                        
      VPHI = WF*POT                                                     
      DTERM = DMULT*VPHI*RT                                             
      IF ( L2SW ) DTERM = DTERM*RT*RT                                   
      DINT = DINT + DTERM                                               
C                                                                       
C        PRINT IF REQUIRED.                                             
C                                                                       
 338  GO TO JUMPP, ( 340, 350 )                                         
 340  VPHI = VPHI/RT                                                    
      WRITE ( IUNIT, 343 ) RT, WF, POT, VPHI, DINT                      
 343  FORMAT ( 1X, F8.3, 4G12.4 )                                       
C                                                                       
C        STORE THE WAVE FUNCTION AND POTENTIAL (IF IPOTYP = 6)          
C        OR JUST THE POTENTIAL.                                         
C                                                                       
 350  GO TO JUMPW, ( 360, 370 )                                         
 360  ARRAY1(I) = WF                                                    
      ARRAY2(I) = POT                                                   
      GO TO 380                                                         
 370  ARRAY1(I) = POT                                                   
C                                                                       
C        INCREMENT R AND END LOOP.                                      
C                                                                       
 380  RT = RT + STEPSZ                                                  
      X = X + XSTEP                                                     
      Y = Y*YM                                                          
 399  CONTINUE                                                          
C                                                                       
 400  IEND = IEND + 1                                                   
C                                                                       
C     IF WE ARE CALCULATING ONLY THE POTENTIAL, ZERO OUT THE REST       
C     OF IT AND QUIT.                                                   
C                                                                       
      IF ( WFSW ) GO TO 430                                             
      IF ( IEND .GT. NUMPTS ) GO TO 600                                 
      DO 419 I = IEND, NUMPTS                                           
      ARRAY1(I) = 0.0D0                                                     
 419  CONTINUE                                                          
      GO TO 600                                                         
C                                                                       
C     WE ARE CALCULATING WAVE FUNCTIONS.                                
C     FINISH THE INTEGRAL ANALYTICALLY, SINCE THE REST IS JUST AN       
C     EXPONENTIAL (TIMES A POLYNOMIAL IN 1/R FOR THE D STATE).          
C                                                                       
 430  XA = SCALX*ALPHA*RT                                               
      YA = SCALW*YNRMS(IL)*DEXP( -XA )                                  
      G = YA                                                            
      U = 0.5D0*YA**2/(SCALX*ALPHA)                                       
      IF ( L2SW ) GO TO 440                                             
C     S-WAVE                                                            
      ASSIGN 470 TO JUMPL                                               
      GO TO 450                                                         
C     D-WAVE                                                            
 440  T = 1.0D0/XA                                                         
      G = G*( 1.0D0 + T*( 3.0D0 + 3.0D0*T ) ) 
      U = U*( 1.0D0 + T*( 6.0D0 + T*(12.0D0 + T*6.0D0) ) )
      ASSIGN 460 TO JUMPL                                               
C     BOTH S AND D                                                      
 450  WINT = WINT + 0.5D0*G*G*STEPSZ + U                                  
      SAVINT(1) = WINT                                                  
      SAVINT(2) = DINT                                                  
C                                                                       
C     IN THE NEXT LOOP WE CALCULATE FROM THE END OF THE                 
C     TABLES ONWARD, WHERE THE NUCLEAR POTENTIALS ARE NEGLIGIBLE.       
C                                                                       
      IF ( IEND .GT. NUMPTS ) GO TO 500                                 
      XP = SCALX*ALPHA*STEPSZ                                           
      YM = DEXP( -XP )                                                  
      ASSIGN 490 TO JUMPP                                               
      IF ( PBUGSW ) ASSIGN 480 TO JUMPP                                 
C                                                                       
      DO 499 I = IEND, NUMPTS                                           
      WF = YA                                                           
      GO TO JUMPL, ( 460, 470 )                                         
C                                                                       
C        D-STATE MODIFICATIONS                                          
C                                                                       
 460  WF = WF*( 1.0D0 + ( 3.0D0 + 3.0D0/XA )/XA )
      XA = XA + XP                                                      
C                                                                       
 470  GO TO JUMPP, ( 480, 490 )                                         
 480  WRITE ( IUNIT, 343 ) RT, WF                                       
      RT = RT + STEPSZ                                                  
C                                                                       
 490  ARRAY1(I) = WF                                                    
      ARRAY2(I) = 0                                                     
      YA = YA*YM                                                        
C                                                                       
 499  CONTINUE                                                          
C                                                                       
C     THE WAVE FUNCTION CALCULATION IS FINISHED.                        
C     NOW GO BACK AND RENORMALIZE IT SO ITS INTEGRAL IS                 
C     UNITY, AS PTOLEMY EXPECTS.                                        
C                                                                       
 500  T = 1.0D0/DSQRT( WINT )                                              
      DO 509 I = 1, NUMPTS                                              
      ARRAY1(I) = T*ARRAY1(I)                                           
 509  CONTINUE                                                          
C                                                                       
C     FINISH THE OUTPUT IF NECESSARY, AND QUIT.                         
C                                                                       
 600  IF ( .NOT. PBUGSW ) RETURN                                        
      WRITE ( IUNIT, 603 ) WINT, DINT                                   
 603  FORMAT (' WINT =', F10.7, '   DINT =', F11.6 )                   
      NUMOUT = NUMPTS + 2                                               
      REWIND IUNIT                                                      
      RETURN                                                            
      END 
      
      SUBROUTINE SHAPE ( ALIAS, MYINTS, IPOTYP, IREQUE, IRETUR, IUNIT,  
     * NUMOUT,L,J,RSTART,STEPSZ, NUMPTS, ARRAY1, ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER,JBLOCK,ISWTCH, INTRNL, FINTRN,CNSTNT, IWAVBL, 
     * FKANDM, ID, ALLOC, ILLOC, FACFR4, KLOC,LENGTH, NALLOC,NAMLOC)
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
      CHARACTER*4 ID                                                    
C                                                                       
C     SIMPLE LINKULE TO USE A FIXED SHAPE TIMES A DEPTH PARAMETER       
C                                                                       
C     DEPENDING ON THE POTENTIAL TYPE, THE NAME OF THE SHAPE IS         
C     ONE OF REALSHAP, IMAGSHAP, REALSOSH, OR IMAGSOSH.  IF NONE        
C     OF THESE OBJECTS EXISTS, THEN SHAPE IS USED.                      
C     IF THE CORRESPONDING ARRAY REALSCAL, IMAGSCAL, REALSOSC, IMAGSOSC,
C     OR SHAPESCA IS DEFINED THEN IT MUST HAVE TWO ELEMENTS THAT ARE    
C     THE INITIAL AND FINAL COORDINATE VALUES FORWHICH THE POTENTIAL    
C     IS GIVEN.  IN THIS CASE THE ARRAY WILL BE INTERPOLATED IF NECESSAR
C     NECESSARY TO THE REQUIRED GRID.                                   
C     ALTERNATIVELY THE SCALE ARRAY MAY HAVE AS MANY ELEMENTS           
C     AS THE SHAPE ARRAY.  IN THIS CASE IT GIVE THE R VALUES            
C     (POSSIBLY UNEQUALLY SPACED) OF THE ELEMENTS IN SHAPE.             
C     IF THERE IS NO SCALE ARRAY TO BE DEFINED ON THE CORRECT GRID.     
C                                                                       
C     12/8/77 - FIRST VERSION - S.P.                                    
C     3/15/78 - ALLOW R GRID TO BE SPECIFIED - S.P.                     
C     1/16/80 - STANDARD MORTRAN - RPG                                  
C     4/14/82 - EXTRAPOLATE FOR LARGE R - S.P.                          
C                                                                       
      CHARACTER*8 NAMES(5)                                              
      DATA NAMES/'REALSHAP','IMAGSHAP','REALSOSH','IMAGSOSH','SHAPE'/
      CHARACTER*8 SCALNM(5)                                             
      DATA SCALNM/'REALSCAL', 'IMAGSCAL', 'REALSOSC', 'IMAGSOSC', 
     *'SHAPESCA'/                                                             
C                                                                       
      CHARACTER*4 WORKNM(2,2)                                           
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                       
      REAL*8 WORK(20)                                                   
      LOGICAL RSW                                                       
C                                                                       
      SMLNUM = CNSTNT(12)                                               
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      CONS = TEMPVS(IPOTYP)                                             
      IF ( CONS .NE. FINTRN(1) ) GO TO 45                               
      WRITE ( IUNIT, 43 )                                               
 43   FORMAT ( '0**** THE WELL DEPTH PARAMETER (V, VI, VSO, VSOI)', 
     * ' MUST BE DEFINED FOR THE SHAPE LINKULE.' )                          
      IRETUR = -1                                                       
      GO TO 900                                                         
 45   CONS = -CONS                                                      
C                                                                       
      GO TO ( 100, 500, 600 ), IREQUE                                   
C                                                                       
C     LOCATE THE SHAPE                                                  
C                                                                       
 100  II = IPOTYP                                                       
      I = NAMLOC(NAMES(II))                                             
      IF ( I .NE. 0 ) GO TO 150                                         
      II = 5                                                            
      I = NAMLOC(NAMES(II))                                             
      IF ( I .NE. 0 ) GO TO 150                                         
      WRITE ( IUNIT, 123 ) NAMES(IPOTYP)                                
 123  FORMAT ('0**** AN OBJECT WITH THE NAME ', A8,  ' OR SHAPE MUST BE 
     *DEFINED TO USE THE SHAPE LINKULE.' )                              
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
C     IS THE SCALING ARRAY DEFINED                                      
C                                                                       
 150  ISCAL = NAMLOC( SCALNM(II) )                                      
      RMAX = RSTART + (NUMPTS-1)*STEPSZ                                 
      NIN = LENGTH(I)                                                   
      IF ( ISCAL .NE. 0 ) GO TO 200                                     
C                                                                       
C     NO, ALL WE CAN DO IS CHECK THE NUMBER OF POINTS                   
C                                                                       
      IF ( NUMPTS .EQ. NIN ) GO TO 180                                  
      WRITE ( IUNIT, 163 ) RMAX, STEPSZ, NUMPTS, NAMES(II), NIN,  SCALNM
     *(II)                                                              
 163  FORMAT ( '0**** ASYMPTOPIA AND STEPSIZE =', 2G15.5,  ' AND REQUIRE
     *', I6, ' POINTS.' /  6X, 'HOWEVER, ', A8, ' HAS', I6, ' POINTS AND
     * ',  A8, ' IS NOT DEFINED.' )                                     
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 180  WRITE ( IUNIT, 183 ) NAMES(II), RSTART, RMAX, STEPSZ              
 183  FORMAT ( X, A8, ' IS BEING ASSUMED TO BE DEFINED FOR',  G15.5, '
     * =< R =<', G15.5, ' WITH STEPSIZE =', G15.5 )                     
      NUMOUT = 1                                                        
      GO TO 400                                                         
C                                                                       
C     SCALLING ARRAY IS DEFINED, WILL INTERPOLATION BE REQUIRED         
C                                                                       
 200  IF ( LENGTH(ISCAL) .EQ. 2 ) GO TO 210                             
      IF ( LENGTH(ISCAL) .EQ. NIN ) GO TO 220                           
      WRITE ( IUNIT, 203 ) SCALNM(II), NIN, LENGTH(ISCAL)               
 203  FORMAT ( '0**** SCALLING ARRAY ', A8, ' SHOULD HAVE 2 OR',  I4, ' 
     *ELEMENTS,',  ' BUT IT HAS', I6, ' ELEMENTS.' )                    
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
C     SCALLING ARRAY SPECIFIES START, STOP                              
C                                                                       
 210  R1 = ALLOC(KLOC(ISCAL))                                            
      R2 = ALLOC(KLOC(ISCAL)+1)                                          
      RSW = .FALSE.                                                     
      IF ( R1.EQ.RSTART.AND.NIN.EQ.NUMPTS.AND.DABS(R2-RMAX)/STEPSZ
     * .LT. 1.D-5 ) GO TO 400                                         
      GO TO 240                                                         
C                                                                       
C     SCALLING SPECIFIES COMPLETE R-GRID                                
C                                                                       
 220  R1 = ALLOC(KLOC(ISCAL))                                            
      R2 = ALLOC(KLOC(ISCAL)+NIN-1)                                      
      RSW = .TRUE.                                                      
C                                                                       
C     MUST RESCALE                                                      
C                                                                       
 240  IOLD = I                                                          
      WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
      ISIZE = NIN+NUMPTS                                                
      I = NALLOC( NUMPTS, WORKNM(1,1), KLOC )                            
      IWORK = NALLOC( ISIZE, WORKNM(1,2), KLOC )                         
      LXIN = KLOC(IWORK)-1                                               
      LXOUT = LXIN+NIN                                                  
      LYIN = KLOC(IOLD)-1                                                
      LYOUT = KLOC(I)-1                                                  
      IF ( RSW ) LXIN = KLOC(ISCAL) - 1                                  
      IF ( RSW ) GO TO 260                                              
C                                                                       
C     SETUP INPUT X'S                                                   
C                                                                       
      R = R1                                                            
      STEPIN = (R2-R1)/(NIN-1)                                          
      DO 259 N = 1, NIN                                                 
      ALLOC(LXIN+N) = R                                                 
      R = R+STEPIN                                                      
 259  CONTINUE                                                          
C                                                                       
C     GET EXTRAPOLATING FUNCTION                                        
C                                                                       
 260  VSTEP = 0.0D0                                                         
      AVAL = 0.0D0                                                          
      YNM1 = ALLOC(LYIN+NIN-1)                                          
      YN = ALLOC(LYIN+NIN)                                              
      IF ( YN*YNM1 .LE. 0.0D0 ) GO TO 270                                   
      AVAL = YN                                                         
      BVAL = DLOG(YNM1/AVAL) / (R2-ALLOC(LXIN+NIN-1))                   
      VSTEP = DEXP(-BVAL*STEPSZ)                                        
C                                                                       
C     SETUP OUTPUT X'S,  USE                                            
C     F(1) FOR R < R1  AND  A*EXP(-B*R) FOR R > R2 .                    
C                                                                       
 270  R = RSTART                                                        
      NSTRT = 1                                                         
      VAL0 = ALLOC(LYIN+1)                                              
      DO 299 N = 1, NUMPTS                                              
      IF ( R .GE. R1 ) GO TO 280                                        
      ALLOC(LYOUT+N) = VAL0                                             
      NSTRT = N+1                                                       
      GO TO 290                                                         
 280  IF ( R .LE. R2 ) GO TO 285                                        
      AVAL = AVAL*VSTEP                                                 
      IF ( DABS(AVAL) .LT. SMLNUM ) AVAL = 0                            
      ALLOC(LYOUT+N) = AVAL                                             
      GO TO 290                                                         
 285  ALLOC(LXOUT+N) = R                                                
      NLAST = N                                                         
 290  R = R + STEPSZ                                                    
 299  CONTINUE                                                          
      NOUT = NLAST+1-NSTRT                                              
C                                                                       
C     DO THE INTERPOLATION                                              
C                                                                       
      CALL AITKEN ( 5,0.0D0,0.0D0,NIN,ALLOC(LXIN+1), ALLOC(LYIN+1), 
     * NOUT,ALLOC(LXOUT+NSTRT),ALLOC(LYOUT+NSTRT),NFAIL,WORST,WORK)
C                                                                       
C     FREE WORK AREA                                                    
C                                                                       
      KLOC(IWORK) = -KLOC(IWORK)                                          
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = I                                                     
      MYINTS(2) = II                                                    
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503) NAMES(MYINTS(2)), CONS                         
 503  FORMAT ( ' THE POTENTIAL SHAPE STORED IN OBJECT ', A8,  ' IS BEING
     * MULTIPLIED BY', G15.5 )                                          
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
 600  LL = KLOC(MYINTS(1)) - 1                                           
      DO 629 I = 1, NUMPTS                                              
      ARRAY1(I) = CONS*ALLOC(LL+I)                                      
 629  CONTINUE                                                          
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
      END  
      
      SUBROUTINE SPLINE ( ALIAS, MYINTS,IPOTYP,IREQUE,IRETUR,IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ,NUMPTS,ARRAY1,ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER,JBLOCK,ISWTCH, INTRNL, FINTRN, CNSTNT,IWAVBL, 
     * FKANDM,ID, ALLOC,ILLOC,FACFR4, KLOC,LENGTH,NALLOC,NAMLOC)
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(2,10),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
      CHARACTER*4 ID                                                    
C                                                                       
C     SPLINE - CONTRUCT POTENTIAL BY SPLINES                            
C                                                                       
C     THIS LINKULE ASSUMES THAT THE PAIRS OF NUMBERS:                   
C       (PARAM1, PARAM2), (PARAM3, PARAM4), ..., (PARAM19, PARAM20)     
C     ARE TO BE INTERPERTED AS                                          
C       (R1, V1), (R2, V2), ..., (R10, V10)                             
C     AND DEFINE A POTENTIAL V(R).  FEWER THAN 10 POINTS MAY            
C     BE USED BY LEAVING THE REMAINING V'S OR R'S UNDEFINED.            
C     THE R'S DO NOT HAVE TO BE ORDERED.                                
C     THE POTENTIAL IS DEFINED AS                                       
C                                                                       
C       V(R) =  V1    R < R1 ,                                          
C               V'(R)   R1 < R < RLAST                                  
C               0       R > RLAST                                       
C                                                                       
C     WHERE V'(R) IS THE RESULT OF CUBIC SPLINE INTERPOLATION           
C     OF THE LOG OF THE POTENTIAL.                                      
C                                                                       
C     4/6/78 - BASED ON SHAPE - S.P.                                    
C     1/16/80 - STANDARD MORTRAN - RPG                                  
C                                                                       
      CHARACTER*4 WORKNM(2,2)                                           
      DATA WORKNM/'*000', 'SHAP', '*000', 'WORK'/                       
      DIMENSION R(20), V(20), B(20), C(20), D(20)                       
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      UNDEF = FINTRN(1)                                                 
C                                                                       
C     GET THE DEFINED POINTS AND VALUES AND ORDER THEM                  
C                                                                       
      NV = 1                                                            
      IV = 1                                                            
  40  R(1) = PARAM(1,IV)                                                
      V(1) = PARAM(2,IV)                                                
      IF ( R(1) .NE. UNDEF .AND. V(1) .NE. UNDEF ) GO TO 50             
      IV = IV+1                                                         
      IF ( IV .LT. 10 ) GO TO 40                                        
      GO TO 950                                                         
C                                                                       
 50   IV = IV + 1                                                       
      IF ( IV .GT. 10 ) GO TO 90                                        
      IF (PARAM(1,IV).EQ.UNDEF.OR.PARAM(2,IV).EQ.UNDEF) GO TO 50
      DO 59 I = 1, NV                                                   
      IF ( PARAM(1,IV) .LT. R(I) ) GO TO 70                             
 59   CONTINUE                                                          
C                                                                       
C     ADD TO END                                                        
C                                                                       
      NV = NV+1                                                         
      R(NV) = PARAM(1,IV)                                               
      V(NV) = PARAM(2,IV)                                               
      GO TO 50                                                          
C                                                                       
C     INSERT POINT; FIRST MOVE REST RIGHT ONE                           
C                                                                       
 70   DO 74 II = I, NV                                                  
      R(NV+I-II+1) = R(NV+I-II)                                         
      V(NV+I-II+1) = V(NV+I-II)                                         
 74   CONTINUE                                                          
      NV = NV+1                                                         
      R(I) = PARAM(1,IV)                                                
      V(I) = PARAM(2,IV)                                                
      GO TO 50                                                          
C                                                                       
C     ALL DEFINITIONS FOUND AND ORDERED                                 
C                                                                       
 90   IF ( NV .LE. 2 ) GO TO 950                                        
C                                                                       
      GO TO ( 100, 500, 600 ), IREQUE                                   
C                                                                       
C     GENERATE RADIUS ARRAY FOR INTRPC                                  
C                                                                       
 100  WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
      I = NALLOC( NUMPTS, WORKNM(1,1), KLOC )                            
C                                                                       
      RR = RSTART                                                       
      DO 159 II = 1, NUMPTS                                             
      ALLOC(KLOC(I)-1+II) = RR                                           
      RR = RR + STEPSZ                                                  
 159  CONTINUE                                                          
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = I                                                     
      MYINTS(2) = II                                                    
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT, 503) ( R(I), V(I), I = 1, NV )                      
 503  FORMAT ( ' THE POTENTIAL IS DEFINED BY THE POINTS:' /  ( 5( F10.3,
     * G13.5 ) ) )                                                      
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
 600  LL = KLOC(MYINTS(1)) - 1                                           
C                                                                       
C     SETUP OUTPUT X'S, DO NOT EXTRAPOLATE; RATHER USE                  
C     F(1) FOR R < R1  AND  EXPONENTIAL DECAY FOR R > R2 .              
C                                                                       
 270  RR = RSTART                                                       
      NSTRT = 1                                                         
      DO 299 N = 1, NUMPTS                                              
      IF ( RR .GE. R(1) ) GO TO 280                                     
      ARRAY1(N) = -V(1)                                                 
      NSTRT = N+1                                                       
      GO TO 290                                                         
 280  IF ( RR .LE. R(NV) ) GO TO 285                                    
      ARRAY1(N) = 0.0D0
      GO TO 290                                                         
 285  NLAST = N                                                         
 290  RR = RR + STEPSZ                                                  
 299  CONTINUE                                                          
      NOUT = NLAST+1-NSTRT                                              
C                                                                       
C     INTERPOLATE THE LOGS                                              
C                                                                       
      DO 719 I = 1, NV                                                  
      V(I) = DLOG(V(I))                                                 
 719  CONTINUE                                                          
      CALL SPLNCB ( NV, R, V, B, C, D )                                 
      CALL INTRPC (NV,R,V,B,C,D,NOUT,ALLOC(LL+NSTRT),ARRAY1(NSTRT))
      DO 749 I = NSTRT, NLAST                                           
      ARRAY1(I) = -DEXP(ARRAY1(I))                                      
 749  CONTINUE                                                          
C                                                                       
C     FILL IN THE END WITH AN EXPONENTIAL DECAY                         
C                                                                       
      TERM = ARRAY1(NLAST)                                              
      NLAST = NLAST+1                                                   
      IF ( NLAST .GT. NUMPTS ) RETURN                                   
      STEP = DEXP( B(NV)*STEPSZ )                                       
      TERM = TERM*STEP                                                  
      DO 769 I = NLAST, NUMPTS                                          
      ARRAY1(I) = TERM                                                  
      IF ( TERM .GT. -1.D-30 ) RETURN                                   
      TERM = TERM*STEP                                                  
 769  CONTINUE                                                          
C                                                                       
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
C                                                                       
 950  WRITE (IUNIT, 953)                                                
 953  FORMAT ( '0**** POTENTIAL MUST BE DEFINED ON AT LEAST', 
     * ' TWO POINTS FOR SPLINE.' ) 
      GO TO 900                                                         
      END       
      
      SUBROUTINE TWOSHA ( ALIAS, MYINTS, IPOTYP,IREQUE,IRETUR,IUNIT, 
     * NUMOUT,L,J,RSTART,STEPSZ, NUMPTS, ARRAY1,ARRAY2,FLOATK,TEMPVS, 
     * PARAM, INTGER,JBLOCK,ISWTCH,INTRNL, FINTRN,CNSTNT,IWAVBL, 
     * FKANDM,ID, ALLOC,ILLOC,FACFR4,KLOC,LENGTH,NALLOC,NAMLOC)
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      DIMENSION MYINTS(2), ARRAY1(1), ARRAY2(1), FLOATK(200), 
     *TEMPVS(6),PARAM(20),INTGER(50),JBLOCK(12),ISWTCH(30),INTRNL(74), 
     *FINTRN(37),CNSTNT(13),IWAVBL(240),FKANDM(26),ALLOC(1),ILLOC(1), 
     *KLOC(1),LENGTH(1)
      INTEGER FACFR4                                                    
      CHARACTER*4 ID                                                    
C                                                                       
C     TWOSHAPE - USE TWO FIXED SHAPES TIMES DEPTH PARAMETERS            
C                                                                       
C     DEPENDING ON THE POTENTIAL TYPE, THE NAMES OF THE SHAPES ARE      
C     ONE OF (REALSHAP, REALADD), (IMAGSHAP, IMAGADD),                  
C     (REALSOSH, REALSOAD), OR (IMAGSOSH, IMAGSOAD).                    
C     IF NONE OF THESE OBJECTS EXISTS, THEN (SHAPE, SHAPEADD) IS USED.  
C     IF THE CORRESPONDING ARRAY REALSCAL, IMAGSCAL, REALSOSC, IMAGSOSC,
C     OR SHAPESCA IS DEFINED THEN IT MUST HAVE TWO ELEMENTS THAT ARE    
C     THE INITIAL AND FINAL COORDINATE VALUES FORWHICH THE SHAPES ARE   
C     GIVEN.  IN THIS CASE THE ARRAYS WILL BE INTERPOLATED IF NECESSARY 
C     NECESSARY TO THE REQUIRED GRID.                                   
C     NOTE THAT BOTH SHAPES MUST BE GIVEN ON THE SAME GRID.             
C     ALTERNATIVELY THE SCALE ARRAY MAY HAVE AS MANY ELEMENTS           
C     AS THE SHAPE ARRAY.  IN THIS CASE IT GIVE THE R VALUES            
C     (POSSIBLY UNEQUALLY SPACED) OF THE ELEMENTS IN SHAPE.             
C     IF THERE IS NO SCALE ARRAY TO BE DEFINED ON THE CORRECT GRID.     
C     THE POTENTIAL USED IS                                             
C                                                                       
C   -V*REALSHAP - PARAM1*REALADD ,                                      
C       -VI*IMAGSHAP - PARAM2*IMAGADD ,                                 
C      -VSO*REALSOSH - PARAM3*REALSOAD , OR                             
C     -VSOI*IMAGSOSH - PARAM4*IMAGSOAD .                                
C                                                                       
C     12/8/77 - FIRST VERSION - S.P.                                    
C     3/15/78 - ALLOW R GRID TO BE SPECIFIED - S.P.                     
C     3/25/78 - TSOSHAPE MADE FROM SHAPE - S.P.                         
C     1/16/80 - STANDARD MORTRAN - RPG                                  
C                                                                       
      CHARACTER*8 NAMES(5)                                              
      DATA NAMES/'REALSHAP','IMAGSHAP','REALSOSH','IMAGSOSH','SHAPE'/
      CHARACTER*8 NAMES2(5)                                             
      DATA NAMES2/'REALADD', 'IMAGADD', 'REALSOAD', 'IMAGSOAD', 
     *'SHAPEADD'/                                                               
      CHARACTER*8 SCALNM(5)                                             
      DATA SCALNM/'REALSCAL', 'IMAGSCAL', 'REALSOSC', 'IMAGSOSC', 
     *'SHAPESCA'/                                                             
C                                                                       
      CHARACTER*4 WORKNM(2,3)                                           
      DATA WORKNM/'*000', 'SHAP', '*000', 'SHP2' , '*000', 'WORK'/      
      REAL*8 WORK(20)                                                   
      LOGICAL RSW                                                       
C                                                                       
      NUMOUT = 0                                                        
      IRETUR = 0                                                        
      CONS = TEMPVS(IPOTYP)                                             
      CONS2 = PARAM(IPOTYP)                                             
      IF ( CONS .NE. FINTRN(1) .AND. CONS2 .NE. FINTRN(1) ) GO TO 45    
      WRITE ( IUNIT, 43 )                                               
 43   FORMAT ( '0**** THE WELL DEPTH PARAMETER (V, VI, VSO, VSOI)', 
     * ' AND CORRESPONDINGLY PARAM1 - PARAM4',  
     * ' MUST BE DEFINED FOR THE TWOSHAPE LINKULE.' )
      IRETUR = -1                                                       
      GO TO 900                                                         
 45   CONS = -CONS                                                      
      CONS2 = -CONS2                                                    
C                                                                       
      GO TO ( 100, 500, 600 ), IREQUE                                   
C                                                                       
C     LOCATE THE SHAPE                                                  
C                                                                       
 100  II = IPOTYP                                                       
      I = NAMLOC(NAMES(II))                                             
      IF ( I .NE. 0 ) GO TO 130                                         
      II = 5                                                            
      I = NAMLOC(NAMES(II))                                             
      IF ( I .NE. 0 ) GO TO 130                                         
      WRITE ( IUNIT, 123 ) NAMES(IPOTYP)                                
 123  FORMAT ('0**** AN OBJECT WITH THE NAME ', A8,  ' OR TWOSHAPE MUST 
     *BE DEFINED TO USE THE SHAPE LINKULE.' )                           
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
C     NOW LOCATE SECOND OBJECT                                          
C                                                                       
 130  I2 = NAMLOC(NAMES2(II))                                           
      IF (I2 .NE. 0 ) GO TO 140                                         
      WRITE ( IUNIT, 133 ) NAMES2(II)                                   
 133  FORMAT ( '0**** AN OBJECT WITH THE NAME ', A8,  ' MUST BE DEFINED 
     *FOR THE TWOSHAPE LINKULE.' )                                      
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
C     CHECK THE TWO ARRAYS FOR COMPATABILITY                            
C                                                                       
 140  NIN = LENGTH(I)                                                   
      IF ( NIN .EQ. LENGTH(I2) ) GO TO 150                              
      WRITE ( IUNIT, 143 ) NAMES(II), NAMES2(II)                        
 143  FORMAT ( '0**** ', A8, ' AND ', A8, ' MUST BE DEFINED ON THE SAME
     * GRID.' )                                                         
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
C     IS THE SCALING ARRAY DEFINED                                      
C                                                                       
 150  ISCAL = NAMLOC( SCALNM(II) )                                      
      RMAX = RSTART + (NUMPTS-1)*STEPSZ                                 
      IF ( ISCAL .NE. 0 ) GO TO 200                                     
C                                                                       
C     NO, ALL WE CAN DO IS CHECK THE NUMBER OF POINTS                   
C                                                                       
      IF ( NUMPTS .EQ. NIN ) GO TO 180                                  
      WRITE (IUNIT,163) RMAX,STEPSZ,NUMPTS,NAMES(II),NIN,SCALNM(II)
 163  FORMAT ( '0**** ASYMPTOPIA AND STEPSIZE =', 2G15.5,  ' AND REQUIRE
     *', I6, ' POINTS.' /  6X, 'HOWEVER, ', A8, ' HAS', I6, ' POINTS AND
     * ',  A8, ' IS NOT DEFINED.' )                                     
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
 180  WRITE ( IUNIT, 183 ) NAMES(II), NAMES2(II), RSTART, RMAX, STEPSZ  
 183  FORMAT ( X, A8, ' AND ', A8, ' ARE BEING ASSUMED TO BE DEFINED',
     * ' FOR', G15.5, ' =< R =<', G15.5, ' WITH STEPSIZE =', G15.5 )    
      NUMOUT = 1                                                        
      GO TO 400                                                         
C                                                                       
C     SCALLING ARRAY IS DEFINED, WILL INTERPOLATION BE REQUIRED         
C                                                                       
 200  IF ( LENGTH(ISCAL) .EQ. 2 ) GO TO 210                             
      IF ( LENGTH(ISCAL) .EQ. NIN ) GO TO 220                           
      WRITE ( IUNIT, 203 ) SCALNM(II), NIN, LENGTH(ISCAL)               
 203  FORMAT ( '0**** SCALLING ARRAY ', A8, ' SHOULD HAVE 2 OR',  I4, ' 
     *ELEMENTS,',  ' BUT IT HAS', I6, ' ELEMENTS.' )                    
      IRETUR = -1                                                       
      GO TO 900                                                         
C                                                                       
C     SCALLING ARRAY SPECIFIES START, STOP                              
C                                                                       
 210  R1 = ALLOC(KLOC(ISCAL))                                            
      R2 = ALLOC(KLOC(ISCAL)+1)                                          
      RSW = .FALSE.                                                     
      IF ( R1.EQ.RSTART .AND. NIN.EQ.NUMPTS .AND. 
     * DABS(R2-RMAX)/STEPSZ.LT.1.D-5 ) GO TO 400
      GO TO 240                                                         
C                                                                       
C     SCALLING SPECIFIES COMPLETE R-GRID                                
C                                                                       
 220  R1 = ALLOC(KLOC(ISCAL))                                            
      R2 = ALLOC(KLOC(ISCAL)+NIN-1)                                      
      RSW = .TRUE.                                                      
C                                                                       
C     MUST RESCALE                                                      
C                                                                       
 240  IOLD = I                                                          
      IOLD2 = I2                                                        
      WORKNM(1,1) = ID                                                  
      WORKNM(1,2) = ID                                                  
      WORKNM(1,3) = ID                                                  
      ISIZE = NIN+NUMPTS                                                
      I = NALLOC( NUMPTS, WORKNM(1,1), KLOC )                            
      I2 = NALLOC( NUMPTS, WORKNM(1,2), KLOC )                           
      IWORK = NALLOC( ISIZE, WORKNM(1,3), KLOC )                         
      LXIN = KLOC(IWORK)-1                                               
      LXOUT = LXIN+NIN                                                  
      LYIN = KLOC(IOLD)-1                                                
      LYOUT = KLOC(I)-1                                                  
      LYIN2 = KLOC(IOLD2)-1                                              
      LYOUT2 = KLOC(I2)-1                                                
      IF ( RSW ) LXIN = KLOC(ISCAL) - 1                                  
      IF ( RSW ) GO TO 270                                              
C                                                                       
C     SETUP INPUT X'S                                                   
C                                                                       
      R = R1                                                            
      STEPIN = (R2-R1)/(NIN-1)                                          
      DO 269 N = 1, NIN                                                 
      ALLOC(LXIN+N) = R                                                 
      R = R+STEPIN                                                      
 269  CONTINUE                                                          
C                                                                       
C     SETUP OUTPUT X'S, DO NOT EXTRAPOLATE; RATHER USE                  
C     F(1) FOR R < R1  AND  0 FOR R > R2 .                              
C                                                                       
 270  R = RSTART                                                        
      NSTRT = 1                                                         
      VAL0 = ALLOC(LYIN+1)                                              
      VAL02 = ALLOC(LYIN2+1)                                            
      DO 299 N = 1, NUMPTS                                              
      IF ( R .GE. R1 ) GO TO 280                                        
      ALLOC(LYOUT+N) = VAL0                                             
      ALLOC(LYOUT2+N) = VAL02                                           
      NSTRT = N+1                                                       
      GO TO 290                                                         
 280  IF ( R .LE. R2 ) GO TO 285                                        
      ALLOC(LYOUT+N) = 0                                                
      ALLOC(LYOUT2+N) = 0                                               
      GO TO 290                                                         
 285  ALLOC(LXOUT+N) = R                                                
      NLAST = N                                                         
 290  R = R + STEPSZ                                                    
 299  CONTINUE                                                          
      NOUT = NLAST+1-NSTRT                                              
C                                                                       
C     DO THE INTERPOLATION                                              
C                                                                       
      CALL AITKEN (5,0.0D0,0.0D0,NIN,ALLOC(LXIN+1),ALLOC(LYIN+1), 
     * NOUT,ALLOC(LXOUT+NSTRT),ALLOC(LYOUT+NSTRT),NFAIL,WORST,WORK)
C                                                                       
      CALL AITKEN ( 5, 0.0D0, 0.0D0, NIN, ALLOC(LXIN+1),ALLOC(LYIN2+1),
     * NOUT,ALLOC(LXOUT+NSTRT),ALLOC(LYOUT2+NSTRT),NFAIL,WORST,WORK)
C                                                                       
C     FREE WORK AREA                                                    
C                                                                       
      KLOC(IWORK) = -KLOC(IWORK)                                          
C                                                                       
C     SAVE ADDRESS OF SHAPE ARRAY                                       
C                                                                       
 400  MYINTS(1) = 256*I + I2                                            
      MYINTS(2) = II                                                    
      IF ( NUMOUT .NE. 0 ) GO TO 900                                    
      RETURN                                                            
C                                                                       
C     PRINT IT OUT                                                      
C                                                                       
 500  WRITE (IUNIT,503) CONS,NAMES(MYINTS(2)),CONS2,NAMES2(MYINTS(2))
 503  FORMAT ( ' THE POTENTIAL IS', T20, G15.5,  ' TIMES THE SHAPE IN ',
     * A8 /  T16, 'PLUS', G15.5, ' TIMES THE SHAPE IN ', A8 )           
      GO TO 900                                                         
C                                                                       
C     DO IT                                                             
C                                                                       
 600  I = MYINTS(1)/256                                                 
      I2 = MYINTS(1) - 256*I                                            
      LL = KLOC(I) - 1                                                   
      LL2 = KLOC(I2) - 1                                                 
      DO 629 I = 1, NUMPTS                                              
      ARRAY1(I) = CONS*ALLOC(LL+I) + CONS2*ALLOC(LL2+I)                 
 629  CONTINUE                                                          
      RETURN                                                            
C                                                                       
 900  NUMOUT = 1                                                        
      REWIND IUNIT                                                      
      RETURN                                                            
      END     
      
C--- SUBROUTINE TO PLOT THEORECTICAL CROSS-SECTIONS ON A LINEPRINTER
C---		     J S WINFIELD  MARCH 1981
C---
	SUBROUTINE LPPLOT(ANG,XSECT,NUMPLT)
C
	CHARACTER *120 XAXIS,LINE
	CHARACTER *9 DATEN,TIMEN
	CHARACTER*1 REACT, HEADER
	COMMON /HEDCOM/ REACT(45), HEADER(65)

	DIMENSION ANG(1),XSECT(1)
	XBIN(XXY) = (ALOG10(XXY)-FLOAT(MINPWR))*10.*FKNTRL+1.5
C--
	CALL DATE(DATEN)
	CALL TIME(TIMEN)
C
C** FIND MAXPWR & MINPWR
C
	FMIN=30.
	FMAX=-30.
	DO 60 I=1,NUMPLT
	IF(XSECT(I).EQ.0.0 ) XSECT(I) = 1.E-30
	TEMP=ALOG10(XSECT(I))
	IF(TEMP.GT.FMAX)  FMAX=TEMP
	IF(TEMP.LT.FMIN)  FMIN=TEMP
  60	CONTINUE
	MAXPWR=IFIX(FMAX)+1
	IF(FMAX.LE.0.) MAXPWR=FMAX
	MINPWR=FMIN
	IF(FMIN.LT.0.) MINPWR=IFIX(FMIN)-1
	NCYCL=MAXPWR-MINPWR
	CLOSE(UNIT=1)
	FKNTRL=10./FLOAT(NCYCL)
	WRITE(6,80) REACT, HEADER
  80  FORMAT ( X, T56, 'P T O L E M Y' / 10X,'COMPUTATION OF CROSS S
     1ECTIONS',  T78, 'SOPHISTICATED PROGRAMS CONTAIN SUBTLE ERRORS'/
     2X, 45A1, 9X, 65A1 / )
	IF(MINPWR.GT.-10) WRITE(6,90)MINPWR,MAXPWR
	IF(MINPWR.LE.-10) WRITE(6,100)MINPWR,MAXPWR
  90	FORMAT(2X,'THETA',4X,'SIGMA',3X,
	1 '10**'I2,39X,13Hmillibarns/sr,39X,'10**'I2,/)
 100	FORMAT(2X,'THETA',4X,'SIGMA',3X,'10**'I3,90X,'10**'I3,/)
C
C**  SET UP INITIAL LINE ...
C
	DO 110 J=1,101
 110	XAXIS(J:J)='-'
	DO 120 J=MINPWR,MAXPWR
	KY=XBIN(10.**J)
 120	XAXIS(KY:KY)='I'
	KY=XBIN(XSECT(1))
	XAXIS(KY:KY)='*'
	WRITE(6,130)ANG(1),XSECT(1),XAXIS
 130	FORMAT(F7.2,1PE11.3,2X,A101)
C
C** NOW CLEAR LINE ( EXCEPT X-AXIS ) AND START PLOTTING !
C
	XAXIS(KY:KY)='-'
	DO 140 J=2,100
 140	LINE(J:J)=' '
	LINE(1:1)='I'
	LINE(101:101)='I'
	DO 160 J=2,NUMPLT
	Y=XSECT(J)
	KY = XBIN(Y)
	LINE(KY:KY)='*'
	WRITE(6,130)ANG(J),XSECT(J),LINE
 150	FORMAT(20X,A101)
 160	LINE(KY:KY)=' '
	WRITE(6,150)XAXIS
	WRITE(6,165)DATEN,TIMEN
 165	FORMAT(//8X'DATE ',A9,'  TIME ',A9)
	END
