C ******************************************
C     VERSION OF JANUARY 26, 1996      
C ******************************************
C
      SUBROUTINE DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
     &          IP1(NM1),IPHES(N)
      LOGICAL CALHES
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      GOTO (1,2,55,55,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
C -----------------------------------------------------------
C
   1  CONTINUE
C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      DO J=1,N
         DO  I=1,N
            E1(I,J)=-FJAC(I,J)
         END DO
         E1(J,J)=E1(J,J)+FAC1
      END DO
      CALL DEC (N,LDE1,E1,IP1,IER)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
      RETURN
C
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE DECOMR
C
C ***********************************************************
C
      SUBROUTINE DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &            M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),
     &          E2R(LDE1,NM1),E2I(LDE1,NM1),IP2(NM1)
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      GOTO (1,2,55,55,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
C -----------------------------------------------------------
C
   1  CONTINUE
C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      DO J=1,N
         DO I=1,N
            E2R(I,J)=-FJAC(I,J)
            E2I(I,J)=0.D0
         END DO
         E2R(J,J)=E2R(J,J)+ALPHN
         E2I(J,J)=BETAN
      END DO
      CALL DECC (N,LDE1,E2R,E2I,IP2,IER)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
      RETURN
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE DECOMC
C
C ***********************************************************
C
      SUBROUTINE SLVRAR(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          M1,M2,NM1,FAC1,E1,LDE1,Z1,F1,IP1,IPHES,IER,IJOB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
     &          IP1(NM1),IPHES(N),Z1(N),F1(N)
        common/datafft/ data(256*256*4)
        COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      GOTO (1,2,55,55,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
C -----------------------------------------------------------
C
   1  CONTINUE
C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      DO I=1,N
         Z1(I)=Z1(I)-F1(I)*FAC1
      END DO
      CALL SOL (N,LDE1,E1,Z1,IP1)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
      DO I=1,N
         Z1(I)=Z1(I)-F1(I)*FAC1
      END DO
C      CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
      DO I=1,N
         DATA(2*I-1)=z1(I)
         DATA(2*I)=0.d0
      END DO
         N2=2*NS
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=FAC1+4.0D0*ALF*NSSQ
         DO I=1,NS
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NS
               TERM=(TI-TCOS(J))*NSSQ
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATA(IJ-1)=DATA(IJ-1)/TERM
               DATA(IJ)=DATA(IJ)/TERM
               DATA(IJ2-1)=DATA(IJ2-1)/TERM
               DATA(IJ2)=DATA(IJ2)/TERM
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         z1(I)=DATA(2*I-1)
      END DO
      RETURN
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE SLVRAR
C
C ***********************************************************
C
      SUBROUTINE SLVRAI(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,Z2,Z3,
     &          F2,F3,CONT,IP2,IPHES,IER,IJOB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),
     &          IP2(NM1),IPHES(N),Z2(N),Z3(N),F2(N),F3(N)
      DIMENSION E2R(LDE1,NM1),E2I(LDE1,NM1)
        common/datafft/ data(256*256*4)
        COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      GOTO (1,2,55,55,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
C -----------------------------------------------------------
C
   1  CONTINUE
C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      DO I=1,N
         S2=-F2(I)
         S3=-F3(I)
         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
      END DO
      CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
      DO I=1,N
         S2=-F2(I)
         S3=-F3(I)
         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
      END DO
C      CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
      DO I=1,N
         DATA(2*I-1)=z2(I)
         DATA(2*I)=z3(I)
      END DO
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=ALPHN+4.0D0*ALF*NSSQ
         N2=2*NS
         DO I=1,NS
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NS
               TERMR=(TI-TCOS(J))*NSSQ
               TERMI=BETAN*NSSQ
               TNORM=TERMR**2+TERMI**2
               TERMR=TERMR/TNORM
               TERMI=-TERMI/TNORM
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATASA=DATA(IJ-1)*TERMR-DATA(IJ)*TERMI
               DATA(IJ)=DATA(IJ)*TERMR+DATA(IJ-1)*TERMI
               DATA(IJ-1)=DATASA
               DATASA=DATA(IJ2-1)*TERMR-DATA(IJ2)*TERMI
               DATA(IJ2)=DATA(IJ2)*TERMR+DATA(IJ2-1)*TERMI
               DATA(IJ2-1)=DATASA
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         z2(I)=DATA(2*I-1)
         z3(I)=DATA(2*I)
      END DO
      RETURN
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE SLVRAI
C
C ***********************************************************
C
      SUBROUTINE SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3,
     &          F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
     &          E2R(LDE1,NM1),E2I(LDE1,NM1),IP1(NM1),IP2(NM1),
     &          IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),F3(N)
        common/datafft/ data(256*256*4)
        COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      GOTO (1,2,55,55,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
C -----------------------------------------------------------
C
   1  CONTINUE
C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      DO I=1,N
         S2=-F2(I)
         S3=-F3(I)
         Z1(I)=Z1(I)-F1(I)*FAC1
         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
      END DO
      CALL SOL (N,LDE1,E1,Z1,IP1)
      CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
      DO I=1,N
         S2=-F2(I)
         S3=-F3(I)
         Z1(I)=Z1(I)-F1(I)*FAC1
         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
      END DO
C      CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
C      CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
      DO I=1,N
         DATA(2*I-1)=z1(I)
         DATA(2*I)=0.d0
      END DO
         N2=2*NS
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=FAC1+4.0D0*ALF*NSSQ
         DO I=1,NS
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NS
               TERM=(TI-TCOS(J))*NSSQ
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATA(IJ-1)=DATA(IJ-1)/TERM
               DATA(IJ)=DATA(IJ)/TERM
               DATA(IJ2-1)=DATA(IJ2-1)/TERM
               DATA(IJ2)=DATA(IJ2)/TERM
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         z1(I)=DATA(2*I-1)
      END DO
      DO I=1,N
         DATA(2*I-1)=z2(I)
         DATA(2*I)=z3(I)
      END DO
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=ALPHN+4.0D0*ALF*NSSQ
         DO I=1,NS
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NS
               TERMR=(TI-TCOS(J))*NSSQ
               TERMI=BETAN*NSSQ
               TNORM=TERMR**2+TERMI**2
               TERMR=TERMR/TNORM
               TERMI=-TERMI/TNORM
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATASA=DATA(IJ-1)*TERMR-DATA(IJ)*TERMI
               DATA(IJ)=DATA(IJ)*TERMR+DATA(IJ-1)*TERMI
               DATA(IJ-1)=DATASA
               DATASA=DATA(IJ2-1)*TERMR-DATA(IJ2)*TERMI
               DATA(IJ2)=DATA(IJ2)*TERMR+DATA(IJ2-1)*TERMI
               DATA(IJ2-1)=DATASA
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         z2(I)=DATA(2*I-1)
         z3(I)=DATA(2*I)
      END DO
      RETURN
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE SLVRAD
C
C ***********************************************************
C
      SUBROUTINE ESTRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,
     &          E1,LDE1,Z1,Z2,Z3,CONT,F1,F2,IP1,IPHES,SCAL,ERR,
     &          FIRST,REJECT,FAC1,RPAR,IPAR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1),
     &     SCAL(N),IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),Y0(N),Y(N)
      DIMENSION CONT(N),RPAR(1),IPAR(1)
      LOGICAL FIRST,REJECT
        common/datafft/ data(256*256*4)
        COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
      HEE1=DD1/H
      HEE2=DD2/H
      HEE3=DD3/H
      GOTO (1,2,55,55,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
   1  CONTINUE
C ------  B=IDENTITY, JACOBIAN A FULL MATRIX
      DO  I=1,N 
         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
         CONT(I)=F2(I)+Y0(I)
      END DO
      CALL SOL (N,LDE1,E1,CONT,IP1) 
      GOTO 77
C
   2  CONTINUE
C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX
      DO I=1,N 
         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
         CONT(I)=F2(I)+Y0(I)
      END DO
      DO I=1,N
         DATA(2*I-1)=CONT(I)
         DATA(2*I)=0.d0
      END DO
         N2=2*NS
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=FAC1+4.0D0*ALF*NSSQ
         DO I=1,NS
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NS
               TERM=(TI-TCOS(J))*NSSQ
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATA(IJ-1)=DATA(IJ-1)/TERM
               DATA(IJ)=DATA(IJ)/TERM
               DATA(IJ2-1)=DATA(IJ2-1)/TERM
               DATA(IJ2)=DATA(IJ2)/TERM
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         CONT(I)=DATA(2*I-1)
      END DO
      GOTO 77
C
C --------------------------------------
C
  77  CONTINUE
      ERR=0.D0
      DO  I=1,N
         ERR=ERR+(CONT(I)/SCAL(I))**2
      END DO
      ERR=MAX(SQRT(ERR/N),1.D-10)
C
      IF (ERR.LT.1.D0) RETURN
      IF (FIRST.OR.REJECT) THEN
          DO I=1,N
             CONT(I)=Y(I)+CONT(I)
          END DO
          CALL FCN(N,X,CONT,F1,RPAR,IPAR)
          NFCN=NFCN+1
          DO I=1,N
             CONT(I)=F1(I)+F2(I)
          END DO
          GOTO (31,32,31,32,31,32,55,55,55,55,31,32,31,32,31), IJOB
C ------ FULL MATRIX OPTION
  31      CONTINUE
          CALL SOL(N,LDE1,E1,CONT,IP1) 
          GOTO 88
C ------ BANDED MATRIX OPTION
 32      CONTINUE
      DO I=1,N
         DATA(2*I-1)=CONT(I)
         DATA(2*I)=0.d0
      END DO
         N2=2*NS
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=FAC1+4.0D0*ALF*NSSQ
         DO I=1,NS
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NS
               TERM=(TI-TCOS(J))*NSSQ
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATA(IJ-1)=DATA(IJ-1)/TERM
               DATA(IJ)=DATA(IJ)/TERM
               DATA(IJ2-1)=DATA(IJ2-1)/TERM
               DATA(IJ2)=DATA(IJ2)/TERM
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         CONT(I)=DATA(2*I-1)
      END DO
         GOTO 88
C -----------------------------------
   88     CONTINUE
          ERR=0.D0 
          DO I=1,N
             ERR=ERR+(CONT(I)/SCAL(I))**2
          END DO
          ERR=MAX(SQRT(ERR/N),1.D-10)
       END IF
       RETURN
C -----------------------------------------------------------
  55   CONTINUE
       RETURN
       END
C
C     END OF SUBROUTINE ESTRAD
C
C ***********************************************************
C
      SUBROUTINE ESTRAV(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          H,DD,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,NS,NNS,
     &          E1,LDE1,ZZ,CONT,FF,IP1,IPHES,SCAL,ERR,
     &          FIRST,REJECT,FAC1,RPAR,IPAR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1),
     &     SCAL(N),IPHES(N),ZZ(NNS),FF(NNS),Y0(N),Y(N)
      DIMENSION DD(NS),CONT(N),RPAR(1),IPAR(1)
      LOGICAL FIRST,REJECT
        common/datafft/ data(256*256*4)
        COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NP,NSSQ
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
      GOTO (1,2,55,55,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
   1  CONTINUE
C ------  B=IDENTITY, JACOBIAN A FULL MATRIX
      DO  I=1,N
         SUM=0.D0
         DO K=1,NS
            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
         END DO
         FF(I+N)=SUM/H
         CONT(I)=FF(I+N)+Y0(I)
      END DO
      CALL SOL (N,LDE1,E1,CONT,IP1) 
      GOTO 77
C
   2  CONTINUE
C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX
      DO  I=1,N
         SUM=0.D0
         DO K=1,NS
            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
         END DO
         FF(I+N)=SUM/H
         CONT(I)=FF(I+N)+Y0(I)
      END DO
      DO I=1,N
         DATA(2*I-1)=CONT(I)
         DATA(2*I)=0.d0
      END DO
         N2=2*NP
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=FAC1+4.0D0*ALF*NSSQ
         DO I=1,NP
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NP
               TERM=(TI-TCOS(J))*NSSQ
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATA(IJ-1)=DATA(IJ-1)/TERM
               DATA(IJ)=DATA(IJ)/TERM
               DATA(IJ2-1)=DATA(IJ2-1)/TERM
               DATA(IJ2)=DATA(IJ2)/TERM
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         CONT(I)=DATA(2*I-1)
      END DO
      GOTO 77
C
C --------------------------------------
C
  77  CONTINUE
      ERR=0.D0
      DO  I=1,N
         ERR=ERR+(CONT(I)/SCAL(I))**2
      END DO
      ERR=MAX(SQRT(ERR/N),1.D-10)
C
      IF (ERR.LT.1.D0) RETURN
      IF (FIRST.OR.REJECT) THEN
          DO I=1,N
             CONT(I)=Y(I)+CONT(I)
          END DO
          CALL FCN(N,X,CONT,FF,RPAR,IPAR)
          NFCN=NFCN+1
          DO I=1,N
             CONT(I)=FF(I)+FF(I+N)
          END DO
          GOTO (31,32,31,32,31,32,55,55,55,55,31,32,31,32,31), IJOB
C ------ FULL MATRIX OPTION
 31      CONTINUE
         CALL SOL (N,LDE1,E1,CONT,IP1) 
          GOTO 88
C ------ BANDED MATRIX OPTION
 32      CONTINUE
      DO I=1,N
         DATA(2*I-1)=CONT(I)
         DATA(2*I)=0.d0
      END DO
         N2=2*NP
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=FAC1+4.0D0*ALF*NSSQ
         DO I=1,NP
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NP
               TERM=(TI-TCOS(J))*NSSQ
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATA(IJ-1)=DATA(IJ-1)/TERM
               DATA(IJ)=DATA(IJ)/TERM
               DATA(IJ2-1)=DATA(IJ2-1)/TERM
               DATA(IJ2)=DATA(IJ2)/TERM
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         CONT(I)=DATA(2*I-1)
      END DO
         GOTO 88
C -----------------------------------
  88      CONTINUE
          ERR=0.D0 
          DO I=1,N
             ERR=ERR+(CONT(I)/SCAL(I))**2
          END DO
          ERR=MAX(SQRT(ERR/N),1.D-10)
       END IF
       RETURN
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
       END
C
C     END OF SUBROUTINE ESTRAV
C
C ***********************************************************
C
      SUBROUTINE SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          M1,M2,NM1,FAC1,E,LDE,IP,DY,AK,FX,YNEW,HD,IJOB,STAGE1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1),
     &          IP(NM1),DY(N),AK(N),FX(N),YNEW(N)
      LOGICAL STAGE1
        common/datafft/ data(256*256*4)
        COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      IF (HD.EQ.0.D0) THEN
         DO  I=1,N
           AK(I)=DY(I)
         END DO
      ELSE
         DO I=1,N
            AK(I)=DY(I)+HD*FX(I)
         END DO
      END IF
C
      GOTO (1,2,55,55,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
C -----------------------------------------------------------
C
   1  CONTINUE
C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      IF (STAGE1) THEN
         DO I=1,N
            AK(I)=AK(I)+YNEW(I)
         END DO
      END IF
      CALL SOL (N,LDE,E,AK,IP)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
C      IF (STAGE1) THEN
C         DO I=1,N
C            AK(I)=AK(I)+YNEW(I)
C         END DO
C      END IF
C      CALL SOLB (N,LDE,E,MLE,MUE,AK,IP)
      IF (STAGE1) THEN
         DO I=1,N
            AK(I)=AK(I)+YNEW(I)
         END DO
      END IF
      DO I=1,N
         DATA(2*I-1)=AK(I)
         DATA(2*I)=0.d0
      END DO
         N2=2*NS
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=FAC1+4.0D0*ALF*NSSQ
         DO I=1,NS
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NS
               TERM=(TI-TCOS(J))*NSSQ
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATA(IJ-1)=DATA(IJ-1)/TERM
               DATA(IJ)=DATA(IJ)/TERM
               DATA(IJ2-1)=DATA(IJ2-1)/TERM
               DATA(IJ2)=DATA(IJ2)/TERM
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         AK(I)=DATA(2*I-1)
      END DO
      RETURN
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE SLVROD
C
C
C ***********************************************************
C
      SUBROUTINE SLVSEU(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          M1,M2,NM1,FAC1,E,LDE,IP,IPHES,DEL,IJOB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1),DEL(N)
      DIMENSION IP(NM1),IPHES(N)
        common/datafft/ data(256*256*4)
        COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      GOTO (1,2,1,2,55,55,55,55,55,55,55,55,55,55,55), IJOB
C
C -----------------------------------------------------------
C
   1  CONTINUE
C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      CALL SOL (N,LDE,E,DEL,IP)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
      DO I=1,N
         DATA(2*I-1)=DEL(I)
         DATA(2*I)=0.d0
      END DO
         N2=2*NS
         CALL FOURN (DATA,NF,NDIM,-1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1)
         FAC1N=FAC1+4.0D0*ALF*NSSQ
         DO I=1,NS
            TI=FAC1N-TCOS(I)
            II=I*2
            DO J=1,NS
               TERM=(TI-TCOS(J))*NSSQ
               IJ=II+(J-1)*N2
               IJ2=2*NSSQ+IJ
               DATA(IJ-1)=DATA(IJ-1)/TERM
               DATA(IJ)=DATA(IJ)/TERM
               DATA(IJ2-1)=DATA(IJ2-1)/TERM
               DATA(IJ2)=DATA(IJ2)/TERM
            END DO
         END DO
         CALL FOURN (DATA,NF,NDIM,1)
         CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1)
      DO I=1,N
         DEL(I)=DATA(2*I-1)
      END DO
      RETURN
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE SLVSEU
C
C
      SUBROUTINE FOURN (DATA,NN,NDIM,ISIGN)
C --- explication in "Numerical Recipies", page 451
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION NN(NDIM),DATA(*)
      NTOT=1
      DO IDIM=1,NDIM
         NTOT=NTOT*NN(IDIM)
      END DO
      NPREV=1
      DO IDIM=1,NDIM
         N=NN(IDIM)
         NREM=NTOT/(N*NPREV)
         IP1=2*NPREV
         IP2=IP1*N
         IP3=IP2*NREM
         I2REV=1
         DO I2=1,IP2,IP1
            IF (I2.LT.I2REV) THEN
               DO I1=I2,I2+IP1-2,2
                  DO I3=I1,IP3,IP2
                     I3REV=I2REV+I3-I2
                     TEMPR=DATA(I3)
                     TEMPI=DATA(I3+1)
                     DATA(I3)=DATA(I3REV)
                     DATA(I3+1)=DATA(I3REV+1)
                     DATA(I3REV)=TEMPR
                     DATA(I3REV+1)=TEMPI
                  END DO
               END DO
            END IF
            IBIT=IP2/2
  1         IF ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) THEN
               I2REV=I2REV-IBIT
               IBIT=IBIT/2
               GO TO 1
            END IF
            I2REV=I2REV+IBIT
         END DO
         IFP1=IP1
  2      IF (IFP1.LT.IP2) THEN
            IFP2=2*IFP1
            THETA=ISIGN*6.2831853071795864769252D0/(IFP2/IP1)
            WPR=-2.D0*SIN(0.5D0*THETA)**2
            WPI=SIN(THETA)
            WR=1.0D0
            WI=0.0D0
            DO I3=1,IFP1,IP1
               DO I1=I3,I3+IP1-2,2
                  DO I2=I1,IP3,IFP2
                     K1=I2
                     K2=K1+IFP1
                     TEMPR=WR*DATA(K2)-WI*DATA(K2+1)
                     TEMPI=WR*DATA(K2+1)+WI*DATA(K2)
                     DATA(K2)=DATA(K1)-TEMPR
                     DATA(K2+1)=DATA(K1+1)-TEMPI
                     DATA(K1)=DATA(K1)+TEMPR
                     DATA(K1+1)=DATA(K1+1)+TEMPI
                  END DO
               END DO
               WTEMP=WR
               WR=WR*WPR-WI*WPI+WR
               WI=WI*WPR+WTEMP*WPI+WI
            END DO
            IFP1=IFP2
            GO TO 2
         END IF
         NPREV=N*NPREV
      END DO
      RETURN
      END
