!   D02MWF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d02mwfe_mod

!      D02MWF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: iset = 1, neq = 5, nin = 5,      &
                                              nout = 6
       INTEGER, PARAMETER                  :: licom = 50 + neq
    CONTAINS
       SUBROUTINE res(neq,t,y,ydot,r,ires,iuser,ruser)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: r(neq)
          REAL (KIND=nag_wp), INTENT (INOUT)  :: ruser(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq), ydot(neq)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. Executable Statements ..
          r(1) = y(3) - ydot(1)
          r(2) = y(4) - ydot(2)
          r(3) = -y(5)*y(1) - ydot(3)
          r(4) = -y(5)*y(2) - 1.0_nag_wp - ydot(4)
          r(5) = y(3)**2 + y(4)**2 - y(5) - y(2)
          RETURN
       END SUBROUTINE res

       SUBROUTINE jac(neq,t,y,ydot,pd,cj,iuser,ruser)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: cj, t
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: pd(*), ruser(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq), ydot(neq)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. Executable Statements ..
          pd(1) = -cj
          pd(3) = -y(5)
          pd(7) = -cj
          pd(9) = -y(5)
          pd(10) = -1.0_nag_wp
          pd(11) = 1.0_nag_wp
          pd(13) = -cj
          pd(15) = 2.0_nag_wp*y(3)
          pd(17) = 1.0_nag_wp
          pd(19) = -cj
          pd(20) = 2.0_nag_wp*y(4)
          pd(23) = -y(1)
          pd(24) = -y(2)
          pd(25) = -1.0_nag_wp
          RETURN
       END SUBROUTINE jac
    END MODULE d02mwfe_mod

    PROGRAM d02mwfe

!      D02MWF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02mwf, d02nef, nag_wp, x04abf
       USE d02mwfe_mod, ONLY : iset, jac, licom, neq, nin, nout, res
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: g1, g2, h0, hmax, t, tout
       INTEGER                             :: i, ifail, ijac, itask, itol,     &
                                              lcom, maxord, nadv
       CHARACTER (8)                       :: jceval
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: atol(:), com(:), rtol(:), y(:),  &
                                              ydot(:)
       REAL (KIND=nag_wp)                  :: ruser(1)
       INTEGER, ALLOCATABLE                :: icom(:)
       INTEGER                             :: iuser(1)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02MWF Example Program Results'
       WRITE (nout,*)
!      Skip heading in data file
       READ (nin,*)
!      neq: number of differential-algebraic equations
       READ (nin,*) maxord
       lcom = 40 + (maxord+4)*neq + neq**2
       ALLOCATE (atol(neq),com(lcom),rtol(neq),y(neq),ydot(neq),icom(licom))
       nadv = nout
       CALL x04abf(iset,nadv)
       READ (nin,*) t, tout
       READ (nin,*) itol, itask
       READ (nin,*) rtol(1:neq)
       READ (nin,*) atol(1:neq)
!      Set initial values
       READ (nin,*) y(1:neq)
       READ (nin,*) hmax, h0
       READ (nin,*) ijac
       IF (ijac==1) THEN
          jceval = 'Analytic'
       ELSE
          jceval = 'Numeric'
       END IF

!      ifail: behaviour on error exit   
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
       ifail = 0
       CALL d02mwf(neq,maxord,jceval,hmax,h0,itol,icom,licom,com,lcom,ifail)

       WRITE (nout,99995) (i,i=1,neq)
       WRITE (nout,99998) t, y(1:neq)
       ydot(1:neq) = 0.0_nag_wp

LOOP:  DO
          CALL d02nef(neq,t,tout,y,ydot,rtol,atol,itask,res,jac,icom,com,lcom, &
             iuser,ruser,ifail)

          WRITE (nout,99998) t, y(1:neq)
          WRITE (nout,99999) itask
          IF ((itask>=0) .AND. (itask<=3)) THEN
             IF (t>=tout) THEN
                g1 = y(1)**2 + y(2)**2 - 1.0_nag_wp
                g2 = y(1)*y(3) + y(2)*y(4)
                WRITE (nout,99997) g1
                WRITE (nout,99996) g2
                EXIT LOOP
             END IF
          ELSE
             EXIT LOOP
          END IF
       END DO LOOP

99999  FORMAT (/' D02NEF returned with ITASK = ',I4/)
99998  FORMAT (1X,F7.4,2X,5(F11.6))
99997  FORMAT (1X,'The position-level constraint G1 = ',E12.4)
99996  FORMAT (1X,'The velocity-level constraint G2 = ',E12.4)
99995  FORMAT (/1X,'  t ',3X,5('       y(',I1,')'))
    END PROGRAM d02mwfe