integrate Subroutine

  subroutine integrate(this,time)
    !* Author: Christopher MacMackin
    !  Date: April 2016
    !
    ! Integrates the cryosphere forward until the specified `time` is
    ! reached.
    !
    class(cryosphere), intent(inout) :: this
    real(r8), intent(in)             :: time
      !! The time to which to integrate the cryosphere
    class(glacier), dimension(:), allocatable :: old_glaciers
    logical :: success, past_fail
    real(r8) :: t, old_t, dt
    real(r8), allocatable, dimension(:) :: sub_state
past_fail = .false.
    if (time <= this%time) then
      call logger%warning('cryosphere%integrate','Request made to '// &
                          'integrate cryosphere to earlier time '//  &
                          'than present state. No action taken.')
      return
    end if

    ! Normally the plume should be solved at the end of the previous
    ! iteration, but in the first iteration obviously there hasn't
    ! been a chance for this to happen yet.
    if (this%first_integration) then
      ! As I am integrating only semi-implicitly and solving the plume
      ! for the current (rather than future) state, I think I should
      ! pass the current time. I only *think* that this is correct,
      ! however.
      call this%sub_ice%solve(this%ice%ice_thickness(), this%ice%ice_density(), &
                              this%ice%ice_temperature(), this%time, success)
      this%first_integration = .false.
      if (.not. success) then
        call logger%fatal('cryosphere%integrate','Failed to solve plume '// &
                          'with initial ice configuration. Writing '//      &
                          'cryosphere state to file "'//hdf_crash_file//'".')
        call this%write_data(hdf_crash_file)
        error stop
      end if
    end if
    allocate(old_glaciers(1), mold=this%ice)

    old_t = this%time
    dt = this%time_step()
    if (dt < 0.5_r8*(time - this%time)) then
      t = old_t + dt
    else if (t + dt > time) then
      t = time
    else
      t = 0.5_r8*(time + old_t)
    end if
    do while (t <= time)
      old_glaciers(1) = this%ice
      call this%ice%integrate(old_glaciers, this%sub_ice%basal_melt(), &
                              this%sub_ice%basal_drag_parameter(),     &
                              this%sub_ice%water_density(), t, success)
        if (past_fail) then
           call this%write_data('isoft_post_failure.h5')
           past_fail = .false.
        end if
      if (.not. success) then
        this%ice = old_glaciers(1)
        if (this%dt_factor > this%min_dt_factor) then
          call logger%warning('cryosphere%integrate','Failure in nonlinear '// &
                              'solver. Reducing time step and trying again.')
          call this%reduce_time_step()
          dt = this%time_step()
          if (dt < 0.5_r8*(time - old_t)) then
            t = old_t + dt
          else if (old_t + dt > time) then
            t = time
          else
            t = 0.5_r8*(time + old_t)
          end if
          cycle
        else
          call logger%fatal('cryosphere%integrate','Failed to integrate '//  &
                            'glacier to time '//trim(str(t))//'! Writing '// &
                            'cryosphere state to file "'//hdf_crash_file//'".')
          call this%write_data(hdf_crash_file)
          iplvl = 2
          t = min(old_t + this%ice%time_step(), time, 0.5_r8*(time + old_t))
          call this%ice%integrate(old_glaciers, this%sub_ice%basal_melt(), &
                                  this%sub_ice%basal_drag_parameter(),     &
                                  this%sub_ice%water_density(), t, success)
          error stop
        end if
      end if

      ! Solve the plume so that it is ready for use in the next step of
      ! the time integration.
      sub_state = this%sub_ice%state_vector()
      call this%sub_ice%solve(this%ice%ice_thickness(), this%ice%ice_density(), &
                              this%ice%ice_temperature(), t, success)

      if (success) then
        call logger%trivia('cryosphere%integrate','Successfully integrated '// &
                           'cryosphere to time '//trim(str(t)))
      else
          call this%write_data('isoft_failed_state.h5')
        this%ice = old_glaciers(1)
        call this%sub_ice%update(sub_state, this%ice%ice_thickness())
        if (this%dt_factor > this%min_dt_factor) then
          call logger%warning('cryosphere%integrate','Failure in plume '// &
                              'solver. Reducing time step and trying again.')
          call this%reduce_time_step()
          dt = this%time_step()
          if (dt < 0.5_r8*(time - old_t)) then
            t = old_t + dt
          else if (old_t + dt > time) then
            t = time
          else
            t = 0.5_r8*(time + old_t)
          end if
          call this%write_data('isoft_pre_failure.h5')
          past_fail = .true.
          cycle
        else
          call logger%fatal('cryosphere%integrate','Failed to solve plume '//   &
                            'at time '//trim(str(t))//'! Writing cryosphere '// &
                            'state to file "'//hdf_crash_file//'".')
          call this%write_data(hdf_crash_file)
          error stop
        end if
      end if

      call this%increase_time_step()
      if (t >= time) exit
      old_t = t
      dt = this%time_step()
      if (dt < 0.5_r8*(time - t)) then
        t = t + dt
      else if (t + dt > time) then
        t = time
      else
        t = 0.5_r8*(time + t)
      end if
      this%time = old_t
    end do

    this%time = time
    call logger%info('cryosphere%integrate','Successfully integrated '// &
                     'cryosphere to time '//trim(str(t)))
  end subroutine integrate