kochergin1987_rate Function

  function kochergin1987_rate(this, velocity, thickness, depth, density_diff, time) &
                                                result(entrainment)
    !* Author: Christopher MacMackin
    !  Date: Feburary 2018
    !
    ! $$e = \frac{c_L^2}{S_m}\sqrt{|\vec{U}|^2+\frac{g'D}{S_m}}.$$
    ! Here, \(c_L\) is an entrainment coefficient, \(\vec{U}\) is the
    ! velocity of the plume, \(g'\) is the reduced gravity, and
    ! \(S_m\) is the turbulent Schmidt number. The Schmidt number is a
    ! function of the Richardson number \(Ri = g'D/|\vec{U}|^2\):
    ! $$ S_m = \frac{Ri}{0.0725(Ri + 0.186 -
    ! \sqrt{Ri^2 - 0.316Ri + 0.0346})}. $$
    !
    class(kochergin1987_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(scalar_field), pointer               :: Ri, Sm
    call velocity%guard_temp(); call thickness%guard_temp() 
    call depth%guard_temp(); call density_diff%guard_temp()
    call depth%allocate_scalar_field(entrainment)
    call entrainment%unset_temp()
    call depth%allocate_scalar_field(Ri)
    call depth%allocate_scalar_field(Sm)
    call Ri%guard_temp(); call Sm%guard_temp()
    entrainment = velocity%norm() ! Have to awkwardly split this operation to prevent ICE
    Ri = this%delta*density_diff*thickness/(entrainment**2)
    Sm = Ri/(0.0725_r8*(Ri + 0.186_r8 - sqrt(Ri**2 - 0.316_r8*Ri + 0.0346_r8)))
    entrainment = this%coefficient*entrainment/Sm * sqrt(1._r8 + Ri/Sm)
    call velocity%clean_temp(); call thickness%clean_temp()
    call depth%clean_temp(); call density_diff%clean_temp()
    call Ri%clean_temp(); call Sm%clean_temp()
    call entrainment%set_temp()
  end function kochergin1987_rate