      PROGRAM PLSADR
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
C
C* Begin Prologue PLSADR
C
C ---------------------------------------------------------------------
C
C* Purpose           Driver testing MEXX with
C                    example pendulum (PL) in sparse (S) matrix mode
C                    using the following options (A)
C                     - global solution output into file PLSADR.GLO
C                     - dense output of the position values at 60
C                       equidistant time intervals
C                     - stop at a root of switching function plswi.f
C                     - time monitor on
C
C* File              plsadr.f
C* Version           1.1
C* Latest Change     96/07/19 (1.6)
C* Modification history
C    1.1             Correct workspace formula.
C                    Avoid free format because of compiler dependency.
C* Library           CodeLib
C* Code              Fortran 77
C                    Double Precision
C* Environment       Standard version for FORTRAN77 environments on
C                    PCs, workstations, and hosts
C
C* Subroutines supplied
C
C     PLSIN          Provide initial values
C     PLSF           Specify the problem
C     PLDOUT         Output interpolated solution values
C     PLSWI          Switching function to stop integration in case of
C                    special solution values
C      
C  Problem dimensions for a maximum of 30 bodies.
C  The actual number of bodies will be requested from standard input.
C
      PARAMETER (
     $     NBODY  = 30,
     $     NPDIM  = 3*NBODY,
     $     NVDIM  = NPDIM,
     $     NUDIM  = 1,
     $     NGDIM  = 2*NBODY,
     $     NLDIM  = NGDIM,
     $     NPVUD  = NPDIM*2 + NUDIM,
     $     NZGDIM = 8*NBODY - 4,
     $     NSWDIM = NBODY - 1
     $     )
C      INTEGER     NL1, NU1, NG1, NGL
C      PARAMETER ( NL1 = MAX (1,NLDIM) ,
C     $            NU1 = MAX (1,NUDIM) ,
C     $            NG1 = MAX (1,NGDIM) ,
C     $            NGL = MAX (1,NGDIM,NLDIM) )
C
C  Workspace dimension
C
      PARAMETER (
     $     NZA    = 2*NZGDIM + NPDIM,
     $     LIWK   = NPDIM + 4*NVDIM + NLDIM + NUDIM + 60 +
     $              12*(NVDIM+NLDIM) + 2*NZGDIM + 5.1*NZA +
     $              5*NSWDIM,
     $     LD     = NVDIM + NLDIM+1 + 2*(2*NZGDIM+NVDIM),
     $     LX     = 12*(NPDIM+3*NVDIM+2*(NLDIM+1)+NUDIM+1),
     $     LO     = (NPDIM+2*NVDIM+NUDIM+NLDIM)*156,
     $     LRWK   = LD + NPDIM**2 + NZGDIM + LX + LO + 6*NSWDIM +
     $              (6+NGDIM+NLDIM)*NPDIM + 9*NVDIM + 6*(NUDIM+1) +
     $               6*(NLDIM+1) + NGDIM+1 + 50
     $     )
C
      EXTERNAL PLSIN, PLSF, DUOUT, PLDOUT, PLSWI
C
      INTEGER MXJOB(150), IWK(LIWK)
      DOUBLE PRECISION RWK(LRWK)
      DIMENSION P(NPDIM), V(NVDIM), RTOL(NPVUD), ATOL(NPVUD)
      DIMENSION A(NVDIM), RLAM(NLDIM), U(NUDIM+1)
C
      COMMON / PEND / PMASS(NBODY), PLEN(NBODY), NPEND
C
      CHARACTER*30 PROTXT
      DATA PROTXT /'Pendulum'/
C
      CHARACTER*67 CMODUL
      DATA CMODUL /'@(#) plsadr.f 1.6 from 96/07/19 (mexx 1.2)
     $     '/
C
      DATA MXJOB /150*0/, RWK /LRWK*0.0D0/, IWK /LIWK*0/
C
C -- Machine constants
      LUSIN = 5
      LUSOUT = 6
C
C -- Start message
      WRITE (LUSOUT,9000) ' '
      WRITE (LUSOUT,9000)
     $     ' ===================================================='
      WRITE (LUSOUT,9000) CMODUL
      WRITE (LUSOUT,9000)
     $     ' ===================================================='
      WRITE (LUSOUT,9000) ' '
