tspadvmodule Module Reference

Data Types
type	tspadvtype

Functions/Subroutines
subroutine, public	adv_cr (advobj, name_model, inunit, iout, fmi, eqnsclfac)
	@ brief Create a new ADV object More...
subroutine	adv_df (this, adv_options)
	Define ADV object. More...
subroutine	adv_ar (this, dis, ibound)
	Allocate and read method for package. More...
subroutine	adv_dt (this, dtmax, msg, thetam)
	Calculate maximum time step length. More...
subroutine	adv_fc (this, nodes, matrix_sln, idxglo, cnew, rhs)
	Fill coefficient method for ADV package. More...
subroutine	advtvd (this, n, cnew, rhs)
	Calculate TVD. More...
real(dp) function	advqtvd (this, n, m, iposnm, cnew)
	Calculate TVD. More...
subroutine	adv_cq (this, cnew, flowja)
	Calculate advection contribution to flowja. More...
subroutine	advtvd_bd (this, cnew, flowja)
	Add TVD contribution to flowja. More...
subroutine	adv_da (this)
	Deallocate memory. More...
subroutine	allocate_scalars (this)
	Allocate scalars specific to the streamflow energy transport (SFE) package. More...
subroutine	read_options (this)
	Read options. More...
real(dp) function	adv_weight (this, iadvwt, ipos, n, m, qnm)
	@ brief Advection weight More...

Function/Subroutine Documentation

adv_ar()

subroutine tspadvmodule::adv_ar ( class(tspadvtype)  this, class(disbasetype), intent(in), pointer  dis, integer(i4b), dimension(:), intent(in), pointer, contiguous  ibound  )

private

Method to allocate and read static data for the ADV package.

Definition at line 110 of file tsp-adv.f90.

    ! -- modules
    ! -- dummy
    class(TspAdvType) :: this
    class(DisBaseType), pointer, intent(in) :: dis
    integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibound
    ! -- local
    ! -- formats
    !
    ! -- adv pointers to arguments that were passed in
    this%dis => dis
    this%ibound => ibound

adv_cq()

subroutine tspadvmodule::adv_cq	(	class(tspadvtype)	this,
		real(dp), dimension(:), intent(in)	cnew,
		real(dp), dimension(:), intent(inout)	flowja
	)

Definition at line 329 of file tsp-adv.f90.

    ! -- modules
    ! -- dummy
    class(TspAdvType) :: this
    real(DP), intent(in), dimension(:) :: cnew
    real(DP), intent(inout), dimension(:) :: flowja
    ! -- local
    integer(I4B) :: nodes
    integer(I4B) :: n, m, idiag, ipos
    real(DP) :: omega, qnm
    !
    ! -- Calculate advection and add to flowja. qnm is the volumetric flow
    !    rate and has dimensions of L^/T.
    nodes = this%dis%nodes
    do n = 1, nodes
      if (this%ibound(n) == 0) cycle
      idiag = this%dis%con%ia(n)
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        m = this%dis%con%ja(ipos)
        if (this%ibound(m) == 0) cycle
        qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
        omega = this%adv_weight(this%iadvwt, ipos, n, m, qnm)
        flowja(ipos) = flowja(ipos) + qnm * omega * cnew(n) + &
                       qnm * (done - omega) * cnew(m)
      end do
    end do
    !
    ! -- TVD
    if (this%iadvwt == 2) call this%advtvd_bd(cnew, flowja)

adv_cr()

subroutine, public tspadvmodule::adv_cr	(	type(tspadvtype), pointer	advobj,
		character(len=*), intent(in)	name_model,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout,
		type(tspfmitype), intent(in), target	fmi,
		real(dp), intent(in), pointer	eqnsclfac
	)

Create a new ADV package

Parameters

[in]

eqnsclfac

governing equation scale factor

