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