gwfstomodule Module Reference

This module contains the storage package methods. More...

Data Types
type	gwfstotype

Functions/Subroutines
subroutine, public	sto_cr (stoobj, name_model, inunit, iout)
	@ brief Create a new package object More...
subroutine	sto_ar (this, dis, ibound)
	@ brief Allocate and read method for package More...
subroutine	sto_rp (this)
	@ brief Read and prepare method for package More...
subroutine	sto_ad (this)
	@ brief Advance the package More...
subroutine	sto_fc (this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
	@ brief Fill A and right-hand side for the package More...
subroutine	sto_fn (this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
	@ brief Fill Newton-Raphson terms in A and right-hand side for the package More...
subroutine	sto_cq (this, flowja, hnew, hold)
	@ brief Calculate flows for package More...
subroutine	sto_bd (this, isuppress_output, model_budget)
	@ brief Model budget calculation for package More...
subroutine	sto_save_model_flows (this, icbcfl, icbcun)
	@ brief Save model flows for package More...
subroutine	sto_da (this)
	@ brief Deallocate package memory More...
subroutine	allocate_scalars (this)
	@ brief Allocate scalars More...
subroutine	allocate_arrays (this, nodes)
	@ brief Allocate package arrays More...
subroutine	read_options (this)
	@ brief Read options for package More...
subroutine	read_data (this)
	@ brief Read data for package More...
subroutine	save_old_ss_sy (this)
	@ brief Save old storage property values More...

Variables
character(len=lenbudtxt), dimension(2)	budtxt = [' STO-SS', ' STO-SY']

Detailed Description

This module contains the methods used to add the effects of storage on the groundwater flow equation. The contribution of specific storage and specific yield can be represented.

Function/Subroutine Documentation

allocate_arrays()

subroutine gwfstomodule::allocate_arrays	(	class(gwfstotype), target	this,
		integer(i4b), intent(in)	nodes
	)

Allocate and initialize STO package arrays.

Parameters

	this	GwfStoType object
[in]	nodes	active model nodes

Definition at line 805 of file gwf-sto.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(GwfStoType), target :: this !< GwfStoType object
    integer(I4B), intent(in) :: nodes !< active model nodes
    ! -- local variables
    integer(I4B) :: n
    !
    ! -- Allocate arrays
    call mem_allocate(this%iconvert, nodes, 'ICONVERT', this%memoryPath)
    call mem_allocate(this%ss, nodes, 'SS', this%memoryPath)
    call mem_allocate(this%sy, nodes, 'SY', this%memoryPath)
    call mem_allocate(this%strgss, nodes, 'STRGSS', this%memoryPath)
    call mem_allocate(this%strgsy, nodes, 'STRGSY', this%memoryPath)
    !
    ! -- Initialize arrays
    this%iss = 0
    do n = 1, nodes
      this%iconvert(n) = 1
      this%ss(n) = dzero
      this%sy(n) = dzero
      this%strgss(n) = dzero
      this%strgsy(n) = dzero
      if (this%integratechanges /= 0) then
        this%oldss(n) = dzero
        if (this%iusesy /= 0) then
          this%oldsy(n) = dzero
        end if
      end if
    end do
    !
    ! -- return
    return

allocate_scalars()

subroutine gwfstomodule::allocate_scalars

(

class(gwfstotype)

this

)

Allocate and initialize scalars for the STO package. The base numerical package allocate scalars method is also called.

Parameters

this	GwfStoType object

Definition at line 768 of file gwf-sto.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    !
    ! -- allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- allocate scalars
    call mem_allocate(this%istor_coef, 'ISTOR_COEF', this%memoryPath)
    call mem_allocate(this%iconf_ss, 'ICONF_SS', this%memoryPath)
    call mem_allocate(this%iorig_ss, 'IORIG_SS', this%memoryPath)
    call mem_allocate(this%iusesy, 'IUSESY', this%memoryPath)
    call mem_allocate(this%satomega, 'SATOMEGA', this%memoryPath)
    call mem_allocate(this%integratechanges, 'INTEGRATECHANGES', &
                      this%memoryPath)
    call mem_allocate(this%intvs, 'INTVS', this%memoryPath)
    !
    ! -- initialize scalars
    this%istor_coef = 0
    this%iconf_ss = 0
    this%iorig_ss = 0
    this%iusesy = 0
    this%satomega = dzero
    this%integratechanges = 0
    this%intvs = 0
    !
    ! -- return
    return

read_data()

subroutine gwfstomodule::read_data

(

class(gwfstotype)

this

)

Read griddata block for STO package.

Parameters

this	GwfStoType object

Definition at line 944 of file gwf-sto.f90.

    ! -- modules
    ! -- dummy variables
    class(GwfStotype) :: this !< GwfStoType object
    ! -- local variables
    character(len=LINELENGTH) :: keyword
    character(len=:), allocatable :: line
    character(len=LINELENGTH) :: cellstr
    integer(I4B) :: istart, istop, lloc, ierr
    logical :: isfound, endOfBlock
    logical :: readiconv
    logical :: readss
    logical :: readsy
    logical :: isconv
    character(len=24), dimension(4) :: aname
    integer(I4B) :: n
    ! -- formats
    !data
    data aname(1)/'                ICONVERT'/
    data aname(2)/'        SPECIFIC STORAGE'/
    data aname(3)/'          SPECIFIC YIELD'/
    data aname(4)/'     STORAGE COEFFICIENT'/
    !
    ! -- initialize
    isfound = .false.
    readiconv = .false.
    readss = .false.
    readsy = .false.
    isconv = .false.
    !
    ! -- get stodata block
    call this%parser%GetBlock('GRIDDATA', isfound, ierr)
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING GRIDDATA'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        call this%parser%GetRemainingLine(line)
        lloc = 1
        select case (keyword)
        case ('ICONVERT')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%iconvert, &
                                        aname(1))
          readiconv = .true.
        case ('SS')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%ss, &
                                        aname(2))
          readss = .true.
        case ('SY')
          call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, &
                                        this%parser%iuactive, this%sy, &
                                        aname(3))
          readsy = .true.
        case default
          write (errmsg, '(a,a)') 'Unknown GRIDDATA tag: ', &
            trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END PROCESSING GRIDDATA'
    else
      write (errmsg, '(a)') 'Required GRIDDATA block not found.'
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- Check for ICONVERT
    if (.not. readiconv) then
      write (errmsg, '(a, a, a)') 'Error in GRIDDATA block: ', &
        trim(adjustl(aname(1))), ' not found.'
      call store_error(errmsg)
    else
      isconv = .false.
      do n = 1, this%dis%nodes
        if (this%iconvert(n) /= 0) then
          isconv = .true.
          this%iusesy = 1
          exit
        end if
      end do
    end if
    !
    ! -- Check for SS
    if (.not. readss) then
      write (errmsg, '(a, a, a)') 'Error in GRIDDATA block: ', &
        trim(adjustl(aname(2))), ' not found.'
      call store_error(errmsg)
    end if
    !
    ! -- Check for SY
    if (.not. readsy .and. isconv) then
      write (errmsg, '(a, a, a)') 'Error in GRIDDATA block: ', &
        trim(adjustl(aname(3))), ' not found.'
      call store_error(errmsg)
    end if
    !
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- Check SS and SY for negative values
    do n = 1, this%dis%nodes
      if (this%ss(n) < dzero) then
        call this%dis%noder_to_string(n, cellstr)
        write (errmsg, '(a,2(1x,a),1x,g0,1x,a)') &
          'Error in SS DATA: SS value in cell', trim(adjustl(cellstr)), &
          'is less than zero (', this%ss(n), ').'
        call store_error(errmsg)
      end if
      if (readsy) then
        if (this%sy(n) < dzero) then
          call this%dis%noder_to_string(n, cellstr)
          write (errmsg, '(a,2(1x,a),1x,g0,1x,a)') &
            'Error in SY DATA: SY value in cell', trim(adjustl(cellstr)), &
            'is less than zero (', this%sy(n), ').'
          call store_error(errmsg)
        end if
      end if
    end do
    !
    ! -- return
    return

