C
C     ..................................................................
C
C        SUBROUTINE JELF
C
C        PURPOSE
C           COMPUTES THE THREE JACOBIAN ELLIPTIC FUNCTIONS SN, CN, DN.
C
C        USAGE
C           CALL JELF(SN,CN,DN,X,SCK)
C
C        DESCRIPTION OF PARAMETERS
C           SN    - RESULT VALUE SN(X)
C           CN    - RESULT VALUE CN(X)
C           DN    - RESULT VALUE DN(X)
C           X     - ARGUMENT OF JACOBIAN ELLIPTIC FUNCTIONS
C           SCK   - SQUARE OF COMPLEMENTARY MODULUS
C
C        REMARKS
C           NONE
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           DEFINITION
C           X=INTEGRAL(1/SQRT((1-T*T)*(1-(K*T)**2)), SUMMED OVER
C           T FROM 0 TO SN), WHERE K=SQRT(1-SCK).
C           SN*SN + CN*CN = 1
C           (K*SN)**2 + DN**2 = 1.
C           EVALUATION
C           CALCULATION IS DONE USING THE PROCESS OF THE ARITHMETIC
C           GEOMETRIC MEAN TOGETHER WITH GAUSS DESCENDING TRANSFORMATION
C           BEFORE INVERSION OF THE INTEGRAL TAKES PLACE.
C           REFERENCE
C           R. BULIRSCH, NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS AND
C                  ELLIPTIC FUNCTIOMS.
C                  HANDBOOK SERIES OF SPECIAL FUNCTIONS
C                  NUMERISCHE MATHEMATIK VOL. 7, 1965, PP. 78-90.
C
C     ..................................................................
C
      SUBROUTINE JELF(SN,CN,DN,X,SCK)
C
C
      DIMENSION ARI(12),GEO(12)
C     TEST MODULUS
      CM=SCK
      Y=X
      IF(SCK)3,1,4
    1 D=EXP(X)
      A=1./D
      B=A+D
      CN=2./B
      DN=CN
      SN=TANH(X)
C        DEGENERATE CASE SCK=0 GIVES RESULTS
C           CN X = DN X = 1/COSH X
C           SN X = TANH X
    2 RETURN
C        JACOBIS MODULUS TRANSFORMATION
    3 D=1.-SCK
      CM=-SCK/D
      D=SQRT(D)
      Y=D*X
    4 A=1.
      DN=1.
      DO 6 I=1,12
      L=I
      ARI(I)=A
      CM=SQRT(CM)
      GEO(I)=CM
      C=(A+CM)*.5
      IF(ABS(A-CM)-1.E-4*A)7,7,5
    5 CM=A*CM
    6 A=C
C
C     START BACKWARD RECURSION
    7 Y=C*Y
      SN=SIN(Y)
      CN=COS(Y)
      IF(SN)8,13,8
    8 A=CN/SN
      C=A*C
      DO 9 I=1,L
      K=L-I+1
      B=ARI(K)
      A=C*A
      C=DN*C
      DN=(GEO(K)+A)/(B+A)
    9 A=C/B
      A=1./SQRT(C*C+1.)
      IF(SN)10,11,11
   10 SN=-A
      GOTO 12
   11 SN=A
   12 CN=C*SN
   13 IF(SCK)14,2,2
   14 A=DN
      DN=CN
      CN=A
      SN=SN/D
      RETURN
      END