jenkins1991_rate Function

  function jenkins1991_rate(this, velocity, thickness, depth, density_diff, time) &
                                                result(entrainment)
    !* Author: Christopher MacMackin
    !  Date: October 2016
    !
    ! $$e = E_0 |\vec{U}\sin(\theta) \simeq E_0|\vec{U}||\nabla b|$$
    ! Here, \(E_0\) is a coefficient typically taken to be 0.036 (the
    ! default value), \(\vec{U}\) is the velocity of the plume, \(\theta\)
    ! is the angle of slope of the ice shelf base, and \(b\) is the
    ! basal depth of the ice shelf.
    !
    ! @Warning
    ! The calculation must be performed as
    ! ```fortran
    ! this%coefficient * depth%d_dx(1) * velocity%norm()
    ! ```
    ! with the variables in a different order than how the equation is
    ! usually formulated. If they are in the correct order then
    ! `gfortran` expects the result to be a `vector_field`. It
    ! is not clear whether this is due to a bug in `gfortran` or in
    ! `factual`.
    !
    class(jenkins1991_entrainment), intent(in) :: this
    class(vector_field), intent(in)            :: velocity
      !! The velocity field of the plume into which fluid is being 
      !! entrained.
    class(scalar_field), intent(in)            :: thickness
      !! The thickness of the plume into which fluid is being
      !! entrained
    class(scalar_field), intent(in)            :: depth
      !! The depth of the upper surface of the plume into which
      !! fluid is being entrained
    class(scalar_field), intent(in)            :: density_diff
      !! The difference between the ambient density and the density of
      !! the plume into which the ambient fluid is being entrained.
    real(r8), intent(in), optional             :: time
      !! The time at which the entrainment is being calculated. If not
      !! present then assumed to be same as previous value passed.
    class(scalar_field), pointer               :: entrainment
      !! The value of the entrainment
    class(vector_field), pointer               :: tmp
    call velocity%guard_temp(); call thickness%guard_temp()
    call depth%guard_temp(); call density_diff%guard_temp()
    call depth%allocate_scalar_field(entrainment)
    call depth%allocate_vector_field(tmp)
    entrainment = velocity%norm()
    call entrainment%unset_temp()
    tmp = .grad. depth
    entrainment = tmp%norm() ! Needed due to ICE when try to put all on one line. TODO: Create minimal example and submit bug report.
    entrainment = this%coefficient * entrainment * velocity%norm()
    call velocity%clean_temp(); call thickness%clean_temp()
    call depth%clean_temp(); call density_diff%clean_temp()
    call entrainment%set_temp()
  end function jenkins1991_rate