C
C -- Initialize problem
      NP = NPDIM
      NG = NGDIM
      NU = NUDIM
      NV=NP
      NL=0
      NBLK=0
      NMRC=0
      NZGMAX=0
      NZFMAX=0
      NSWMAX=0
C
      CALL PLSIN (NP, NV, NL, NG, NU, NBLK, NMRC, NZGMAX,
     $     NZFMAX, NSWMAX, T0, P, V, U, A, RLAM, PROTXT)
C
C
C
C  Set options for integration
C
      T=T0
      TEND = 3.D0
      H = 1.D-3
      ITOL = 0
      TOL = 1.D-3
      RTOL(1) = TOL
      ATOL(1) = TOL
C
C  Select options for sparse linear algebra
C
      MMODE = 1
      JOB = 10
      IF (JOB .EQ. 10) THEN
         MGPFL = 1
      ELSE
         MGPFL = 0
      ENDIF
C     
      MXJOB(2) = JOB
      MXJOB(3) = MMODE
      MXJOB(4) = NBLK 
      MXJOB(5) = NMRC
      MXJOB(6) = MGPFL
      MXJOB(7) = NZGMAX
      MXJOB(8) = NZFMAX
C
C  Select I/O options
C     - General print output
      MXJOB(11) = 2
      MXJOB(12) = LUSOUT
C     - Global solution
      MXJOB(15) = 1
      MXJOB(16) = 20
      OPEN (UNIT=20,FILE='PLSADR.GLO')
C     - Time monitor
      MXJOB(19) = 1
      MXJOB(20) = LUSOUT
C
C  Select options for dense output
C
      MDOUT = 1
      NDOUT = 59
      MDIND = 1
C
      MXJOB(31) = MDOUT
      MXJOB(32) = NDOUT
      MXJOB(33) = MDIND
      RWK(51) = 0.1D0
C
      IF (MDIND .GT. 0) THEN
C 
C        Components for dense output and global solution output
C
         NDP = NP
         IWK(51) = NDP
C
         INDP = 60
         INDV = INDP + NP
         INDU = INDV + NV
         INDA = INDU + NU
         INDL = INDA + NV
C
         DO 1000 I=1,NDP
            IWK(INDP + I) = I
 1000    CONTINUE
C
      ENDIF
C
C  Set options for root finder
C
      MSWIT = 1
C
C MSWIT, NSWMAX, NSUB
      IF (NSWMAX .GT. 0) THEN
         MXJOB(35) = MSWIT
         MXJOB(36) = NSWMAX
         MXJOB(37) = 0
      ENDIF
C
C  Document options
C
      WRITE (LUSOUT,9000) ' Settings chosen.         '
      WRITE (LUSOUT,9100) ' MOUT, MDOUT, MDIND:      ',
     $     MOUT, MDOUT, MDIND
      WRITE (LUSOUT,9200) ' Delta-t-max:             ',
     $     RWK(51)
C     
      WRITE (LUSOUT,9100) ' MMODE, JOB:',
     $     MMODE, JOB
C
      WRITE (LUSOUT,9200) ' RTOL, ATOL:              ',
     $     RTOL(1), ATOL(1)
C
      IDID = 0
C
      CALL MEXX (NP, NV, NL, NG, NU,
     $     PLSF, T, TEND, P, V, U, A, RLAM,
     $     ITOL, RTOL, ATOL,
     $     H, MXJOB, IDID, LIWK, IWK, LRWK, RWK,
     $     DUOUT, PLDOUT, PLSWI)
C     
      IF (IDID .EQ. 0) THEN
         WRITE (LUSOUT,9000) ' o.k.'
      ELSE IF (IDID .EQ. 1) THEN
         WRITE (LUSOUT,9200) ' Root found: Integration stopped at', T
         LISW = 60 + NP +3*NV + NU
         DO 1010 I=1,NSWMAX
            IF (IWK(LISW+I) .NE. 0) THEN
               WRITE (LUSOUT,9100) '             in component', I
            ENDIF
 1010    CONTINUE
      ELSE
         WRITE (LUSOUT,9100) ' Error return code from MEXX:', IDID
      ENDIF
