C**************************** D P L *********************************
      FUNCTION DPL(n,x,y,m,t)
C
C        Version 3.0  - 15/1/97       Author: A. SALIN
C
C        Calculation of a function by interpolation. Double precision.
C        See comments in SUBROUTINE DSPLIN
C
C      IMPLICIT REAL*8 (a-h,o-z)
      IMPLICIT REAL (a-h,o-z)
c
      DIMENSION x(n),y(n)
      REAL m(n)
C      REAL*8 m(n)
      LOGICAL order
c
      COMMON/KODSPL/kod,k,iw,iligne
c
    1 FORMAT(' DPL: extrapolation - ',1PE15.8,' <',E15.8)
    2 FORMAT(' DPL: extrapolation - ',1PE15.8,' >',E15.8)
c
C        do i=1,n
C        write(*,*)  x(i),y(i),m(i)
c        enddo
      order=x(2).GT.x(1)
      IF(k.LE.1.OR.k.GT.n) THEN
	klo=0
	khi=n+1
	GO TO 100
      END IF
C
      inc=1
      klo=k-1
      IF(t.GT.x(klo).EQV.order) THEN
 10     khi=klo+inc
	IF(khi.GT.n) THEN
	  khi=n+1
	ELSE IF(t.GT.x(khi).EQV.order) THEN
	  klo=khi
	  inc=inc+inc
	  GO TO 10
	END IF
      ELSE
	khi=klo
 20     klo=khi-inc
	IF(klo.LT.1) THEN
	  klo=0
	ELSE IF(t.LT.x(klo).EQV.order) THEN
	  khi=klo
	  inc=inc+inc
	  GO TO 20
	END IF
      END IF
C
 100  CONTINUE
      IF(khi-klo.GT.1) THEN
	km=(khi+klo)/2
	IF(t.GT.x(km).EQV.order) THEN
	  klo=km
	ELSE
	  khi=km
	END IF
	GO TO 100
      END IF
	k=klo+1
C
      IF(t.EQ.x(1)) k=2
      IF(k.LE.1) THEN
	e=x(2)-x(1)
	DPL=((y(2)-y(1))/e-m(2)*e/6.D0)*(t-x(1))+y(1)
	kod=1
	IF(iw.GT.0.AND.iligne.GT.0) THEN
	  IF(order) THEN
	    WRITE(iw,1) t,x(1)
	  ELSE
	    WRITE(iw,2) t,x(1)
	  END IF
	  iligne=iligne-1
	END IF
      ELSE IF(k.GT.n) THEN
	e=x(n)-x(n-1)
	DPL=((y(n)-y(n-1))/e+m(n-1)*e/6.D0)*(t-x(n))+y(n)
	kod=2
	IF(iw.GT.0.AND.iligne.GT.0) THEN
	  IF(order) THEN
	    WRITE(iw,2) t,x(n)
	  ELSE
	    WRITE(iw,1) t,x(n)
	  END IF
	  iligne=iligne-1
	END IF
      ELSE
	f=x(k)-t
	g=t-x(k-1)
	e=f+g
	DPL=(-g*f*(m(k-1)*(f+e)+m(k)*(g+e))+6.D0*(g*y(k)+f*y(k-1)))/
     1      (6.D0*e)
	kod=0
      END IF
      END FUNCTION
C**************************** D P L P *******************************
      FUNCTION DPLP(n,x,y,m,t)
C
C        Version 3.0  - 17/1/97       Author: A. SALIN
C
C        Calculation of the derivative of a function by interpolation.
C        Double precision. See comments in SUBROUTINE DSPLIN.
C        If t is outside the interval over which the function is
C        defined, the derivative is extrapolated linearly from its
C        value for the first two (last two) points of the interval.
C
C      IMPLICIT REAL*8 (a-h,o-z)
      IMPLICIT REAL (a-h,o-z)
c
      DIMENSION x(n),y(n)
C      REAL*8 m(n)
      REAL m(n)
      LOGICAL order
c
      COMMON/KODSPL/kod,k,iw,iligne
c
    1 FORMAT(' DPLP: extrapolation - ',1PE15.8,' <',E15.8)
    2 FORMAT(' DPLP: extrapolation - ',1PE15.8,' >',E15.8)
c
      order=x(2).GT.x(1)
      IF(k.LE.1.OR.k.GT.n) THEN
        klo=0
        khi=n+1
        GO TO 100
      END IF
C
      inc=1
      klo=k-1
      IF(t.GT.x(klo).EQV.order) THEN
 10     khi=klo+inc
        IF(khi.GT.n) THEN
          khi=n+1
        ELSE IF(t.GT.x(khi).EQV.order) THEN
          klo=khi
          inc=inc+inc
          GO TO 10
        END IF
      ELSE
        khi=klo
 20     klo=khi-inc
        IF(klo.LT.1) THEN
          klo=0
        ELSE IF(t.LT.x(klo).EQV.order) THEN
          khi=klo
          inc=inc+inc
          GO TO 20
        END IF
      END IF
C
 100  CONTINUE
      IF(khi-klo.GT.1) THEN
        km=(khi+klo)/2
        IF(t.GT.x(km).EQV.order) THEN
          klo=km
        ELSE
          khi=km
        END IF
        GO TO 100
      END IF
        k=klo+1
