Ponded water/alternative solutions - object-oriented

src/biogeophys/SoilHydrologyMod.F90

@@ -58,6 +58,13 @@ module SoilHydrologyMod
   private :: QflxH2osfcSurf      ! Compute qflx_h2osfc_surf
   private :: QflxH2osfcDrain     ! Compute qflx_h2osfc_drain
 
+  type, extends(flux_type) :: drainPond_type
+     real(r8) :: smoothScale
+     real(r8) :: drainMax
+   contains
+     procedure :: getFlux => drainPondFlux
+  end type drainPond_type
+
   !-----------------------------------------------------------------------
   real(r8), private :: baseflow_scalar = 1.e-2_r8
 
@@ -585,7 +592,7 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
      logical           , intent(inout) :: truncate_h2osfc_to_zero( bounds%begc: ) ! whether h2osfc should be truncated to 0 to correct for roundoff errors, in order to maintain bit-for-bit the same answers as the old code
      !
      ! !LOCAL VARIABLES:
-     procedure(fluxTemplate), pointer :: funcName ! function name
+     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
@@ -628,10 +635,16 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
 
               ! implicit Euler solution with operator splitting
            case(ixImplicitEuler)
-              funcName=>drainPond
-              call implicitEuler(funcName, dtime, h2osfc(c), h2osfc1, err, message)
+              ! 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), h2osfc1, err, message)
               if(err/=0) call endrun(subname // ':: '//trim(message))
-              call drainPond(h2osfc1, qflx_h2osfc_drain(c))
+              call drainPond_inst%getFlux(h2osfc1, qflx_h2osfc_drain(c))
 
               ! analytical solution with operator splitting
            case(ixAnalytical)
@@ -649,24 +662,22 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
         end if
      end do
 
-   contains
-
-     ! model-specific flux routine
-     subroutine drainPond(storage, flux, dfdx)
-       ! dummy variables
-       real(r8), intent(in)           :: storage    ! storage
-       real(r8), intent(out)          :: flux       ! drainage flux
-       real(r8), intent(out),optional :: dfdx       ! derivative
-       ! local variables
-       real(r8)                       :: arg        ! temporary argument
-       ! procedure starts here
-       arg  = exp(-storage/smoothScale)
-       flux = -drainMax*(1._r8 - arg)
-       if(present(dfdx)) dfdx = -drainMax*arg/smoothScale
-     end subroutine drainPond
-
    end subroutine QflxH2osfcDrain
 
+   subroutine drainPondFlux(this, x, f, dfdx)
+     ! dummy variables
+     class(drainPond_type), intent(in) :: this
+     real(r8), intent(in)           :: x    ! storage
+     real(r8), intent(out)          :: f    ! drainage flux
+     real(r8), intent(out),optional :: dfdx ! derivative
+     ! local variables
+     real(r8)                       :: arg        ! temporary argument
+     ! procedure starts here
+     arg  = exp(-x/this%smoothScale)
+     f    = -this%drainMax*(1._r8 - arg)
+     if(present(dfdx)) dfdx = -this%drainMax*arg/this%smoothScale
+   end subroutine drainPondFlux
+
 
    !-----------------------------------------------------------------------
    subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, &

src/utils/computeFluxMod.F90

@@ -5,24 +5,21 @@ module computeFluxMod
   use shr_kind_mod      , only : r8 => shr_kind_r8
   implicit none
 
+  type, abstract :: flux_type
+   contains
+     procedure(getFlux_interface), deferred :: getFlux
+  end type flux_type
+
   abstract interface
-     subroutine fluxTemplate(x,f,dfdx)
+     subroutine getFlux_interface(this,x,f,dfdx)
        use shr_kind_mod      , only : r8 => shr_kind_r8
+       import :: flux_type
        implicit none
+       class(flux_type), intent(in) :: this
        real(r8), intent(in)            :: x    ! state
        real(r8), intent(out)           :: f    ! f(x)
        real(r8), intent(out), optional :: dfdx ! derivative
-     end subroutine fluxTemplate
+     end subroutine getFlux_interface
   end interface
 
-contains
-
-  subroutine getFlux(func, x, f, dfdx)
-    procedure(fluxTemplate), pointer :: func
-    real(r8), intent(in)           :: x    ! state
-    real(r8), intent(out)          :: f    ! f(x)
-    real(r8), intent(out),optional :: dfdx ! derivative
-    call func(x,f,dfdx)
-  end subroutine getFlux
-
 end module computeFluxMod

src/utils/implicitEulerMod.F90

@@ -1,14 +1,14 @@
 module implicitEulerMod
 
   use shr_kind_mod      , only : r8 => shr_kind_r8
-  use computeFluxMod
+  use computeFluxMod    , only : flux_type
   implicit none
 
 contains
 
-  subroutine implicitEuler(funcName, dt, xInit, xNew, err, message)
+  subroutine implicitEuler(flux_inst, dt, xInit, xNew, err, message)
     ! dummy variables
-    procedure(fluxTemplate), pointer :: funcName
+    class(flux_type), intent(in) :: flux_inst
     real(r8), intent(in)     :: dt            ! time step
     real(r8), intent(in)     :: xInit         ! initial state
     real(r8), intent(out)    :: xNew          ! updated state
@@ -31,7 +31,7 @@ subroutine implicitEuler(funcName, dt, xInit, xNew, err, message)
     ! iterate
     do iter=1,maxiter
        ! get the flux and derivative
-       call getFlux(funcName, xNew, flux, dfdx)
+       call flux_inst%getFlux(xNew, flux, dfdx)
        ! get the residual and the update
        xRes = xNew - (xInit + dt*flux)
        delX = -xRes/(1._r8 - dt*dfdx)