! This program demonstrates how to use the solution queue
! trimming option which should be used for problems with
! otherwise generae large solution queues. (It also
! illustrates how to use the interpolation option.)
! Refer to the header comments at the beginning of the
! main program.

MODULE define_DDEs

  IMPLICIT NONE
  
  INTEGER, PARAMETER :: NEQN=1,NLAGS=1

  ! Physical parameters assigned in the main program:
  DOUBLE PRECISION :: a, b, bp, ap, bm, am, kappa
  ! Pi:
  INTEGER, PARAMETER :: DP = KIND (1.d0)
  REAL(DP), PARAMETER :: Pi = 3.141592653589793238462643_DP

  ! Needed in subroutine TRIM_GET below:
  DOUBLE PRECISION, ALLOCATABLE :: USER_TINT(:), USER_YINT(:,:)
  INTEGER USER_NINT,USER_START
  
CONTAINS

  SUBROUTINE DDES(T,Y,Z,DY)

    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY
    DOUBLE PRECISION, DIMENSION(NEQN,NLAGS) :: Z
     
     INTENT(IN)  :: T,Y,Z
     INTENT(OUT) :: DY

     DOUBLE PRECISION :: Zp, Zm, A_of_Z

     Zp = bp/kappa/ap*(ap-1D0)
     Zm = -bm/kappa/am*(am-1D0)

     IF (Z(1,1)>Zp) THEN
          A_of_Z=bp+bp/ap*(DTANH(kappa*ap/bp*(Z(1,1)-Zp))-1D0)
     ELSEIF (Z(1,1)<Zm) THEN
          A_of_Z=-bm-bm/am*(DTANH(kappa*am/bm*(Zm-Z(1,1)))-1D0)
     ELSE
          A_of_Z=kappa*Z(1,1)
     ENDIF
     
     !DY(1) = 2D0*Z(1,1)/(1D0 + Z(1,1)**9.65D0) - Y(1)
     DY(1) = -a*A_of_Z + b*DCOS(2D0*Pi*T)
     !DY(1) = -a*Z(1,1)

    RETURN
  END SUBROUTINE DDES
!______________________________________________________________

  SUBROUTINE TRIM_GET

! ***Caution*** You should not make any changes to this subroutine.

! TRIM_GET gets the solution information saved by TRIM_SAVE in the
! files 'dropped_tqueue.dat' and 'dropped_queue.dat' with unit numbers
! UNUM655=655 and UNUM656=656, respectively. It then interpolates the
! values in USER_TINT(:). The interpolated solution is stored in the
! array USER_YINT(:,:).