Here is the call graph for this function:

read_options()

subroutine gwfstomodule::read_options

(

class(gwfstotype)

this

)

Read options block for STO package.

Parameters

this	GwfStoType object

Definition at line 846 of file gwf-sto.f90.

    ! -- modules
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    ! -- local variables
    character(len=LINELENGTH) :: keyword, fname
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    ! -- formats
    character(len=*), parameter :: fmtisvflow = &
      "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
      &WHENEVER ICBCFL IS NOT ZERO.')"
    character(len=*), parameter :: fmtflow = &
      &"(4x, 'FLOWS WILL BE SAVED TO FILE: ', a, /4x, 'OPENED ON UNIT: ', I7)"
    character(len=*), parameter :: fmtorigss = &
      "(4X,'ORIGINAL_SPECIFIC_STORAGE OPTION:',/, &
      &1X,'The original specific storage formulation will be used')"
    character(len=*), parameter :: fmtstoc = &
      "(4X,'STORAGECOEFFICIENT OPTION:',/, &
      &1X,'Read storage coefficient rather than specific storage')"
    character(len=*), parameter :: fmtconfss = &
      "(4X,'SS_CONFINED_ONLY OPTION:',/, &
      &1X,'Specific storage changes only occur under confined conditions')"
    !
    ! -- get options block
    call this%parser%GetBlock('OPTIONS', isfound, ierr, &
                              supportopenclose=.true., blockrequired=.false.)
    !
    ! -- parse options block if detected
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING STORAGE OPTIONS'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('SAVE_FLOWS')
          this%ipakcb = -1
          write (this%iout, fmtisvflow)
        case ('STORAGECOEFFICIENT')
          this%istor_coef = 1
          write (this%iout, fmtstoc)
        case ('SS_CONFINED_ONLY')
          this%iconf_ss = 1
          this%iorig_ss = 0
          write (this%iout, fmtconfss)
        case ('TVS6')
          if (this%intvs /= 0) then
            errmsg = 'Multiple TVS6 keywords detected in OPTIONS block.'// &
                     ' Only one TVS6 entry allowed.'
            call store_error(errmsg, terminate=.true.)
          end if
          call this%parser%GetStringCaps(keyword)
          if (trim(adjustl(keyword)) /= 'FILEIN') then
            errmsg = 'TVS6 keyword must be followed by "FILEIN" '// &
                     'then by filename.'
            call store_error(errmsg, terminate=.true.)
          end if
          call this%parser%GetString(fname)
          this%intvs = getunit()
          call openfile(this%intvs, this%iout, fname, 'TVS')
          call tvs_cr(this%tvs, this%name_model, this%intvs, this%iout)
          !
          ! -- right now these are options that are only available in the
          !    development version and are not included in the documentation.
          !    These options are only available when IDEVELOPMODE in
          !    constants module is set to 1
        case ('DEV_ORIGINAL_SPECIFIC_STORAGE')
          this%iorig_ss = 1
          write (this%iout, fmtorigss)
        case ('DEV_OLDSTORAGEFORMULATION')
          call this%parser%DevOpt()
          this%iconf_ss = 1
          this%iorig_ss = 0
          write (this%iout, fmtconfss)
        case default
          write (errmsg, '(a,a)') 'Unknown STO option: ', &
            trim(keyword)
          call store_error(errmsg, terminate=.true.)
        end select
      end do
      write (this%iout, '(1x,a)') 'END OF STORAGE OPTIONS'
    end if
    !
    ! -- set omega value used for saturation calculations
    if (this%inewton > 0) then
      this%satomega = dem6
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