C
      IF(t.EQ.x(1)) k=2
      IF(k.LE.1) THEN
        e=x(2)-x(1)
        g=t-x(1)
        f=x(2)-t
        DPLP=(y(2)-y(1))/e+(m(2)*(2*g-f)+m(1)*(g-2*f))/6.D0
        kod=1
        IF(iw.GT.0.AND.iligne.GT.0) THEN
          IF(order) THEN
            WRITE(iw,1) t,x(1)
          ELSE
            WRITE(iw,2) t,x(1)
          END IF
          iligne=iligne-1
        END IF
      ELSE IF(k.GT.n) THEN
        e=x(n)-x(n-1)
        g=t-x(n-1)
        f=x(n)-t
        DPLP=(y(n)-y(n-1))/e+(m(n)*(2*g-f)+m(n-1)*(g-2*f))/6.D0
        kod=2
        IF(iw.GT.0.AND.iligne.GT.0) THEN
          IF(order) THEN
            WRITE(iw,2) t,x(n)
          ELSE
            WRITE(iw,1) t,x(n)
          END IF
          iligne=iligne-1
        END IF
      ELSE
        f=x(k)-t
        g=t-x(k-1)
        e=f+g
        e2=e*e
        DPLP=(-m(k-1)*(3.D0*f*f-e2)+m(k)*(3.D0*g*g-e2)
     1       +6.D0*(y(k)-y(k-1)))/(6.D0*e)
        kod=0
      END IF
      END FUNCTION
C**************************** D S P L I N *************************** 
C        Interpolation package. Given a set of n couples x(i),y(i), 
C     DSPLIN defines a cubic spline function F(x) such that: 
C        1)- F(x(i)) = y(i) 
C        2)- F(x), F'(x), F''(x) are continuous in the interval 
C     [x(1),x(n)]. 
C        Required subroutine: TRIDIA. 
C 
      BLOCK DATA BLKSPL 
      COMMON/KODSPL/kod,k,iw,iligne 
      DATA iw,iligne,k/6,100,2/ 
      END 
C******************************************************************** 
C 
      SUBROUTINE DSPLIN(n,x,y,cm,cm1,ic1,cmn,icn,alpha,beta,b) 
