      FUNCTION SQUANK ( A, BIG, ERROR, FIFTH, RUM, NO, FUN )

c*********************************************************************72
c
C  SQUANK stands for Simpson Quadrature Used Adaptively, Noise Killed.
C
c  Discussion:
c
c    A line which read "FX5 = FUN ( X3 )" has been corrected 
c    to "FX5 = FUN ( X5 )"
c    Thanks to Ken Ebisawa, 02 January 2019.
c
c  Modified:
c
c    02 January 2019
c
C  CALLING PROGRAM REQUIRES:
C
C    A      THE LOWER LIMIT OF INTEGRATION;
C    BIG    THE UPPER LIMIT;
C    ERROR  THE REQUIRED TOLERANCE (ABSOLUTE ERROR);
C    FUN    THE EXTERNAL FUNCTION TO BE INTEGRATED.
C
C  OUTPUT:
C
C    SQUANK THE FIFTH ORDER RESULT = THIRD + FIFTH;
C    FIFTH  THE FIFTH ORDER ADJUSTMENT TERM;
C    RUM    THE CLAIMED TOLERANCE (ADJUSTED FOR ROUNDOFF ERROR).
C    NO     THE NUMBER OF FUNCTION EVALUATIONS REQUIRED.
C
C  NOTES ON USE.  
c
c    (1) DISCONTINUOUS FUNCTIONS.
C
C    THIS ROUTINE IS BASED ON DEGREE 3 AND DEGREE 5 LOCAL POLYNOMIAL APPROX-
C    IMATION.  CONSEQUENTLY, IT SHOULD NOT BE USED WITH FUNCTIONS WHICH HAVE
C    DISCONTINUITIES IN THE FOURTH OR LOWER DERIVATIVES WITHIN THE INTERVAL OF
C    INTEGRATION.  IF THERE ARE SUCH DISCONTINUITIES, THIS ROUTINE WILL TAKE
C    THIS TO BE EVIDENCE OF ROUND OFF ERROR IN FUNCTION VALUES, AND WILL ADJUST
C    THE TOLERANCE.
C
C    IF THE LOCATIONS OF SUCH DISCONTINUITIES ARE KNOWN, THE ROUTINE MAY
C    BE USED SEPARATELY FOR EACH INTERVAL BETWEEN CONSECUTIVE DISCONTINUITIES.
C    WHILE, LIKE ALL SUCH ROUTINES, IT DISLIKES DISCONTINUITIES, IT CAN HANDLE
C    THEM IF THEY ARE LOCATED AT THE END POINTS OF THE INTEGRATION INTERVAL.
C
C    (2) FUNCTIONS WITH HIGH-FREQUENCY OSCILLATIONS.
C
C    THE ROUTINE WILL RETURN UNRELIABLE RESULTS FOR FUNCTIONS LIKE G(X)
C    TIMES COS(100*X).  IF THE HIGHEST PERIOD LIKELY TO BE ENCOUNTERED IS
C    KNOWN, THE INTERVAL SHOULD BE SUBDIVIDED IN SUCH A WAY THAT,
C      (ONE) THERE ARE NOT MORE THAN THREE PERIODS PER INTERVAL, AND
C      (TWO), THE PERIOD *P* DIVIDED BY THE SUBINTERVAL, *(B-A)*, IS NOT A
C    SIMPLE FRACTION *N/M* WITH N OR M LESS THAN 9.
C
C    (3)  INTERVAL SUB-DIVISION.
C
C    THE FAILURES DESCRIBED ABOVE ARE GENERALLY WORSE FOR *SQUANK* THAN FOR
C    OTHER ROUTINES BECAUSE *SQUANK* TAKES THE INCONVENIENT BEHAVIOR AS AN
C    INDICATION OF ROUND OFF ERROR.  IN GENERAL, SUB-DIVISION OF THE INTERVAL IS
C    ADVOCATED.  ESSENTIALLY, THE USER CARRIES OUT, UNDER DRIVING PROGRAM
C    CONTROL, A SEQUENCE OF CALCULATIONS WHICH SHOULD HAVE BEEN CARRIED OUT
C    IN THE SUBROUTINE IN ANY CASE.  IN THIS WAY, HE PREVENTS CHANCE LOW ORDER
C    FALSE CONVERGENCE AT VIRTUALLY NO ADDITIONAL COST.  NOTE THAT THE SUM OF
C    THE PARAMETERS *ERROR* FOR THE SUB-INTERVALS SHOULD CORRESPOND TO THE
C    VALUE REQUIRED FOR THE WHOLE INTERVAL.
C
C  NIM NUMBERING SYSTEM AND LOGIC.
C
C    THE INTERVAL (A,B) IS DEFINED NIM = 1, LEVEL = 0.
C    THE INTERVAL NIM = N, LEVEL = L, IS BISECTED, IF NECESSARY, INTO
C    TWO INTERVALS, NIM = 2*N AND NIM = 2*N+1, BOTH AT LEVEL = L+1.
C    IF INTERVAL NIM = N, LEVEL = L DOES NOT CONVERGE, THE NEXT INTERVAL
C    CONSIDERED IS NIM = 2*N, LEVEL = L+1.
C    IF INTERVAL NIM = N, LEVEL = L DOES CONVERGE, THE NEXT INTERVAL CONSIDERED
C    IS NIM = M(R)+1, LEVEL = L-R, WHERE M(R) IS THE FIRST
C    EVEN MEMBER OF THE SEQUENCE M(Q) = N, M(S+1) = (M(S)-1)/2.  IF THIS
C    GIVES LEVEL = 0, THE CALCULATION IS COMPLETE.
C
C  SCALING TO AVOID EXCESSIVE DIVISION BY TWO.
C
C    THE INTERVAL (X1,X5) IS OF LENGTH H = X5-X1.  THE POINTS X1, X2, X3, X4, X5
C    ARE THE POINTS OF QUARTERSECTION OF THIS INTERVAL AND FX1, FX2, FX3, FX4, 
C    FX5 ARE THE CORRESPONDING FUNCTION VALUES.
C    EST IS APPROXIMATION TO (6.0/H)*INTEGRAL(X1,X5).
C    (EST1+EST2) IS APPROXIMATION TO (12.0/H)*INTEGRAL(X1,X5).
C    SUM IS APPROXIMATION TO (12.0)*INTEGRAL(A,X1).
C
C  STORAGE.
C
C    X3ST(L) = 0.5*(X5ST(L) + X1).  THUS, X3ST(L) COULD BE RECALCULATED
C    AT EACH STAGE TO AVOID STORAGE.  ESTST(L) IS SAME IN THIS RESPECT.
C    THE RESULTS OF ABOVE RECALCULATIONS ARE IDENTICAL MACHINE NUMBERS.
C
C    X5ST(L) = X1+(B-A)*(2**(-L)).  THIS COULD ALSO BE RECALCULATED.  BUT
C    IN THIS CASE, CALCULATION IS EXCESSIVE AND THERE IS A POSSIBILITY OF
C    ROUND OFF ERROR ARISING BECAUSE THE SAME POINT IS BEING CALCULATED
C    IN TWO OR MORE DIFFERENT WAYS.
C
C  AVOIDANCE OF ROUND OFF ERROR TROUBLE.
C
C    IF INTERVAL DOES NOT CONVERGE, FOLLOWING INTERVAL SHOULD HAVE ADIFF
C    VALUE APPROXIMATELY EQUAL TO (1/16) TIMES PREVIOUS ADIFF VALUE, CALLED
C    ADIFF1 IN THE CODE.  THERE IS A THEOREM WHICH STATES THAT, UNLESS THE
C    FOURTH DERIVATIVE OF FUN(X) VANISHES IN THE PREVIOUS INTERVAL, ADIFF
C    IS LESS THAN OR EQUAL TO ADIFF1.  IF THIS DOES NOT HAPPEN, IT IS TAKEN
C    TO BE AN INDICATION OF POSSIBLE ROUND OFF LEVEL.  IN THIS CASE,
C    UNLESS LEV IS LESS THAN FIVE, THE CURRENT TOLERANCE LEVEL, CEPS,
C    IS APPROPRIATELY ADJUSTED.  HOWEVER CEPS IS RESET AS AND WHEN
C    APPEARS THAT IT SHOULD BE ADJUSTED EITHER UP OR DOWN.  IT IS REDUCED
C    IF CONVERGENCE OCCURS WITH A NONZERO ADIFF STRICTLY LESS THAN
C    0.25*CEPS.  AN INVOLVED SECTION OF CODING GUARDS TO SOME EXTENT AGAINST
C    AN UNDOUBLE PRECISIONISTIC VALUE ARISING AS A RESULT OF A ZERO IN THE FOURTH
C    DERIVATIVE.  A FACTOR EFACT IS CALCULATED WHICH ADJUSTS THE CLAIMED
C    TOLERANCE TO TAKE INTO ACCOUNT THESE ALTERATIONS IN THE TOLERANCE LEVEL.
C
C    THE ROUTINE ENTERS THESE INVOLVED SECTIONS OF CODING ONLY IF ROUND
C    OFF ERROR APPEARS TO BE PRESENT.  IN A NORMAL (ROUND OFF ERROR FREE)
C    RUN, THESE SECTIONS ARE SKIPPED AT A COST OF A SINGLE COMPARISON PER
C    ITERATION (TWO FUNCTION EVALUATIONS).
C
C  ARBITRARY CONSTANTS.
C
C    THE FOLLOWING CONSTANTS HAVE BEEN ASSIGNED IN THE LIGHT OF EXPERIENCE
C    WITH NO THEORETICAL JUSTIFICATION.
C
C    (1) NO CONVERGENCE IS ALLOWED AT LEVEL **ZERO**.  THIS MEANS THAT THE
C    ROUTINE IS CONSTRAINED TO BASE THE RESULT ON AT LEAST 9 FUNCTION VALUES.
C
C    (2) NO UPWARD ADJUSTMENT OF THE TOLERANCE LEVEL IS CONSIDERED AT LEVELS
C    LOWER THAN LEVEL **FIVE**.  THE OINT SPACING IS THEN (BIG-A)/128.0.
C
C    (3) PHYSICAL LIMIT.  HIGHEST LEVEL ALLOWED IS LEVEL = **THIRTY**.  HERE
C    CONVERGENCE IS ASSIGNED WHETHER OR NOT THE INTERVAL HAS CONVERGED.  THE
C    POINT SPACING IS THEN ABOUT (BIG-A)*2.0*10**(-10).
C
C    (4) UPWARD ADJUSTMENT OF TOLERANCE LEVEL IS LIMITED IN GENERAL TO
C    A FACTOR OF **2.0** OR LESS.
C
C    (5) DOWNWARD ADJUSTMENT OF TOLERANCE LEVEL IS INHIBITED IN GENERAL UNLESS
C    BY A FACTOR GREATER THAN **4.0**.
C
C  SOME NOTATION.
C
C    SUM AND SIM ARE RUNNING SUMS.  INCREASED AT STAGE EIGHT, THEY ARE
C    RESPECTIVELY 12.0 * ( THIRD ORDER APPROXIMATION TO THE INTEGRAL )
C    AND -180.0 * ( FIFTH ORDER ADJUSTMENT TO THE INTEGRAL).
C
C    CEPSF IS THE REQUIRED (SCALED) TOLERANCE.
C
C    CEPS IS THE RUNNING VALUE OF THE ADJUSTED TOLERANCE.
C
C    QCEPS = 0.25 * CEPS.
C
C    LEVTAG = -1 OR 0, 2, 3, INDICATES WHETHER TOLERANCE IS NOT OR IS
C    CURRENTLY ADJUSTED.  (SEE COMMENT IN STAGE SEVEN).
C
C    EFACT IS RUNNING SUM CORRESPONDING TO 180.0 * RUN.
C
C    FACERR = 15.0 OR 1.0, DEPENDING ON WHETHER TOLERANCE IS OR IS NOT
C    CURRENTLY ADJUSTED.  IF IT IS, THERE IS NO JUSTIFICATION FOR THE
C    FIFTH ORDER ADJUSTMENT AND ACCURACY IS NOT EXPECTED TO BE (1/15) TIMES
C    DIFFERENCE OF APPROXIMATIONS.  FACERR = 15.0 REMOVES THE BUILT-IN 15.0
C    FACTOR FOR CALCULATION OF EFACT.
C
C    EPMACH, THE MACHINE ACCURACY PARAMETER.  THE ROUND OFF ERROR GUARD DOES
C    NOT REQUIRE THIS NUMBER.  IT IS MACHINE INDEPENDENT.  THIS IS ONLY
C    USED TO HELP IN AN INITIAL GUESS IN STAGE TWO IF THE VALUE OF ERROR
C    HAPPENS TO BE ZERO.  ANY NON-ZERO NUMBER MAY BE USED INSTEAD, WITH A
C    VERY SMALL PENALTY IN NUMBER OF FUNCTION EVALUATIONS IF A COMPLETELY
C    UNREASONABLE NUMBER IS USED.
C
      IMPLICIT NONE

      DOUBLE PRECISION A
      DOUBLE PRECISION ADIFF
      DOUBLE PRECISION ADIFF1
      DOUBLE PRECISION BIG
      DOUBLE PRECISION CEPS
      DOUBLE PRECISION CEPSF
      DOUBLE PRECISION CEPST
      DOUBLE PRECISION CRIT
      DOUBLE PRECISION DIFF
      DOUBLE PRECISION EFACT
      DOUBLE PRECISION EPMACH
      DOUBLE PRECISION ERROR
      DOUBLE PRECISION EST
      DOUBLE PRECISION EST1
      DOUBLE PRECISION EST2
      DOUBLE PRECISION ESTST(30)
      DOUBLE PRECISION FACERR
      DOUBLE PRECISION FIFTH
      DOUBLE PRECISION FUN
      DOUBLE PRECISION FX1
      DOUBLE PRECISION FX2
      DOUBLE PRECISION FX3
      DOUBLE PRECISION FX3ST(30)
      DOUBLE PRECISION FX4
      DOUBLE PRECISION FX5
      DOUBLE PRECISION FX5ST(30)
      INTEGER LEV
      INTEGER LEVTAG
      INTEGER NIM
      INTEGER NO
      INTEGER NOM
      INTEGER NUM
      DOUBLE PRECISION PREDIF(30)
      DOUBLE PRECISION QCEPS
      DOUBLE PRECISION RUM
      DOUBLE PRECISION SIM
      DOUBLE PRECISION SQUANK
      DOUBLE PRECISION SUM
      DOUBLE PRECISION THIRD
      DOUBLE PRECISION X1
      DOUBLE PRECISION X2
      DOUBLE PRECISION X3
      DOUBLE PRECISION X3ST(30)
      DOUBLE PRECISION X4
      DOUBLE PRECISION X5
      DOUBLE PRECISION X5ST(30)
      DOUBLE PRECISION XZERO

      EPMACH = 0.0000000000075D+00