save_old_ss_sy()

subroutine gwfstomodule::save_old_ss_sy

(

class(gwfstotype)

this

)

Save ss and sy values from the previous time step for use with storage integration when integratechanges is non-zero.

Parameters

this	GwfStoType object

Definition at line 1078 of file gwf-sto.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    ! -- local variables
    integer(I4B) :: n
    !
    ! -- Allocate TVS arrays if needed
    if (.not. associated(this%oldss)) then
      call mem_allocate(this%oldss, this%dis%nodes, 'OLDSS', this%memoryPath)
    end if
    if (this%iusesy == 1 .and. .not. associated(this%oldsy)) then
      call mem_allocate(this%oldsy, this%dis%nodes, 'OLDSY', this%memoryPath)
    end if
    !
    ! -- Save current specific storage
    do n = 1, this%dis%nodes
      this%oldss(n) = this%ss(n)
    end do
    !
    ! -- Save current specific yield, if used
    if (this%iusesy == 1) then
      do n = 1, this%dis%nodes
        this%oldsy(n) = this%sy(n)
      end do
    end if
    !
    ! -- Return
    return

sto_ad()

subroutine gwfstomodule::sto_ad

(

class(gwfstotype)

this

)

Advance data in the STO package.

Parameters

this	GwfStoType object

Definition at line 256 of file gwf-sto.f90.

    ! -- modules
    use tdismodule, only: kstp
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    !
    ! -- Store ss and sy values from end of last time step if needed
    if (this%integratechanges /= 0 .and. kstp > 1) then
      call this%save_old_ss_sy()
    end if
    !
    ! -- TVS
    if (this%intvs /= 0) then
      call this%tvs%ad()
    end if
    !
    ! -- return
    return

sto_ar()

subroutine gwfstomodule::sto_ar	(	class(gwfstotype)	this,
		class(disbasetype), intent(in), pointer	dis,
		integer(i4b), dimension(:), pointer, contiguous	ibound
	)

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

Parameters

	this	GwfStoType object
[in]	dis	model discretization object
	ibound	model ibound array

