!   D02MWF Example Program Text
!   Mark 28.4 Release. NAG Copyright 2022.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: jac, res
!     .. Parameters ..
      Integer, Parameter, Public       :: iset = 1, neq = 5, nin = 5, nout = 6
      Integer, Parameter, Public       :: licom = 50 + neq
    Contains
      Subroutine res(neq,t,y,ydot,r,ires,iuser,ruser)

!       .. 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)

!       .. 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 d02mwfe_mod, Only: iset, jac, licom, neq, nin, nout, res
      Use nag_library, Only: d02mwf, d02nef, nag_wp, x04abf
!     .. 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