! Usage:
! To use this subroutine, the interpolation arrays USER_TINT(1:USER_NINT)
! and USER_TINT(1:USER_NINT,1:NEQN) must be allocated in the main program
! where USER_NINT, also defined in the main program, is the number of
! desired interpolation points. The interpolation times
! USER_TINT(1:USER_NINT) must be increasing and spanned by the integration
! limits (TSPAN(1) and TSPAN(2) for the problem. The flag USER_START must
! be initialized to a value of 1 in the main program.
! Include the following in the header for your module:
!    DOUBLE PRECISION, ALLOCATABLE :: USER_TINT(:), USER_YINT(:,:)
!    INTEGER USER_NINT,USER_START
! Include the following in your main program:
!    USER_NINT = ...
!    ALLOCATE(USER_TINT(USER_NINT),USER_YINT(USER_NINT,NEQN),STAT=IER)
!    CALL MY_CHECK_STAT(IER,1)
!    USER_TINT(1:USER_INT) = ,,,
!    USER_START = 1

!    Get access to coefficient subroutine DDE_WSET.
     USE DDE_SOLVER_M

     IMPLICIT NONE

! .. Scalar Arguments ..
     INTEGER :: IDROPM1,IER,ITINT,ISTART,J,INEW,INDEXO,IBEGIN,MYMETHOD
     LOGICAL POLYS,POLYD,DEBUG
     DOUBLE PRECISION TQLAST,TVAL,TO,TN,DELINT,C
! ..
! .. Local Arrays ..
     DOUBLE PRECISION, ALLOCATABLE :: TQUEUE(:), QUEUE(:,:)
     DOUBLE PRECISION W(0:9),WD(0:9)
! ..
! .. Parameters (must match those used in DDE_SOLVER_M) ..
     INTEGER, PARAMETER :: UNUM655 = 655, UNUM656 = 656
! ..
!    OPEN the output files written by DDE_SOLVER_M.
     OPEN(UNIT=UNUM655,FILE='dropped_tqueue.dat',FORM="UNFORMATTED",&
          STATUS="OLD",POSITION="REWIND")
     OPEN(UNIT=UNUM656,FILE='dropped_queue.dat',FORM="UNFORMATTED",&
          STATUS="OLD",POSITION="REWIND")

!    Set parameters needed by the coefficient subroutine DDE_WSET.
     POLYS = .TRUE.
     POLYD = .FALSE.
     MYMETHOD = 4
     DEBUG = .FALSE.

!    Caution:
!    The READs must match the WRITEs in DDE_SOLVER_M.

!    Get the last time in processed previous queue.
     READ(UNUM655) TQLAST

!    Get the number of times in this queue:
     READ(UNUM655) IDROPM1
     IF (DEBUG) THEN
        PRINT *, ' In USER_TRIM_GET, IDROPM1 = ', IDROPM1
     END IF

!    Allocate, read and process the latest discarded queue data.
     ALLOCATE (TQUEUE(0:IDROPM1),QUEUE(NEQN,10*IDROPM1),STAT=IER)
     CALL MY_CHECK_STAT(IER,2)
     TQUEUE(0) = TQLAST
     READ(UNUM655) TQUEUE(1:IDROPM1)
     READ(UNUM656) QUEUE(1:NEQN,1:10*IDROPM1)

!    Perform the interpolation for the user times spanned by the
!    current solution queue.

!    Set the index of the start bracket check:
     ISTART = 0

     DO ITINT = USER_START, USER_NINT

        IF (ITINT < USER_NINT) THEN
           IF (USER_TINT(ITINT) >= USER_TINT(ITINT+1)) THEN
              PRINT *, ' USER_TINT is not increasing in TRIM_GET.'
              STOP
           END IF
        END IF

        TVAL = USER_TINT(ITINT)

!       Check that TVAL is spanned by the current solution queue.
        IF (TVAL < TQUEUE(0)) THEN
           PRINT *, ' Illegal interpolation requested in TRIM_GET.'
           PRINT *, ' tval/tqueue(0): ', tval,tqueue(0)
           STOP
        END IF

!       Stop interpolating when TVAL is beyond the present queue.
        IF (TVAL > TQUEUE(IDROPM1)) GOTO 20

!       Bracket TVAL in the present queue.
        INEW = IDROPM1

!       Linear search; remember last bracketing interval.

!       Caution:
!       TQUEUE is assumed to be increasing and spanned by the
!       integration limits supplied to DDE_SOLVER_M.

          IF (TVAL<=TQUEUE(ISTART)) THEN

            DO J = ISTART, 0, -1
              IF (TVAL>=TQUEUE(J)) THEN
!               TQUEUE(J) <= TVAL <= TQUEUE(J+1)
                ISTART = J
                INDEXO = J + 1
                TO = TQUEUE(J)
                TN = TQUEUE(INDEXO)
                GOTO 10
              END IF
            END DO

          ELSE

!           TVAL > TQUEUE(ISTART)
            DO J = ISTART + 1, INEW
              IF (TVAL<=TQUEUE(J)) THEN
!               TQUEUE(J-1) <= TVAL <= TQUEUE(J)
                ISTART = J - 1
                INDEXO = J
                TO = TQUEUE(ISTART)
                TN = TQUEUE(INDEXO)
                GOTO 10
              END IF
            END DO

          END IF

!       The bracketing search ends here.
10      CONTINUE

!       Start here in USER_TINT on the next pass.
        USER_START = ITINT + 1

!       The data to be interpolated corresponds to TO = TOLD, TN = TNEW.
!       INDEXO is the pointer for the corresponding slug of data. The
!       data for TOLD to be interpolated begins in queue location
!       (1,IBEGIN+1).
        IBEGIN = (INDEXO-1)*10

!       Define the step size corresponding to the interpolation data.
        DELINT = TN - TO

!       Determine the value of C at which the solution is to be
!       interpolated.
        IF (ABS(DELINT)<=0.0D0) THEN
          C = 0.0D0
        ELSE
          C = (TVAL-TO)/DELINT
        END IF

!       Calculate the interpolation coefficients.
!       (No derivative interpolation yet)
        CALL DDE_WSET(C,W,WD,POLYS,POLYD,MYMETHOD)

!       Interpolate the solution.
        USER_YINT(ITINT,1:NEQN) = QUEUE(1:NEQN,IBEGIN+1) + DELINT * &
          (W(0)*QUEUE(1:NEQN,IBEGIN+2)+W(2)*QUEUE(1:NEQN,IBEGIN+3)+ &
           W(3)*QUEUE(1:NEQN,IBEGIN+4)+W(4)*QUEUE(1:NEQN,IBEGIN+5)+ &
           W(5)*QUEUE(1:NEQN,IBEGIN+6)+W(6)*QUEUE(1:NEQN,IBEGIN+7)+ &
           W(7)*QUEUE(1:NEQN,IBEGIN+8)+W(8)*QUEUE(1:NEQN,IBEGIN+9)+ &
           W(9)*QUEUE(1:NEQN,IBEGIN+10))

     END DO

20   CONTINUE

!    Caution:
!    Do not exit this subroutine without doing the following two things.

!    CLOSE the output files written by DDE_SOLVER_M.
     CLOSE(UNUM655)
     CLOSE(UNUM656)

!    Get rid of the local allocatable arrays.
     DEALLOCATE(TQUEUE,QUEUE)
     CALL MY_CHECK_STAT(IER,3)

    RETURN
  END SUBROUTINE TRIM_GET
!______________________________________________________________

  SUBROUTINE MY_CHECK_STAT(IER,CALLED_FROM)

! .. Scalar Arguments ..
        INTEGER, INTENT (IN) :: CALLED_FROM, IER
! ..
  IF (IER/=0) THEN
     PRINT *, ' A storage allocation error occurred at code location ', &
     CALLED_FROM
     STOP
   END IF

   RETURN
   END SUBROUTINE MY_CHECK_STAT

END MODULE define_DDEs

!******************************************************************

PROGRAM ilya

!
! The DDE is defined in the module define_DDEs.  The problem
! is solved here with DDE_SOLVER and its output written to a
! file. An auxilary function imports the data into Matlab and
! plots it.

  USE define_DDEs
  USE DDE_SOLVER_M

  IMPLICIT NONE

  ! The quantities
  !
  !   NEQN = number of equations
  !   NLAGS = number of delays
  !
  ! are defined in the module define_DDEs as PARAMETERs so 
  ! they can be used for dimensioning arrays here. They are 
  ! passed to the solver in the array NVAR.

  INTEGER, DIMENSION(2) :: NVAR = (/NEQN,NLAGS/)

! Constant delay:
  !DOUBLE PRECISION, PARAMETER :: LAG=2D0
  !DOUBLE PRECISION, DIMENSION(NLAGS) :: DELAY=(/ LAG /)
  !DOUBLE PRECISION :: LAG=2D0
  DOUBLE PRECISION, DIMENSION(NLAGS) :: DELAY

! Constant history:
  DOUBLE PRECISION, DIMENSION(NEQN) :: HISTORY= (/ 1D0 /)
  TYPE(DDE_SOL) :: SOL
  TYPE(DDE_OPTS) :: OPTS
  TYPE(DDE_INT) :: Y,YLAG

! Output and interpolation:
  DOUBLE PRECISION :: T0, TFINAL, Tmin, Length, Sp, Sm, S, Var
  INTEGER, PARAMETER :: NOUT=10**5
  DOUBLE PRECISION, DIMENSION(NOUT) :: T,YYT
  DOUBLE PRECISION, ALLOCATABLE :: MY_SAVE(:,:)
  INTEGER IER,ITOL,NTOL,I,J,JMAX
  DOUBLE PRECISION TOL,OVERRUN,MAXDELAY
  LOGICAL QUEUE_TRIM,USER_INT
  DOUBLE PRECISION :: UROUND = EPSILON(1D0)
! Timing parameters:
  REAL DDTIME,DDTIME1,DDTIME2

  ! Files to which various output will be written:
  OPEN(UNIT=6, FILE='out.txt')  
  OPEN(UNIT=15, FILE='stat.txt')
  OPEN(UNIT=26, FILE='mysol.txt')
  
  ! Storage to save the results for each tolerance:
  ALLOCATE(MY_SAVE(NOUT,4),STAT=IER)
  CALL MY_CHECK_STAT(IER,4)

  ! Do we wish to do solution queue trimming?
  
  ! TRUE gives the solution with trimming. In this case,
  ! subroutine TRIM_GET will be called by the solver just
  ! before any trim is done to perform user interpolation.
    QUEUE_TRIM = .TRUE.
    
  ! FALSE gives the solution without trimming. In this
  ! case, the DDE_VAL interpolation option may then be
  ! following the call todde_solver.)
  ! Uncomment the following statement if you do not wish
  ! to do solution queue trimming.
  ! QUEUE_TRIM = .FALSE.

  ! If solution queue trimming is to be done, do we
  ! wish to first process the trimmed information
  ! using subroutine TRIM_GRT for interpolation
  ! purposes?
  !   (ignored if QUEUE_TRIM is FALSE)
  USER_INT = .TRUE.
  
  ! Problem parameters:
  a=1D0
  b=1D0
  bp=1D0
  bm=1D0
  ap=1D0
  am=1D0
  kappa=1D2
  
  ! Output parameters:
  T0=0D0
  TFINAL=1D3
  Length=5D2
  Tmin=TFINAL-Length

  ! Times at which the interpolated solution is desired:
  DO I = 1, NOUT
     T(I)=Tmin+real(I-1)*(Length/real(NOUT-1))
  END DO
  
  ! Interpolation arrays used in subroutine TRIM_GET:
  IF (QUEUE_TRIM .AND. USER_INT) THEN
     USER_NINT = NOUT
     ALLOCATE(USER_TINT(USER_NINT),USER_YINT(USER_NINT,NEQN),STAT=IER)
     CALL MY_CHECK_STAT(IER,5)
     USER_TINT(1:NOUT) = T(1:NOUT)
  END IF
  
