glacier_integrate Subroutine

  subroutine glacier_integrate(this, old_states, basal_melt, basal_drag, &
                               water_density, time, success)
    !* Author: Chris MacMackin
    !  Date: November 2016
    !
    ! Integrates the glacier's state to `time`. This is done using the
    ! NITSOL package of iterative Krylov solvers. If a different
    ! algorithm for the integration is desired, then this method may
    ! be overridden in the concrete implementations of the glacier
    ! type.
    !
    class(glacier), intent(inout)            :: this
    class(glacier), dimension(:), intent(in) :: old_states
      !! Previous states of the glacier, with the most recent one
      !! first.
    class(scalar_field), intent(in)          :: basal_melt
      !! The melt rate that the bottom of the glacier experiences
      !! during this time step.
    class(scalar_field), intent(in)          :: basal_drag
      !! A paramter, e.g. coefficient of friction, needed to calculate
      !! the drag on basal surface of the glacier.
    real(r8), intent(in)                     :: water_density
      !! The density of the water below the glacier.
    real(r8), intent(in)                     :: time
      !! The time to which the glacier should be integrated
    logical, intent(out)                     :: success
      !! True if the integration is successful, false otherwise

    logical                                   :: first_call
    integer, save                             :: nval, kdmax = 20
    real(r8), dimension(:), allocatable       :: state
    integer, dimension(10)                    :: input
    integer, dimension(6)                     :: info
    real(r8), dimension(:), allocatable, save :: work
    real(r8), dimension(1)                    :: real_param
    integer, dimension(1)                     :: int_param
    integer                                   :: flag

    call basal_melt%guard_temp(); call basal_drag%guard_temp()
    first_call = .true.
    nval = this%data_size()
    if (allocated(work)) then
      if (size(work) < nval*(kdmax+5) + kdmax*(kdmax+3)) then
        deallocate(work)
        allocate(work(nval*(kdmax+5) + kdmax*(kdmax+3)))
      end if
    else
      allocate(work(nval*(kdmax+5) + kdmax*(kdmax+3)))
    end if
    state = this%state_vector()
    call this%set_time(time)
    input = 0
    input(4) = kdmax
    input(5) = 1
    input(9) = -1
    input(10) = 3
    etafixed = 0.3_r8

#ifdef DEBUG
    call logger%debug('glacier%integrate','Calling NITSOL (nonlinear solver)')
#endif
    call nitsol(nval, state, nitsol_residual, nitsol_precondition, &
                1.e-10_r8*nval, 1.e-10_r8*nval, input, info, work, &
                real_param, int_param, flag, ddot, dnrm2)
    call this%update(state)
