!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File:  ct30ker.f90
! Kernel for the control interface.
!
! 06 Feb 2010: Current version.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

module ct30ker
  use precision,    only : ip, rp
  use ctUsrModule,  only : ctProb
  use ctWorkModule, only : ctWork, ctGrid, c0grid, c0gtrsh, c0wtrsh
  use ctGlobal,     only : nY, nU, nP, nC, nPhs, objL, ytype, ctype, &
                           nnY, nlY, nnC, nlC, neJ, neG, &
                           iDisc, iIntf, iRefn, iCprt, iRefL, refTol
  use ct11prnt,     only : c1prnt
  use ct16refn,     only : c1refn
  use ct25opt,      only : c2solv
  use ct29set,      only : c2setup
  implicit none

  private
  public  :: c3ker
  private :: c3err

contains

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine c3ker ( Start, prob, x, hs, odecon, algcon, varbds, &
                     INFO, mincw, miniw, minrw, nS, nInf, sInf, &
                     cu, lencu, iu, leniu, ru, lenru, &
                     cw, lencw, iw, leniw, rw, lenrw )
    integer(ip),  intent(in)    :: Start, lencu, leniu, lenru, &
                                   lencw, leniw, lenrw
    integer(ip),  intent(inout) :: nS, iu(leniu), iw(leniw)
    real(rp),     intent(inout) :: rw(lenrw), ru(lenru)
    character(8), intent(inout) :: cw(lencw), cu(lencu)
    integer(ip),  intent(inout), pointer :: hs(:)
    real(rp),     intent(inout), pointer :: x(:)
    integer(ip),  intent(out)   :: INFO, mincw, miniw, minrw, nInf
    real(rp),     intent(out)   :: sInf
    external             :: odecon, algcon, varbds
    type(ctProb), target :: prob

    !===========================================================================
    ! On entry, global option pointers have been set.
    ! 1. Sets up global problem pointers.
    ! 2. Create the uniform grid based on user data and set up the workspace.
    ! 3. Call c2solv to solve problem.
    ! 4. Print solution.
    ! 5. Check for refinement.
    !
    ! 26 Dec 2008: First version.
    ! 09 Feb 2010: v5.
    !===========================================================================
    integer(ip)    :: nItn, refLoop, nAdded, i
    real(rp)       :: r, infBnd
    character(80)  :: str
    type(ctGrid)   :: grid
    type(ctWork)   :: work

    ! Error checks
    call c3err ( iIntf, prob, iw, leniw )

    ! If refining, set to TR method.
    if ( iRefn /= 0 ) &
         iDisc = 0

    nItn    = 0
    refLoop = 0
    if ( iRefL < 0 ) &   ! Refine until no nodes are added.
         refLoop = 1

    infBnd = rw(70)
    if ( infBnd <= 0.0 ) infBnd = 1.0d+20

    ! Setup global information.
    nY    => prob%nY
    nU    => prob%nU
    nP    => prob%nP
    nC    => prob%nC
    nPhs  => prob%nPhs
    objL  => prob%objL
    ytype => prob%ytype

    if ( nC > 0 ) then
       ctype => prob%ctype
    else
       allocate ( ctype(1,nPhs) )
       ctype = 0
    end if

    allocate ( nnY(nPhs), nlY(nPhs), nnC(nPhs), nlC(nPhs) )
    nnY = 0
    nlY = 0
    nnC = 0
    nlC = 0

    nlY = sum ( ytype, dim=1 )
    nnY = nY - nlY

    if ( nC > 0 ) then
       nlC = sum ( ctype, dim=1 )
       nnC = nC - nlC
    end if

    ! Derivatives in each phase.
    if ( iIntf == 0 ) then
       neJ => prob%neJ

       if ( nC > 0 ) then
          neG => prob%neG
       else
          allocate ( neG(nPhs) )
          neG = 0
       end if

    else
       allocate ( neJ(nPhs), neG(nPhs) )
       neJ = nY*( nY + nU + nP )
       neG = nC*( nY + nU + nP )
    end if


    ! Set the (uniform) grid.
    call c0grid ( iDisc, grid, nPhs, prob%phsPt, prob%npInt )

    ! Set the workspace.
    call c2setup ( iDisc, iIntf, nPhs, nY, nU, nP, nC, objL, ytype, ctype, &
                   nnY, nlY, nnC, nlC, neJ, neG, x, hs, &
                   grid, work, odecon, algcon, varbds, infBnd, &
                   cu, lencu, iu, leniu, ru, lenru, iw, leniw )

    ! Solve the problem.
    call c2solv ( Start, prob%probName, x, hs, work, grid, odecon, algcon, &
                   INFO, mincw, miniw, minrw, nS, nInf, sInf, &
                   cu, lencu, iu, leniu, ru, lenru, &
                   cw, lencw, iw, leniw, rw, lenrw )
    nItn = nItn + 1

    ! Print solution
    if ( iCPrt /= 0 ) &
         call c1prnt ( prob, grid, x, work, iw, leniw )

    r = INFO / 10.0
    if ( iRefn == 0 .or. floor(r) /= 0 ) &
         go to 900


    !---------------------------------------------------------------------------
    ! Adaptive Refinement
    ! sncRefn will create a new grid (replacing the old one) and calculate a
    ! new initial point for the refined problem.

