MethodModel.f90

Go to the documentation of this file.

module methodmodelmodule
  use kindmodule, only: dp, i4b
  use constantsmodule, only: dzero, done
  use methodmodule, only: methodtype, level_model
  use particlemodule, only: particletype
  use celldefnmodule, only: celldefntype, saturation_dry, &
                            saturation_watertable, &
                            saturation_saturated
  use particleeventmodule, only: particleeventtype
  use mathutilmodule, only: is_close
  use errorutilmodule, only: pstop
 
  private
  public :: methodmodeltype
 
  type, abstract, extends(methodtype) :: methodmodeltype
  contains
    procedure, public :: assess
    procedure, public :: get_level
    ! cell load utilities
    procedure :: cap_cell_wt_flow
    procedure :: load_cell_no_exit_face
    procedure :: load_cell_saturation_status
  end type methodmodeltype
 
contains
 
  !> @brief Check particle reporting/termination status
  subroutine assess(this, particle, cell_defn, tmax)
    ! dummy
    class(methodmodeltype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    type(celldefntype), pointer, intent(inout) :: cell_defn
    real(DP), intent(in) :: tmax
    ! noop
  end subroutine assess
 
  !> @brief Get the model method's level.
  function get_level(this) result(level)
    class(methodmodeltype), intent(in) :: this
    integer(I4B) :: level
    level = level_model
  end function get_level
 
  !> @brief Prevent non-boundary upwards flow at the water table.
  !!
  !! Unless the top face is an assigned boundary with outflow,
  !! cells which contain a water table should not have upward
  !! flow through the top (i.e., the water table). Prevent it
  !! by capping the top face flow at zero in these conditions.
  !!
  !! Assumes cell properties and flows are already loaded.
  !<
  subroutine cap_cell_wt_flow(this, defn)
    class(methodmodeltype), intent(inout) :: this
    type(celldefntype), pointer, intent(inout) :: defn
    ! local
    integer(I4B) :: itopface
 
    ! If the cell contains a water table that is not an
    ! assigned boundary face with upflow, cap flow at 0.
    itopface = this%fmi%max_faces ! fmi's lateral face indices are not closed
    if (this%fmi%is_boundary_face(defn%icell, itopface)) return
    if (defn%isatstat == saturation_watertable) then
      itopface = defn%npolyverts + 3 ! cell defn's lateral face indices are closed
      defn%faceflow(itopface) = max(dzero, defn%faceflow(itopface))
    end if
 
  end subroutine cap_cell_wt_flow
 
  !> @brief Set flag indicating if the cell has any faces with outflow.
  !! Assumes cell properties and flows are already loaded.
  subroutine load_cell_no_exit_face(this, defn)
    ! dummy
    class(methodmodeltype), intent(inout) :: this
    type(celldefntype), pointer, intent(inout) :: defn
    ! local
    integer(I4B) :: m, nfaces
 
    defn%inoexitface = 1
    nfaces = defn%npolyverts + 3
    do m = 1, nfaces
      if (defn%faceflow(m) < dzero) defn%inoexitface = 0
    end do
 
  end subroutine load_cell_no_exit_face
 
  !> @brief Set the saturation status of a cell.
  !! See the status enumeration in CellDefnModule.
  subroutine load_cell_saturation_status(this, defn)
    ! dummy
    class(methodmodeltype), intent(inout) :: this
    type(celldefntype), pointer, intent(inout) :: defn
    ! local
    integer(I4B) :: ic, ictopnbr, itopface, idiag, ipos
 
    ic = defn%icell
    defn%isatstat = saturation_saturated
 
    ! dry?
    if (this%fmi%ibdgwfsat0(ic) == 0) then
      defn%isatstat = saturation_dry
      return
    end if
 
    ! partially saturated?
    if (this%fmi%gwfsat(ic) < done) then
      defn%isatstat = saturation_watertable
      return
    end if
 
    ! no top neighbor?
    itopface = defn%npolyverts + 3 ! cell defn's lateral face indices are closed
    if (defn%facenbr(itopface) == 0) then
      defn%isatstat = saturation_saturated
      return
    end if
 
    ! dry top neighbor?
    idiag = this%fmi%dis%con%ia(ic)
    ipos = idiag + defn%facenbr(itopface)
    ictopnbr = this%fmi%dis%con%ja(ipos)
    if (this%fmi%ibdgwfsat0(ictopnbr) == 0) then
      defn%isatstat = saturation_watertable
      return
    end if
 
  end subroutine load_cell_saturation_status
 
end module methodmodelmodule