Definition at line 107 of file gwf-sto.f90.

    ! -- modules
    use memorymanagermodule, only: mem_setptr
    use memoryhelpermodule, only: create_mem_path
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    class(DisBaseType), pointer, intent(in) :: dis !< model discretization object
    integer(I4B), dimension(:), pointer, contiguous :: ibound !< model ibound array
    ! -- local variables
    ! -- formats
    character(len=*), parameter :: fmtsto = &
      "(1x,/1x,'STO -- STORAGE PACKAGE, VERSION 1, 5/19/2014', &
      &' INPUT READ FROM UNIT ', i0, //)"
    !
    ! --print a message identifying the storage package.
    write (this%iout, fmtsto) this%inunit
    !
    ! -- store pointers to arguments that were passed in
    this%dis => dis
    this%ibound => ibound
    !
    ! -- set pointer to gwf iss
    call mem_setptr(this%iss, 'ISS', create_mem_path(this%name_model))
    !
    ! -- Allocate arrays
    call this%allocate_arrays(dis%nodes)
    !!
    !! -- Register side effect handlers
    !call this%register_handlers()
    !
    ! -- Read storage options
    call this%read_options()
    !
    ! -- read the data block
    call this%read_data()
    !
    ! -- TVS
    if (this%intvs /= 0) then
      call this%tvs%ar(this%dis)
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sto_bd()

subroutine gwfstomodule::sto_bd	(	class(gwfstotype)	this,
		integer(i4b), intent(in)	isuppress_output,
		type(budgettype), intent(inout)	model_budget
	)

Budget calculation for the STO package components. Components include specific storage and specific yield storage.

Parameters

	this	GwfStoType object
[in]	isuppress_output	flag to suppress model output
[in,out]	model_budget	model budget object

Definition at line 636 of file gwf-sto.f90.

    ! -- modules
    use tdismodule, only: delt
    use budgetmodule, only: budgettype, rate_accumulator
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    integer(I4B), intent(in) :: isuppress_output !< flag to suppress model output
    type(BudgetType), intent(inout) :: model_budget !< model budget object
    ! -- local variables
    real(DP) :: rin
    real(DP) :: rout
    !
    ! -- Add confined storage rates to model budget
    call rate_accumulator(this%strgss, rin, rout)
    call model_budget%addentry(rin, rout, delt, budtxt(1), &
                               isuppress_output, '         STORAGE')
    !
    ! -- Add unconfined storage rates to model budget
    if (this%iusesy == 1) then
      call rate_accumulator(this%strgsy, rin, rout)
      call model_budget%addentry(rin, rout, delt, budtxt(2), &
                                 isuppress_output, '         STORAGE')
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sto_cq()

subroutine gwfstomodule::sto_cq	(	class(gwfstotype)	this,
		real(dp), dimension(:), intent(inout), contiguous	flowja,
		real(dp), dimension(:), intent(in), contiguous	hnew,
		real(dp), dimension(:), intent(in), contiguous	hold
	)

Flow calculation for the STO package components. Components include specific storage and specific yield storage.

Parameters

	this	GwfStoType object
[in,out]	flowja	connection flows
[in]	hnew	current head
[in]	hold	previous head

