C./ ADD NAME=SSU
C     SUBROUTINE TO TEST ACCURACY OF ROOTS OF CUBIC EQUATION,
C
C     SUBROUTINE ROOTST ( K, KAPPA, H )
C     IMPLICIT REAL*8 (A-H,O-Z)
C     REAL*8 K, K2, K3, KAPPA
C
C     K2 = K*K
C     K3 = K2*K
C     T1 = DABS(2.*K3)
C     T2 = H*K2
C     DT = DABS(T1-DABS(KAPPA-T2))
C     T = T1
C     IF (DABS(KAPPA).GT.T) T = DABS(KAPPA)
C     IF (DABS(T2).GT.T) T = DABS(T2)
C     DT = DT/T
C     IF (DT.LT..0001) RETURN
C     PRINT 1000, KAPPA, H, K, DT
C1000 FORMAT (' ERROR IN ROOT, KAPPA = ',E12.4,' H = ',E12.4,'
C    * K = ',E12.4,' DT = ',E12.4 )
C     END
C
C     SUBROUTINE TO CALCULATE THE COMPLETE ELLIPTIC INTEGRAL
C     OF THE THIRD KIND,
C
      SUBROUTINE CEL3( CE, KP, M, FP )
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 KC, M, M0, KP
      KC = KP
      P = FP
      IF (KC*P.EQ.0.D0) GO TO 10
      KC = DABS(KC)
      E = KC
      M0 = 1.D0
      IF (P.LE.0.D0) GO TO 2
      C = 1.D0
      P = DSQRT(P)
      D = 1.D0/P
      GO TO 3
 2    G = 1.D0 -P
      F = KC*KC - P
      P = DSQRT(F/G)
      D = -M/(G*P)
      C = 0.D0
 3    F = C
      C = C + D/P
      G = E/P
      D = 2.D0 * ( F*G + D )
      P = G + P
      G = M0
      M0 = KC + M0
      IF ( DABS(G-KC).LE.1.D-4*G ) GO TO 4
      KC = 2.D0 * DSQRT(E)
      E = KC * M0
      GO TO 3
 4    CE = 1.5707963D0 * (C*M0+D)/(M0*(M0+P))
      RETURN
 10   CE = (1.D18)**2
      RETURN
      END
C
      SUBROUTINE SS( E, DE, S, DELTA, EC, BET )
C     SIGNALS
C     DLE = .0 AND DE NOT= 0, NORMAL,
C     DLE > .0, COMPUTES COEFFICIENT OF LN(E-EL) NEAR E = EL,
C     DLE < .0, COMPUTES COEFFICIENT OF (E-EL)**1/3 FOR E NEAR
C     EL AND EL = EMID,
C     DE = .0, GIVES END CORRECTIONS, S FOR E < OR E > EMAX,
C
      IMPLICIT REAL*8 (A-H,O-Z)
      real*8 darcos
      DIMENSION S(3), X(3)
      REAL*8 KAPPA, KAPRM, K(3), KL, NUM(3), K2, KC, KK, KC2
      COMMON KAPPA, KAPRM, H, Q, RQ, RP, H2, R1
      COMMON/A/NA, NS, KL, ISIG
      COMMON/C/EPS, DLE
      COMMON/D/S12, S11 /PI/PI
      PI = DATAN(1.D0)*4.
      DO 1 N = 1, 3
 1    S(N) = 0.D0
      DDE = .1*DE
      E2 = EC - DELTA/3.D0
      E3 = EC + DELTA/3.D0
      EMIN = .0
      IF (E2.LT.0.D0) EMIN = E2
      EMAX = .0
      IF (E3.GT.0.D0) EMAX = E3
      EMID = E2
      IF (EMID.EQ.EMIN) EMID = .0
      IF (EMID.EQ.EMAX) EMID = E3
      IF (DABS(E).LT.DDE) E = .0
      IF (DABS(E-E2).LT.DDE) E = E2
      IF (DABS(E-E3).LT.DDE) E = E3
      B2 = BET**.66666666D0
      B4 = B2**2
      RP = B4 + 2.D0/B2
      RQ = 18.D0/RP
      IF (DLE.GT.0.) E = E-EPS
      H = RP*E-2.D0*EC/B2
      H2 = 6.D0*E - 4.D0*EC
      Q = B4*(2.D0*H**2/3.D0)-BET**2*E**2-2.*((E-EC)**2+DELTA**2/9.D0)
      Q = 3.D0*Q/(RP*B4)
      R1 = 2.D0*H*(2.*BET**2+1.D0)/(RP*B2)
      KAPPA = E*(E-E2)*(E-E3)
      KAPRM = E*(E-E2) +E*(E-E3) +(E-E2)*(E-E3)
      XP = KAPRM/RP
      Y = ( Q + RQ*XP )
      YP = ( R1 + 18.D0*H2/RP**2)
      P23 = 2.D0*3.1415927D0/3.D0
      IF (DLE.GT.0.D0) GO TO 50
C     IF (DLE.EQ.0.D0) GO TO 7
C     EL = EMID
C     E = EL
C     XA = ((EMID-EMIN)*(EMAX-EMID)*DE/2.D0)**.66666666D0
C     A = DSQRT(3.D0)*XA
C     RA = DSQRT(A)
C     AK = DSQRT((1.D0-DSQRT(3.D0)/2.)/2.)
C     AP = A+XA-XP
C     CALL CEL1(KK, AK, IER)
C     FJ1 = -KK/RA
C     FJ2 = FJ1/(XP*RP)
C     FJ3 = FJ2/(XP*RP)
C     EK = .0
C     GO TO 30
C
C     NORMAL CASE
 7    G = H/3.D0
      IF (DE.EQ.0.) GO TO 40
      IF (DABS(H**3).LT.DABS(KAPPA*1.E-10)) GO TO 4
C     THREE ROOT REGION
      C = 2.D0* KAPPA*(3.D0/H)**3 -1.D0
      IF (C*C.GE.1.) GO TO 10
      PHI = DARCOS(C)/3.D0
      K(1) = DABS( G*(DCOS(PHI) -.5D0 ))
C
      K(2) = DABS( G*(DCOS(PHI-P23) -.5D0 ))
C
      K(3) = DABS( G*(DCOS(PHI+P23) -.5D0 ))
C
      IF (ISIG.GT.1) PRINT 1000, (K(L), L = 1, 3 )
 1000 FORMAT (' K S = ', 3E12.4 )
      DO 3 L = 1,3
 3    X(L) = K(L)**2
      IF(DABS(E).LT.DDE.OR.DABS(E-E2).LT.DDE.OR.DABS(E-E3).LT.DDE)
     *  GO TO 18
      K2 = (X(2)-X(1))/(X(3)-X(1))
      AK = DSQRT(K2)
      ALPHA2 = (X(2)-X(1))/(XP-X(1))
      KC2 = (X(3)-X(2))/(X(3)-X(1))
      KC = DSQRT(KC2)
      FP = (XP-X(2))/(XP-X(1))
      CALL CEL1( KK, AK, IER )
      CALL CEL2( EK, AK, 1.D0, KC2, IER )
      CALL CEL3( PIK, KC, K2, FP )
      W = ( KAPPA - H*XP )**2 - 4.*XP**3
      RTW = DSQRT(W)
      WP = -2.*H2*(6.D0*XP**2-XP*H**2+H*KAPPA)/RP
      RAC = DSQRT(X(3)-X(1))
      X1P = X(1) -XP
      X2P = X(2) -XP
      X3P = X(3) -XP
      FJ1 = -2.*KK/RAC
      FJ2 = 2.*PIK/(RP*RAC*X1P) - PI/(RP*RTW)
      FJ3 = -2.*RAC*EK/(W*RP**2)+.5D0*KK/(RAC*X1P*X2P*RP**2)
     * -(6.D0*XP**2-XP*H**2+H*KAPPA)*PIK/(RAC*W*X1P*RP**2)
      FJ3 = 2.*FJ3+PI*(6.D0*XP**2-XP*H**2+H*KAPPA)/(W*RTW*RP**2)
      S(1) = YP*FJ2-H2*Y*FJ3
      S(2) = 18.D0*FJ1/RP**2-Y*FJ2
      S(3) = 2.*(X3P*KK/RAC-RAC*EK)*RP**2/3.D0
      GO TO 20
C
C     ONE ROOT REGION
 4    K(3) = (DABS(KAPPA/2.))**.33333333D0
C
      GO TO 11
 10   TEMP = ( DABS(C) + DSQRT(C*C-1.D0))**.33333333D0
      TEMP = DSIGN(TEMP,C)
      K(3) = DABS( .5D0*( TEMP -1.D0 +1.D0/TEMP )*G)
C
 11   K(1) = -1.
      IF (ISIG.GT.1) PRINT 2000, K(3)
 2000 FORMAT (' K(3) = ',E12.4 )
 18   XA = K(3)**2
      A2 = (6.D0*XA**2-XA*H**2+H*KAPPA)/2.D0
      A = DSQRT(A2)
      B1 = H**2/8. -XA/2.
      K2 = ( A +B1 -XA )/(2.*A)
      KC2 = (A -B1 +XA )/(2.*A)
      KC = DSQRT(KC2)
      AP = A+XA-XP
      AM = A-XA+XP
C     IF (DABS(AM/A).LT..001) AM = XP+XA*(DSQRT(3.D0+H/K(3))-1.)
      XAP = XA-XP
      ALPHA = AM/AP
      FP = AP**2/(4.*A*XAP)
      RT = DSQRT(K2+AM**2/(4.*A*XAP))
      RA = DSQRT(A)
      IF (K2.GT.0) GO TO 19
      AK = .0D0
      KK = PI/2.D0
      EK = KK
      PIK = KK/DSQRT(FP)
      GO TO 17
 19   AK = DSQRT(K2)
      CALL CEL1 ( KK, AK, IER )
      CALL CEL2 ( EK, AK, 1.D0, KC2, IER )
      CALL CEL3 ( PIK, KC, K2, FP )
 17   FJ1 = -KK/RA
      AL2M = -(4.*A*XAP)/AP**2
      A2OM = -AM**2/(4.*A*XAP)
      AKK = KC2*ALPHA**2 +K2
      IF (DABS(ALPHA).LT..01) GO TO 25
      FJ2 = (AP*PIK/(2.*XAP) -KK)/(RA*AM*RP)
      FJ3 = FJ2/(AM*RP) -(A*PIK/(2.*XAP**2) -4.*A**2*(A2OM*EK
     * -AKK*KK/AL2M +(A2OM**2-K2)*PIK)/(ALPHA*AP**3*AKK))
     * /(RA*AM*RP**2)
      GO TO 30
 25   FJ2 = KK/(2.*RA*RP*XAP)
      IF (K2.GE..0001) TEMP = (KK-EK)/K2
      IF (K2.LT..0001) TEMP = PI/4.
      FJ3 = (.5D0*TEMP-KK)/(2.*RA*(RP*XAP)**2)
 30   S(1) = YP*FJ2 -H2*Y*FJ3
      S(2) = 18.D0*FJ1/RP**2 -Y*FJ2
      S(3) = (AP*KK/RA -2.*RA*EK)*RP**2/3.D0
 5    IF (DABS(E-EMID).GT.DDE.OR.DLE.LT.0.) GO TO 20
C
C     CORRECT FOR DISCONTINUITY AT E = EMID
      NUM(2) = -Q/KAPRM**2
      NUM(1) = R1/KAPRM**2 +H2*NUM(2)/KAPRM
      NUM(3) = RP/3.D0
      G = 1.5707963D0*DABS(KAPRM)/(2.*K(3))
      S11 = G*NUM(1)
      S12 = G*NUM(2)
      DO 8 L = 1, 3
 8    S(L) = S(L) + G * NUM(L)
      GO TO 20
C
C     CALCULATE OUTSIDE S = S11, S12
 40   RT = DSQRT((KAPPA-H*XP)**2-4.*XP**3)
      S12 = -PI*Y/(RP*RT)
      IF (ISIG.GT.3) PRINT 5000, E, S12
 5000 FORMAT (' E = ',E12.4,'  S12 = ',E13.6)
      GO TO 20
C
C     CALCULATE COEFFICIENT OF LN TERM,E NEAR EL
 50   XC = (H/3.D0)**2
      P = KAPRM - RP*XC
      TEMP = -(DABS(P)/DSQRT(3.*XC))
C     TEMP = TEMP*((2.+3.*EPS/DE)*ALOG(DABS(H*P*(DE+EPS)/(81.*XC**2)))
C    * +(2.-3.*EPS/DE)*ALOG(DABS(H*F*(DE-EPS)/(81.*XC**2))) -6.)
      P2 = P**2
      S(3) = TEMP*RP/3.D0
      S(2) = -TEMP*(Q+RQ*XC)/P2
      S(1) = TEMP*R1/P2 +H2*S(2)/P
      E = E + EPS
C
C     COMMON TERMINATION
 20   IF (ISIG.GT.3) PRINT 3000, E, ( S(L), L= 1,3 )
 3000 FORMAT (' E = ',E12.4,' S S = ',3E14.8 )
C     IF (IABS(ISIG).GT.1) PRINT 4000, FJ1, FJ2, FJ3
 4000 FORMAT (' FJ S ',3E14.8 )
      RETURN
      END
C
      REAL*8 FUNCTION DARCOS(X)
      REAL*8 X,PI
      COMMON/PI/PI
      IF (X**2-.0001) 1, 1, 2
 1    DARCOS = PI/2. - X - X**3/6.
      RETURN
 2    DARCOS = DATAN( DSQRT(1.D0-X**2)/X)
      IF (X.LT.0.) DARCOS = DARCOS + PI
      RETURN
      END
C
C  .............................................................................
C      SUBROUTINE CEL1
C
C      PURPOSE
C         CALOULATE COMPLETE ELLIPTIC INTEGRAL OF FIRST KIND
C
C      USAGE
C         CALL CEL1(RES,AK,IER)
C
C      DESCRIPTION OF PARAMETERS
C         REC   - RESULT VALUE
C         AK    - MODULUS (INPUT)
C         IER   - RESULTANT ERROR CODE WHERE
C                 IER=0  NO ERROR
C                 IER=1  AK NOT IN RENGE -1 TO +1
C
C      REMARKS
C         THE RESULT IS SET TO 1,E75 IF ABS(AK) GE 1
C         THE RESULT IS SET TO 1,0E38 IF ABS(AK) GE 1
C         FOR MODULUS AK AND COMFLEMENTARY MODULUS CK,
C         EQUATION AK*AK*CK*CK=1,0 IS USED
C         AK MUST BE IN THE RANGE -1 TO +1
C
C      SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C         NONE
C         SPMAX
C
C      METHOD
C         DEFINITION
C         CEL1(AK)=INTEGRAL(1/SQRT((1+T*T)*(1+(CKKT)**2)), SUMMED
C         OVER T FROM 0 TO INFINITY),
C         EQUIVALENT ARE THE DEFINITIONS
C         CEL1(AK)=INTEGRAL(1/(COS(T)SQRT(1+(CK*TAN(T))**2)),SUMMED
C         OVER T FROM 0 TO PI/2),
C         CEL1(AK)=INTEGRAL(1/SQRT(1-(AK*SIN(T))**2),SUMMED OVER I
C         FROM 0 TO PI/2), WHERE K=SQRT(1,-CK*CK),
C         EVALUATION
C         LANDENS TRANSFORMATION IS USED FOR CALCULATION,
C         REFERENCE
C         R,BULIRSCH, 'NUMERICAL CALOULATION OF ELLIPTIC INTEGRALS
C         AND ELLIPTIC FUNCTIONS',HANDBOOK SERIES SPECIAL FUNCTIONS,
C         NUMERISCHE MATHEMATIK VOL, 7, 1965, PP, 78-90,
C
C   ............................................................................
C
      SUBROUTINE CEL1(RES,AK,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      IER=0
      ARI=2.
      GEO=(0.5D0-AK)+0.5D0
      GEO=GEO+GEO*AK
      RES=0.5D0
      IF(GEO)1,2,4
    1 IER=1
C   2 RES=1.E75
    2 RES = SPMAX(1)
      RETURN
    3 GEO=GEO*AARI
    4 GEO=DSQRT(GEO)
      GEO=GEO+GEO
      AARI=ARI
      ARI=ARI+GEO
      RES=RES+RES
      IF(GEO/AARI-0.9999)3,5,5
    5 RES=RES/ARI*6.283185D0
      RETURN
      END
C
C     ..........................................................................
C
C        SUBROUTINE CEL2
C
C        PURPOSE
C           COMPUTE THE GENERALIZED COMPLETE ELLIPTIC INTEGRAL OF
C           SECOND KIND,
C
C        USAGE
C           CALL CEL2(RES,AK,A,B,IER)
C        DESCRIPTION OF PARAMETERS
C           REC   - RESULT VALUE
C           AK    - MODULUS (INRUT)
C           A     - SONSTANT TERM IN NUMERATOR
C           B     - FACTOR OF QUADRATIC TERM IN NUMERATOR
C           IER   - RESULTANT ERROR CODE WHERE
C                   IER=0  NO ERROR
C                   IER=1  AK NOT IN RENGE -1 TO +1
C
C        REMARKS
C           FOR ABS(AK) GE 1 THE RESULT IS SET TO 1,E75 IF B IS
C           POSITIVE, TO -1,E75 IF B IS NEGATIVE,
C           FOR ABS(AK) GE 1 THE RESULT IS SET TO 1,0E38 IF B IS
C           POSITIVE, TO -1,0E38 IF B IS NEGATIVE,
C           SPECIAL CASES ARE
C           K(K) OBTAINED WITH A = 1, B = 1
C           E(K) OBTAINED WITH A = 1, B = CK*CK WHERE CK IS
C           COMPLEMENTARY MODULUS,
C           B(K) OBTAINED WITH A = 1, B = 0
C           D(K) OBTAINED WITH A = 0, B = 1
C           WHERE K, E, B, D DEFINE SPECIAL CASES OF THE GENERALIZED
C           COMPLETE ELLIPTIC INTEGRAL OF SEOOND KIND IN THE USUAL
C           NOTATION, AND THE ARGUMENT K OF THESE FUNCTIONS MEANS
C           THE MODULS,
C
C        SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C           SPMAX
C
C        METHOD
C           DEFINITION
C           RES=INTEGRAL((A+B*T*T)/(SQRT((1+T*T)*(1+(CK*T)**2))*(1+T*T))
C        SUMMED OVER T FROM 0 TO INFINITY),
C        EVALUATION
C        LANDENS TRANSFORMATION IS USED FOR CALCULATION,
C        REFERENCE
C        R,BULIRCH, 'NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS
C        AND ELLIPTIC FUNCTIONS' , HANDBOOK SERIES SPECIAL FUNCTIONS,
C        NUMERISOHE MATHEMATIK VOL, 7, 1965, PP, 78-90,
C
C     ..........................................................................
C
      SUBROUTINE CEL2(RES,AK,A,B,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      IER=0
      ARI=2.
      GEO=(0.5D0-AK)+0.5D0
      GEO=GEO*AK + GEO
      RES=A
      A1=A+B
      B0=B+B
      IF(GEO)1,2,6
    1 IER=1
    2 IF(B)3,8,4
C   3 RES=-1.E75
    3 RES=-SPMAX(1)
      RETURN
C   4 RES=1.E75
    4 RES=SPMAX(1)
      RETURN
    5 GEO=GEO*AARI
    6 GEO=DSQRT(GEO)
      GEO=GEO+GEO
      AARI=ARI
      ARI=ARI+GEO
      B0=B0+RES*GEO
      RES=A1
      B0=B0+B0
      A1=B0/ARI+A1
      IF(GEO/AARI-0.9999)5,7,7
    7 RES=A1/ARI
      RES=RES+0.5707963D0*RES
    8 RETURN
      END
      REAL*8 FUNCTION SPMAX(K)
      SPMAX=1.0D38
      RETURN
      END
