glens_ice_viscosity Function

  function glens_ice_viscosity(this, velocity, temperature, time) &
                                                    result(viscosity)
    !* Author: Christopher MacMackin
    !  Date: April 2017
    !
    ! Calculates the viscosity of ice using Glen's flow law. See the
    ! documentation of the [[glens_law_viscosity]] object for a
    ! description of this parameterisation.
    !
    class(glens_law_viscosity), intent(in) :: this
    class(vector_field), intent(in)        :: velocity
      !! The velocity field of the ice for which the velocity is
      !! being calculated
    real(r8), intent(in)                   :: temperature
      !! The temperature of the ice for which viscosity is being
      !! calculated.
    real(r8), intent(in), optional         :: time
      !! The time at which the viscosity is being calculated. If not
      !! present then assumed to be same as previous value passed.
    class(scalar_field), pointer           :: viscosity
      !! The value of the viscosity
    call velocity%guard_temp()
    if (velocity%dimensions() > 1) then
       error stop ('Multidimensional case not implemented')
       ! This should all work, but there seems to be yet another
       ! compiler bug. At present I don't need the multi-dimensional
       ! case anyway.
       
!      viscosity => velocity%component_d_dx(1, 1)**2 + &
!                   velocity%component_d_dx(2, 2)**2 + &
!                   0.25_r8*(velocity%component_d_dx(1, 2) + &
!                            velocity%component_d_dx(2, 1)**2)
!      viscosity => sqrt(viscosity)
    else
      viscosity => abs(velocity%component_d_dx(1, 1, 1))
    end if
    viscosity => 0.5_r8*this%b_val*viscosity**(1._r8/this%index - 1._r8)
    call velocity%clean_temp()
    call viscosity%set_temp()
  end function glens_ice_viscosity