C
C  STAGE ONE
C  INITIALIZE ALL QUANTITIES REQUIRED FOR CENTRAL CALCULATION (STAGE 3).
C
      SUM = 0.0D+00
      SIM = 0.0D+00
      CEPSF = 180.0D+00 * ERROR / ( BIG - A )
      CEPS = CEPSF
      ADIFF = 0.0D+00
      LEVTAG = -1
      FACERR = 1.0D+00
      XZERO = A
      EFACT = 0.0D+00
      NIM = 1
      LEV = 0
C
C  FIRST INTERVAL.
C
      X1 = A
      X5 = BIG
      X3 = 0.5D+00 * ( A + BIG )
      FX1 = FUN ( X1 )
      FX3 = FUN ( X3 )
c
c  This line has been corrected from "FX5 = FUN ( X3 )"
c  Thanks to Ken Ebisawa, 02 January 2019.
c
      FX5 = FUN ( X5 )
      NO = 3
      EST = FX1 + FX5 + 4.0D+00 * FX3
C
C  STAGE TWO
C  SET A STARTING VALUE FOR TOLERANCE IN CASE THAT CEPSF = 0.0.
C
      IF ( CEPSF ) 295, 205, 295

205   continue

      LEVTAG = 0
      FACERR = 15.0D+00
      CEPS = EPMACH * ABS ( FX1 )
      IF ( FX1 ) 295, 210, 295