! PARAMETER VARIATION CYCLE

! The parameter variation cycle begins here:
  JMAX = 11
  DO J = 1, JMAX
  PRINT *, ' J = ', J

     ! The error tolerance loop begins here:  
     NTOL = 1
     !NTOL = 4
     DO ITOL = 1, NTOL
        TOL = 1.0D0 / 10.0D0**(4+2*ITOL)
        ! Smallest allowable nonzero relative error tolerance:
        TOL = MAX(TOL,2.0D0*UROUND+1.0D-12)
        PRINT *, ' Error tolerance = ', TOL
        ! Delay
         DELAY(1) = Pi/2D0/kappa*10D0**REAL((J-1D0)/(JMAX-1D0))
        !DELAY(1) = 2D0
        MAXDELAY = MAXVAL(DELAY(1:NLAGS))

        ! Initialize the TRIM_GET search index:
        USER_START = 1

        ! DDE solver options:
        IF (QUEUE_TRIM) THEN
           ! If trimming the solution queue ...
           OPTS = DDE_SET(RE=TOL,AE=TOL,MAX_STEPS=10**5, &
                          MAX_DELAY=MAXDELAY,TRIM_FREQUENCY=1000)
        ELSE
           ! If not trimming the solution queue ...
           OPTS = DDE_SET(RE=TOL,AE=TOL,INTERPOLATION=.TRUE.,MAX_STEPS=10**5)
        END IF
        
        CALL CPU_TIME(DDTIME1)        
        ! DDE solution:
        IF (QUEUE_TRIM .AND. USER_INT) THEN
           ! If proceesing the trimmed results in TRIM_GET in order to do
           ! interpolation:
           SOL = DDE_SOLVER(NVAR,DDES,DELAY,HISTORY, TSPAN=(/ T0,TFINAL /), &
                            OPTIONS=OPTS,USER_TRIM_GET=TRIM_GET)
        ELSE
           ! If not processing the trimmed results (so, no interpolation) ...
           SOL = DDE_SOLVER(NVAR,DDES,DELAY,HISTORY, TSPAN=(/ T0,TFINAL /), &
                            OPTIONS=OPTS)
        END IF
        CALL CPU_TIME(DDTIME2)
        DDTIME = DDTIME2 - DDTIME1
        
        ! Was the solver successful?
        
        IF (SOL%FLAG == 0) THEN

           ! Integration statistics:
           CALL PRINT_STATS(SOL)

           ! Integrated solution:
           WRITE(UNIT=26,FMT='(2(F10.5))') (SOL%T(I),SOL%Y(I,1),I=1,SOL%NPTS)
           
           IF (QUEUE_TRIM) THEN
              IF (USER_INT) THEN
                ! Analyze the interpolated values calculated in TRIM_GET:
                 T(1:NOUT) = USER_TINT(1:NOUT)
                 MY_SAVE(1:NOUT,ITOL) = USER_YINT(1:NOUT,1)
                 YYT(1:NOUT) = USER_YINT(1:NOUT,1)
                 ! Release the intermediate work arrays and the option
                 ! and solution structures to prevent memory leak:
                 IF (ITOL < NTOL) THEN  
                    CALL RELEASE_ARRAYS(SOL,OPTS)
                 END IF
              ELSE
                 PRINT *, ' This test program currently USER_INT to be'
                 PRINT *, ' TRUE if TRIM_QUEUE is TRUE. Interpolated'
                 PRINT *, ' values are not available otherwise.'
                 STOP
              END IF
           ELSE
              ! Analyze the interpolated values calculated using the
              ! INTERPOLATION option:
              CALL CPU_TIME(DDTIME1)   
              Y = DDE_VAL(T,SOL)
              CALL CPU_TIME(DDTIME2)
              DDTIME = DDTIME + (DDTIME2 - DDTIME1)
              MY_SAVE(1:NOUT,ITOL) = Y%YT(1:NOUT,1)
              YYT(1:NOUT) = Y%YT(1:NOUT,1)
              ! Release the intermediate work arrays, the option and
              ! solution structures,and the INTERPOLATION structure
              ! to prevent memory leak:
              IF (ITOL < NTOL) THEN  
                 CALL RELEASE_ARRAYS(SOL,OPTS)
                 CALL RELEASE_INT(Y)
              END IF
           END IF
           
           ! Statistics;
           Sp=SUM(YYT,MASK=YYT.GT.0D0)/REAL(NOUT)
           Sm=SUM(YYT,MASK=YYT.LT.0D0)/REAL(NOUT)
           S=SUM(YYT)/REAL(NOUT)
           Var=SUM((YYT-S)**2)/REAL(NOUT)
           WRITE(UNIT=15,FMT='(5(1X,F10.5))') DELAY(1),Sp,Sm,S,Var

           ! Execution time for this pass:
           PRINT *, ' Execution time = ', DDTIME
           PRINT *, ' _____________________________________________ '           
         ELSE
           
           PRINT *, ' Abnormal return from DDE_SOLVER with FLAG = ', &
                      SOL%FLAG
           STOP
           
        END IF ! End of the (SOL%FLAG == 0) check

     END DO ! End of ITOL loop
     
  END DO ! End of JMAX loop
  
  ! Write the last interpolated solution to a file for subsequent
  ! plotting in Matlab:
  DO I = 1, NOUT
     WRITE(UNIT=6,FMT='(2F20.10)') T(I),YYT(I)
  END DO
        
  ! Calculate the estimated maximum error overrun for each tolerance:
  IF (NTOL > 1) THEN
     DO ITOL = 1, NTOL-1
        TOL = 1.0D0 / 10.0D0**(4+2*ITOL)
        OVERRUN = 0.0D0
        DO I = 1, NOUT
           OVERRUN = MAX(OVERRUN,ABS(MY_SAVE(I,ITOL)-MY_SAVE(I,NTOL)) / &
                         (TOL*(1.0D0+ABS(MY_SAVE(I,NTOL)))))
        END DO
        PRINT *, ' ITOL = ', ITOL, ' Estimated error overrun = ', OVERRUN
     END DO  
  END IF

  ! Clean up:
  DEALLOCATE(MY_SAVE,STAT=IER)
  CALL MY_CHECK_STAT(IER,6)    
  IF (QUEUE_TRIM) THEN
     CALL RELEASE_ARRAYS(SOL,OPTS)
     IF (USER_INT) THEN
        DEALLOCATE(USER_TINT,USER_YINT,STAT=IER)
        CALL MY_CHECK_STAT(IER,7)
     END IF
  ELSE
  END IF
  
  STOP
 END PROGRAM ilya