!   D02BJ_T2W_F Example Program Text
!   Mark 28.3 Release. NAG Copyright 2022.

    Module d02bj_t2w_fe_mod

!     Data for D02BJ_T2W_F example program

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: cos, nagad_t2w_w_rtype, tan, Operator (+),      &
                               Operator (/), Operator (*), Operator (**)
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: fcn, g
!     .. Parameters ..
      Integer, Parameter, Public       :: n = 3, nin = 5, nout = 6
!     n: number of differential equations
    Contains
      Subroutine fcn(ad_handle,x,y,f,iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_t2w_w_rtype), Intent (In) :: x
!       .. Array Arguments ..
        Type (nagad_t2w_w_rtype), Intent (Inout) :: f(*), ruser(*)
        Type (nagad_t2w_w_rtype), Intent (In) :: y(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Type (nagad_t2w_w_rtype)       :: alpha, beta
!       .. Executable Statements ..
        alpha = ruser(1)
        beta = ruser(2)
        f(1) = tan(y(3))
        f(2) = alpha*tan(y(3))/y(2) + beta*y(2)/cos(y(3))
        f(3) = alpha/y(2)**2
        Return
      End Subroutine fcn
      Subroutine g(ad_handle,x,y,retval,iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_t2w_w_rtype), Intent (Out) :: retval
        Type (nagad_t2w_w_rtype), Intent (In) :: x
!       .. Array Arguments ..
        Type (nagad_t2w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_t2w_w_rtype), Intent (In) :: y(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Executable Statements ..
        retval = y(1)
        Return
      End Subroutine g
    End Module d02bj_t2w_fe_mod
    Program d02bj_t2w_fe

!     D02BJ_T2W_F Example Main Program

!     .. Use Statements ..
      Use d02bj_t2w_fe_mod, Only: fcn, g, n, nin, nout
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: d02bj_t2w_f, d02bj_t2w_x, nagad_t2w_w_rtype,    &
                               x10aa_t2w_f, x10ab_t2w_f, Assignment (=)
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: alpha = -0.032E0_nag_wp
      Real (Kind=nag_wp), Parameter    :: beta = -0.02E0_nag_wp
!     .. Local Scalars ..
      Type (c_ptr)                     :: ad_handle
      Type (nagad_t2w_w_rtype)         :: tol, x, xend_ad
      Real (Kind=nag_wp)               :: da, db, xend, xinit
      Integer                          :: i, ifail, iw, kinit
!     .. Local Arrays ..
      Type (nagad_t2w_w_rtype)         :: ruser(4)
      Type (nagad_t2w_w_rtype), Allocatable :: w(:), y(:), y_in(:)
      Real (Kind=nag_wp), Allocatable  :: yinit(:)
      Integer                          :: iuser(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'D02BJ_T2W_F Example Program Results'

      iw = 20*n
      Allocate (w(iw),y(n),yinit(n),y_in(n))
!     Skip heading in data file
      Read (nin,*)
!     xinit: initial x value,        xend: final x value.
!     yinit: initial solution values
      Read (nin,*) xinit, xend
      Read (nin,*) yinit(1:n)
      Read (nin,*) kinit

      Write (nout,99996) 'no intermediate output, root-finding'
      tol = 1.0E-5_nag_wp
      Write (nout,*)
      Write (nout,99999) ' Calculation with TOL =', tol%value%value

!     Create AD configuration data object
      ifail = 0
      Call x10aa_t2w_f(ad_handle,ifail)

      y_in(1:n) = yinit(1:n)
      ruser(1) = alpha
      ruser(2) = beta
      ruser(3) = (xend-xinit)/real(kinit+1,kind=nag_wp)
      ruser(4) = xend

      Do i = 1, 2
        ruser(i)%value%tangent = 1.0_nag_wp
        ruser(i)%tangent%value = 1.0_nag_wp
        x = xinit
        xend_ad = xend
        y(1:n) = y_in(1:n)

        ifail = 0
        Call d02bj_t2w_f(ad_handle,x,xend_ad,n,y,fcn,tol,'Default',            &
          d02bj_t2w_x,g,w,iuser,ruser,ifail)

        ruser(i)%tangent%value = 0.0_nag_wp
        If (i==1) Then
          da = x%tangent%tangent
        Else
          db = x%tangent%tangent
        End If
        ruser(i)%value%tangent = 0.0_nag_wp
      End Do
      Do i = 1, 3
        x = xinit
        xend_ad = xend
        y(1:n) = y_in(1:n)

        y(i)%value%tangent = 1.0_nag_wp
        y(i)%tangent%value = 1.0_nag_wp

        ifail = 0
        Call d02bj_t2w_f(ad_handle,x,xend_ad,n,y,fcn,tol,'Default',            &
          d02bj_t2w_x,g,w,iuser,ruser,ifail)

        yinit(i) = x%tangent%tangent

        y(i)%tangent%value = 0.0_nag_wp
        y(i)%value%tangent = 0.0_nag_wp
      End Do

      Write (nout,99998) '  Root of Y(1) = 0.0 at', x%value%value
      Write (nout,99997) '  Solution is', (y(i)%value%value,i=1,n)
      Write (nout,*)

99999 Format (1X,A,E8.1)
99998 Format (1X,A,F7.3)
99997 Format (1X,A,3F13.4)
99996 Format (1X,'Case : ',A)

      Write (nout,*)
      Write (nout,*) ' Derivatives calculated: Second order tangents'
      Write (nout,*) ' Computational mode    : algorithmic'

      Write (nout,*)
      Write (nout,*) ' Derivatives: (hit point w.r.t. parameters)'
      Write (nout,99995) '     d^2x/dg^2      =', -da
      Write (nout,99995) '     d^2x/ddrag^2   =', -db
      Write (nout,99995) '     d^2x/dheight^2 =', yinit(1)
      Write (nout,99995) '     d^2x/dvel^2    =', yinit(2)
      Write (nout,99995) '     d^2x/dangle^2  =', yinit(3)
99995 Format (1X,A,1X,E12.5)

!     Remove computational data object
      ifail = 0
      Call x10ab_t2w_f(ad_handle,ifail)

    End Program d02bj_t2w_fe