210   continue

      CEPS = EPMACH * ABS ( FX3 )
      LEVTAG = 3
      IF ( FX3 ) 295, 215, 295

215   continue

      CEPS = EPMACH * ABS ( FX5 )
      IF ( FX5 ) 295, 220, 295

220   continue

      CEPS = EPMACH

295   continue

      QCEPS = 0.25D+00 * CEPS
C
C  INITIALIZING COMPLETE.
C
C  STAGE THREE
C  CENTRAL CALCULATION.
C  REQUIRED X1, X3, X5, FX1, FX3, FX5, EST, ADIFF.
C
300   CONTINUE

      X2 = 0.5D+00 * ( X1 + X3 )
      X4 = 0.5D+00 * ( X3 + X5 )
      FX2 = FUN ( X2 )
      FX4 = FUN ( X4 )
      NO = NO + 2
      EST1 = FX1 + 4.0D+00 * FX2 + FX3
      EST2 = FX3 + 4.0D+00 * FX4 + FX5
      ADIFF1 = ADIFF
      DIFF = EST + EST - EST1 - EST2
      IF ( LEV - 30 ) 305, 800, 800

305   continue

      ADIFF = ABS ( DIFF )
      CRIT = ADIFF - CEPS
      IF ( CRIT ) 700, 700, 400