Definition at line 49 of file tsp-adv.f90.

    ! -- dummy
    type(TspAdvType), pointer :: advobj
    character(len=*), intent(in) :: name_model
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    type(TspFmiType), intent(in), target :: fmi
    real(DP), intent(in), pointer :: eqnsclfac !< governing equation scale factor
    !
    ! -- Create the object
    allocate (advobj)
    !
    ! -- create name and memory path
    call advobj%set_names(1, name_model, 'ADV', 'ADV')
    !
    ! -- Allocate scalars
    call advobj%allocate_scalars()
    !
    ! -- Set variables
    advobj%inunit = inunit
    advobj%iout = iout
    advobj%fmi => fmi
    advobj%eqnsclfac => eqnsclfac

Here is the caller graph for this function:

adv_da()

subroutine tspadvmodule::adv_da ( class(tspadvtype)  this )

private

Definition at line 387 of file tsp-adv.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(TspAdvType) :: this
    !
    ! -- Deallocate arrays if package was active
    if (this%inunit > 0) then
    end if
    !
    ! -- nullify pointers
    this%ibound => null()
    !
    ! -- Scalars
    call mem_deallocate(this%iadvwt)
    call mem_deallocate(this%ats_percel)
    !
    ! -- deallocate parent
    call this%NumericalPackageType%da()

adv_df()

subroutine tspadvmodule::adv_df ( class(tspadvtype)  this, type(tspadvoptionstype), intent(in), optional  adv_options  )

private

Define the ADV package

Parameters

[in]

adv_options

the optional options, for when not constructing from file

Definition at line 78 of file tsp-adv.f90.

    ! -- dummy
    class(TspAdvType) :: this
    type(TspAdvOptionsType), optional, intent(in) :: adv_options !< the optional options, for when not constructing from file
    ! -- local
    character(len=*), parameter :: fmtadv = &
      "(1x,/1x,'ADV-- ADVECTION PACKAGE, VERSION 1, 8/25/2017', &
      &' INPUT READ FROM UNIT ', i0, //)"
    !
    ! -- Read or set advection options
    if (.not. present(adv_options)) then
      !
      ! -- Initialize block parser (adv has no define, so it's
      ! not done until here)
      call this%parser%Initialize(this%inunit, this%iout)
      !
      ! --print a message identifying the advection package.
      write (this%iout, fmtadv) this%inunit
      !
      ! --read options from file
      call this%read_options()
    else
      !
      ! --set options from input arg
      this%iadvwt = adv_options%iAdvScheme
    end if

adv_dt()

subroutine tspadvmodule::adv_dt ( class(tspadvtype)  this, real(dp), intent(out)  dtmax, character(len=*), intent(inout)  msg, real(dp), dimension(:), intent(in)  thetam  )

private

Return the largest time step that meets stability constraints

Parameters

	this	this instance
[out]	dtmax	maximum allowable dt subject to stability constraint
[in,out]	msg	package/cell dt constraint message
[in]	thetam	porosity

Definition at line 128 of file tsp-adv.f90.

    ! dummy
    class(TspAdvType) :: this !< this instance
    real(DP), intent(out) :: dtmax !< maximum allowable dt subject to stability constraint
    character(len=*), intent(inout) :: msg !< package/cell dt constraint message
    real(DP), dimension(:), intent(in) :: thetam !< porosity
    ! local
    integer(I4B) :: n
    integer(I4B) :: m
    integer(I4B) :: ipos
    integer(I4B) :: nrmax
    character(len=LINELENGTH) :: cellstr
    real(DP) :: dt
    real(DP) :: flowmax
    real(DP) :: flowsumpos
    real(DP) :: flowsumneg
    real(DP) :: flownm
    real(DP) :: cell_volume
    dtmax = dnodata
    nrmax = 0
    msg = ''
 
    ! If ats_percel not specified by user, then return without making
    ! the courant time step calculation
    if (this%ats_percel == dnodata) then
      return
    end if
 
    ! Calculate time step lengths based on stability constraint for each cell
    ! and store the smallest one
    do n = 1, this%dis%nodes
      if (this%ibound(n) == 0) cycle
      flowsumneg = dzero
      flowsumpos = dzero
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        if (this%dis%con%mask(ipos) == 0) cycle
        m = this%dis%con%ja(ipos)
        if (this%ibound(m) == 0) cycle
        flownm = this%fmi%gwfflowja(ipos)
        if (flownm < dzero) then
          flowsumneg = flowsumneg - flownm
        else
          flowsumpos = flowsumpos + flownm
        end if
      end do
      flowmax = max(flowsumneg, flowsumpos)
      if (flowmax < dprec) cycle
      cell_volume = this%dis%get_cell_volume(n, this%dis%top(n))
      dt = cell_volume * this%fmi%gwfsat(n) * thetam(n) / flowmax
      dt = dt * this%ats_percel
      if (dt < dtmax) then
        dtmax = dt
        nrmax = n
      end if
    end do
    if (nrmax > 0) then
      call this%dis%noder_to_string(nrmax, cellstr)
      write (msg, *) adjustl(trim(this%memoryPath))//'-'//trim(cellstr)
    end if

