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

    Module d02bj_a1w_fe_mod

!     Data for D02BJ_A1W_F example program

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: cos, nagad_a1w_w_rtype, tan, Assignment (=),    &
                               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_a1w_w_rtype), Intent (In) :: x
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Inout) :: f(*), ruser(*)
        Type (nagad_a1w_w_rtype), Intent (In) :: y(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Type (nagad_a1w_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_a1w_w_rtype), Intent (Out) :: retval
        Type (nagad_a1w_w_rtype), Intent (In) :: x
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_a1w_w_rtype), Intent (In) :: y(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Executable Statements ..
        retval = y(1)
        Return
      End Subroutine g
    End Module d02bj_a1w_fe_mod
    Program d02bj_a1w_fe

!     D02BJ_A1W_F Example Main Program

!     .. Use Statements ..
      Use d02bj_a1w_fe_mod, Only: fcn, g, n, nin, nout
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: d02bj_a1w_f, d02bj_a1w_x,                       &
                               nagad_a1w_get_derivative,                       &
                               nagad_a1w_ir_interpret_adjoint_sparse,          &
                               nagad_a1w_ir_register_variable,                 &
                               nagad_a1w_ir_remove, nagad_a1w_set_derivative,  &
                               nagad_a1w_w_rtype, x10aa_a1w_f, x10ab_a1w_f,    &
                               x10za_a1w_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_a1w_w_rtype)         :: tol, x, xend_ad
      Real (Kind=nag_wp)               :: da, db, xend, xinit
      Integer                          :: i, ifail, iw, kinit
!     .. Local Arrays ..
      Type (nagad_a1w_w_rtype)         :: ruser(4)
      Type (nagad_a1w_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_A1W_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

      x = xinit
      xend_ad = xend


!     Create AD tape
      Call x10za_a1w_f

!     Create AD configuration data object
      ifail = 0
      Call x10aa_a1w_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

!     Register variables to differentiate w.r.t.
      Call nagad_a1w_ir_register_variable(ruser(1:2))
      Call nagad_a1w_ir_register_variable(y_in)

      y(1:n) = y_in(1:n)

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

      Write (nout,99998) '  Root of Y(1) = 0.0 at', x%value
      Write (nout,99997) '  Solution is', (y(i)%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: First order adjoints'
      Write (nout,*) ' Computational mode    : algorithmic'

!     Setup evaluation of derivatives via adjoints
      Call nagad_a1w_set_derivative(x,1.0_nag_wp)
      ifail = 0
      Call nagad_a1w_ir_interpret_adjoint_sparse(ifail)

      da = nagad_a1w_get_derivative(ruser(1))
      db = nagad_a1w_get_derivative(ruser(2))
      yinit = nagad_a1w_get_derivative(y_in)

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

!     Remove computational data object and tape
      Call x10ab_a1w_f(ad_handle,ifail)
      Call nagad_a1w_ir_remove

    End Program d02bj_a1w_fe