C
C  END OF CENTRAL LOOP.
C  NEXT STAGE IS STAGE FOUR IN CASE OF NO NATURAL CONVERGENCE.
C  NEXT STAGE IS STAGE SEVEN IN CASE OF NATURAL CONVERGENCE.
C
C  STAGE FOUR.
C  NO NATURAL CONVERGENCE.  A COMPLEXT SEQUENCE OF INSTRUCTIONS
C  FOLLOWS WHICH ASSIGNS CONVERGENCE AND/OR ALTERS TOLERANCE
C  LEVEL IN UPWARD DIRECTION IF THERE ARE INDICATIONS OF ROUND OFF
C  ERROR.
C
400   CONTINUE

      IF ( ADIFF1 - ADIFF ) 410, 410, 500
C
C  IN A NORMAL RUN WITH NO ROUND OFF ERROR PROBLEM, ADIFF1 IS GREATER THAN
C  ADIFF, AND THE REST OF STAGE FOUR IS OMITTED.
C
410   continue

      IF ( LEV - 5 ) 500, 415, 415

415   continue

      EFACT = EFACT + CEPS * ( X1 - XZERO ) * FACERR
      XZERO = X1
      FACERR = 15.0D+00
C
C  THE REST OF STAGE FOUR DEALS WITH UPWARD ADJUSTMENT OF TOLERANCE (CEPS)
C  BECAUSE OF SUSPECTED ROUND OFF ERROR TROUBLE.
C
      IF ( ADIFF - 2.0D+00 * CEPS ) 420, 420, 425
