1      subroutine compel(dk, dellk)
2c!purpose
3c  calculate complete elliptic integral of first kind
4c!calling sequence
5c     subroutine compel(dk, dellk)
6c     double precision dk,dellk
7c!
8      double precision dpi, domi,dellk
9      double precision de, dgeo, dk, dri, dari, dtest
10      external slamch, dlamch
11      double precision dlamch,flma
12      real slamch
13c
14      data de /1.0d+0/
15      flma=2.0d+0**(int(slamch('l'))-2)
16      dpi=4.0d+0*atan(1.0d+0)
17      domi=2.0d+0*dlamch('p')
18      dgeo = de - dk*dk
19      if (dgeo .le. 0) then
20         goto 10
21      else
22         goto 20
23      endif
24  10  dellk = flma
25      return
26  20  dgeo = sqrt(dgeo)
27      dri = de
28  30  dari = dri
29      dtest = dari*domi
30      dri = dgeo + dri
31      if((dari-dgeo-dtest).le.0.0d+0) goto 50
32      dgeo = sqrt(dari*dgeo)
33      dri = 0.50d+0*dri
34      go to 30
35  50  dellk = dpi/dri
36      return
37      end
38