C
      WRITE (LUSOUT,9100) ' '
      WRITE (LUSOUT,9100) 'C==================================='
      WRITE (LUSOUT,9100) ' '
      WRITE (LUSOUT,9200) ' TOL       ', RTOL(1)
      WRITE (LUSOUT,9100) ' '
      WRITE (LUSOUT,9100) ' Statistics:  NSTEP, NACCPT, NREJCT'
      WRITE (LUSOUT,9100) '           ', MXJOB(51), MXJOB(52), MXJOB(53)
      WRITE (LUSOUT,9100) ' '
      WRITE (LUSOUT,9100) '             NFPROB,   NFCN,   NGCN'
      WRITE (LUSOUT,9100) '           ', MXJOB(54), MXJOB(55), MXJOB(56)
      WRITE (LUSOUT,9100) ' '
      WRITE (LUSOUT,9100) '              NMULT,   NDEC,   NSOL'
      WRITE (LUSOUT,9100) '           ', MXJOB(57), MXJOB(58), MXJOB(59)
      WRITE (LUSOUT,9100) ' '
      WRITE (LUSOUT,9100) 'C==================================='
C      
C
      WRITE (LUSOUT,9000) ' '
      WRITE (LUSOUT,9000) ' End of example ', PROTXT
      WRITE (LUSOUT,9000) ' Bye'
C
      STOP
C
 9000 FORMAT (2A)
 9100 FORMAT (A, 4I8)
 9200 FORMAT (A, 2(1PD8.1))
C
C  end of program PLSADR
C
      END
C
      SUBROUTINE PLSF (NP, NV, NL, NG, NU, LDG,
     $     T, P, V, U, RLAM, AM, GP, F,
     $     PDOT, UDOT, G, GI, FL, QFLAG,
     $     NBLK, NMRC, NZGMAX, NZFMAX,
     $     IROWG, JCOLG, NZGACT, IPCG,
     $     IROWF, JCOLF, NZFACT, IPCF,
     $     IFAIL)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C  Input parameters:
      INTEGER NP, NV, NL, NG, NU, LDG
      LOGICAL QFLAG(9)
      DOUBLE PRECISION T, P(NP), V(NV), U(NU), RLAM(NL)
      INTEGER NBLK, NMRC, NZGMAX, NZFMAX
C  Output parameters:
      DOUBLE PRECISION AM(NBLK*NMRC,NMRC), GP(NZGMAX), F(NV), PDOT(NP),
     $     UDOT(NU), G(NG), GI(NG), FL(NZFMAX)
      INTEGER IROWG(NZGMAX), JCOLG(NZGMAX), NZGACT, IPCG, IROWF(NZFMAX),
     $     JCOLF(NZFMAX), NZFACT, IPCF, IFAIL
C
C* File              plsf.f
C* Version           1.1
C* Latest Change     96/05/03 (1.5)
C* Modification history
C    1.1             New free format starts text output with ' '.
C     
      PARAMETER (NBMAX = 30)
      COMMON / PEND / PMASS(NBMAX), PLEN(NBMAX), NPEND
C
      SAVE
C
      DATA GRAV/10./
C
C  Machine constant (logical unit stdout)
      LUSOUT = 6
C
      IFAIL = 0
C
      IF (QFLAG(1)) THEN
         DO 1010 J=1,NMRC
            DO 1000 I=1,NMRC*NBLK
               AM(I,J) = 0.D0
 1000       CONTINUE
 1010    CONTINUE
         DO 1020 I=1,NPEND
            L = 3*(I-1) + 1
            AM(1+(L-1)*NMRC,1) = PMASS(I)
            AM(1+  L  *NMRC,1) = PMASS(I)
            AM(1+(L+1)*NMRC,1) = PLEN(I)**2*PMASS(I)/3.D0
 1020    CONTINUE
      ENDIF
C     
      IF (QFLAG(2) .OR. QFLAG(3)) THEN
         LL = 0
         DO 1030 I=1,NPEND
            IROW = (I-1)*2 + 1
            ICOL = (I-1)*3 + 1
            LL = LL + 1
            GP(LL) = -1.D0
            IROWG(LL) = IROW
            JCOLG(LL) = ICOL