C
C  SMALL JUMP IN CEPS, ASSIGN CONVERGENCE.
C
420   continue

      CEPS = ADIFF
      LEVTAG = 0
      GO TO 780

425   continue

      IF ( ADIFF1 - ADIFF ) 435, 430, 435
C
C  LARGE JUMP IN CEPS.
C
430   continue

      CEPS = ADIFF
      GO TO 445
C
C  FACTOR TWO JUMP IN CEPS.
C
435   continue

      CEPS = 2.0D+00 * CEPS
      IF ( LEVTAG - 3 ) 440, 445, 445

440   continue

      LEVTAG = 2

445   continue

      QCEPS = 0.25D+00 * CEPS
C
C  STAGE FIVE
C  NO ACTUAL CONVERGENCE.  STORE RIGHT HAND ELEMENTS.
C
500   CONTINUE

      NIM = 2 * NIM
      LEV = LEV + 1
      ESTST(LEV) = EST2
      X3ST(LEV) = X4
      X5ST(LEV) = X5
      FX3ST(LEV) = FX4
      FX5ST(LEV) = FX5
      PREDIF(LEV) = ADIFF
C
C  STAGE SIX
C  SET UP QUANTITIES FOR CENTRAL CALCULATION.
C
C  READY TO GO AHEAD AT LEVEL LOWER WITH LEFT HAND ELEMENTS.
C  X1 AND FX1 ARE THE SAME AS BEFORE.
C
      X5 = X3
      X3 = X2
      FX5 = FX3
      FX3 = FX2
      EST = EST1
      GO TO 300
C
C  STAGE SEVEN.
C  NATURAL CONVERGENCE IN PREVIOUS INTERVAL.  THE FOLLOWING COMPLEX SEQUENCE
C  CHECKS PRIMARILY THAT TOLERANCE LEVEL IS NOT TOO HIGH.  UNDER CERTAIN
C  CIRCUMSTANCES, NONCONVERGENCE IS ASSIGNED AND/OR TOLERANCE LEVEL
C  IS RESET.
C
700   CONTINUE
C
C  CHECK THAT IT WAS NOT LEVEL ZERO INTERVAL.  IF SO, ASSIGN NONCONVERGENCE.
C
      IF ( LEV ) 400, 400, 705
