alternative numerical solutions for drainage of ponded water

src/biogeophys/SoilHydrologyMod.F90

@@ -28,6 +28,10 @@ module SoilHydrologyMod
   use LandunitType      , only : lun                
   use ColumnType        , only : col                
   use PatchType         , only : patch                
+
+  use computeFluxMod                            ! template for the flux module
+  use implicitEulerMod  , only : implicitEuler  ! implicit Euler solver for a single state
+
   !
   ! !PUBLIC TYPES:
   implicit none
@@ -581,8 +585,16 @@ 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:
-     integer :: fc, c
-     real(r8) :: dtime         ! land model time step (sec)
+     procedure(fluxTemplate), pointer :: funcName ! function name
+     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
+     integer(i4b)        :: err                   ! error code
+     character(len=128)  :: message               ! error message
 
      character(len=*), parameter :: subname = 'QflxH2osfcDrain'
      !-----------------------------------------------------------------------
@@ -603,9 +615,31 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, &
            truncate_h2osfc_to_zero(c) = .true.
         else
 
+           ! define maximum drainage
+           ! NOTE: check fraction
+           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
+              case(ixImplicitEuler)
+                 funcName=>drainPond
+                 call implicitEuler(funcName, dtime, h2osfc(c), h2osfc1, err, message)
+                 if(err/=0) call endrun(subname // ':: '//trim(message))
+                 call drainPond(h2osfc1, 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 = -alog(1._r8 - uFunc)*smoothScale
+                 qflx_h2osfc_drain(c)= (h2osfc1 - h2osfc(c))/dtime
 
-           qflx_h2osfc_drain(c)=min(frac_h2osfc(c)*qinmax(c),h2osfc(c)/dtime)
            if(h2osfcflag==0) then
               ! ensure no h2osfc
               qflx_h2osfc_drain(c)= max(0._r8,h2osfc(c)/dtime)
@@ -614,6 +648,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 = -inputRate*(1._dp - arg)
+      if(present(dfdx)) dfdx = -inputRate*arg/smoothScale
+      end subroutine drainPond
+
    end subroutine QflxH2osfcDrain
 
 

src/utils/computeFluxMod.F90

@@ -0,0 +1,28 @@
+module computeFluxMod
+
+ ! provides interface for the flux modules
+
+ USE nrtype
+ implicit none
+
+ abstract interface
+ subroutine fluxTemplate(x,f,dfdx)
+  USE nrtype
+  implicit none
+  real(dp), intent(in)            :: x    ! state
+  real(dp), intent(out)           :: f    ! f(x)
+  real(dp), intent(out), optional :: dfdx ! derivative
+ end subroutine fluxTemplate
+ end interface
+
+ contains
+
+  subroutine getFlux(func, x, f, dfdx)
+   procedure(fluxTemplate), pointer :: func
+   real(dp), intent(in)           :: x    ! state
+   real(dp), intent(out)          :: f    ! f(x)
+   real(dp), intent(out),optional :: dfdx ! derivative
+   call func(x,f,dfdx)
+  end subroutine getFlux
+
+end module computeFluxMod

src/utils/implicitEulerMod.F90

@@ -0,0 +1,49 @@
+module implicitEulerMod
+
+ USE nrtype
+ USE computeFluxMod
+ implicit none
+
+ contains
+
+ subroutine implicitEuler(funcName, dt, xInit, xNew, err, message)
+ ! dummy variables
+ procedure(fluxTemplate), pointer :: funcName
+ real(dp), intent(in)     :: dt            ! time step
+ real(dp), intent(in)     :: xInit         ! initial state
+ real(dp), intent(out)    :: xNew          ! updated state
+ integer(i4b),intent(out) :: err           ! error code
+ character(*),intent(out) :: message       ! error message
+ ! local variables
+ integer(i4b)             :: iter          ! iteration index
+ integer(i4b),parameter   :: maxiter=20    ! maximum number of iterations
+ real(dp), parameter      :: xTol=1.e-8_dp ! convergence tolerance
+ real(dp)                 :: flux          ! flux
+ real(dp)                 :: dfdx          ! derivative
+ real(dp)                 :: xRes          ! residual error
+ real(dp)                 :: delX          ! state update
+ ! initialiuze routine
+ err=0; message='implicitEuler/'
+
+ ! initialize
+ xNew = xInit
+
+ ! iterate
+ do iter=1,maxiter
+  ! get the flux and derivative
+  call getFlux(funcName, xNew, flux, dfdx)
+  ! get the residual and the update
+  xRes = xNew - (xInit + dt*flux)
+  delX = -xRes/(1._dp - dt*dfdx)
+  ! update and check for convergence
+  xNew = xNew + delX
+  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
+  endif
+ end do
+ end subroutine implicitEuler
+
+end module implicitEulerMod