C     
            LL = LL + 1
            GP(LL) = PLEN(I)*COS(P(ICOL+2))
            IROWG(LL) = IROW
            JCOLG(LL) = ICOL + 2
C
            LL = LL + 1
            GP(LL) = -1.D0
            IROWG(LL) = IROW + 1
            JCOLG(LL) = ICOL + 1
C
            LL = LL + 1
            GP(LL) = PLEN(I)*SIN(P(ICOL+2))
            IROWG(LL) = IROW + 1
            JCOLG(LL) = ICOL + 2
C
            IF (I .GT. 1) THEN
C
               LL = LL + 1
               GP(LL) = 1.D0
               IROWG(LL) = IROW
               JCOLG(LL) = ICOL - 3
C     
               LL = LL + 1
               GP(LL) = PLEN(I-1)*COS(P(ICOL-1))
               IROWG(LL) = IROW
               JCOLG(LL) = ICOL - 1
C
               LL = LL + 1
               GP(LL) = 1.D0 
               IROWG(LL) = IROW + 1
               JCOLG(LL) = ICOL - 2
C
               LL = LL + 1
               GP(LL) = PLEN(I-1)*SIN(P(ICOL-1))
               IROWG(LL) = IROW + 1
               JCOLG(LL) = ICOL - 1
            ENDIF
 1030    CONTINUE
         NZGACT = LL
         IPCG = 0
      ENDIF         
      IF (QFLAG(4)) THEN
         DO 1040 I=1,NPEND
            L = 3*(I-1) + 1
            F(L) = 0.D0
            F(L+1) = -GRAV*PMASS(I)  
            F(L+2) = 0.D0
 1040    CONTINUE
      ENDIF
      IF (QFLAG(5)) THEN 
         DO 1050 I=1,NV
            PDOT(I)=V(I)
 1050    CONTINUE
      ENDIF
      IF (QFLAG(7)) THEN
         G(1) = -P(1) + PLEN(1)*SIN(P(3))
         G(2) = -P(2) - PLEN(1)*COS(P(3))
         DO 1060 I=2,NPEND
            L2 = 2*(I - 1) + 1
            L3 = 3*(I-1) + 1
            G(L2) = P(L3-3) + PLEN(I-1)*SIN(P(L3-1))
     $           -P(L3) + PLEN(I)*SIN(P(L3+2))
            G(L2 + 1) = P(L3-2) - PLEN(I-1)*COS(P(L3-1))
     $           -P(L3+1) - PLEN(I)*COS(P(L3+2))
 1060    CONTINUE
      ENDIF
      IF (QFLAG(8)) THEN
         DO 1070 I=1,NG
            GI(I) = 0.D0
 1070    CONTINUE  
      ENDIF
      IF (QFLAG(9)) THEN
         NZFACT = 0
         IPCF = 0
         WRITE (LUSOUT,*) 'PLSF (FL not available)'
      ENDIF      
C
      RETURN
      END
C
      SUBROUTINE PLSIN (NP, NV, NL, NG, NU, NBLK, NMRC, NZGMAX,
     $     NZFMAX, NSWMAX, T0, P, V, U, A, RLAM, PROTXT)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      INTEGER NP, NV, NL, NG, NU, NBLK, NMRC, NZGMAX, NZFMAX
      DIMENSION P(*), V(*), U(*), A(*), RLAM(*)
      CHARACTER*(*) PROTXT
C
C **********************************************************************
C
C   Pendulum (variable number of bodies < 31)
C
C
C   Parameter:  masses                  pmass(i) 
C               lengths                 plen(i)  
C
C
C **********************************************************************
C
C* File              plsin.f
C* Version           1.1
C* Latest Change     96/07/19 (1.4)
C* Modification history
C    1.1             Correct initial values of RLAM.
C                    New free format starts text output with ' '.
C
C -------------------------------------------------------
C
      PARAMETER (NBMAX = 30)
      COMMON / PEND / PMASS(NBMAX), PLEN(NBMAX), NPEND
C
      SAVE