C
C  LEVTAG = -1, CEPS = CEPSF, ITS ORIGINAL VALUE.
C  LEVTAG =  0, CEPS IS GREATER THAN CEPSF, REGULAR SITUATION.
C  LEVTAG =  2, CEPS IS GREATER THAN CEPSF.  CEPS PREVIOUSLY ASKED FOR A BIG
C               JUMP, BUT DID NOT GET ONE.
C  LEVTAG =  3, CEPS IS GREATER THAN CEPSF.  CEPS PREVIOUSLY HAD A BIG JUMP.
C
705   continue

      IF ( LEVTAG ) 800, 710, 710
C
C  IN A NORMAL RUN WITH NO ROUND OFF ERROR PROBLEM, LEVTAG = -1, AND THE
C  REST OF STAGE SEVEN IS OMITTED.
C
710   continue

      CEPST = 15.0D+00 * CEPS
C
C  CEPST HERE IS FACERR * CURRENT VALUE OF CEPS.
C
      IF ( CRIT ) 715, 800, 800

715   continue

      IF ( LEVTAG - 2 ) 720, 740, 750
C
C  LEVTAG = 0.
C
720   continue

      IF ( ADIFF ) 800, 800, 725

725   continue

      IF ( ADIFF - QCEPS ) 730, 800, 800

730   continue

      IF ( ADIFF - CEPSF ) 770, 770, 735

735   continue

      LEVTAG = 0
      CEPS = ADIFF
      EFACT = EFACT + CEPST * ( X1 - XZERO )
      XZERO = X1
      GO TO 445
C
C  LEVTAG = 2.
C
740   continue

      LEVTAG = 0
      IF ( ADIFF ) 765, 765, 725
C
C  LEVTAG = 3.
C
750   continue

      LEVTAG = 0
      IF ( ADIFF ) 775, 775, 730

765   continue

      CEPS = ADIFF1
      GO TO 775

770   continue

      LEVTAG = -1
      FACERR = 1.0D+00
      CEPS = CEPSF

775   continue

      EFACT = EFACT + CEPST * ( X1 - XZERO )
      XZERO = X1

780   CONTINUE

      QCEPS = 0.25D+00 * CEPS
C
C  STAGE EIGHT.
C  ACTUAL CONVERGENCE IN PREVIOUS INTERVAL.  INCREMENTS ADDED INTO
C  RUNNING SUMS.
C
C  ADD INTO SUM AND SIM
C
800   CONTINUE

      SUM = SUM + ( EST1 + EST2 ) * ( X5 - X1 )
      IF ( LEVTAG ) 805, 810, 810
C
C  WE ADD INTO SIM ONLY IF WE ARE CLEAR OF ROUND OFF LEVEL.
C
805   continue

      SIM = SIM + DIFF * ( X5 - X1 )

810   CONTINUE
C
C  STAGE NINE.
C  SORT OUT WHICH LEVEL TO GO TO.  THIS INVOLVES NIM NUMBERING SYSTEM
C  DESCRIBED BEFORE STAGE ONE.
C
905   continue

      NUM = NIM / 2
      NOM = NIM - 2 * NUM
      IF ( NOM ) 910, 915, 910

910   continue

      NIM = NUM
      LEV = LEV - 1
      GO TO 905

915   continue

      NIM = NIM + 1
C
C  NEW LEVEL IS SET.  IF LEV = 0, WE HAVE FINISHED.
C
      IF ( LEV ) 1100, 1100, 1000
C
C  STAGE TEN
C  SET UP QUANTITIES FOR CENTRAL CALCULATION.
C
1000  CONTINUE

      X1 = X5
      FX1 = FX5
      X3 = X3ST(LEV)
      X5 = X5ST(LEV)
      FX3 = FX3ST(LEV)
      FX5 = FX5ST(LEV)
      EST = ESTST(LEV)
      ADIFF = PREDIF(LEV)
      GO TO 300
C
C  STAGE ELEVEN.
C  CALCULATION NOW COMPLETE.  FINALIZE.
C
1100  CONTINUE

      EFACT = EFACT + CEPS * ( BIG - XZERO ) * FACERR
      RUM = EFACT / 180.0D+00
      THIRD = SUM / 12.0D+00
      FIFTH = - SIM / 180.0D+00
      SQUANK = THIRD + FIFTH

      RETURN
      END
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end