100 call c1refn ( iIntf, refTol, nPhs, nY, nU, nP, neJ, work%n, x, hs, grid, &
                  odecon, nAdded, cu, lencu, iu, leniu, ru, lenru )

    ! Print a header...
    str = ''
    call snPRNT (13, str, iw, leniw)
    do i = 1, 80
       str = str(1:i) // '*'
    end do
    call snPRNT (13, str, iw, leniw)
    write(str,*) '   Refinement Run: ', nAdded, ' nodes added.'
    call snPRNT ( 3, str, iw, leniw)

    if ( nAdded == 0 ) &
         go to 900

    ! Set up the optimal control problem.
    call c0wtrsh ( work )
    call c2setup ( iDisc, iIntf, nPhs, nY, nU, nP, nC, objL, ytype, ctype, &
                   nnY, nlY, nnC, nlC, neJ, neG, x, hs, &
                   grid, work, odecon, algcon, varbds, infBnd, &
                   cu, lencu, iu, leniu, ru, lenru, iw, leniw )

    ! Solve the problem.
    call c2solv ( 0_ip, prob%probName, x, hs, work, grid, odecon, algcon, &
                  INFO, mincw, miniw, minrw, nS, nInf, sInf, &
                  cu, lencu, iu, leniu, ru, lenru, &
                  cw, lencw, iw, leniw, rw, lenrw )

    nItn = nItn + 1

    ! Print solution
    if ( iCPrt /= 0 ) &
         call c1prnt ( prob, grid, x, work, iw, leniw )

    ! Error exit?
    r = INFO / 10.0
    if ( floor(r) /= 0 ) &
         go to 900

    if ( nItn <= iRefL .or. refLoop == 1 ) &
         go to 100


    ! Clean up.
900 call c0wtrsh ( work )
    call c0gtrsh ( grid )

    deallocate ( nnY, nlY, nnC, nlC )

    if ( nC == 0 ) &
         deallocate ( ctype )

    if ( iIntf /= 0 ) then
       deallocate ( neJ, neG )
    else
       if ( nC == 0 ) &
            deallocate ( neG )
    end if

  end subroutine c3ker

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine c3err ( iIntf, prob, iw, leniw )
    integer(ip), intent(in)    :: iIntf, leniw
    integer(ip), intent(inout) :: iw(leniw)
    type(ctProb) :: prob

    !===========================================================================
    ! Very simple error checks on the ctProb structure.
    !
    ! 21 Aug 2008: First version of errCheck.
    ! 24 Dec 2008: Updated for v4 of snctrl.
    ! 09 Feb 2010: v5.
    !===========================================================================

    if ( prob%nY < 0 ) go to 900
    if ( prob%nU < 0 ) go to 900
    if ( prob%nC < 0 ) go to 900
    if ( prob%nP < 0 ) go to 900
    if ( prob%nPhs < 0 ) go to 900
    if ( .not. associated(prob%objL) ) go to 900
    if ( .not. associated(prob%ytype) ) go to 900

    if ( prob%nC > 0 ) then
       if ( .not. associated(prob%ctype) ) go to 900
    end if
    if ( .not. associated(prob%npInt) ) go to 900
    if ( .not. associated(prob%phsPt) ) go to 900

    if ( iIntf == 0 ) then
       if ( .not. associated(prob%neJ) ) go to 900
       if ( prob%nC > 0 ) then
          if ( .not. associated(prob%neG) ) go to 900
       end if
    end if

    return

900 call snPRNT (13, 'snctrl: STOPPING because of input error.', iw, leniw)
    stop

  end subroutine c3err

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

end module ct30ker