C
      CHARACTER*67 CMODUL
      DATA CMODUL /'@(#) plsin.f 1.4 from 96/07/19 (mexx 1.2)
     $     '/
C
      PROTXT = 'Pendulum A (sparse)'
C
C  Machine constants (logical units stdin and stdout)
      LUSIN = 5
      LUSOUT = 6
C
      WRITE (LUSOUT,*) 'Start example ', PROTXT
      WRITE (LUSOUT,*) ' '
C
C  Dimensions
C  Get problem parameters
C
      WRITE (LUSOUT,*) 'Enter number of bodies (default 8):'
C      READ (LUSIN,*, ERR=1000, END=1000) NPEND
C      GOTO 1010
 1000 NPEND = 8
 1010 WRITE (LUSOUT,*) 'Given:', NPEND
      IF (NPEND .LT. 1 .OR. NPEND .GT. NBMAX) THEN
         WRITE (LUSOUT,*)
     $        'Number of bodies is out of range (default value chosen)'
         NPEND = 8
      ENDIF
      WRITE (LUSOUT,*) ' '
C
      NP = 3*NPEND
      NV = NP
      NG = 2*NPEND
      NL = NG
      NU = 0
      NMRC = 1
      NBLK = NP
      NZGMAX = 8*NPEND - 4
      NZFMAX = 0
      NSWMAX = NPEND - 1
C
C  Model parameters
      T0 = 0.D0
      PMASSH = 1.0D0/DBLE(NPEND)
      PLENH = 0.5D0/DBLE(NPEND)
      DO 1020 I=1,NPEND
         PMASS(I) = PMASSH
         PLEN(I) = PLENH  
 1020 CONTINUE
C     
C  Initial values
      PIHALF = ASIN(1.D0)
      P(1) = PLEN(1)
      P(2) = 0.D0 
      P(3) = PIHALF  
C     
      V(1) = 0.D0
      V(2) = 0.D0
      V(3) = 0.D0
C     
      RLAM(1) = 0.D0
      RLAM(2) = 0.D0
C     
      DO 1030 I=2,NPEND
         L = (I-1)*3 + 1
         P(L) = P(L-3) + PLEN(I-1) + PLEN(I)
         P(L+1) = P(2) 
         P(L+2) = P(3)
         V(L) = 0.D0
         V(L+1) = 0.D0
         V(L+2) = 0.D0
         L = (I-1)*2 + 1
         RLAM(L) = 0.D0
         RLAM(L+1) = 0.D0
 1030 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE PLDOUT (NDP, NDV, NDU, NDA, NDL, T,
     $                   DP, DV, DU, DA, DL, INDP,
     $                   INDV, INDU, INDA, INDL, INFO, IFAIL)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      INTEGER NDP, NDV, NDU, NDA, NDL
      DOUBLE PRECISION T, DP(*), DV(*), DU(*), DA(*), DL(*)
      INTEGER INDP(*), INDV(*), INDU(*), INDA(*), INDL(*), INFO(50),
     $     IFAIL 
