anchor.diff:          integer(ip), parameter :: global_nx(2) = [50, 30], local_nx(2) = [23, 90]
   
           type(domain_type), allocatable :: domains(:)
  -        integer(ip) :: i, j
  +        integer(ip) :: i, j, provide_function
           real(dp) :: err_tol, err, time, dt
           type(forcing_patch_type), pointer :: domain_forcing_context
           integer(ip), parameter :: nd = 2
  +        procedure(forcing_time), pointer :: forcing_time_function
   
           err_tol = spacing(10.0_dp)*10
   
  -        allocate(domains(nd))
  -
  -        ! Setup the main test domain
  -        call setup_domain_for_test(domains(1), global_lw, global_nx, global_ll)
  -        ! Backup the initial condition for later usage
  -        domains(1)%backup_U = domains(1)%U
  -
  -        ! Setup another test domain
  -        call setup_domain_for_test(domains(2), local_lw, local_nx, local_ll)
  -        domains(2)%backup_U = domains(2)%U
  -
  -        ! Suppose we just evolved from 0.0 to 10.0
  -        time = 10.0_dp
  -        dt = 10.0_dp
  -
  -        !
  -        ! Tests of add_forcing_work_to_U_over_time
  -        !
  -
  -        ! Test 1 -- time-step exactly equals range of forcing
  -        do j = 1, nd
  -            call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  -                ! unpack the forcing_context, which is stored as a c_ptr in the domain
  -            call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  -                time, dt, start_time = 0.0_dp, end_time = 10.0_dp)
  -            err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work)))
  -            if(err < err_tol) then
  -                print*, 'PASS'
  -            else
  -                print*, 'FAIL', err, err_tol, __LINE__
  +        do provide_function = 1, 2
  +            ! Loop to test with/without forcing function
  +
  +            if(provide_function == 1) then
  +                ! First time we set forcing_time_function => NULL,
  +                forcing_time_function => NULL()
  +            else if(provide_function == 2) then
  +                ! Second time we set it to identity_function
  +                ! While the result should be the same, different code paths are exercised
  +                forcing_time_function => identity_function
               end if
  -            domains(j)%U = domains(j)%backup_U
  +
  +            allocate(domains(nd))
  +
  +            ! Setup the main test domain
  +            call setup_domain_for_test(domains(1), global_lw, global_nx, global_ll)
  +            ! Backup the initial condition for later usage
  +            domains(1)%backup_U = domains(1)%U
  +
  +            ! Setup another test domain
  +            call setup_domain_for_test(domains(2), local_lw, local_nx, local_ll)
  +            domains(2)%backup_U = domains(2)%U
  +
  +            ! Suppose we just evolved from 0.0 to 10.0
  +            time = 10.0_dp
  +            dt = 10.0_dp
  +
  +            !
  +            ! Tests of add_forcing_work_to_U_over_time
  +            !
  +
  +            ! Test 1 -- time-step exactly equals range of forcing, { i.e. dt == (end_time - start_time) }
  +            do j = 1, nd
  +                call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  +                    ! unpack the forcing_context, which is stored as a c_ptr in the domain
  +
  +                call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  +                    time, dt, start_time = 0.0_dp, end_time = 10.0_dp, forcing_time_function=forcing_time_function)
  +                err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work)))
  +                if(err < err_tol) then
  +                    print*, 'PASS'
  +                else
  +                    print*, 'FAIL', err, err_tol, __LINE__
  +                end if
  +                domains(j)%U = domains(j)%backup_U
  +            end do
  +
  +
  +            ! Test 2 -- time-step starts before range of forcing {i.e. time-dt < start_time }
  +            do j = 1, nd
  +                call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  +                call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  +                    time, dt, start_time = 4.0_dp, end_time = 10.0_dp, forcing_time_function=forcing_time_function)
  +                err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work)))
  +                if(err < err_tol ) then
  +                    print*, 'PASS'
  +                else
  +                    print*, 'FAIL', err, err_tol, __LINE__
  +                end if
  +                domains(j)%U = domains(j)%backup_U
  +            end do
  +
  +            ! Test 3 -- time-step starts after range of forcing {i.e. time - dt > start_time}
  +            do j = 1, nd
  +                call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  +                call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  +                    time, dt, start_time = -4.0_dp, end_time = 10.0_dp, forcing_time_function=forcing_time_function)
  +                err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/14.0_dp)))
  +                if(err < err_tol) then
  +                    print*, 'PASS'
  +                else
  +                    print*, 'FAIL', err, err_tol, __LINE__, &
  +                       minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/14.0_dp)), &
  +                       maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/14.0_dp))
  +                end if
  +                domains(j)%U = domains(j)%backup_U
  +            end do
  +
  +            ! Test 4 -- time-step fully inside the forcing time
  +            do j = 1, nd
  +                call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  +                call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  +                    time, dt, start_time = -15.0_dp, end_time = 15.0_dp, forcing_time_function=forcing_time_function)
  +                err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/30.0_dp)))
  +                if(err < err_tol) then
  +                    print*, 'PASS'
  +                else
  +                    print*, 'FAIL', err, err_tol, __LINE__, &
  +                       minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/30.0_dp)), &
  +                       maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/30.0_dp))
  +                end if
  +                domains(j)%U = domains(j)%backup_U
  +            end do
  +
  +            ! Test 5 -- time-step overshoots the forcing end time
  +            do j = 1, nd
  +                call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  +                call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work,&
  +                    time, dt, start_time = -15.0_dp, end_time = 5.0_dp, forcing_time_function=forcing_time_function)
  +                err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 5.0_dp/20.0_dp)))
  +                if(err < err_tol) then
  +                    print*, 'PASS'
  +                else
  +                    print*, 'FAIL', err, err_tol, __LINE__, &
  +                       minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 5.0_dp/20.0_dp)), &
  +                       maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 5.0_dp/20.0_dp))
  +                end if
  +                domains(j)%U = domains(j)%backup_U
  +            end do
  +
  +            ! Test 6 -- time-step doesn't overlap with the forcing
  +            do j = 1, nd
  +                call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  +                call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  +                    time, dt, start_time = -15.0_dp, end_time = -5.0_dp, forcing_time_function=forcing_time_function)
  +                err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 0.0_dp)))
  +                if(err == 0.0_dp) then
  +                    print*, 'PASS'
  +                else
  +                    print*, 'FAIL', err, err_tol, __LINE__, &
  +                       minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 0.0_dp)), &
  +                       maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 0.0_dp))
  +                end if
  +                domains(j)%U = domains(j)%backup_U
  +            end do
  +
  +            ! Test 7 -- a timestepping type test, where we should end up adding the full contribution of forcing_work
  +            ! Use some erratic numbers -- it shouldn't matter so long as we fully step through start_time to end_time
  +            do j = 1, nd
  +                call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  +                time = -10.01242_dp
  +                dt = 1.4444_dp
  +                do i = 1, 30
  +                    ! When the loop is complete we will have covered the time
  +                    call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  +                        time, dt, start_time = -1.332_dp, end_time = 5.15430_dp, forcing_time_function=forcing_time_function)
  +                    time = time + dt
  +                end do
  +                err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 1.0_dp)))
  +                if(err < err_tol) then
  +                    print*, 'PASS'
  +                else
  +                    print*, 'FAIL', err, err_tol, __LINE__, &
  +                       minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work)), &
  +                       maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work))
  +                end if
  +                domains(j)%U = domains(j)%backup_U
  +            end do
  +
  +            do j = 1, nd
  +                call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  +                deallocate(domain_forcing_context%forcing_work)
  +                deallocate(domain_forcing_context)
  +            end do
  +
  +            deallocate(domains)
  +
           end do
   
  -
  -        ! Test 2 -- time-step starts before range of forcing
  -        do j = 1, nd
  -            call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  -            call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  -                time, dt, start_time = 4.0_dp, end_time = 10.0_dp)
  -            err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work)))
  -            if(err < err_tol ) then
  -                print*, 'PASS'
  -            else
  -                print*, 'FAIL', err, err_tol, __LINE__
  -            end if
  -            domains(j)%U = domains(j)%backup_U
  -        end do
  -
  -        ! Test 3 -- time-step starts after range of forcing
  -        do j = 1, nd
  -            call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  -            call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  -                time, dt, start_time = -4.0_dp, end_time = 10.0_dp)
  -            err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/14.0_dp)))
  -            if(err < err_tol) then
  -                print*, 'PASS'
  -            else
  -                print*, 'FAIL', err, err_tol, __LINE__, &
  -                   minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/14.0_dp)), &
  -                   maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/14.0_dp))
  -            end if
  -            domains(j)%U = domains(j)%backup_U
  -        end do
  -
  -        ! Test 4 -- time-step fully inside the forcing time
  -        do j = 1, nd
  -            call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  -            call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  -                time, dt, start_time = -15.0_dp, end_time = 15.0_dp)
  -            err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/30.0_dp)))
  -            if(err < err_tol) then
  -                print*, 'PASS'
  -            else
  -                print*, 'FAIL', err, err_tol, __LINE__, &
  -                   minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/30.0_dp)), &
  -                   maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 10.0_dp/30.0_dp))
  -            end if
  -            domains(j)%U = domains(j)%backup_U
  -        end do
  -
  -        ! Test 5 -- time-step overshoots the forcing end time
  -        do j = 1, nd
  -            call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  -            call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work,&
  -                time, dt, start_time = -15.0_dp, end_time = 5.0_dp)
  -            err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 5.0_dp/20.0_dp)))
  -            if(err < err_tol) then
  -                print*, 'PASS'
  -            else
  -                print*, 'FAIL', err, err_tol, __LINE__, &
  -                   minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 5.0_dp/20.0_dp)), &
  -                   maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 5.0_dp/20.0_dp))
  -            end if
  -            domains(j)%U = domains(j)%backup_U
  -        end do
  -
  -        ! Test 6 -- time-step doesn't overlap with the forcing
  -        do j = 1, nd
  -            call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  -            call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  -                time, dt, start_time = -15.0_dp, end_time = -5.0_dp)
  -            err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 0.0_dp)))
  -            if(err == 0.0_dp) then
  -                print*, 'PASS'
  -            else
  -                print*, 'FAIL', err, err_tol, __LINE__, &
  -                   minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 0.0_dp)), &
  -                   maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 0.0_dp))
  -            end if
  -            domains(j)%U = domains(j)%backup_U
  -        end do
  -
  -        ! Test 7 -- a timestepping type test, where we should end up adding the full contribution of forcing_work
  -        ! Use some erratic numbers -- it shouldn't matter so long as we fully step through start_time to end_time
  -        do j = 1, nd
  -            call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  -            time = -10.01242_dp
  -            dt = 1.4444_dp
  -            do i = 1, 30
  -                ! When the loop is complete we will have covered the time
  -                call add_forcing_work_to_U_over_time(domains(j)%U, domain_forcing_context%forcing_work, &
  -                    time, dt, start_time = -1.332_dp, end_time = 5.15430_dp)
  -                time = time + dt
  -            end do
  -            err = maxval(abs(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work * 1.0_dp)))
  -            if(err < err_tol) then
  -                print*, 'PASS'
  -            else
  -                print*, 'FAIL', err, err_tol, __LINE__, &
  -                   minval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work)), &
  -                   maxval(domains(j)%U - (domains(j)%backup_U + domain_forcing_context%forcing_work))
  -            end if
  -            domains(j)%U = domains(j)%backup_U
  -        end do
  -
  -        do j = 1, nd
  -            call c_f_pointer(domains(j)%forcing_context_cptr, domain_forcing_context)
  -            deallocate(domain_forcing_context%forcing_work)
  -            deallocate(domain_forcing_context)
  -        end do
  -
       end subroutine
   
   end module
  
positive.code:                 call generic_stop
              end if
  
              ! !This if-statement must always be .FALSE. given the outer if
              !if(start_time > end_time) then
              !    write(log_output_unit, *) &
              !        'Error: start_time must be > end_time in add_forcing_work_to_U_over_time'
              !    call generic_stop
              !end if
  
              ! Determine what fraction of forcing_work should be added to U
              if(start_time == end_time) then
                  update_fraction = 1.0_dp
              else
                  ! The time-step might not be completely between start_time and end_time.
                  ! This is the fraction of the forcing we actually need to apply
                  local_dt = min(time, end_time) - max((time-dt), start_time)
                  update_fraction = local_dt/(end_time - start_time)
              end if
              !print*, update_fraction
  
  
negative.code: 
      end subroutine
  
      function identity_function(x) result(p)
          ! Useful for testing
          real(dp), intent(in) :: x
          real(dp) :: p
          p = x
      end function
  
      subroutine test_forcing_mod
          !! Test the forcing subroutines by setting up 2 domains, and checking the forcing works
  
  