anchor.diff:  
       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
   
  
positive.code:             !! Update the domain for a time-step from (time-dt, time). Note "time" is the END of the timestep
          real(dp), intent(in) :: start_time, end_time
              !! Apply the entire forcing at a steady rate between start_time and end_time
  
          real(dp) :: update_fraction, local_dt
          integer(ip) :: i, j, k
  
          if( ((time - dt) <= end_time) .and. (time >= start_time)) then
              ! The forcing is only applied between start_time and end_time - so
              ! we need to make sure the current time-step overlapped that period.
              ! (Note we are at the end of the timestep, and time was already updated).
  
negative.code:         integer(ip), parameter :: global_nx(2) = [50, 30], local_nx(2) = [23, 90]
  
          type(domain_type), allocatable :: domains(:)
          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
  
          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
  
              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
  
      end subroutine
  
  end module
  