ice_shelf_integrate_layers Subroutine

  subroutine ice_shelf_integrate_layers(this, old_states, time, success)
    !* Author: Chris MacMackin
    !  Date: September 2018
    !
    ! Integrate the Taylor coefficients representing the vertical
    ! structure of internal reflectors forward to the specified
    ! time. This is done using an implicit method, with the resulting
    ! linear system solved using GMRES.
    !
    class(ice_shelf), intent(inout)          :: this
    class(glacier), dimension(:), intent(in) :: old_states
      !! Previous states of the ice_shelf, with the most recent one
      !! first.
    real(r8), intent(in)                     :: time
      !! The time to which the ice_shelf should be integrated
    logical, intent(out)                     :: success
      !! True if the integration is successful, false otherwise

    integer :: n, flag
    real(r8), dimension(:), allocatable :: solution
    real(r8) :: dt, resid
    type(jacobian_block) :: precond_block

    if (.not. allocated(this%kappa)) return
    select type(old_states)
    class is(ice_shelf)
      dt = time - old_states(1)%time
      do n = 1, size(this%kappa)
        solution = this%kappa(n)%raw()
        precond_block = jacobian_block(dt*this%velocity%component(1), 1, &
                                       boundary_locs=[this%thickness_size], &
                                       boundary_types=[dirichlet], &
                                       coef=real(-n, r8)) + 1._r8
        call gmres_solve(solution, operator, old_states(1)%kappa(n)%raw(), resid, &
                         flag, precond=preconditioner, krylov_dim=40)
        call this%kappa(n)%set_from_raw(solution)
        if (flag /= 0) then
          call logger%error('ice_shelf%integrate_layers','GMRES solver '// &
                            'returned with error code '//str(flag))
          success = .false.
          return
        end if
        call this%kappa(n)%set_from_raw(solution)
      end do
    class default
      call logger%fatal('ice_shelf%integrate_layers','Type other than `ice_shelf` '// &
                        'passed to `ice_shelf` object as a previous state.')
      error stop
    end select

  contains

    function operator(v, xcur, rhs, rpar, ipar, success)
      real(r8), dimension(:), intent(in)    :: v
        !! The vector to be multiplied
      real(r8), dimension(:), intent(in)    :: xcur
        !! Array containing the current estimate of the independent
        !! variables in the linear system. This may not be needed, but
        !! is provided just in case.
      real(r8), dimension(:), intent(in)    :: rhs
        !! Array containing the right hand side of the linear
        !! system. This may not be needed, but is provided just in
        !! case.
      real(r8), dimension(*), intent(inout) :: rpar
        !! Parameter/work array 
      integer, dimension(*), intent(inout)  :: ipar
        !! Parameter/work array
      logical, intent(out)                  :: success
        !! Indicates whether operation was completed succesfully
      real(r8), dimension(size(xcur))       :: operator
        !! Result of the operation

      type(cheb1d_scalar_field) :: kappa
      class(scalar_field), pointer :: tmp
      
      call kappa%assign_meta_data(this%kappa(n))
      call kappa%set_from_raw(v)
      tmp => dt*this%velocity .dot. (.grad. kappa)
      kappa = (1._r8 - n*dt*(.div.this%velocity))*kappa + tmp
      operator = kappa%raw()
      operator(this%thickness_size) = v(this%thickness_size)
      success = .true.
    end function operator


    function preconditioner(v, xcur, rhs, rpar, ipar, success)
      real(r8), dimension(:), intent(in)    :: v
        !! The vector to be multiplied
      real(r8), dimension(:), intent(in)    :: xcur
        !! Array containing the current estimate of the independent
        !! variables in the linear system. This may not be needed, but
        !! is provided just in case.
      real(r8), dimension(:), intent(in)    :: rhs
        !! Array containing the right hand side of the linear
        !! system. This may not be needed, but is provided just in
        !! case.
      real(r8), dimension(*), intent(inout) :: rpar
        !! Parameter/work array 
      integer, dimension(*), intent(inout)  :: ipar
        !! Parameter/work array
      logical, intent(out)                  :: success
        !! Indicates whether operation was completed succesfully
      real(r8), dimension(size(xcur))       :: preconditioner
        !! Result of the operation

      type(cheb1d_scalar_field) :: tmp
      
      call tmp%assign_meta_data(this%kappa(n))
      call tmp%set_from_raw(v)
      tmp = precond_block%solve_for(tmp)
      preconditioner = tmp%raw()
      success = .true.
    end function preconditioner
  end subroutine ice_shelf_integrate_layers