MODULE define_DDEs

! Autocatalytic oscillations (Borelli & Coleman, p. 260)

  IMPLICIT NONE

! Define the number of ODEs.
  INTEGER, PARAMETER :: NEQN=2,NLAGSD=0
  
CONTAINS
    
  SUBROUTINE ODES(T,Y,DY)
! Define the derivatives for ODEAVG.
    DOUBLE PRECISION :: T,W
    DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY
    INTENT(IN) :: T,Y
    INTENT(OUT) :: DY
    W = 500.0D0*EXP(-0.002*T)
    DY(1) = 0.002D0*W - 0.08D0*Y(1) - Y(1)*Y(2)*Y(2)
    DY(2) =     -Y(2) + 0.08D0*Y(1) + Y(1)*Y(2)*Y(2)
    RETURN
  END SUBROUTINE ODES

  SUBROUTINE DDES(T,Y,Z,DY)
! Define the derivatives for pure ODE check solution.
    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY
    DOUBLE PRECISION, DIMENSION(NEQN,:) :: Z
    INTENT(IN) :: T,Y,Z
    INTENT(OUT) :: DY
    CALL ODES(T,Y,DY)
    RETURN
  END SUBROUTINE DDES

END MODULE define_DDEs

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

PROGRAM autocatalyticf90

! The DDEs are defined in the module define_DDEs. The problem is
! solved here with ODEAVG/DDE_SOLVER and its output is written
! to a file autocatalyticf90.dat. The auxilary function
! autocatalyticf90.m imports the data into Matlab and plots it.

  USE define_DDEs
  USE DDE_SOLVER_M

  IMPLICIT NONE

  ! The quantities passed to the solver in NVAR are
  !   NEQN  = number of ODEs
  !   AVG_OPT = averaging option
  ! NEQN is defined in the module define_DDEs as a
  ! PARAMETER so it can be used for dimensioning
  ! arrays here.

  INTEGER, DIMENSION(2) :: NVAR
  INTEGER, DIMENSION(2) :: NVARD = (/NEQN,NLAGSD/)
  DOUBLE PRECISION, DIMENSION(NEQN) :: Y0
  DOUBLE PRECISION, DIMENSION(1) :: DELTA,BETA
  DOUBLE PRECISION, DIMENSION(NEQN) :: ABSERR,RELERR

  TYPE(DDE_SOL) :: SOL,SOLODE
  !   SOL%NPTS         -- NPTS,number of output points.
  !   SOL%T(NPTS)      -- values of independent variable, T.  
  !   SOL%Y(NPTS,NEQN) -- values of dependent variable, Y,
  !                       corresponding to values of SOL%T. 
  TYPE(DDE_OPTS) :: OPTS,OPTODE
  DOUBLE PRECISION :: T0,TF,XINIT,YINIT
  INTEGER :: AVG_OPT(2),I,J,MAVG,MKIND
  
! Open the output files used by the autocatalyticf90.m plotting script.
  OPEN(UNIT=6, FILE='autocatalyticf90.dat')
  OPEN(UNIT=7, FILE='autocatalyticf90d.dat')
  OPEN(UNIT=8, FILE='autocatalyticf90ode.dat')
  OPEN(UNIT=9, FILE='autocatalyticf90temp.dat')
  
! Define the averaging option to be used.
  PRINT *,' Please enter 1 or 2, the number of averages to be'
  PRINT *,' computed, followed by <return>.'
  READ  *, MAVG
! PRINT *,' Please enter 1 for solution averaging or 2 for'
! PRINT *,' RMS amplitude averaging, followed by <return>.'
! READ  *, MKIND
  PRINT *,' '
  MKIND = 1
  AVG_OPT(1) = MAVG
  AVG_OPT(2) = MKIND

  AVG_OPT(1) = MAX(AVG_OPT(1),1)
  AVG_OPT(1) = MIN(AVG_OPT(1),2)
  AVG_OPT(2) = MAX(AVG_OPT(2),1)
  AVG_OPT(2) = MIN(AVG_OPT(2),2)

! Averaging parameter:
  DELTA(1) = 25.0D0

! Save DELTA and AVG_OPT for use in the Matlab plotting
! script autocatalyticf90.m.
  WRITE(7,*) DELTA(1), AVG_OPT(1)

! Define the integration parameters.
  NVAR = (/NEQN,AVG_OPT(1)/)
  T0 = 0.0D0
  XINIT = 0.0D0
  YINIT = 0.0D0
  TF = 1000.0D0
  Y0 = (/XINIT,YINIT/)

! Supply the integration options for the pure ODE check solution.
  ABSERR(1:NEQN) = 1.0D-10
  RELERR(1:NEQN) = 1.0D-4
  ABSERR(1:NEQN) = 1.0D-6
  RELERR(1:NEQN) = 1.0D-3
  BETA(1) = 0.0D0
  OPTODE = DDE_SET(RE_VECTOR=RELERR,AE_VECTOR=ABSERR)
  SOLODE = DDE_SOLVER(NVARD,DDES,BETA,Y0,TSPAN=(/T0,TF/),OPTIONS=OPTODE) 
  IF (SOLODE%FLAG == 0) THEN
!    Print integration statistics.
     CALL PRINT_STATS(SOLODE)
!    If the solver was successful, write the solution
!    to a file for subsequent plotting in Matlab.
     DO I = 1, SOLODE%NPTS
        WRITE(UNIT=8,FMT='(3D12.4)') SOLODE%T(I),(SOLODE%Y(I,J),J=1,NEQN)
     END DO
     PRINT *,' Normal return from DDE_SOLVER with results'
     PRINT *,' written to the file autocatalyticf90ode.dat.'
     PRINT *,' '
  ELSE
     PRINT *,' Abnormal return from DDE_SOLVER with FLAG = ',&
     SOLODE%FLAG
     STOP
  END IF

! Supply the integration options for the ODEAVG solution.
  ABSERR(1:NEQN) = 1.0D-10
  RELERR(1:NEQN) = 1.0D-4
  ABSERR(1:NEQN) = 1.0D-6
  RELERR(1:NEQN) = 1.0D-3
  OPTS = DDE_SET(RE_VECTOR=RELERR,AE_VECTOR=ABSERR)

! Perform the integration.
  SOL = ODEAVG(AVG_OPT,DELTA,NVAR,ODES,Y0,TSPAN=(/T0,TF/),OPTIONS=OPTS) 
  IF (SOL%FLAG == 0) THEN
!    Print integration statistics.
     CALL PRINT_STATS(SOL)   
!    If the solver was successful, write the solution
!    to a file for subsequent plotting in Matlab.
     DO I = 1, SOL%NPTS
        WRITE(UNIT=6,FMT='(3D12.4)') SOL%T(I),(SOL%Y(I,J),J=1,NEQN)
     END DO
     PRINT *,' Normal return from DDE_SOLVER with results'
     PRINT *,' written to the file .dat.'
     PRINT *,' '
     PRINT *,' These results can be accessed in Matlab'
     PRINT *,' and plotted by'
     PRINT *,' '
     PRINT *,' >> [t,y] = autocatalyticf90.m'
     PRINT *,' '
  ELSE
     PRINT *,' Abnormal return from DDE_SOLVER with FLAG = ', SOL%FLAG
     STOP
  END IF

  STOP
END PROGRAM autocatalyticf90