adv_fc()

subroutine tspadvmodule::adv_fc ( class(tspadvtype)  this, integer, intent(in)  nodes, class(matrixbasetype), pointer  matrix_sln, integer(i4b), dimension(:), intent(in)  idxglo, real(dp), dimension(:), intent(in)  cnew, real(dp), dimension(:), intent(inout)  rhs  )

private

Method to calculate coefficients and fill amat and rhs.

Definition at line 193 of file tsp-adv.f90.

    ! -- modules
    ! -- dummy
    class(TspAdvType) :: this
    integer, intent(in) :: nodes
    class(MatrixBaseType), pointer :: matrix_sln
    integer(I4B), intent(in), dimension(:) :: idxglo
    real(DP), intent(in), dimension(:) :: cnew
    real(DP), dimension(:), intent(inout) :: rhs
    ! -- local
    integer(I4B) :: n, m, idiag, ipos
    real(DP) :: omega, qnm
    !
    ! -- Calculate advection terms and add to solution rhs and hcof.  qnm
    !    is the volumetric flow rate and has dimensions of L^/T.
    do n = 1, nodes
      if (this%ibound(n) == 0) cycle
      idiag = this%dis%con%ia(n)
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        if (this%dis%con%mask(ipos) == 0) cycle
        m = this%dis%con%ja(ipos)
        if (this%ibound(m) == 0) cycle
        qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
        omega = this%adv_weight(this%iadvwt, ipos, n, m, qnm)
        call matrix_sln%add_value_pos(idxglo(ipos), qnm * (done - omega))
        call matrix_sln%add_value_pos(idxglo(idiag), qnm * omega)
      end do
    end do
    !
    ! -- TVD
    if (this%iadvwt == 2) then
      do n = 1, nodes
        if (this%ibound(n) == 0) cycle
        call this%advtvd(n, cnew, rhs)
      end do
    end if

adv_weight()

real(dp) function tspadvmodule::adv_weight	(	class(tspadvtype)	this,
		integer, intent(in)	iadvwt,
		integer, intent(in)	ipos,
		integer, intent(in)	n,
		integer, intent(in)	m,
		real(dp), intent(in)	qnm
	)

Calculate the advection weight

Definition at line 504 of file tsp-adv.f90.

    ! -- return
    real(DP) :: omega
    ! -- dummy
    class(TspAdvType) :: this
    integer, intent(in) :: iadvwt
    integer, intent(in) :: ipos
    integer, intent(in) :: n
    integer, intent(in) :: m
    real(DP), intent(in) :: qnm
    ! -- local
    real(DP) :: lnm, lmn
 
    select case (iadvwt)
    case (1)
      ! -- calculate weight based on distances between nodes and the shared
      !    face of the connection
      if (this%dis%con%ihc(this%dis%con%jas(ipos)) == 0) then
        ! -- vertical connection; assume cell is fully saturated
        lnm = dhalf * (this%dis%top(n) - this%dis%bot(n))
        lmn = dhalf * (this%dis%top(m) - this%dis%bot(m))
      else
        ! -- horizontal connection
        lnm = this%dis%con%cl1(this%dis%con%jas(ipos))
        lmn = this%dis%con%cl2(this%dis%con%jas(ipos))
      end if
      omega = lmn / (lnm + lmn)
    case (0, 2)
      ! -- use upstream weighting for upstream and tvd schemes
      if (qnm > dzero) then
        omega = dzero
      else
        omega = done
      end if
    end select