Definition at line 508 of file gwf-sto.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    real(DP), dimension(:), contiguous, intent(inout) :: flowja !< connection flows
    real(DP), dimension(:), contiguous, intent(in) :: hnew !< current head
    real(DP), dimension(:), contiguous, intent(in) :: hold !< previous head
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: rate
    real(DP) :: tled
    real(DP) :: sc1
    real(DP) :: sc2
    real(DP) :: rho1
    real(DP) :: rho2
    real(DP) :: sc1old
    real(DP) :: sc2old
    real(DP) :: rho1old
    real(DP) :: rho2old
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: snold
    real(DP) :: snnew
    real(DP) :: aterm
    real(DP) :: rhsterm
    !
    ! -- initialize strg arrays
    do n = 1, this%dis%nodes
      this%strgss(n) = dzero
      this%strgsy(n) = dzero
    end do
    !
    ! -- Set strt to zero or calculate terms if not steady-state stress period
    if (this%iss == 0) then
      !
      ! -- set variables
      tled = done / delt
      !
      ! -- Calculate storage change
      do n = 1, this%dis%nodes
        if (this%ibound(n) <= 0) cycle
        ! -- aquifer elevations and thickness
        tp = this%dis%top(n)
        bt = this%dis%bot(n)
        !
        ! -- aquifer saturation
        if (this%iconvert(n) == 0) then
          snold = done
          snnew = done
        else
          snold = squadraticsaturation(tp, bt, hold(n), this%satomega)
          snnew = squadraticsaturation(tp, bt, hnew(n), this%satomega)
        end if
        !
        ! -- primary storage coefficient
        sc1 = sscapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
        rho1 = sc1 * tled
        !
        if (this%integratechanges /= 0) then
          ! -- Integration of storage changes (e.g. when using TVS):
          !    separate the old (start of time step) and new (end of time step)
          !    primary storage capacities
          sc1old = sscapacity(this%istor_coef, tp, bt, this%dis%area(n), &
                              this%oldss(n))
          rho1old = sc1old * tled
        else
          ! -- No integration of storage changes: old and new values are
          !    identical => normal MF6 storage formulation
          rho1old = rho1
        end if
        !
        ! -- calculate specific storage terms and rate
        call ssterms(this%iconvert(n), this%iorig_ss, this%iconf_ss, tp, bt, &
                     rho1, rho1old, snnew, snold, hnew(n), hold(n), &
                     aterm, rhsterm, rate)
        !
        ! -- save rate
        this%strgss(n) = rate
        !
        ! -- add storage term to flowja
        idiag = this%dis%con%ia(n)
        flowja(idiag) = flowja(idiag) + rate
        !
        ! -- specific yield
        rate = dzero
        if (this%iconvert(n) /= 0) then
          !
          ! -- secondary storage coefficient
          sc2 = sycapacity(this%dis%area(n), this%sy(n))
          rho2 = sc2 * tled
          !
          if (this%integratechanges /= 0) then
            ! -- Integration of storage changes (e.g. when using TVS):
            !    separate the old (start of time step) and new (end of time
            !    step) secondary storage capacities
            sc2old = sycapacity(this%dis%area(n), this%oldsy(n))
            rho2old = sc2old * tled
          else
            ! -- No integration of storage changes: old and new values are
            !    identical => normal MF6 storage formulation
            rho2old = rho2
          end if
          !
          ! -- calculate specific yield storage terms and rate
          call syterms(tp, bt, rho2, rho2old, snnew, snold, &
                       aterm, rhsterm, rate)
 
        end if
        this%strgsy(n) = rate
        !
        ! -- add storage term to flowja
        idiag = this%dis%con%ia(n)
        flowja(idiag) = flowja(idiag) + rate
      end do
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sto_cr()

subroutine, public gwfstomodule::sto_cr	(	type(gwfstotype), pointer	stoobj,
		character(len=*), intent(in)	name_model,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout
	)

Create a new storage (STO) object

Parameters

	stoobj	GwfStoType object
[in]	name_model	name of model
[in]	inunit	package input file unit
[in]	iout	model listing file unit

Definition at line 75 of file gwf-sto.f90.

    ! -- dummy variables
    type(GwfStoType), pointer :: stoobj !< GwfStoType object
    character(len=*), intent(in) :: name_model !< name of model
    integer(I4B), intent(in) :: inunit !< package input file unit
    integer(I4B), intent(in) :: iout !< model listing file unit
    !
    ! -- Create the object
    allocate (stoobj)
    !
    ! -- create name and memory path
    call stoobj%set_names(1, name_model, 'STO', 'STO')
    !
    ! -- Allocate scalars
    call stoobj%allocate_scalars()
    !
    ! -- Set variables
    stoobj%inunit = inunit
    stoobj%iout = iout
    !
    ! -- Initialize block parser
    call stoobj%parser%Initialize(stoobj%inunit, stoobj%iout)
    !
    ! -- return
    return

Here is the caller graph for this function:

sto_da()

subroutine gwfstomodule::sto_da

(

class(gwfstotype)

this

)

Deallocate STO package scalars and arrays.

Parameters

this	GwfStoType object

Definition at line 717 of file gwf-sto.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    !
    ! -- TVS
    if (this%intvs /= 0) then
      call this%tvs%da()
      deallocate (this%tvs)
    end if
    !
    ! -- Deallocate arrays if package is active
    if (this%inunit > 0) then
      call mem_deallocate(this%iconvert)
      call mem_deallocate(this%ss)
      call mem_deallocate(this%sy)
      call mem_deallocate(this%strgss)
      call mem_deallocate(this%strgsy)
      !
      ! -- deallocate TVS arrays
      if (associated(this%oldss)) then
        call mem_deallocate(this%oldss)
      end if
      if (associated(this%oldsy)) then
        call mem_deallocate(this%oldsy)
      end if
    end if
    !
    ! -- Deallocate scalars
    call mem_deallocate(this%istor_coef)
    call mem_deallocate(this%iconf_ss)
    call mem_deallocate(this%iorig_ss)
    call mem_deallocate(this%iusesy)
    call mem_deallocate(this%satomega)
    call mem_deallocate(this%integratechanges)
    call mem_deallocate(this%intvs)
    !
    ! -- deallocate parent
    call this%NumericalPackageType%da()
    !
    ! -- return
    return