C
C
C* File              pldout.f
C* Version           1.1
C* Latest Change     96/04/30 (1.5)
C* Modification history
C    1.1             Avoid free format because of compiler dependency.
C      
C* Input arguments:
C
C   NDP       Int   Number of dense output components of position 
C                   vector P. Dimension of arrays DP and INDP. The
C                   current value is either the user prescibed value
C                   IWK(51) (if MDIND=INFO(33)=1 was set) or NP, the
C                   dimension of the position vector P (set by MEXX if
C                   INFO(33)=0 was set).
C
C   NDV       Int   as NDP but now for the velocity vector V.
C
C   NDU       Int   as NDP but now for the vector of dynamic 
C                   variables U.
C
C   NDA       Int   as NDP but now for the acceleration vector A.
C
C   NDL       Int   as NDP but now for the vector of Lagrange 
C                   multipliers.
C
C   T         Dble  Current time
C
C   DP        Dble  the dense output position vector. 
C
C   DV        Dble  the dense output velocity vector. 
C
C   DU        Dble  the dense output dynamic variable vector. 
C
C   DA        Dble  the dense output acceleration vector. 
C
C   DL        Dble  the dense output Lagrange multiplier vector. 
C
C   INDP(NDP) Int   indicator array of length NP. List of components of
C                   the position vecor p which was selected for dense
C                   output. INDP(I)=J means that component J of the
C                   vector P is available as component I of vector DP.
C                   (I.e. DP(I)=P(INDP(I))
C                   Note, that MEXX automatically provides INDP(I)=I
C                   if INFO(33)=0 was set.
C
C   INDV(NDV) Int   as INDP but now for the velocity vector V.
C
C   INDU(NDU) Int   as INDP but now the vector of dynamic variables U.
C
C   INDU(NDA) Int   as INDP but now for the acceleration vector A.
C
C   INDL(NDL) Int   as INDP but now the vector of Lagrange multipliers.
C  
C   INFO(50)  Int   Info array
C                   INFO(1) .. INFO(40) see MXJOB(1) .. MXJOB(40) 
C                                           in MEXX
C                   INFO(41): Type of call: 
C                             = 0: the first call
C                             = 1: an intermediate call
C                             = 2: the last call at TFIN
C                             = 3: at a root of switching function
C                                  and integration continues
C                             = 4: at a root of switching function
C                                  and integration will terminate
C                                  (this is the last call)
C
C                   INFO(42): number of current integration step
C
C                   INFO(43): interpolation indicator
C                             = 0: current t is an integration point
C                             = 1: the solution values are interpolated
C                                  values
C
C                   INFO(44): call counter
C
C   IFAIL     Int   Error indicator. IFAIL=0 on input
C   
C   Do NOT alter any of the input arguments.
C
C
C* Output arguments:
C
C   IFAIL     Int  Error indicator
C                  If IFAIL.LT.0  on return, MEXX will terminate
C 
C----------------------------------------------------------------------
C
C
      INTEGER LUNIT
C
      CHARACTER*11 FNAME
C
      PARAMETER ( LUNIT = 11 , FNAME = 'PLSADR.DNO' )
C
      SAVE
C
      IFAIL = 0
C
      IF(INFO(41).EQ.0) THEN
         OPEN (UNIT=LUNIT,FILE=FNAME)
         WRITE(LUNIT,'(A)') 
     $    'C Dense Output Data from MEXX by routine pldout'
         WRITE(LUNIT,'(A)') 
     $    'C  the user selected dense output options are:'
         WRITE(LUNIT,'(A)') 'C    MDOUT     NDOUT     MDIND     MMODE'
         WRITE(LUNIT,'(3I10)') INFO(31),INFO(32),INFO(33),INFO(3)
         WRITE(LUNIT,'(A)')  'C  the dense output dimensions are:'
         WRITE(LUNIT,'(A)') 'C NDP  NDV  NDU  NDA  NDL'
         WRITE(LUNIT,'(5I5)') NDP,NDV,NDU,NDA,NDL 
         WRITE(LUNIT,'(A)') 'C  the selected components of P are:'
         WRITE(LUNIT,'(10I7)') (INDP(I),I=1,NDP)   
         WRITE(LUNIT,'(A)') 'C  the selected components of V are:'
         WRITE(LUNIT,'(10I7)') (INDV(I),I=1,NDV)   
         WRITE(LUNIT,'(A)') 'C  the selected components of U are:'
         WRITE(LUNIT,'(10I7)') (INDU(I),I=1,NDU)   
         WRITE(LUNIT,'(A)') 'C  the selected components of A are:'
         WRITE(LUNIT,'(10I7)') (INDA(I),I=1,NDV)   
         WRITE(LUNIT,'(A)') 'C  the selected components of LAMBDA are:'
         WRITE(LUNIT,'(10I7)') (INDL(I),I=1,NDL)   
      ENDIF
C
      WRITE(LUNIT,'(A)') 'C  No                   T Type Ipol'
      WRITE(LUNIT,'(I5,1PD20.10,I5,I5)') INFO(44),T,INFO(41),INFO(43)
      WRITE(LUNIT,'(A)') 'C  DP:'
      WRITE(LUNIT,'(4(1PD20.10))') (DP(I),I=1,NDP)   
      WRITE(LUNIT,'(A)') 'C  DV:'
      WRITE(LUNIT,'(4(1PD20.10))') (DV(I),I=1,NDV)   
      WRITE(LUNIT,'(A)') 'C  DU:'
      WRITE(LUNIT,'(4(1PD20.10))') (DU(I),I=1,NDU)   
      WRITE(LUNIT,'(A)') 'C  DA:'
      WRITE(LUNIT,'(4(1PD20.10))') (DA(I),I=1,NDA)   
      WRITE(LUNIT,'(A)') 'C  DL:'
      WRITE(LUNIT,'(4(1PD20.10))') (DL(I),I=1,NDL)   
C
      IF(INFO(41).EQ.2) THEN
         WRITE(LUNIT,'(A)') 
     $    'C MEXX terminates at final point of integration'
      ENDIF
      IF(INFO(41).EQ.4) THEN
         WRITE(LUNIT,'(A)') 
     $    'C MEXX terminates at a root of the switching functions'
      ENDIF
      IF(INFO(41).EQ.2 .OR. INFO(41).EQ.4) THEN
         CLOSE (UNIT=LUNIT)
      ENDIF
C
C
C
C
      RETURN 
C
C end of subroutine PLDOUT
C
      END
C
      SUBROUTINE PLSWI (NP, NV, NL, NU, T, P, V, U, A, RLAM, NSWIF, G)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION P(NP), V(NV), U(*), A(NV), RLAM(NL), G(*)
C
C* File              plswi.f
C* Version           1.0
C* Latest Change     93/09/15 (1.1)
C      
      DIMENSION DELPHI(101)
C
      SAVE
      PI=2.D0*ASIN(1.D0)
C
      NSWIF=0
      DO 1000 I=6,NP,3
         DELPHI(I)=P(I)-P(I-3)
         NSWIF=NSWIF+1
         G(NSWIF)=DELPHI(I)**2-PI**2
 1000 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE DUOUT (NP, NV, NU, NL, T, P, V, U, A, RL, INFO, IRTRN)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      INTEGER INFO(50), NP, NV, NU
      INTEGER IRTRN
      DOUBLE PRECISION T
      DIMENSION P(NP), V(NV), U(NU), A(NV), RL(NL)
C
C* Begin Prologue DUOUT
C
C* Purpose           Dummy routine instead of solution output for MEXX
C
C* File              duout.f
C* Version           1.0
C* Latest Change     96/03/05 (1.2)
C
C  
C* Input arguments:
C
C   NP        Int   Length of the position vector P
C
C   NV        Int   Length of the velocity vector V
C
C   NU        Int   Length of the vector of dynamic variables U
C
C   NA        Int   Length of the acceleration vector A
C
C   NL        Int   Length of the vector of Lagrange multipliers
C
C   T         Dble  Current time
C
C   P(NP)     Dble  The position vector. 
C
C   V(NV)     Dble  The velocity vector. 
C
C   U(NU)     Dble  The dynamic variable vector. 
C
C   A(NV      Dble  the dense output acceleration vector. 
C
C   RL(NL)    Dble  the dense output Lagrange multiplier vector. 
C  
C   INFO(50)  Int   Info array
C                   INFO(1) .. INFO(40) see MXJOB(1) .. MXJOB(40) 
C                                           in MEXX
C                   INFO(41): Type of call: 
C                             = 0: the first call
C                             = 1: an intermediate call
C                             = 2: the last call at TFIN
C                             = 3: at a root of switching function
C                                  and integration continues
C                             = 4: at a root of switching function
C                                  and integration will terminate
C                                  (this is the last call)
C
C                   INFO(42): number of current integration step
C
C                   INFO(43): interpolation indicator
C                             = 0: current t is an integration point
C                             = 1: the solution values are interpolated
C                                  values
C
C                   INFO(44): call counter
C   
C   Do NOT alter any of the input arguments
C
C* Output argument:
C
C   IRTRN     Int  Stop indicator
C                  If IRTRN .NE. 0 on return, MEXX will terminate
C 
C
      LUSOUT = 6
C
      WRITE (LUSOUT,*) ' +++ Warning in DUOUT'
      WRITE (LUSOUT,*)
     $     '     This dummy routine should never have been called.'
C
      IRTRN = 0
C
      RETURN
      END