advqtvd()

real(dp) function tspadvmodule::advqtvd ( class(tspadvtype)  this, integer(i4b), intent(in)  n, integer(i4b), intent(in)  m, integer(i4b), intent(in)  iposnm, real(dp), dimension(:), intent(in)  cnew  )

private

Use explicit scheme to calculate the advective component of transport. TVD is an acronym for Total-Variation Diminishing

Definition at line 264 of file tsp-adv.f90.

    ! -- modules
    use constantsmodule, only: dprec
    ! -- return
    real(DP) :: qtvd
    ! -- dummy
    class(TspAdvType) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: m
    integer(I4B), intent(in) :: iposnm
    real(DP), dimension(:), intent(in) :: cnew
    ! -- local
    integer(I4B) :: ipos, isympos, iup, idn, i2up, j
    real(DP) :: qnm, qmax, qupj, elupdn, elup2up
    real(DP) :: smooth, cdiff, alimiter
    !
    ! -- initialize
    qtvd = dzero
    !
    ! -- Find upstream node
    isympos = this%dis%con%jas(iposnm)
    qnm = this%fmi%gwfflowja(iposnm)
    if (qnm > dzero) then
      ! -- positive flow into n means m is upstream
      iup = m
      idn = n
    else
      iup = n
      idn = m
    end if
    elupdn = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos)
    !
    ! -- Find second node upstream to iup
    i2up = 0
    qmax = dzero
    do ipos = this%dis%con%ia(iup) + 1, this%dis%con%ia(iup + 1) - 1
      j = this%dis%con%ja(ipos)
      if (this%ibound(j) == 0) cycle
      qupj = this%fmi%gwfflowja(ipos)
      isympos = this%dis%con%jas(ipos)
      if (qupj > qmax) then
        qmax = qupj
        i2up = j
        elup2up = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos)
      end if
    end do
    !
    ! -- Calculate flux limiting term
    if (i2up > 0) then
      smooth = dzero
      cdiff = abs(cnew(idn) - cnew(iup))
      if (cdiff > dprec) then
        smooth = (cnew(iup) - cnew(i2up)) / elup2up * &
                 elupdn / (cnew(idn) - cnew(iup))
      end if
      if (smooth > dzero) then
        alimiter = dtwo * smooth / (done + smooth)
        qtvd = dhalf * alimiter * qnm * (cnew(idn) - cnew(iup))
        qtvd = qtvd * this%eqnsclfac
      end if
    end if

advtvd()

subroutine tspadvmodule::advtvd ( class(tspadvtype)  this, integer(i4b), intent(in)  n, real(dp), dimension(:), intent(in)  cnew, real(dp), dimension(:), intent(inout)  rhs  )

private

Use explicit scheme to calculate the advective component of transport. TVD is an acronym for Total-Variation Diminishing

Definition at line 236 of file tsp-adv.f90.

    ! -- modules
    ! -- dummy
    class(TspAdvType) :: this
    integer(I4B), intent(in) :: n
    real(DP), dimension(:), intent(in) :: cnew
    real(DP), dimension(:), intent(inout) :: rhs
    ! -- local
    real(DP) :: qtvd
    integer(I4B) :: m, ipos
    !
    ! -- Loop through each n connection.  This will
    do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
      if (this%dis%con%mask(ipos) == 0) cycle
      m = this%dis%con%ja(ipos)
      if (m > n .and. this%ibound(m) /= 0) then
        qtvd = this%advqtvd(n, m, ipos, cnew)
        rhs(n) = rhs(n) - qtvd
        rhs(m) = rhs(m) + qtvd
      end if
    end do