sto_fc()

subroutine gwfstomodule::sto_fc	(	class(gwfstotype)	this,
		integer(i4b), intent(in)	kiter,
		real(dp), dimension(:), intent(in)	hold,
		real(dp), dimension(:), intent(in)	hnew,
		class(matrixbasetype), pointer	matrix_sln,
		integer(i4b), dimension(:), intent(in)	idxglo,
		real(dp), dimension(:), intent(inout)	rhs
	)

Fill the coefficient matrix and right-hand side with the STO package terms.

Parameters

	this	GwfStoType object
[in]	kiter	outer iteration number
[in]	hold	previous heads
[in]	hnew	current heads
	matrix_sln	A matrix
[in]	idxglo	global index model to solution
[in,out]	rhs	right-hand side

Definition at line 281 of file gwf-sto.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    integer(I4B), intent(in) :: kiter !< outer iteration number
    real(DP), intent(in), dimension(:) :: hold !< previous heads
    real(DP), intent(in), dimension(:) :: hnew !< current heads
    class(MatrixBaseType), pointer :: matrix_sln !< A matrix
    integer(I4B), intent(in), dimension(:) :: idxglo !< global index model to solution
    real(DP), intent(inout), dimension(:) :: rhs !< right-hand side
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: tled
    real(DP) :: sc1
    real(DP) :: sc2
    real(DP) :: rho1
    real(DP) :: rho2
    real(DP) :: sc1old
    real(DP) :: sc2old
    real(DP) :: rho1old
    real(DP) :: rho2old
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: snold
    real(DP) :: snnew
    real(DP) :: aterm
    real(DP) :: rhsterm
    ! -- formats
    character(len=*), parameter :: fmtsperror = &
      &"('DETECTED TIME STEP LENGTH OF ZERO.  GWF STORAGE PACKAGE CANNOT BE ', &
      &'USED UNLESS DELT IS NON-ZERO.')"
    !
    ! -- test if steady-state stress period
    if (this%iss /= 0) return
    !
    ! -- Ensure time step length is not zero
    if (delt == dzero) then
      write (errmsg, fmtsperror)
      call store_error(errmsg, terminate=.true.)
    end if
    !
    ! -- set variables
    tled = done / delt
    !
    ! -- loop through and calculate storage contribution to hcof and rhs
    do n = 1, this%dis%nodes
      idiag = this%dis%con%ia(n)
      if (this%ibound(n) < 1) cycle
      !
      ! -- aquifer elevations and thickness
      tp = this%dis%top(n)
      bt = this%dis%bot(n)
      !
      ! -- aquifer saturation
      if (this%iconvert(n) == 0) then
        snold = done
        snnew = done
      else
        snold = squadraticsaturation(tp, bt, hold(n), this%satomega)
        snnew = squadraticsaturation(tp, bt, hnew(n), this%satomega)
      end if
      !
      ! -- storage coefficients
      sc1 = sscapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
      rho1 = sc1 * tled
      !
      if (this%integratechanges /= 0) then
        ! -- Integration of storage changes (e.g. when using TVS):
        !    separate the old (start of time step) and new (end of time step)
        !    primary storage capacities
        sc1old = sscapacity(this%istor_coef, tp, bt, this%dis%area(n), &
                            this%oldss(n))
        rho1old = sc1old * tled
      else
        ! -- No integration of storage changes: old and new values are
        !    identical => normal MF6 storage formulation
        rho1old = rho1
      end if
      !
      ! -- calculate specific storage terms
      call ssterms(this%iconvert(n), this%iorig_ss, this%iconf_ss, tp, bt, &
                   rho1, rho1old, snnew, snold, hnew(n), hold(n), &
                   aterm, rhsterm)
      !
      ! -- add specific storage terms to amat and rhs
      call matrix_sln%add_value_pos(idxglo(idiag), aterm)
      rhs(n) = rhs(n) + rhsterm
      !
      ! -- specific yield
      if (this%iconvert(n) /= 0) then
        rhsterm = dzero
        !
        ! -- secondary storage coefficient
        sc2 = sycapacity(this%dis%area(n), this%sy(n))
        rho2 = sc2 * tled
        !
        if (this%integratechanges /= 0) then
          ! -- Integration of storage changes (e.g. when using TVS):
          !    separate the old (start of time step) and new (end of time step)
          !    secondary storage capacities
          sc2old = sycapacity(this%dis%area(n), this%oldsy(n))
          rho2old = sc2old * tled
        else
          ! -- No integration of storage changes: old and new values are
          !    identical => normal MF6 storage formulation
          rho2old = rho2
        end if
        !
        ! -- calculate specific storage terms
        call syterms(tp, bt, rho2, rho2old, snnew, snold, &
                     aterm, rhsterm)
