Move different drainPond calculation options to their own methods

Main point is to decrease the clutter in the main QflxH2osfcDrain
routine - trying to keep the focus more on the science than on the
alternative numerical solutions.

src/biogeophys/SoilHydrologyMod.F90

@@ -59,10 +59,15 @@ module SoilHydrologyMod
   private :: QflxH2osfcDrain     ! Compute qflx_h2osfc_drain
 
   type, extends(computeFlux_type) :: drainPond_type
-     real(r8) :: smoothScale
-     real(r8) :: drainMax
+     private
+     real(r8), public :: drainMax ! maximum drainage rate of ponded water
+     real(r8), public :: dtime    ! model time step
+     real(r8) :: smoothScale = 0.05_r8  ! smoothing scale
    contains
-     procedure :: getFlux => drainPondFlux
+     procedure :: getFlux => drainPondFlux  ! required method to get the current flux estimate
+     procedure :: drainPondExplicitEuler    ! compute the drainPond flux using an explicit euler method
+     procedure :: drainPondImplicitEuler    ! compute the drainPond flux using an implicit euler method
+     procedure :: drainPondAnalytical       ! compute the drainPond flux using an analytical solution
   end type drainPond_type
 
   !-----------------------------------------------------------------------
@@ -595,11 +600,6 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
      type(drainPond_type) :: drainPond_inst
      integer             :: fc, c               
      real(r8)            :: dtime                 ! land model time step (sec)
-     real(r8)            :: h2osfc1               ! h2osfc at the end of the time step, given qflx_h2osfc_drain only
-     real(r8)            :: drainMax              ! maximum drainage rate of ponded water
-     real(r8)            :: const                 ! constant in analytical integral
-     real(r8)            :: uFunc                 ! analytical integral
-     real(r8), parameter :: smoothScale=0.05_r8   ! smoothing scale
 
      character(len=*), parameter :: subname = 'QflxH2osfcDrain'
      !-----------------------------------------------------------------------
@@ -611,6 +611,7 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
      SHR_ASSERT_ALL((ubound(truncate_h2osfc_to_zero) == (/bounds%endc/)), errMsg(sourcefile, __LINE__))
 
      dtime = get_step_size()
+     drainPond_inst%dtime = dtime
 
      do fc = 1, num_hydrologyc
         c = filter_hydrologyc(fc)
@@ -622,32 +623,18 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
 
            ! define maximum drainage
            ! NOTE: check fraction
-           drainMax = frac_h2osfc(c)*qinmax(c)
+           drainPond_inst%drainMax = frac_h2osfc(c)*qinmax(c)
 
            ! switch between different numerical solutions
            select case(ixSolution)
 
-              ! constrained Explicit Euler solution with operator splitting (original)
            case(ixExplicitEuler)
-              qflx_h2osfc_drain(c)=min(drainMax,h2osfc(c)/dtime)
-
-              ! implicit Euler solution with operator splitting
+              call drainPond_inst%drainPondExplicitEuler(h2osfc(c), qflx_h2osfc_drain(c))
            case(ixImplicitEuler)
-              ! The next two lines look like extra overhead, but in a real implementation
-              ! (1) smoothScale could probably be local to the drainPond type; (2)
-              ! drainMax would probably also be local to the drainPond type, and if not,
-              ! it would be an array that is just set once outside a loop, rather than
-              ! being set for each loop iteration.
-              drainPond_inst%smoothScale = smoothScale
-              drainPond_inst%drainMax = drainMax
-              call implicitEuler(drainPond_inst, dtime, h2osfc(c), qflx_h2osfc_drain(c))
-
+              call drainPond_inst%drainPondImplicitEuler(h2osfc(c), qflx_h2osfc_drain(c))
               ! analytical solution with operator splitting
            case(ixAnalytical)
-              const   = 1._r8 - 1._r8/(1._r8 - exp(-h2osfc(c)/smoothScale))
-              uFunc   = -1._r8/(const*exp(drainMax*dtime/smoothScale) - 1._r8)
-              h2osfc1 = -log(1._r8 - uFunc)*smoothScale
-              qflx_h2osfc_drain(c)= (h2osfc1 - h2osfc(c))/dtime
+              call drainPond_inst%drainPondAnalytical(h2osfc(c), qflx_h2osfc_drain(c))
            end select
 
            if(h2osfcflag==0) then
@@ -674,8 +661,43 @@ subroutine drainPondFlux(this, x, f, dfdx)
      if(present(dfdx)) dfdx = -this%drainMax*arg/this%smoothScale
    end subroutine drainPondFlux
 
+   subroutine drainPondExplicitEuler(this, h2osfc, qflx_h2osfc_drain)
+     ! Compute the drainPond flux using an explicit euler method
+     class(drainPond_type), intent(in) :: this
+     real(r8) , intent(in)  :: h2osfc            ! drainPond storage
+     real(r8) , intent(out) :: qflx_h2osfc_drain ! drainage flux
 
-   !-----------------------------------------------------------------------
+     ! constrained Explicit Euler solution with operator splitting (original)
+     qflx_h2osfc_drain = min(this%drainMax, h2osfc/this%dtime)
+   end subroutine drainPondExplicitEuler
+
+   subroutine drainPondImplicitEuler(this, h2osfc, qflx_h2osfc_drain)
+     ! Compute the drainPond flux using an implicit euler method
+     class(drainPond_type), intent(in) :: this
+     real(r8) , intent(in)  :: h2osfc            ! drainPond storage
+     real(r8) , intent(out) :: qflx_h2osfc_drain ! drainage flux
+
+     ! implicit Euler solution with operator splitting
+     call implicitEuler(this, this%dtime, h2osfc, qflx_h2osfc_drain)
+   end subroutine drainPondImplicitEuler
+
+   subroutine drainPondAnalytical(this, h2osfc, qflx_h2osfc_drain)
+     ! Compute the drainPond flux using an analytical solution
+     class(drainPond_type), intent(in) :: this
+     real(r8) , intent(in)  :: h2osfc            ! drainPond storage
+     real(r8) , intent(out) :: qflx_h2osfc_drain ! drainage flux
+
+     real(r8) :: h2osfc1 ! h2osfc at the end of the time step, given qflx_h2osfc_drain only
+     real(r8) :: const   ! constant in analytical integral
+     real(r8) :: uFunc   ! analytical integral
+
+     const   = 1._r8 - 1._r8/(1._r8 - exp(-h2osfc/this%smoothScale))
+     uFunc   = -1._r8/(const*exp(this%drainMax*this%dtime/this%smoothScale) - 1._r8)
+     h2osfc1 = -log(1._r8 - uFunc)*this%smoothScale
+     qflx_h2osfc_drain= (h2osfc1 - h2osfc)/this%dtime
+   end subroutine drainPondAnalytical
+
+     !-----------------------------------------------------------------------
    subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, &
         num_urbanc, filter_urbanc, &
         waterflux_inst, soilhydrology_inst, saturated_excess_runoff_inst, waterstate_inst)