C 
C        Author: A. Salin   Version 3.2  20/1/86 - 10/6/97      
C 
C     ****SUBROUTINE DSPLIN: calculates the vector F''(x(i)) which 
C     defines the cubic spline function F(x). 
C          N= number of pivots x(i). 
C          X= array of abcissae. Should be stored by increasing or     
C     decreasing order. 
C          Y= array of values of Y(i). 
C          CM= array of F''(x(i)) with dimension equal to that of X 
C     and Y.       
C          IC1,CM1,ICN,CMN: define the condition at x(1) and x(n) 
C     If IC1=0, function is the same in first two intervals. 
C     If IC1=1, second derivative at x(1) given in CM1 
C     If IC1=2, first derivative at x(1) given in CM1 
C     IF IC1>2, zero second derivative at x(1) 
C          Similar definitions for ICN and CMN around x(n) 
C          ALPHA, BETA, B are working arrays of dimension N at least. 
C     For cyclic spline use the program TRISPL. 
C 
C     ****DPLCOF: the interpolated function at R is: 
C       F=C(1)*CM(K-1)+C(2)*CM(K)+C(3)*Y(K-1)+C(4)*Y(K) 
C       where C is obtained from: 
C       CALL DPLCOF(n,x,r,c,k) 
C       The vector C should be of dimension 4. 
C 
C     ****Functions DPL,DPLP,DPLP2: calculate respectively the function 
C     F(x), its first or second derivative:  
C             N= number of pivots 
C             X,Y,M= same as X,Y,CM in DSPLIN. 
C             T= value for which the function (or its derivative) 
C     must be determined. 
C     WARNING: when T is outside the interval [X(1),X(N)], F(X)  
C     *******  (resp. F', F") is determined by linear extrapolation 
C              using the value of F (resp. F', F") for the first two   
C              (last two) points in the interval. 
C 
C     ****Subroutine DINITI et DINTSP: calculate the integral of the 
C     function F(x) from X(1) to T. The subroutine DINITI should first 
C     be called before the first CALL DINTSP concerning a given 
C     function F(x). DINITI calculates the array (CI) of values of the 
C     integral of F(x) from X(1) to all X(i). 
C 
C        Arguments of DINITI: 
C             N,X,Y,CM: as in subroutine DSPLIN. 
C             CI: array of dimension at least N.   
C        Arguments of DINTSP: 
C             N,X,Y,CM,CI: as in DINITI 
C             T: value of the upper limit of the integration (X(1).LT.T. 
C     LE.X(N)). 
C             B: value of the integral. 
C       The calculations in DINTSP are nearly as rapid as for the deter 
C     mination of the function F(x) by DPL. 
C 
C     Parameters of COMMON/KODSPL/: 
C     ----------------------------------------------------- 
C       -KOD= after execution of one subroutine of the package, KOD 
C     takes the value: 
C              *0: no error 
C              *-1: values of X(i) are not stored in correct order. 
C              *1: T outside the interval [X(1),X(N)]. Linear extrapo- 
C     lation done using X(1) and X(2). 
C              *2: T outside the interval [X(1),X(N)]. Linear extrapo- 
C     lation done using X(n-1) and X(n). 
C       -K= after the execution of DPL, DPLP or DPLP2, K is such that 
C     T lies in the interval [X(k-1),X(k)]. 
C       -IW= error messages are printed on unit IW unless IW.LE.0 
C     Default: 6. 
C       -ILIGNE: when T is outside the interval [X(1),X(n)], DPL, DPLP 
C     and DPLP2 print a message if ILIGNE.GT.0. The initial value of 
C     ILIGNE (defined by DATA) is 100. For every extrapolation, the 
C     value is decreased by 1. 
C 
C     *************************************************************** 
C
C     IMPLICIT REAL*8 (a-h,o-z) 
      IMPLICIT REAL (a-h,o-z) 
      PARAMETER (six=6.D0,douze=12.D0,zero=0.D0) 
C 
      DIMENSION x(n),y(n),cm(n) 
      DIMENSION alpha(n),beta(n),b(n) 
C      
      COMMON/KODSPL/kod,k,iw,iligne 
C 
   91 FORMAT('Error in the data for DSPLIN',/,'Check abcissae:', 
     1       1PD15.8,1X,' and ',D15.8,/,'Program stopped') 
C 
C                                           Contour condition at X(1) 
C        do i=1,n
C        write(*,*)  x(i),y(i)
C        enddo
      IF(ic1.EQ.0) THEN 
        s1=zero 
        fac1=(x(2)-x(1))*(1.D0+(x(2)-x(1))/(x(3)-x(2)))/six 
        fab1=-(x(2)-x(1))**2/(six*(x(3)-x(2))) 
      ELSE IF(ic1.EQ.1) THEN 
        s1=-cm1*(x(2)-x(1))/six 
        fac1=zero 
        fab1=zero 
      ELSE IF(ic1.EQ.2) THEN 
        s1=0.5D0*(cm1-(y(2)-y(1))/(x(2)-x(1))) 
        fac1=-(x(2)-x(1))/douze 
        fab1=zero 
      ELSE 
C        WRITE(iw,*) 'DSPLIN: unknown condition - use default' 
        s1=zero 
        fac1=zero 
        fab1=zero 
      END IF 
C                                           Contour condition at X(n) 
      IF(icn.EQ.0) THEN 
        sn=zero 
        facn=(x(n)-x(n-1))*(1.D0+(x(n)-x(n-1))/(x(n-1)-x(n-2)))/six 
        fabn=-(x(n)-x(n-1))**2/(six*(x(n-1)-x(n-2))) 
      ELSE IF(icn.EQ.1) THEN 
        sn=-cmn*(x(n)-x(n-1))/six 
        facn=zero 
        fabn=zero 
      ELSE IF(icn.EQ.2) THEN 
        sn=-0.5D0*(cmn-(y(n)-y(n-1))/(x(n)-x(n-1))) 
        facn=-(x(n)-x(n-1))/douze 
        fabn=zero 
      ELSE 
C        WRITE(iw,*) 'DSPLIN: unknown condition - use default' 
        sn=zero 
        facn=zero 
        fabn=zero 
      END IF 
C 
      kod=0 
      c=x(2)-x(1) 
      e=(y(2)-y(1))/c 
      DO i=3,n 
        i2=i-2 
        a=c 
        c=x(i)-x(i-1) 
        IF(a*c.LE.zero) THEN    
          kod=-1 
          IF(iw.GT.0) WRITE(iw,91) x(i-1),x(i) 
          STOP 
          END IF 
        alpha(i2)=(a+c)/3.D0 
        beta(i2)=c/six 
        cm(i2)=beta(i2) 
        d=e 
        e=(y(i)-y(i-1))/c 
        b(i2)=e-d 
      END DO 
C 
      b(1)=b(1)+s1 
      b(n-2)=b(n-2)+sn 
      alpha(1)=alpha(1)+fac1 
      alpha(n-2)=alpha(n-2)+facn 
      cm(n-3)=cm(n-3)+fabn 
      beta(1)=beta(1)+fab1 
C                                            Solve tridiagonal system 
      CALL TRIDIA(alpha,cm,beta,b,cm(2),n-2) 
C      
      IF(ic1.EQ.0) THEN 
        cm(1)=cm(2)*(1.D0+(x(2)-x(1))/(x(3)-x(2))) 
     1        -cm(3)*(x(2)-x(1))/(x(3)-x(2)) 
      ELSE IF(ic1.EQ.1) THEN 
        cm(1)=cm1 
      ELSE IF(ic1.EQ.2) THEN 
        cm(1)=-six*s1/(x(2)-x(1))-cm(2)/2.D0 
      ELSE 
        cm(1)=zero 
      END IF 
C       
      IF(icn.EQ.0) THEN 
        cm(n)=cm(n-1)*(1.D0+(x(n)-x(n-1))/(x(n-1)-x(n-2))) 
     1        -cm(n-2)*(x(n)-x(n-1))/(x(n-1)-x(n-2)) 
      ELSE IF(icn.EQ.1) THEN 
        cm(n)=cmn 
      ELSE IF(icn.EQ.2) THEN 
        cm(n)=-six*sn/(x(n)-x(n-1))-cm(n-1)/2.D0 
      ELSE 
        cm(n)=zero 
      END IF 
C 
      k=2 
      RETURN 
      END SUBROUTINE
c Incluye todas las subrutinas que utilizan las programas de dinamica de
c superficies
c
C**************************** D I F S Y *******************************
C                                                             
C       Traduction du programme ALGOL de Bulirsch and Stoer, Numerische
C     Mathematik 8(1966)1.
C       DIFSY calcule un pas d'integration H inferieur ou egal a H0 pour
C     une equation differentielle ordinaire du premier ordre de forme:
C     dZ/dX=F(X,Z)
C     dont le membre de droite doit etre donne par un sous-programme
C     commencant par:
C     SUBROUTINE F(X,Z,DZ)
C       Le programme prend pour pas H le premier des nombres H0, H0/2,
C     H0/4 pour lequel 9 etapes d'extrapolation suffisent a obtenir le
C     resultat avec la precision requise. Si H est different de H0, le
C     parametre FIN prend la valeur .TRUE.
C       X et Y sont les valeurs initiales. Apres execution de DIFSY, les
C     valeurs initiales de X et Y sont remplacees par X+H et Y(X+H). H0
C     peut egalement etre modifie (changement de pas automatique). La
C     nouvelle valeur de H0 est la meilleure valeur estimee pour le pas
C     suivant.
C       S et la constante EPS commandent la precision des valeurs
C     calculees. Le calcul est arrete lorsque deux valeurs successives
C     de Y different au plus de EPS*S. EPS ne doit pas etre inferieur a
C     10**(-D+3) ou D est le nombre de chiffres decimaux dans la machine
C     utilisee.
C         Il est conseille de faire S=0. pour le premier pas
C     d'integration. A la sortie de DIFSY, S est devenu S=MAX(S,Y(XSI))
C     pour XSI compris dans l'intervalle (X,X+H).
C     *
      SUBROUTINE DIFSY(x,y,h0,eps,s,f,fin)
      INTEGER r,sr
      DIMENSION dt(7),d(7),yg(8),yh(8)
C      REAL*8 a,d,fc,g,h0,s,x,eps
C      REAL*8 y,dz,ya,yl,ym,u,yh,yg,dt,v,ta,c,b1,b,dy
      REAL a,d,fc,g,h0,s,x,eps
      REAL y,dz,ya,yl,ym,u,yh,yg,dt,v,ta,c,b1,b,dy
      LOGICAL konv,bo,bh,fin
C     *
      CALL F(x,y,dz)
      bh=.FALSE.
      fin=.FALSE.
      ya=y
   20 CONTINUE
      a=h0+x
      fc=1.5D0
      bo=.FALSE.
      m=1
      r=2
      sr=3
      jjp1=0
C
C                                         Extrapolations successives.
      DO 170 jp1=1,10
        IF(bo) THEN
          d(2)=16.D0/9.D0
          d(4)=64.D0/9.D0
          d(6)=256.D0/9.D0
        ELSE
          d(2)=9.D0/4.D0
          d(4)=9.D0
          d(6)=36.D0
        END IF
        konv=.FALSE.
        IF(jp1.GT.3) konv=.TRUE.
        IF(jp1.GT.7) THEN
          lp1=7
          fc=0.6D0*fc
          d(7)=64.D0
        ELSE
          lp1=jp1
           d(LP1)=FLOAT(m*m)
C          d(LP1)=DFLOAT(m*m)
        END IF
        m=2*m
C        g=h0/DFLOAT(m)
        g=h0/FLOAT(m)
        b=2.D0*g
        IF(bh.AND.jp1.LT.9) THEN
          ym=yh(jp1)
          yl=yg(jp1)
        ELSE
          kk=(m-2)/2
C
C                                                       Calcul de T(H)
          m=m-1
          yl=ya
          ym=ya+g*dz
          DO 100 k=1,m
            CALL F(x+FLOAT(k)*g,ym,dy)
C            CALL F(x+DFLOAT(k)*g,ym,dy)
            u=yl+b*dy
            yl=ym
            ym=u
C            u=DABS(u)
            u=ABS(u)
            IF(u.GT.s)  s=u
            IF(k.NE.kk.OR.k.EQ.2) GO TO 100
            jjp1=1+jjp1
            yh(jjp1)=ym
            yg(jjp1)=yl
  100     CONTINUE
        END IF
        CALL F(a,ym,dy)
        v=dt(1)
        ta=(ym+yl+g*dy)/2.D0
C
C                                 Calcul de la valeur extrapolee T(0).
        c=ta
        dt(1)=ta
        IF(lp1.GE.2) THEN
          DO 150 kp1=2,lp1
            b1=d(kp1)*v
            b=b1-c
            u=v
            IF(b.NE.0.D0) THEN
              b=(c-v)/b
              u=c*b
              c=b1*b
            END IF
            v=dt(kp1)
            dt(kp1)=u
            ta=u+ta
  150     CONTINUE
        END IF
C        IF(DABS(y-ta).GT.eps*s) konv=.FALSE.
        IF(ABS(y-ta).GT.eps*s) konv=.FALSE.
        y=ta
        IF(konv) GO TO 180
        d(3)=4.D0
        d(5)=16.D0
        bo=.NOT.bo
        m=r
        r=sr
        sr=m*2
  170 CONTINUE
C
C          Plus de 9 iterations etaient necessaires. On recommence en
C          divisant le pas par 2.
        bh=.NOT.bh
        fin=.TRUE.
        h0=h0/2.D0
        GO TO 20
  180 h0=fc*h0
      x=a
      RETURN
      END SUBROUTINE
c
C**************************** D I F S Y 2 ***************************
      SUBROUTINE DIFSY2(n,x,y,h0,eps,s,f,fin)
C 
C       Version 0.1     29/1/97
C
C       Traduction du programme ALGOL de Bulirsch and Stoer, 
C     Numerische Mathematik 8(1966)1. Double precision.
C       DIFSY2 calcule un pas d'integration H inferieur ou egal a H0 
C     pour un systeme de N equations differentielles ordinaires du 
C     premier ordre de forme:
C     dY/dX=F(X,Y)
C     dont le membre de droite doit etre donne par un sous-programme
C     commencant par:
C     SUBROUTINE F(x,y,dy)
C     REAL*8 y(n),dy(n)
C     et dont le nom est passe comme argument (f) dans la sequence
C     d'appel de DIFSY2.
C       Le programme prend pour pas H le premier des nombres H0, H0/2,
C     H0/4 pour lequel 9 etapes d''extrapolation suffisent a obtenir le
C     resultat avec la precision requise. Si H est different de H0, le
C     parametre FIN prend la valeur .TRUE.
C       X et Y sont les valeurs initiales. Apres execution de DIFSY2, 
C     les valeurs initiales de X et Y sont remplacees par X+H et 
C     Y(X+H). H0 peut egalement etre modifie (changement de pas 
C     automatique). La nouvelle valeur de H0 est la meilleure valeur 
C     estimee pour le pas suivant.
C       Le tableau S et la constante EPS commandent la precision des 
C     valeurs calculees. Le calcul est arrete lorsque deux valeurs 
C     successives de Y(i) different au plus de EPS*S(i). EPS ne doit 
C     pas etre inferieur a 10**(-D+3) ou D est le nombre de chiffres 
C     decimaux dans la machine utilisee.
C       Il est conseille de faire S(i)=0. pour le premier pas
C     d''integration. A la sortie de DIFSY2, S(i) est devenu 
C     S(i)=MAX[S(i),Y(i,xsi)] 
C     pour xsi compris dans l''intervalle [X,X+H].
C     *
      PARAMETER (id=2)
      REAL d2,d4,d6,dp2
C      REAL*8 d2,d4,d6,dp2
      PARAMETER(d2=16.D0/9.D0,d4=64.D0/9.D0,d6=256.D0/9.D0,
     1          dp2=9.D0/4.D0)
C
      INTEGER r,sr
      LOGICAL bo,bh,fin,konv
      REAL d(7),dt(id,7),dy(id),dz(id),s(id),y(id),ya(id),yg(8,id),
     1       yh(8,id),yl(id),ym(id)
C      REAL*8 d(7),dt(id,7),dy(id),dz(id),s(id),y(id),ya(id),yg(8,id),
C     1       yh(8,id),yl(id),ym(id)
C      REAL*8 a,b1,b,c,eps,fc,g,h0,ta,u,v,x
      REAL a,b1,b,c,eps,fc,g,h0,ta,u,v,x
C
C      write(*,*) x,y(1),y(2)
      IF(n.GT.id) THEN
        WRITE(*,*) 'ID too small in DIFSY2'
	STOP
      END IF
      	
      CALL F(x,y,dz)
      bh=.FALSE.
      fin=.FALSE.
      DO i=1,n
        ya(i)=y(i)
      END DO
C                                                          Iterations
   20 CONTINUE
      a=h0+x
      fc=1.5D0
      bo=.FALSE.
      m=1
      r=2
      sr=3
      jjp1=0
C                                          Extrapolations successives
      DO 170 jp1=1,10
        IF(bo) THEN    
          d(2)=d2
          d(4)=d4
          d(6)=d6
        ELSE
          d(2)=dp2
          d(4)=9.D0
          d(6)=36.D0
        END IF
        konv=.FALSE.
        IF(jp1.GT.3) konv=.TRUE.
        IF(jp1.GT.7) THEN    
          lp1=7
          fc=0.6D0*fc
          d(7)=64.D0
        ELSE
          lp1=jp1
C          d(lp1)=DFLOAT(m**2)
          d(lp1)=FLOAT(m**2)
        END IF
        m=2*m
C        g=h0/DFLOAT(m)
        g=h0/FLOAT(m)
        b=2.D0*g
        IF(bh.AND.jp1.LT.9) THEN
          DO i=1,n
            ym(i)=yh(jp1,i)
            yl(i)=yg(jp1,i)
          END DO
        ELSE
C                                                      Calcul de T(H)
          kk=(m-2)/2
          m=m-1
          DO  i=1,n
            yl(i)=ya(i)
            ym(i)=ya(i)+g*dz(i)
          END DO
          DO k=1,m
C            CALL F(x+DFLOAT(k)*g,ym,dy)
            CALL F(x+FLOAT(k)*g,ym,dy)
            DO i=1,N
              u=yl(i)+b*dy(i)
              yl(i)=ym(i)
              ym(i)=u
C              u=DABS(u)
              u=ABS(u)
              IF(u.GT.s(i))  s(i)=u
            END DO  
            IF(k.NE.kk.OR.k.EQ.2) GO TO 100
            jjp1=1+jjp1
            DO i=1,n
              yh(jjp1,i)=ym(i)
              yg(jjp1,i)=yl(i)
            END DO  
 100        CONTINUE
          END DO  
        END IF   
        CALL F(a,ym,dy)
        DO 160 i=1,n
          v=dt(i,1)
          ta=(ym(i)+yl(i)+g*dy(i))/2.D0
C                                 Calcul de la valeur extrapolee T(0)
          c=ta
          dt(i,1)=ta
          IF(lp1.GE.2) THEN     
            DO kp1=2,lp1
              b1=d(kp1)*v
              b=b1-c
              u=v
              IF(b.NE.0.D0) THEN     
                b=(c-v)/b
                u=c*b
                c=b1*b
              END IF
              v=dt(i,kp1)
              dt(i,kp1)=u
              ta=u+ta
            END DO  
          END IF
C          IF(DABS(y(i)-ta).GT.eps*s(i)) konv=.FALSE.
          IF(ABS(y(i)-ta).GT.eps*s(i)) konv=.FALSE.
          y(i)=ta
 160    CONTINUE          
        IF(konv) GO TO 180
        d(3)=4.D0
        d(5)=16.D0
        bo=.NOT.bo
        m=r
        r=sr
        sr=m*2
 170  CONTINUE
C
C     Plus de 9 iterations etaient necessaires. On recommence en
C     divisant le pas par 2.
C     
      bh=.NOT.bh
      fin=.TRUE.
      h0=h0/2.D0
      GO TO 20
C 
 180  h0=fc*h0
      x=a
      RETURN
      END SUBROUTINE
c
C**************************** D I F S Y 2 S *************************
      SUBROUTINE DIFSY2S(n,x,y,h0,eps,s,f,switch)
C 
C       Version 0.1     29/1/97
C
C       Traduction du programme ALGOL de Bulirsch and Stoer, 
C     Numerische Mathematik 8(1966)1. Double precision.
C       Cette version inclue un switch logique (SWITCH) determiné par 
C     la routine F qui arrete l'integration quand SWITCH prend la 
C     valeur .TRUE. 
C       DIFSY2S calcule un pas d'integration H inferieur ou egal a H0 
C     pour un systeme de N equations differentielles ordinaires du 
C     premier ordre de forme:
C     dY/dX=F(X,Y)
C     dont le membre de droite doit etre donne par un sous-programme
C     commencant par:
C     SUBROUTINE F(x,y,dy,switch)
C     REAL*8 y(n),dy(n)
C     et dont le nom est passe comme argument (f) dans la sequence
C     d'appel de DIFSY2.
C       Le programme prend pour pas H le premier des nombres H0, H0/2,
C     H0/4 pour lequel 9 etapes d'extrapolation suffisent a obtenir le
C     resultat avec la precision requise.
C       X et Y sont les valeurs initiales. Apres execution de DIFSY2, 
C     les valeurs initiales de X et Y sont remplacees par X+H et 
C     Y(X+H). H0 peut egalement etre modifie (changement de pas 
C     automatique). La nouvelle valeur de H0 est la meilleure valeur 
C     estimee pour le pas suivant.
C       Le tableau S et la constante EPS commandent la precision des 
C     valeurs calculees. Le calcul est arrete lorsque deux valeurs 
C     successives de Y(i) different au plus de EPS*S(i). EPS ne doit 
C     pas etre inferieur a 10**(-D+3) ou D est le nombre de chiffres 
C     decimaux dans la machine utilisee.
C       Il est conseille de faire S(i)=0. pour le premier pas
C     d'integration. A la sortie de DIFSY2, S(i) est devenu 
C     S(i)=MAX[S(i),Y(i,xsi)] 
C     pour xsi compris dans l'intervalle [X,X+H].
C     *
      PARAMETER (id=99)
C      REAL*8 d2,d4,d6,dp2
      REAL d2,d4,d6,dp2
      PARAMETER(d2=16.D0/9.D0,d4=64.D0/9.D0,d6=256.D0/9.D0,
     1          dp2=9.D0/4.D0)
C
      INTEGER r,sr
      LOGICAL bo,bh,konv,switch
      REAL d(7),dt(id,7),dy(id),dz(id),s(id),y(id),ya(id),yg(8,id),
     1       yh(8,id),yl(id),ym(id)
C      REAL*8 d(7),dt(id,7),dy(id),dz(id),s(id),y(id),ya(id),yg(8,id),
C     1       yh(8,id),yl(id),ym(id)
C      REAL*8 a,b1,b,c,eps,fc,g,h0,ta,u,v,x
      REAL a,b1,b,c,eps,fc,g,h0,ta,u,v,x
C
      IF(n.GT.id) THEN
        WRITE(*,*) 'ID too small in DIFSYS'
	STOP
      END IF
      	
      switch=.FALSE.
      CALL F(x,y,dz,switch)
      IF(switch) RETURN
      bh=.FALSE.
      DO i=1,n
        ya(i)=y(i)
      END DO
C                                                          Iterations
   20 CONTINUE
      a=h0+x
      fc=1.5D0
      bo=.FALSE.
      m=1
      r=2
      sr=3
      jjp1=0
C                                          Extrapolations successives
      DO 170 jp1=1,10
        IF(bo) THEN    
          d(2)=d2
          d(4)=d4
          d(6)=d6
        ELSE
          d(2)=dp2
          d(4)=9.D0
          d(6)=36.D0
        END IF
        konv=.FALSE.
        IF(jp1.GT.3) konv=.TRUE.
        IF(jp1.GT.7) THEN    
          lp1=7
          fc=0.6D0*fc
          d(7)=64.D0
        ELSE
          lp1=jp1
C          d(lp1)=DFLOAT(m**2)
          d(lp1)=FLOAT(m**2)
        END IF
        m=2*m
C        g=h0/DFLOAT(m)
        g=h0/FLOAT(m)
        b=2.D0*g
        IF(bh.AND.jp1.LT.9) THEN
          DO i=1,n
            ym(i)=yh(jp1,i)
            yl(i)=yg(jp1,i)
          END DO
        ELSE
C                                                      Calcul de T(H)
          kk=(m-2)/2
          m=m-1
          DO  i=1,n
            yl(i)=ya(i)
            ym(i)=ya(i)+g*dz(i)
          END DO
          DO k=1,m
C            CALL F(x+DFLOAT(k)*g,ym,dy,switch)
            CALL F(x+FLOAT(k)*g,ym,dy,switch)
            IF(switch) RETURN
            DO i=1,N
              u=yl(i)+b*dy(i)
              yl(i)=ym(i)
              ym(i)=u
C              u=DABS(u)
              u=ABS(u)
              IF(u.GT.s(i))  s(i)=u
            END DO  
            IF(k.NE.kk.OR.k.EQ.2) GO TO 100
            jjp1=1+jjp1
            DO i=1,n
              yh(jjp1,i)=ym(i)
              yg(jjp1,i)=yl(i)
            END DO  
 100        CONTINUE
          END DO  
        END IF   
        CALL F(a,ym,dy,switch)
        IF(switch) RETURN
        DO 160 i=1,n
          v=dt(i,1)
          ta=(ym(i)+yl(i)+g*dy(i))/2.D0
C                                 Calcul de la valeur extrapolee T(0)
          c=ta
          dt(i,1)=ta
          IF(lp1.GE.2) THEN     
            DO kp1=2,lp1
              b1=d(kp1)*v
              b=b1-c
              u=v
              IF(b.NE.0.D0) THEN     
                b=(c-v)/b
                u=c*b
                c=b1*b
              END IF
              v=dt(i,kp1)
              dt(i,kp1)=u
              ta=u+ta
            END DO  
          END IF
          IF(ABS(y(i)-ta).GT.eps*s(i)) konv=.FALSE.
C          IF(DABS(y(i)-ta).GT.eps*s(i)) konv=.FALSE.
          y(i)=ta
 160    CONTINUE          
        IF(konv) GO TO 180
        d(3)=4.D0
        d(5)=16.D0
        bo=.NOT.bo
        m=r
        r=sr
        sr=m*2
 170  CONTINUE
C
C     Plus de 9 iterations etaient necessaires. On recommence en
C     divisant le pas par 2.
C     
      bh=.NOT.bh
      h0=h0/2.D0
      GO TO 20
C 
 180  h0=fc*h0
      x=a
      RETURN
      END SUBROUTINE
c
C**************************** S A K E N *****************************
C        SAKEN   VERSION 6.2   31/10/75 - 5/1/00      A. SALIN
C
C        Zero of a function
C
C        Double precision
C
C        Start with an upper (RSUP) and lower (RINF) bound of the
C        root.
C        Firstly, the interval [RINF,RSUP] is decreased such that
C        |RSUP-RINF| < DR. Then the root is found with precision EPS
C        by successive interpolations using Aitken scheme.
C
C        Constraints:
C             The first part assumes that the root is not located at
C             an extremum of the function.
C        The function is given in a FUNCTION.
C
C        RINF: lower bound on R.
C        RSUP: upper bound on R.
C        DR: interval between RINF and RSUP before starting
C        interpolation.
C        EPS: required precision on R.
C        F: name of FUNCTION procedure defining the function.
C        R: set to the root position.
C        KOD:
C                  On input:
C             KOD=0 : do interpolation even if the root is not
C        between RINF and RSUP after decreasing the interval to DR.
C             KOD=-1 : message written if the root is not in the
C        defined interval.
C                  On output:
C             KOD=0 : convergence obtained.
C             KOD=1 : no single root between RINF and RSUP
C             KOD=2 : interpolation has not converged
C        ICONS: unit for error messages
C
      SUBROUTINE SAKEN(rinf,rsup,dr,eps,F,r,kod,icons)
      PARAMETER (nmax=100)
C      REAL*8 rsup,rinf,dr,eps,r1,r2,r3,y1,y2,y3,a,ext,x,r,xx,F,prod
      REAL rsup,rinf,dr,eps,r1,r2,r3,y1,y2,y3,a,ext,x,r,xx,F,prod
      DIMENSION a(nmax),x(nmax)
C
      r1=rinf
      r2=rsup
      y1=F(r1)
      IF(y1.EQ.0.d0) THEN
        r=r1
        kod=0
        RETURN
      END IF
      y2=F(r2)
      IF(y2.EQ.0.d0) THEN
        r=r2
        kod=0
        RETURN
      END IF
C                                                     Reduce interval
   11 CONTINUE
      IF(ABS(r2-r1).GT.ABS(dr))  THEN
C      IF(DABS(r2-r1).GT.DABS(dr))  THEN
        r3=(r1+r2)/2.D0
        y3=F(r3)
        IF(y3.EQ.0.d0) THEN
          r=r3
          kod=0
          RETURN
        END IF
        IF((y3*y1).LE.0.D0) THEN
          y2=y3
          r2=r3
        ELSE
          y1=y3
          r1=r3
        END IF
      GO TO 11
      END IF
C                                              Check location of root
      IF(y1*y2.GT.0.D0) THEN
        IF(kod.EQ.-1) THEN
          WRITE(icons,3) r1,y1,r2,y2
    3 FORMAT('SAKEN: Root not between RINF and RSUP',/,' R1=',
     1       1PD12.6,3X,'Y1=',D12.6,3X,'R2=',D12.6,3X,'Y2=',D12.6)
          kod=1
          RETURN
        ELSE IF(kod.NE.0) THEN
          kod=1
          RETURN
        ELSE
          kod=1
        END IF
      END IF
C                                                       Interpolation
      a(1)=r1
      x(1)=y1
      ext=r1
      x(2)=y2
      r=(r2-r1)/(y2-y1)
      a(2)=r
      prod=-y1
      r=ext+prod*r
C      IF(DABS(r-ext).LE.(DABS(r)*eps)) THEN
      IF(ABS(r-ext).LE.(ABS(r)*eps)) THEN
        IF(kod.NE.1) kod=0
        RETURN
      END IF
      ext=r
      DO n=3,nmax
        nm=n-1
        x(n)=F(r)
        xx=x(n)
        DO i=1,nm
          r=(r-a(i))/(xx-x(i))
        END DO
        a(n)=r
        prod=-prod*x(nm)
        r=ext+prod*r
C      IF(DABS(r-ext).LE.(DABS(r)*eps)) THEN
      IF(ABS(r-ext).LE.(ABS(r)*eps)) THEN
        IF(kod.NE.1) kod=0
        RETURN
      END IF
        ext=r
      END DO
C                                            Convergence not achieved
      WRITE(icons,2) nmax
    2 FORMAT('SAKEN: more than ',i2,' points required')
      kod=2
      RETURN
      END SUBROUTINE
C     *********************** T R I C Y C L E ***********************
      SUBROUTINE TRICYCLE(alpha,beta,gamma,cbot,ctop,b,x,n,z)
C
C     Solution of tridiagonal cyclic system
C     Required subroutine: TRIDIA
C
C         A. SALIN    Version 1    9/8/98
C
C     ALPHA: main diagonal (destroyed in TRICYCLE)
C     BETA: subdiagonal
C     GAMMA: superdiagonal
C     CBOT: element of bottom left corner of matrix
C     CTOP: element of top right corner of matrix
C     B: righthand side (destroyed in TRICYCLE)
C     X: solution
C     N: dimension of system
C     Z: scratch vector of dimension N at least

C      IMPLICIT REAL*8 (a-h,o-z)
      IMPLICIT REAL (a-h,o-z)
      DIMENSION alpha(n),beta(n),gamma(n),b(n),x(n),z(n)
C
      alp=-alpha(1)
      alpha(1)=alpha(1)-alp
      alpha(n)=alpha(n)-cbot*ctop/alp
      DO i=1,n
        z(i)=alpha(i)
      END DO
      CALL TRIDIA(alpha,beta,gamma,b,x,n)
C
      DO i=1,n
        alpha(i)=z(i)
      END DO
      b(1)=alp
      b(n)=cbot
      DO i=2,n-1
        b(i)=0.D0
      END DO
      CALL TRIDIA(alpha,beta,gamma,b,z,n)
      fact=(x(1)+ctop*x(n)/alp)/(1.D0+z(1)+ctop*z(n)/alp)
      DO i=1,n
        x(i)=x(i)-fact*z(i)
      END DO

      END SUBROUTINE
C**************************** T R I D I A *************************** 
      SUBROUTINE TRIDIA(alpha,beta,gamma,b,x,n) 
C 
C        Solution of tridiagonal systems 
C 
C     ALPHA: main diagonal (destroyed in TRIDIA). 
C     BETA: subdiagonal (first element has index 1). 
C     GAMMA: superdiagonal (first element has index 1). 
C     B: right-hand side (destroyed in TRIDIA). 
C     X: solution. 
C     N: dimension of system. 
C 
C      IMPLICIT REAL*8 (a-h,o-z) 
      IMPLICIT REAL (a-h,o-z)
      DIMENSION alpha(n),beta(n),gamma(n),b(n),x(n) 
C 
      DO i=2,n 
        rap = beta(i-1) / alpha(i-1) 
        alpha(i) = alpha(i) - rap*gamma(i-1) 
        b(i) = b(i) - rap*b(i-1) 
      END DO 
      x(n) = b(n) / alpha(n) 
       
      DO j=2,n 
        i = n-j+1 
        x(i) = ( b(i) - gamma(i)*x(i+1) ) / alpha(i) 
      END DO 
      RETURN 
      END SUBROUTINE
C     ************************ T R I S P L ***************************

C     Spline interpolation of a periodic function with period
C     [x(n)-x(1)] - i.e. y(n)=y(1).

C     See general comments in DSPLIN.
C     The vectors ALPHA, BETA, GAMMA, B, BZ are scratch vectors of
C     dimension at list N.
C     Required subroutine: TRIDIA, TRICYCLE
C
C        Author: A. Salin   Version 1  09/06/98
C
C      BLOCK DATA BLKTRI
C      COMMON/KODSPL/kod,k,iw,iligne
C      DATA iw,iligne,k/6,100,2/
C      END
C     ***************************************************************
      SUBROUTINE TRISPL(n,x,y,cm,alpha,beta,gamma,b,bz)

C      IMPLICIT REAL*8 (a-h,o-z)
      IMPLICIT REAL (a-h,o-z)
      PARAMETER (six=6.D0,trois=3.D0,zero=0.D0)
C
      DIMENSION x(n),y(n),cm(n)
      DIMENSION alpha(n),beta(n),gamma(n),b(n),bz(n)
C
      COMMON/KODSPL/kod,k,iw,iligne
C
   91 FORMAT('Error in the data for TRISPL',/,'Check abcissae:',
     1       1PD15.8,1X,' and ',D15.8,/,'Program stopped')
C
      kod=0
      c=x(2)-x(1)
      e=(y(2)-y(1))/c
      DO i=3,n
        i2=i-2
        a=c
        c=x(i)-x(i-1)
        IF(a*c.LE.zero) THEN
          kod=-1
          IF(iw.GT.0) WRITE(iw,91) x(i-1),x(i)
          STOP
        END IF
        alpha(i2)=(a+c)/trois
        beta(i2)=c/six
        gamma(i2)=beta(i2)
        d=e
        e=(y(i)-y(i-1))/c
        b(i2)=e-d
      END DO
C     
      a=c
      c=x(2)-x(1)
      alpha(n-1)=(a+c)/trois
      beta(n-1)=c/six
      gamma(n-1)=beta(n-1)
      d=e
      e=(y(2)-y(1))/c
      b(n-1)=e-d
      cbot=beta(n-1)
      ctop=cbot
C
C                                            Solve tridiagonal system
      CALL TRICYCLE(alpha,beta,gamma,cbot,ctop,b,cm(2),n-1,bz)
C
      cm(1)=cm(n)
C
      k=2
      END SUBROUTINE