!
        ! -- add specific yield terms to amat and rhs
        call matrix_sln%add_value_pos(idxglo(idiag), aterm)
        rhs(n) = rhs(n) + rhsterm
      end if
    end do
    !
    ! -- return
    return

Here is the call graph for this function:

sto_fn()

subroutine gwfstomodule::sto_fn	(	class(gwfstotype)	this,
		integer(i4b), intent(in)	kiter,
		real(dp), dimension(:), intent(in)	hold,
		real(dp), dimension(:), intent(in)	hnew,
		class(matrixbasetype), pointer	matrix_sln,
		integer(i4b), dimension(:), intent(in)	idxglo,
		real(dp), dimension(:), intent(inout)	rhs
	)

Fill the coefficient matrix and right-hand side with STO package with Newton-Raphson terms.

Parameters

	this	GwfStoType object
[in]	kiter	outer iteration number
[in]	hold	previous heads
[in]	hnew	current heads
	matrix_sln	A matrix
[in]	idxglo	global index model to solution
[in,out]	rhs	right-hand side

Definition at line 411 of file gwf-sto.f90.

    ! -- modules
    use tdismodule, only: delt
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    integer(I4B), intent(in) :: kiter !< outer iteration number
    real(DP), intent(in), dimension(:) :: hold !< previous heads
    real(DP), intent(in), dimension(:) :: hnew !< current heads
    class(MatrixBaseType), pointer :: matrix_sln !< A matrix
    integer(I4B), intent(in), dimension(:) :: idxglo !< global index model to solution
    real(DP), intent(inout), dimension(:) :: rhs !< right-hand side
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: idiag
    real(DP) :: tled
    real(DP) :: sc1
    real(DP) :: sc2
    real(DP) :: rho1
    real(DP) :: rho2
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: tthk
    real(DP) :: h
    real(DP) :: snnew
    real(DP) :: derv
    real(DP) :: rterm
    real(DP) :: drterm
    !
    ! -- test if steady-state stress period
    if (this%iss /= 0) return
    !
    ! -- set variables
    tled = done / delt
    !
    ! -- loop through and calculate storage contribution to hcof and rhs
    do n = 1, this%dis%nodes
      idiag = this%dis%con%ia(n)
      if (this%ibound(n) <= 0) cycle
      !
      ! -- aquifer elevations and thickness
      tp = this%dis%top(n)
      bt = this%dis%bot(n)
      tthk = tp - bt
      h = hnew(n)
      !
      ! -- aquifer saturation
      snnew = squadraticsaturation(tp, bt, h)
      !
      ! -- storage coefficients
      sc1 = sscapacity(this%istor_coef, tp, bt, this%dis%area(n), this%ss(n))
      sc2 = sycapacity(this%dis%area(n), this%sy(n))
      rho1 = sc1 * tled
      rho2 = sc2 * tled
      !
      ! -- calculate newton terms for specific storage
      !    and specific yield
      if (this%iconvert(n) /= 0) then
        !
        ! -- calculate saturation derivative
        derv = squadraticsaturationderivative(tp, bt, h)
        !
        ! -- newton terms for specific storage
        if (this%iconf_ss == 0) then
          if (this%iorig_ss == 0) then
            drterm = -rho1 * derv * (h - bt) + rho1 * tthk * snnew * derv
          else
            drterm = -(rho1 * derv * h)
          end if
          call matrix_sln%add_value_pos(idxglo(idiag), drterm)
          rhs(n) = rhs(n) + drterm * h
        end if
        !
        ! -- newton terms for specific yield
        !    only calculated if the current saturation
        !    is less than one
        if (snnew < done) then
          ! -- calculate newton terms for specific yield
          if (snnew > dzero) then
            rterm = -rho2 * tthk * snnew
            drterm = -rho2 * tthk * derv
            call matrix_sln%add_value_pos(idxglo(idiag), drterm + rho2)
            rhs(n) = rhs(n) - rterm + drterm * h + rho2 * bt
          end if
        end if
      end if
    end do
    !
    ! -- return
    return

Here is the call graph for this function:

sto_rp()

subroutine gwfstomodule::sto_rp

(

class(gwfstotype)

this

)

Method to read and prepare stress period data for the STO package.

Parameters

this	GwfStoType object

