!   D02TVF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    Module d02tvfe_mod

!     D02TVF 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                           :: ffun, fjac, gafun, gajac, gbfun,     &
                                          gbjac, guess
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: mmax = 1, neq = 6, nin = 5,          &
                                          nlbc = 3, nout = 6, nrbc = 3
!     .. Local Scalars ..
      Real (Kind=nag_wp), Public, Save :: beta0, eta, lambda, mu
!     .. Local Arrays ..
      Integer, Public, Save            :: m(neq) = (/1,1,1,1,1,1/)
    Contains
      Subroutine ffun(x,y,neq,m,f,iuser,ruser)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(neq)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: y(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: beta
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos
!       .. Executable Statements ..
        beta = beta0*(one+cos(two*x01aaf(beta)*x))
        f(1) = mu - beta*y(1,0)*y(3,0)
        f(2) = beta*y(1,0)*y(3,0) - y(2,0)/lambda
        f(3) = y(2,0)/lambda - y(3,0)/eta
        f(4:6) = zero
        Return
      End Subroutine ffun
      Subroutine fjac(x,y,neq,m,dfdy,iuser,ruser)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: dfdy(neq,neq,0:*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: y(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: beta
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos
!       .. Executable Statements ..
        beta = beta0*(one+cos(two*x01aaf(beta)*x))
        dfdy(1,1,0) = -beta*y(3,0)
        dfdy(1,3,0) = -beta*y(1,0)
        dfdy(2,1,0) = beta*y(3,0)
        dfdy(2,2,0) = -one/lambda
        dfdy(2,3,0) = beta*y(1,0)
        dfdy(3,2,0) = one/lambda
        dfdy(3,3,0) = -one/eta
        Return
      End Subroutine fjac
      Subroutine gafun(ya,neq,m,nlbc,ga,iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: neq, nlbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: ga(nlbc)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: ya(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        ga(1) = ya(1,0) - ya(4,0)
        ga(2) = ya(2,0) - ya(5,0)
        ga(3) = ya(3,0) - ya(6,0)
        Return
      End Subroutine gafun
      Subroutine gbfun(yb,neq,m,nrbc,gb,iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: neq, nrbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: gb(nrbc)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: yb(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        gb(1) = yb(1,0) - yb(4,0)
        gb(2) = yb(2,0) - yb(5,0)
        gb(3) = yb(3,0) - yb(6,0)
        Return
      End Subroutine gbfun
      Subroutine gajac(ya,neq,m,nlbc,dgady,iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: neq, nlbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: dgady(nlbc,neq,0:*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: ya(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        dgady(1,1,0) = one
        dgady(1,4,0) = -one
        dgady(2,2,0) = one
        dgady(2,5,0) = -one
        dgady(3,3,0) = one
        dgady(3,6,0) = -one
        Return
      End Subroutine gajac
      Subroutine gbjac(yb,neq,m,nrbc,dgbdy,iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: neq, nrbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: dgbdy(nrbc,neq,0:*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: yb(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        dgbdy(1,1,0) = one
        dgbdy(1,4,0) = -one
        dgbdy(2,2,0) = one
        dgbdy(2,5,0) = -one
        dgbdy(3,3,0) = one
        dgbdy(3,6,0) = -one
        Return
      End Subroutine gbjac
      Subroutine guess(x,neq,m,y,dym,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: dym(neq)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*), y(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        y(1:3,0) = one
        y(4,0) = y(1,0)
        y(5,0) = y(2,0)
        y(6,0) = y(3,0)
        dym(1:neq) = zero
        Return
      End Subroutine guess
    End Module d02tvfe_mod
    Program d02tvfe

!     D02TVF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02tlf, d02tvf, d02tyf, d02tzf, nag_wp
      Use d02tvfe_mod, Only: beta0, eta, ffun, fjac, gafun, gajac, gbfun,      &
                             gbjac, guess, lambda, m, mmax, mu, neq, nin,      &
                             nlbc, nout, nrbc, one
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: dx, ermx
      Integer                          :: i, iermx, ifail, ijermx, licomm,     &
                                          lrcomm, mxmesh, ncol, nmesh
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: mesh(:), rcomm(:), tols(:), y(:,:)
      Real (Kind=nag_wp)               :: ruser(1)
      Integer, Allocatable             :: icomm(:), ipmesh(:)
      Integer                          :: iuser(2)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'D02TVF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) ncol, nmesh, mxmesh
      Allocate (mesh(mxmesh),tols(neq),y(neq,0:mmax-1),ipmesh(mxmesh))

      Read (nin,*) beta0, eta, lambda, mu
      Read (nin,*) tols(1:neq)

      dx = one/real(nmesh-1,kind=nag_wp)
      mesh(1) = 0.0_nag_wp
      Do i = 2, nmesh - 1
        mesh(i) = mesh(i-1) + dx
      End Do
      mesh(nmesh) = one

      ipmesh(1) = 1
      ipmesh(2:nmesh-1) = 2
      ipmesh(nmesh) = 1

!     Workspace query to get size of rcomm and icomm
      ifail = 0
      Call d02tvf(neq,m,nlbc,nrbc,ncol,tols,mxmesh,nmesh,mesh,ipmesh,ruser,0,  &
        iuser,2,ifail)
      lrcomm = iuser(1)
      licomm = iuser(2)
      Allocate (rcomm(lrcomm),icomm(licomm))

!     Initialize
      ifail = 0
      Call d02tvf(neq,m,nlbc,nrbc,ncol,tols,mxmesh,nmesh,mesh,ipmesh,rcomm,    &
        lrcomm,icomm,licomm,ifail)

!     Solve
      ifail = -1
      Call d02tlf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rcomm,icomm,iuser,   &
        ruser,ifail)

!     Extract mesh.
      ifail = -1
      Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rcomm,icomm,      &
        ifail)

      If (ifail/=1) Then
!       Print mesh statistics
        Write (nout,99999) nmesh, ermx, iermx, ijermx
        Write (nout,99998)(i,ipmesh(i),mesh(i),i=1,nmesh)

!       Print solution on mesh.
        Write (nout,99997)
        Do i = 1, nmesh

          ifail = 0
          Call d02tyf(mesh(i),y,neq,mmax,rcomm,icomm,ifail)
          Write (nout,99996) mesh(i), y(1:3,0)

        End Do
      End If

99999 Format (/,' Used a mesh of ',I4,' points',/,' Maximum error = ',E10.2,   &
        '  in interval ',I4,' for component ',I4,/)
99998 Format (/,' Mesh points:',/,4(I4,'(',I1,')',F7.4))
99997 Format (/,' Computed solution at mesh points',/,'    x       y1   ',     &
        '      y2         y3')
99996 Format (1X,F6.3,1X,3E11.3)
    End Program d02tvfe