Move final flux calculation into implicitEuler

This simplifies the use of implicitEuler. Martyn Clark supported doing
this, since it will be needed in general:
https://github.com/NCAR/clm-ctsm/pull/12#discussion_r137650687

Also, move the endrun clal into implicitEuler, rather than requiring the
caller to check the error code. This simplifies the interface, and works
more smoothly now that the final flux calculation is also in this
routine.

src/biogeophys/SoilHydrologyMod.F90

@@ -600,8 +600,6 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
      real(r8)            :: const                 ! constant in analytical integral
      real(r8)            :: uFunc                 ! analytical integral
      real(r8), parameter :: smoothScale=0.05_r8   ! smoothing scale
-     integer        :: err                   ! error code
-     character(len=128)  :: message               ! error message
 
      character(len=*), parameter :: subname = 'QflxH2osfcDrain'
      !-----------------------------------------------------------------------
@@ -642,9 +640,7 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
               ! being set for each loop iteration.
               drainPond_inst%smoothScale = smoothScale
               drainPond_inst%drainMax = drainMax
-              call implicitEuler(drainPond_inst, dtime, h2osfc(c), h2osfc1, err, message)
-              if(err/=0) call endrun(subname // ':: '//trim(message))
-              call drainPond_inst%getFlux(h2osfc1, qflx_h2osfc_drain(c))
+              call implicitEuler(drainPond_inst, dtime, h2osfc(c), qflx_h2osfc_drain(c))
 
               ! analytical solution with operator splitting
            case(ixAnalytical)

src/utils/computeFluxMod.F90

@@ -6,6 +6,11 @@ module computeFluxMod
   implicit none
 
   type, abstract :: computeFlux_type
+     ! flux_name is just used to generate more helpful error messages. This needs to be
+     ! set in initialization. (If that's awkward to accomplish, then we could have a
+     ! deferred procedure getFluxName that simply returns a character string giving the
+     ! name of this flux.)
+     character(len=128) :: flux_name
    contains
      procedure(getFlux_interface), deferred :: getFlux
   end type computeFlux_type

src/utils/implicitEulerMod.F90

@@ -2,28 +2,25 @@ module implicitEulerMod
 
   use shr_kind_mod      , only : r8 => shr_kind_r8
   use computeFluxMod    , only : computeFlux_type
+  use abortutils        , only : endrun
   implicit none
 
 contains
 
-  subroutine implicitEuler(computeFlux_inst, dt, xInit, xNew, err, message)
+  subroutine implicitEuler(computeFlux_inst, dt, xInit, flux)
     ! dummy variables
     class(computeFlux_type), intent(in) :: computeFlux_inst
     real(r8), intent(in)     :: dt            ! time step
     real(r8), intent(in)     :: xInit         ! initial state
-    real(r8), intent(out)    :: xNew          ! updated state
-    integer,intent(out) :: err           ! error code
-    character(*),intent(out) :: message       ! error message
+    real(r8), intent(out)    :: flux          ! flux
     ! local variables
     integer             :: iter          ! iteration index
     integer,parameter   :: maxiter=20    ! maximum number of iterations
     real(r8), parameter      :: xTol=1.e-8_r8 ! convergence tolerance
-    real(r8)                 :: flux          ! flux
+    real(r8)                 :: xNew          ! updated state
     real(r8)                 :: dfdx          ! derivative
     real(r8)                 :: xRes          ! residual error
     real(r8)                 :: delX          ! state update
-    ! initialiuze routine
-    err=0; message='implicitEuler/'
 
     ! initialize
     xNew = xInit
@@ -40,10 +37,12 @@ subroutine implicitEuler(computeFlux_inst, dt, xInit, xNew, err, message)
        if(abs(delX) < xTol) exit
        ! check for non-convergence
        if(iter==maxiter)then
-          message=trim(message)//'the implicit Euler solution did not converge!'
-          err=20; return
+          call endrun(computeFlux_inst%flux_name // ': the implicit Euler solution did not converge!')
        endif
     end do
+
+    ! Get final flux
+    call computeFlux_inst%getFlux(xNew, flux)
   end subroutine implicitEuler
 
 end module implicitEulerMod