Definition at line 157 of file gwf-sto.f90.

    ! -- modules
    use tdismodule, only: kper, nper
    implicit none
    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    ! -- local variables
    integer(I4B) :: ierr
    logical :: isfound, readss, readsy, endOfBlock
    character(len=16) :: css(0:1)
    character(len=LINELENGTH) :: line, keyword
    ! -- formats
    character(len=*), parameter :: fmtlsp = &
      &"(1X,/1X,'REUSING ',A,' FROM LAST STRESS PERIOD')"
    character(len=*), parameter :: fmtblkerr = &
      &"('Error.  Looking for BEGIN PERIOD iper.  Found ', a, ' instead.')"
    ! -- data
    data css(0)/'       TRANSIENT'/
    data css(1)/'    STEADY-STATE'/
! ------------------------------------------------------------------------------
    !
    ! -- Store ss and sy values from end of last stress period if needed
    if (this%integratechanges /= 0) then
      call this%save_old_ss_sy()
    end if
    !
    ! -- get stress period data
    if (this%ionper < kper) then
      !
      ! -- get period block
      call this%parser%GetBlock('PERIOD', isfound, ierr, &
                                supportopenclose=.true., &
                                blockrequired=.false.)
      if (isfound) then
        !
        ! -- read ionper and check for increasing period numbers
        call this%read_check_ionper()
      else
        !
        ! -- PERIOD block not found
        if (ierr < 0) then
          ! -- End of file found; data applies for remainder of simulation.
          this%ionper = nper + 1
        else
          ! -- Found invalid block
          call this%parser%GetCurrentLine(line)
          write (errmsg, fmtblkerr) adjustl(trim(line))
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
      end if
    end if
    !
    ! -- read data if ionper == kper
    ! are these here to anticipate reading ss,sy per stress period?
    readss = .false.
    readsy = .false.
 
    !stotxt = aname(2)
    if (this%ionper == kper) then
      write (this%iout, '(//,1x,a)') 'PROCESSING STORAGE PERIOD DATA'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('STEADY-STATE')
          this%iss = 1
        case ('TRANSIENT')
          this%iss = 0
        case default
          write (errmsg, '(a,a)') 'Unknown STORAGE data tag: ', &
            trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END PROCESSING STORAGE PERIOD DATA'
      !else
      !  write(this%iout,fmtlsp) 'STORAGE VALUES'
    end if
 
    write (this%iout, '(//1X,A,I0,A,A,/)') &
      'STRESS PERIOD ', kper, ' IS ', trim(adjustl(css(this%iss)))
    !
    ! -- TVS
    if (this%intvs /= 0) then
      call this%tvs%rp()
    end if
    !
    ! -- return
    return

Here is the call graph for this function:

sto_save_model_flows()

subroutine gwfstomodule::sto_save_model_flows	(	class(gwfstotype)	this,
		integer(i4b), intent(in)	icbcfl,
		integer(i4b), intent(in)	icbcun
	)

Save cell-by-cell budget terms for the STO package.

Parameters

	this	GwfStoType object
[in]	icbcfl	flag to output budget data
[in]	icbcun	cell-by-cell file unit number

Definition at line 669 of file gwf-sto.f90.

    ! -- dummy variables
    class(GwfStoType) :: this !< GwfStoType object
    integer(I4B), intent(in) :: icbcfl !< flag to output budget data
    integer(I4B), intent(in) :: icbcun !< cell-by-cell file unit number
    ! -- local variables
    integer(I4B) :: ibinun
    integer(I4B) :: iprint, nvaluesp, nwidthp
    character(len=1) :: cdatafmp = ' ', editdesc = ' '
    real(DP) :: dinact
    !
    ! -- Set unit number for binary output
    if (this%ipakcb < 0) then
      ibinun = icbcun
    elseif (this%ipakcb == 0) then
      ibinun = 0
    else
      ibinun = this%ipakcb
    end if
    if (icbcfl == 0) ibinun = 0
    !
    ! -- Record the storage rates if requested
    if (ibinun /= 0) then
      iprint = 0
      dinact = dzero
      !
      ! -- storage(ss)
      call this%dis%record_array(this%strgss, this%iout, iprint, -ibinun, &
                                 budtxt(1), cdatafmp, nvaluesp, &
                                 nwidthp, editdesc, dinact)
      !
      ! -- storage(sy)
      if (this%iusesy == 1) then
        call this%dis%record_array(this%strgsy, this%iout, iprint, -ibinun, &
                                   budtxt(2), cdatafmp, nvaluesp, &
                                   nwidthp, editdesc, dinact)
      end if
    end if
    !
    ! -- return
    return

Variable Documentation

budtxt

character(len=lenbudtxt), dimension(2) gwfstomodule::budtxt = [' STO-SS', ' STO-SY']

Definition at line 29 of file gwf-sto.f90.

  character(len=LENBUDTXT), dimension(2) :: budtxt = & !< text labels for budget terms
    &['          STO-SS', '          STO-SY']