!!$    if (flag == 6 .and. input(9) > -1) then
!!$      input(9) = -1
!!$      call logger%trivia('glacier%integrate','Backtracking failed in NITSOL '// &
!!$                         'at simulation time '//str(time)//'. Trying again '//  &
!!$                         'without backtracking.')
!!$      call nitsol(nval, state, nitsol_residual, nitsol_precondition, &
!!$                  1.e-7_r8, 1.e-7_r8, input, info, work, real_param, &
!!$                  int_param, flag, ddot, dnrm2)
!!$    end if
#ifdef DEBUG
    call logger%debug('glacier%integrate','NITSOL required '//       &
                      trim(str(info(5)))//' nonlinear iterations '// &
                      'and '//trim(str(info(1)))//' function calls.')
#endif

    select case(flag)
    case(0)
      call logger%trivia('glacier%integrate','Integrated glacier to time '// &
                         trim(str(time)))
      success = .true.
    case(1)
      call logger%error('glacier%integrate','Reached maximum number of'// &
                        ' iterations integrating glacier')
      success = .false.
    !case(5)
    !  call logger%debug('glacier%integrate','Solution diverging. Trying '// &
    !                    'again with backtracking.')
    !  state = old_states(1)%state_vector()
    !  input(9) = 0
    !  call nitsol(nval, state, nitsol_residual, nitsol_precondition, &
    !              1.e-7_r8, 1.e-7_r8, input, info, work, real_param, &
    !              int_param, flag, ddot, dnrm2)
    !  call this%update(state)
    !  if (flag == 0) then
    !    call logger%trivia('glacier%integrate','Integrated glacier to time '// &
    !                       trim(str(time)))
    !    success = .true.
    !  else
    !    call logger%error('glacier%integrate','NITSOL failed when integrating'// &
    !                      ' glacier with error code '//trim(str(flag)))
    !    success = .false.
    !  end if
    case default
      call logger%error('glacier%integrate','NITSOL failed when integrating'// &
                        ' glacier with error code '//trim(str(flag)))
      success = .false.
    end select

    if (success) then
      call this%integrate_layers(old_states, time, success)
    end if

    call basal_melt%clean_temp(); call basal_drag%clean_temp()

  contains
    
    subroutine nitsol_residual(n, xcur, fcur, rpar, ipar, itrmf)
      !! A routine matching the interface expected by NITSOL which
      !! returns the residual for the glacier.
      integer, intent(in)                   :: n
        !! Dimension of the problem
      real(r8), dimension(n), intent(in)    :: xcur
        !! Array of length `n` containing the current \(x\) value
      real(r8), dimension(n), intent(out)   :: fcur
        !! Array of length `n` containing f(xcur) on output
      real(r8), dimension(*), intent(inout) :: rpar
        !! Parameter/work array
      integer, dimension(*), intent(inout)  :: ipar
        !! Parameter/work array
      integer, intent(out)                  :: itrmf
        !! Termination flag. 0 means normal termination, 1 means
        !! failure to produce f(xcur)

      logical :: success
      ! If this is the first call of this routine then the
      ! basal_surface object will already be in the same state as
      ! reflected in xcur
      if (first_call) then
        first_call = .false.
      else
        call this%update(xcur(1:n))
      end if
      call this%solve_velocity(basal_drag, success)
      if (.not. success) then
        itrmf = 1
        return
      end if
      fcur(1:n) = this%residual(old_states,basal_melt,basal_drag,water_density)
      !print*, fcur(1:n)
      itrmf = 0
    end subroutine nitsol_residual

    subroutine nitsol_precondition(n, xcur, fcur, ijob, v, z, rpar, ipar, itrmjv)
      !! A subroutine matching the interface expected by NITSOL, which
      !! acts as a preconditioner.
      integer, intent(in)                   :: n
        ! Dimension of the problem
      real(r8), dimension(n), intent(in)    :: xcur
        ! Array of lenght `n` containing the current $x$ value
      real(r8), dimension(n), intent(in)    :: fcur
        ! Array of lenght `n` containing the current \(f(x)\) value
      integer, intent(in)                   :: ijob
        ! Integer flat indicating which product is desired. 0
        ! indicates \(z = J\vec{v}\). 1 indicates \(z = P^{-1}\vec{v}\).
      real(r8), dimension(n), intent(in)    :: v
        ! An array of length `n` to be multiplied
      real(r8), dimension(n), intent(out)   :: z
        ! An array of length n containing the desired product on
        ! output.
      real(r8), dimension(*), intent(inout) :: rpar
        ! Parameter/work array 
      integer, dimension(*), intent(inout)  :: ipar
        ! Parameter/work array
      integer, intent(out)                  :: itrmjv
        ! Termination flag. 0 indcates normal termination, 1
        ! indicatesfailure to prodce $J\vec{v}$, and 2 indicates
        ! failure to produce \(P^{-1}\vec{v}\)
      if (ijob /= 1) then
        itrmjv = 0
        return
      end if
      z(1:n) = this%precondition(old_states, basal_melt, basal_drag, &
                                 water_density, v(1:n))
      itrmjv = 0
    end subroutine nitsol_precondition

  end subroutine glacier_integrate