advtvd_bd()

subroutine tspadvmodule::advtvd_bd ( class(tspadvtype)  this, real(dp), dimension(:), intent(in)  cnew, real(dp), dimension(:), intent(inout)  flowja  )

private

Definition at line 361 of file tsp-adv.f90.

    ! -- modules
    ! -- dummy
    class(TspAdvType) :: this
    real(DP), dimension(:), intent(in) :: cnew
    real(DP), dimension(:), intent(inout) :: flowja
    ! -- local
    real(DP) :: qtvd, qnm
    integer(I4B) :: nodes, n, m, ipos
    !
    nodes = this%dis%nodes
    do n = 1, nodes
      if (this%ibound(n) == 0) cycle
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        m = this%dis%con%ja(ipos)
        if (this%ibound(m) /= 0) then
          qnm = this%fmi%gwfflowja(ipos)
          qtvd = this%advqtvd(n, m, ipos, cnew)
          flowja(ipos) = flowja(ipos) + qtvd
        end if
      end do
    end do

allocate_scalars()

subroutine tspadvmodule::allocate_scalars

(

class(tspadvtype)

this

)

Definition at line 411 of file tsp-adv.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    ! -- dummy
    class(TspAdvType) :: this
    ! -- local
    !
    ! -- allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- Allocate
    call mem_allocate(this%iadvwt, 'IADVWT', this%memoryPath)
    call mem_allocate(this%ats_percel, 'ATS_PERCEL', this%memoryPath)
    !
    ! -- Initialize
    this%iadvwt = 0
    this%ats_percel = dnodata
    !
    ! -- Advection creates an asymmetric coefficient matrix
    this%iasym = 1

read_options()

subroutine tspadvmodule::read_options

(

class(tspadvtype)

this

)

Read the options block

Definition at line 437 of file tsp-adv.f90.

    ! -- modules
    use constantsmodule, only: linelength
    use simmodule, only: store_error
    ! -- dummy
    class(TspAdvType) :: this
    ! -- local
    character(len=LINELENGTH) :: errmsg, keyword
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    ! -- formats
    character(len=*), parameter :: fmtiadvwt = &
      &"(4x,'ADVECTION WEIGHTING SCHEME HAS BEEN SET TO: ', a)"
    !
    ! -- get options block
    call this%parser%GetBlock('OPTIONS', isfound, ierr, blockrequired=.false., &
                              supportopenclose=.true.)
    !
    ! -- parse options block if detected
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING ADVECTION OPTIONS'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('SCHEME')
          call this%parser%GetStringCaps(keyword)
          select case (keyword)
          case ('UPSTREAM')
            this%iadvwt = 0
            write (this%iout, fmtiadvwt) 'UPSTREAM'
          case ('CENTRAL')
            this%iadvwt = 1
            write (this%iout, fmtiadvwt) 'CENTRAL'
          case ('TVD')
            this%iadvwt = 2
            write (this%iout, fmtiadvwt) 'TVD'
          case default
            write (errmsg, '(a, a)') &
              'Unknown scheme: ', trim(keyword)
            call store_error(errmsg)
            write (errmsg, '(a, a)') &
              'Scheme must be "UPSTREAM", "CENTRAL" or "TVD"'
            call store_error(errmsg)
            call this%parser%StoreErrorUnit()
          end select
        case ('ATS_PERCEL')
          this%ats_percel = this%parser%GetDouble()
          if (this%ats_percel == dzero) this%ats_percel = dnodata
          write (this%iout, '(4x,a,1pg15.6)') &
            'User-specified fractional cell distance for adaptive time &
            &steps: ', this%ats_percel
        case default
          write (errmsg, '(a,a)') 'Unknown ADVECTION option: ', &
            trim(keyword)
          call store_error(errmsg, terminate=.true.)
        end select
      end do
      write (this%iout, '(1x,a)') 'END OF ADVECTION OPTIONS'
    end if

Here is the call graph for this function: