gwtdspmodule Module Reference

Data Types
type	gwtdsptype

Functions/Subroutines
subroutine, public	dsp_cr (dspobj, name_model, input_mempath, inunit, iout, fmi)
	Create a DSP object. More...
subroutine	dsp_df (this, dis, dspOptions)
	Define DSP object. More...
subroutine	dsp_ac (this, moffset, sparse)
	Add connections to DSP. More...
subroutine	dsp_mc (this, moffset, matrix_sln)
	Map DSP connections. More...
subroutine	dsp_ar (this, ibound, thetam)
	Allocate and read method for package. More...
subroutine	dsp_ad (this)
	Advance method for the package. More...
subroutine	dsp_fc (this, kiter, nodes, nja, matrix_sln, idxglo, rhs, cnew)
	Fill coefficient method for package. More...
subroutine	dsp_cq (this, cnew, flowja)
	@ brief Calculate flows for package More...
subroutine	allocate_scalars (this)
	@ brief Allocate scalar variables for package More...
subroutine	allocate_arrays (this, nodes)
	@ brief Allocate arrays for package More...
subroutine	dsp_da (this)
	@ brief Deallocate memory More...
subroutine	log_options (this, found)
	Write user options to list file. More...
subroutine	source_options (this)
	Update simulation mempath options. More...
subroutine	log_griddata (this, found)
	Write dimensions to list file. More...
subroutine	source_griddata (this)
	Update DSP simulation data from input mempath. More...
subroutine	calcdispellipse (this)
	Calculate dispersion coefficients. More...
subroutine	calcdispcoef (this)
	Calculate dispersion coefficients. More...

Function/Subroutine Documentation

allocate_arrays()

subroutine gwtdspmodule::allocate_arrays	(	class(gwtdsptype)	this,
		integer(i4b), intent(in)	nodes
	)

Method to allocate arrays for the package.

Definition at line 409 of file gwt-dsp.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    use constantsmodule, only: dzero
    ! -- dummy
    class(GwtDspType) :: this
    integer(I4B), intent(in) :: nodes
    ! -- local
    !
    ! -- Allocate
    call mem_allocate(this%alh, nodes, 'ALH', trim(this%memoryPath))
    call mem_allocate(this%alv, nodes, 'ALV', trim(this%memoryPath))
    call mem_allocate(this%ath1, nodes, 'ATH1', trim(this%memoryPath))
    call mem_allocate(this%ath2, nodes, 'ATH2', trim(this%memoryPath))
    call mem_allocate(this%atv, nodes, 'ATV', trim(this%memoryPath))
    call mem_allocate(this%diffc, nodes, 'DIFFC', trim(this%memoryPath))
    call mem_allocate(this%d11, nodes, 'D11', trim(this%memoryPath))
    call mem_allocate(this%d22, nodes, 'D22', trim(this%memoryPath))
    call mem_allocate(this%d33, nodes, 'D33', trim(this%memoryPath))
    call mem_allocate(this%angle1, nodes, 'ANGLE1', trim(this%memoryPath))
    call mem_allocate(this%angle2, nodes, 'ANGLE2', trim(this%memoryPath))
    call mem_allocate(this%angle3, nodes, 'ANGLE3', trim(this%memoryPath))
    !
    ! -- Allocate dispersion coefficient array if xt3d not in use
    if (this%ixt3d == 0) then
      call mem_allocate(this%dispcoef, this%dis%njas, 'DISPCOEF', &
                        trim(this%memoryPath))
    else
      call mem_allocate(this%dispcoef, 0, 'DISPCOEF', trim(this%memoryPath))
    end if
    !
    ! -- Return
    return

allocate_scalars()

subroutine gwtdspmodule::allocate_scalars ( class(gwtdsptype)  this )

private

Method to allocate scalar variables for the package.

Definition at line 356 of file gwt-dsp.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    use constantsmodule, only: dzero
    ! -- dummy
    class(GwtDspType) :: this
    ! -- local
    !
    ! -- allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- Allocate
    call mem_allocate(this%idiffc, 'IDIFFC', this%memoryPath)
    call mem_allocate(this%idisp, 'IDISP', this%memoryPath)
    call mem_allocate(this%ialh, 'IALH', this%memoryPath)
    call mem_allocate(this%ialv, 'IALV', this%memoryPath)
    call mem_allocate(this%iath1, 'IATH1', this%memoryPath)
    call mem_allocate(this%iath2, 'IATH2', this%memoryPath)
    call mem_allocate(this%iatv, 'IATV', this%memoryPath)
    call mem_allocate(this%ixt3doff, 'IXT3DOFF', this%memoryPath)
    call mem_allocate(this%ixt3drhs, 'IXT3DRHS', this%memoryPath)
    call mem_allocate(this%ixt3d, 'IXT3D', this%memoryPath)
    call mem_allocate(this%id22, 'ID22', this%memoryPath)
    call mem_allocate(this%id33, 'ID33', this%memoryPath)
    call mem_allocate(this%iangle1, 'IANGLE1', this%memoryPath)
    call mem_allocate(this%iangle2, 'IANGLE2', this%memoryPath)
    call mem_allocate(this%iangle3, 'IANGLE3', this%memoryPath)
    !
    ! -- Initialize
    this%idiffc = 0
    this%idisp = 0
    this%ialh = 0
    this%ialv = 0
    this%iath1 = 0
    this%iath2 = 0
    this%iatv = 0
    this%ixt3doff = 0
    this%ixt3drhs = 0
    this%ixt3d = 0
    this%id22 = 1
    this%id33 = 1
    this%iangle1 = 1
    this%iangle2 = 1
    this%iangle3 = 1
    !
    ! -- Return
    return

calcdispcoef()

subroutine gwtdspmodule::calcdispcoef ( class(gwtdsptype)  this )

private

Definition at line 791 of file gwt-dsp.f90.

    ! -- modules
    use hgeoutilmodule, only: hyeff
    use gwfconductanceutilsmodule, only: staggered_thkfrac
    ! -- dummy
    class(GwtDspType) :: this
    ! -- local
    integer(I4B) :: nodes, n, m, idiag, ipos
    real(DP) :: clnm, clmn, dn, dm
    real(DP) :: vg1, vg2, vg3
    integer(I4B) :: ihc, isympos
    integer(I4B) :: iavgmeth
    real(DP) :: satn, satm, topn, topm, botn, botm
    real(DP) :: hwva, cond, cn, cm, denom
    real(DP) :: anm, amn, thksatn, thksatm
    !
    ! -- set iavgmeth = 1 to use arithmetic averaging for effective dispersion
    iavgmeth = 1
    !
    ! -- Process connections
    nodes = size(this%d11)
    do n = 1, nodes
      if (this%fmi%ibdgwfsat0(n) == 0) cycle
      idiag = this%dis%con%ia(n)
      do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
        !
        ! -- Set m to connected cell
        m = this%dis%con%ja(ipos)
        !
        ! -- skip for lower triangle
        if (m < n) cycle
        isympos = this%dis%con%jas(ipos)
        this%dispcoef(isympos) = dzero
        if (this%fmi%ibdgwfsat0(m) == 0) cycle
        !
        ! -- cell dimensions
        hwva = this%dis%con%hwva(isympos)
        clnm = this%dis%con%cl1(isympos)
        clmn = this%dis%con%cl2(isympos)
        ihc = this%dis%con%ihc(isympos)
        topn = this%dis%top(n)
        topm = this%dis%top(m)
        botn = this%dis%bot(n)
        botm = this%dis%bot(m)
        !
        ! -- flow model information
        satn = this%fmi%gwfsat(n)
        satm = this%fmi%gwfsat(m)
        !
        ! -- Calculate dispersion coefficient for cell n in the direction
        !    normal to the shared n-m face and for cell m in the direction
        !    normal to the shared n-m face.
        call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, ipos)
        dn = hyeff(this%d11(n), this%d22(n), this%d33(n), &
                   this%angle1(n), this%angle2(n), this%angle3(n), &
                   vg1, vg2, vg3, iavgmeth)
        dm = hyeff(this%d11(m), this%d22(m), this%d33(m), &
                   this%angle1(m), this%angle2(m), this%angle3(m), &
                   vg1, vg2, vg3, iavgmeth)
        !
        ! -- Calculate dispersion conductance based on NPF subroutines and the
        !    effective dispersion coefficients dn and dm.
        if (ihc == 0) then
          clnm = satn * (topn - botn) * dhalf
          clmn = satm * (topm - botm) * dhalf
          anm = hwva
          !
          ! -- n is convertible and unsaturated
          if (satn == dzero) then
            anm = dzero
          else if (n > m .and. satn < done) then
            anm = dzero
          end if
          !
          ! -- m is convertible and unsaturated
          if (satm == dzero) then
            anm = dzero
          else if (m > n .and. satm < done) then
            anm = dzero
          end if
          !
          ! -- amn is the same as anm for vertical flow
          amn = anm
          !
        else
          !
          ! -- horizontal conductance
          !
          ! -- handle vertically staggered case
          if (ihc == 2) then
            thksatn = staggered_thkfrac(topn, botn, satn, topm, botm)
            thksatm = staggered_thkfrac(topm, botm, satm, topn, botn)
          else
            thksatn = (topn - botn) * satn
            thksatm = (topm - botm) * satm
          end if
          !
          ! -- calculate the saturated area term
          anm = thksatn * hwva
          amn = thksatm * hwva
          !
          ! -- n or m is unsaturated, so no dispersion
          if (satn == dzero .or. satm == dzero) then
            anm = dzero
            amn = dzero
          end if
          !
        end if
        !
        ! -- calculate conductance using the two half cell conductances
        cn = dzero
        if (clnm > dzero) cn = dn * anm / clnm
        cm = dzero
        if (clmn > dzero) cm = dm * amn / clmn
        denom = cn + cm
        if (denom > dzero) then
          cond = cn * cm / denom
        else
          cond = dzero
        end if
        !
        ! -- Assign the calculated dispersion conductance
        this%dispcoef(isympos) = cond
        !
      end do
    end do
    !
    ! -- Return
    return

Here is the call graph for this function:

calcdispellipse()

subroutine gwtdspmodule::calcdispellipse

(

class(gwtdsptype)

this

)

Definition at line 677 of file gwt-dsp.f90.

    ! -- modules
    ! -- dummy
    class(GwtDspType) :: this
    ! -- local
    integer(I4B) :: nodes, n
    real(DP) :: q, qx, qy, qz
    real(DP) :: alh, alv, ath1, ath2, atv, a
    real(DP) :: al, at1, at2
    real(DP) :: qzoqsquared
    real(DP) :: dstar
    !
    ! -- loop through and calculate dispersion coefficients and angles
    nodes = size(this%d11)
    do n = 1, nodes
      !
      ! -- initialize
      this%d11(n) = dzero
      this%d22(n) = dzero
      this%d33(n) = dzero
      this%angle1(n) = dzero
      this%angle2(n) = dzero
      this%angle3(n) = dzero
      if (this%fmi%ibdgwfsat0(n) == 0) cycle
      !
      ! -- specific discharge
      qx = dzero
      qy = dzero
      qz = dzero
      q = dzero
      qx = this%fmi%gwfspdis(1, n)
      qy = this%fmi%gwfspdis(2, n)
      qz = this%fmi%gwfspdis(3, n)
      q = qx**2 + qy**2 + qz**2
      if (q > dzero) q = sqrt(q)
      !
      ! -- dispersion coefficients
      alh = dzero
      alv = dzero
      ath1 = dzero
      ath2 = dzero
      atv = dzero
      if (this%idisp > 0) then
        alh = this%alh(n)
        alv = this%alv(n)
        ath1 = this%ath1(n)
        ath2 = this%ath2(n)
        atv = this%atv(n)
      end if
      dstar = dzero
      if (this%idiffc > 0) then
        ! -- Multiply diffusion coefficient by mobile porosity, defined
        !    as volume of voids in the mobile domain per aquifer volume
        dstar = this%diffc(n) * this%thetam(n)
      end if
      !
      ! -- Calculate the longitudal and transverse dispersivities
      al = dzero
      at1 = dzero
      at2 = dzero
      if (q > dzero) then
        qzoqsquared = (qz / q)**2
        al = alh * (done - qzoqsquared) + alv * qzoqsquared
        at1 = ath1 * (done - qzoqsquared) + atv * qzoqsquared
        at2 = ath2 * (done - qzoqsquared) + atv * qzoqsquared
      end if
      !
      ! -- Calculate and save the diagonal components of the dispersion tensor
      this%d11(n) = al * q + dstar
      this%d22(n) = at1 * q + dstar
      this%d33(n) = at2 * q + dstar
      !
      ! -- Angles of rotation if velocity based dispersion tensor
      if (this%idisp > 0) then
        !
        ! -- angles of rotation from model coordinates to direction of velocity
        ! qx / q = cos(a1) * cos(a2)
        ! qy / q = sin(a1) * cos(a2)
        ! qz / q = sin(a2)
        !
        ! -- angle3 is zero
        this%angle3(n) = dzero
        !
        ! -- angle2
        a = dzero
        if (q > dzero) a = qz / q
        this%angle2(n) = asin(a)
        !
        ! -- angle1
        a = q * cos(this%angle2(n))
        if (a /= dzero) then
          a = qx / a
        else
          a = dzero
        end if
        !
        ! -- acos(1) not defined, so set to zero if necessary
        if (a <= -done) then
          this%angle1(n) = dpi
        elseif (a >= done) then
          this%angle1(n) = dzero
        else
          this%angle1(n) = acos(a)
        end if
        !
      end if
    end do
    !
    ! -- Return
    return

dsp_ac()

subroutine gwtdspmodule::dsp_ac ( class(gwtdsptype)  this, integer(i4b), intent(in)  moffset, type(sparsematrix), intent(inout)  sparse  )

private

Add connections for extended neighbors to the sparse matrix

Definition at line 170 of file gwt-dsp.f90.

    ! -- modules
    use sparsemodule, only: sparsematrix
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(GwtDspType) :: this
    integer(I4B), intent(in) :: moffset
    type(sparsematrix), intent(inout) :: sparse
    ! -- local
    !
    ! -- Add extended neighbors (neighbors of neighbors)
    if (this%ixt3d > 0) call this%xt3d%xt3d_ac(moffset, sparse)
    !
    ! -- Return
    return

dsp_ad()

subroutine gwtdspmodule::dsp_ad ( class(gwtdsptype)  this )

private

Definition at line 233 of file gwt-dsp.f90.

    ! -- modules
    use tdismodule, only: kstp, kper
    ! -- dummy
    class(GwtDspType) :: this
    ! -- local
    !
    ! -- xt3d
    ! TODO: might consider adding a new mf6 level set pointers method, and
    ! doing this stuff there instead of in the time step loop.
    if (kstp * kper == 1) then
      if (this%ixt3d > 0) then
        call this%xt3d%xt3d_ar(this%fmi%ibdgwfsat0, this%d11, this%id33, &
                               this%d33, this%fmi%gwfsat, this%id22, this%d22, &
                               this%iangle1, this%iangle2, this%iangle3, &
                               this%angle1, this%angle2, this%angle3)
      end if
    end if
    !
    ! -- Fill d11, d22, d33, angle1, angle2, angle3 using specific discharge
    call this%calcdispellipse()
    !
    ! -- Recalculate dispersion coefficients if the flows were updated
    if (this%fmi%iflowsupdated == 1) then
      if (this%ixt3d == 0) then
        call this%calcdispcoef()
      else if (this%ixt3d > 0) then
        call this%xt3d%xt3d_fcpc(this%dis%nodes, .false.)
      end if
    end if
    !
    ! -- Return
    return

dsp_ar()

subroutine gwtdspmodule::dsp_ar	(	class(gwtdsptype)	this,
		integer(i4b), dimension(:), pointer, contiguous	ibound,
		real(dp), dimension(:), pointer, contiguous	thetam
	)

Method to allocate and read static data for the package.

Definition at line 211 of file gwt-dsp.f90.

    ! -- modules
    ! -- dummy
    class(GwtDspType) :: this
    integer(I4B), dimension(:), pointer, contiguous :: ibound
    real(DP), dimension(:), pointer, contiguous :: thetam
    ! -- local
    ! -- formats
    character(len=*), parameter :: fmtdsp = &
      "(1x,/1x,'DSP-- DISPERSION PACKAGE, VERSION 1, 1/24/2018', &
      &' INPUT READ FROM UNIT ', i0, //)"
    !
    ! -- dsp pointers to arguments that were passed in
    this%ibound => ibound
    this%thetam => thetam
    !
    ! -- Return
    return

dsp_cq()

subroutine gwtdspmodule::dsp_cq ( class(gwtdsptype)  this, real(dp), dimension(:), intent(inout)  cnew, real(dp), dimension(:), intent(inout)  flowja  )

private

Method to calculate dispersion contribution to flowja

Definition at line 322 of file gwt-dsp.f90.

    ! -- modules
    ! -- dummy
    class(GwtDspType) :: this
    real(DP), intent(inout), dimension(:) :: cnew
    real(DP), intent(inout), dimension(:) :: flowja
    ! -- local
    integer(I4B) :: n, m, ipos, isympos
    real(DP) :: dnm
    !
    ! -- Calculate dispersion and add to flowja
    if (this%ixt3d > 0) then
      call this%xt3d%xt3d_flowja(cnew, flowja)
    else
      do n = 1, this%dis%nodes
        if (this%fmi%ibdgwfsat0(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%fmi%ibdgwfsat0(m) == 0) cycle
          isympos = this%dis%con%jas(ipos)
          dnm = this%dispcoef(isympos)
          flowja(ipos) = flowja(ipos) + dnm * (cnew(m) - cnew(n))
        end do
      end do
    end if
    !
    ! -- Return
    return

dsp_cr()

subroutine, public gwtdspmodule::dsp_cr	(	type(gwtdsptype), pointer	dspobj,
		character(len=*), intent(in)	name_model,
		character(len=*), intent(in)	input_mempath,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout,
		type(tspfmitype), intent(in), target	fmi
	)

Definition at line 77 of file gwt-dsp.f90.

    ! -- modules
    use kindmodule, only: lgp
    use memorymanagerextmodule, only: mem_set_value
    ! -- dummy
    type(GwtDspType), pointer :: dspobj
    character(len=*), intent(in) :: name_model
    character(len=*), intent(in) :: input_mempath
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    type(TspFmiType), intent(in), target :: fmi
    ! -- locals
    ! -- formats
    character(len=*), parameter :: fmtdsp = &
      "(1x,/1x,'DSP-- DISPERSION PACKAGE, VERSION 1, 1/24/2018', &
      &' INPUT READ FROM MEMPATH: ', A, //)"
    !
    ! -- Create the object
    allocate (dspobj)
    !
    ! -- create name and memory path
    call dspobj%set_names(1, name_model, 'DSP', 'DSP', input_mempath)
    !
    ! -- Allocate scalars
    call dspobj%allocate_scalars()
    !
    ! -- Set variables
    dspobj%inunit = inunit
    dspobj%iout = iout
    dspobj%fmi => fmi
    !
    if (dspobj%inunit > 0) then
      !
      ! -- Print a message identifying the dispersion package.
      if (dspobj%iout > 0) then
        write (dspobj%iout, fmtdsp) input_mempath
      end if
    end if
    !
    ! -- Return
    return

Here is the caller graph for this function:

dsp_da()

subroutine gwtdspmodule::dsp_da

(

class(gwtdsptype)

this

)

Method to deallocate memory for the package.

Definition at line 448 of file gwt-dsp.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    use memorymanagerextmodule, only: memorylist_remove
    use simvariablesmodule, only: idm_context
    ! -- dummy
    class(GwtDspType) :: this
    ! -- local
    !
    ! -- Deallocate input memory
    call memorylist_remove(this%name_model, 'DSP', idm_context)
    !
    ! -- deallocate arrays
    if (this%inunit /= 0) then
      call mem_deallocate(this%alh)
      call mem_deallocate(this%alv, 'ALV', trim(this%memoryPath))
      call mem_deallocate(this%ath1)
      call mem_deallocate(this%ath2, 'ATH2', trim(this%memoryPath))
      call mem_deallocate(this%atv, 'ATV', trim(this%memoryPath))
      call mem_deallocate(this%diffc)
      call mem_deallocate(this%d11)
      call mem_deallocate(this%d22)
      call mem_deallocate(this%d33)
      call mem_deallocate(this%angle1)
      call mem_deallocate(this%angle2)
      call mem_deallocate(this%angle3)
      call mem_deallocate(this%dispcoef)
      if (this%ixt3d > 0) call this%xt3d%xt3d_da()
    end if
    !
    ! -- deallocate objects
    if (this%ixt3d > 0) deallocate (this%xt3d)
    !
    ! -- deallocate scalars
    call mem_deallocate(this%idiffc)
    call mem_deallocate(this%idisp)
    call mem_deallocate(this%ialh)
    call mem_deallocate(this%ialv)
    call mem_deallocate(this%iath1)
    call mem_deallocate(this%iath2)
    call mem_deallocate(this%iatv)
    call mem_deallocate(this%ixt3doff)
    call mem_deallocate(this%ixt3drhs)
    call mem_deallocate(this%ixt3d)
    call mem_deallocate(this%id22)
    call mem_deallocate(this%id33)
    call mem_deallocate(this%iangle1)
    call mem_deallocate(this%iangle2)
    call mem_deallocate(this%iangle3)
    !
    ! -- deallocate variables in NumericalPackageType
    call this%NumericalPackageType%da()
    !
    ! -- nullify pointers
    this%thetam => null()
    !
    ! -- Return
    return

Here is the call graph for this function:

dsp_df()

subroutine gwtdspmodule::dsp_df	(	class(gwtdsptype)	this,
		class(disbasetype), pointer	dis,
		type(gwtdspoptionstype), intent(in), optional	dspOptions
	)

Parameters

[in]

dspoptions

the optional DSP options, used when not creating DSP from file

Definition at line 122 of file gwt-dsp.f90.

    ! -- modules
    ! -- dummy
    class(GwtDspType) :: this
    class(DisBaseType), pointer :: dis
    type(GwtDspOptionsType), optional, intent(in) :: dspOptions !< the optional DSP options, used when not
                                                                !! creating DSP from file
    ! -- local
    !
    ! -- Store pointer to dis
    this%dis => dis
    !
    !
    ! -- set default xt3d representation to on and lhs
    this%ixt3d = 1
    !
    ! -- Read dispersion options
    if (present(dspoptions)) then
      this%ixt3d = dspoptions%ixt3d
      !
      ! -- Allocate only, grid data will not be read from file
      call this%allocate_arrays(this%dis%nodes)
    else
      !
      ! -- Source options
      call this%source_options()
      call this%allocate_arrays(this%dis%nodes)
      !
      ! -- Source dispersion data
      call this%source_griddata()
    end if
    !
    ! -- xt3d create
    if (this%ixt3d > 0) then
      call xt3d_cr(this%xt3d, this%name_model, this%inunit, this%iout, &
                   ldispopt=.true.)
      this%xt3d%ixt3d = this%ixt3d
      call this%xt3d%xt3d_df(dis)
    end if
    !
    ! -- Return
    return

Here is the call graph for this function:

dsp_fc()

subroutine gwtdspmodule::dsp_fc	(	class(gwtdsptype)	this,
		integer(i4b)	kiter,
		integer(i4b), intent(in)	nodes,
		integer(i4b), intent(in)	nja,
		class(matrixbasetype), pointer	matrix_sln,
		integer(i4b), dimension(nja), intent(in)	idxglo,
		real(dp), dimension(nodes), intent(inout)	rhs,
		real(dp), dimension(nodes), intent(inout)	cnew
	)

Method to calculate and fill coefficients for the package.

Definition at line 272 of file gwt-dsp.f90.

    ! -- modules
    ! -- dummy
    class(GwtDspType) :: this
    integer(I4B) :: kiter
    integer(I4B), intent(in) :: nodes
    integer(I4B), intent(in) :: nja
    class(MatrixBaseType), pointer :: matrix_sln
    integer(I4B), intent(in), dimension(nja) :: idxglo
    real(DP), intent(inout), dimension(nodes) :: rhs
    real(DP), intent(inout), dimension(nodes) :: cnew
    ! -- local
    integer(I4B) :: n, m, idiag, idiagm, ipos, isympos, isymcon
    real(DP) :: dnm
    !
    if (this%ixt3d > 0) then
      call this%xt3d%xt3d_fc(kiter, matrix_sln, idxglo, rhs, cnew)
    else
      do n = 1, nodes
        if (this%fmi%ibdgwfsat0(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 (m < n) cycle
          if (this%fmi%ibdgwfsat0(m) == 0) cycle
          isympos = this%dis%con%jas(ipos)
          dnm = this%dispcoef(isympos)
          !
          ! -- Contribution to row n
          call matrix_sln%add_value_pos(idxglo(ipos), dnm)
          call matrix_sln%add_value_pos(idxglo(idiag), -dnm)
          !
          ! -- Contribution to row m
          idiagm = this%dis%con%ia(m)
          isymcon = this%dis%con%isym(ipos)
          call matrix_sln%add_value_pos(idxglo(isymcon), dnm)
          call matrix_sln%add_value_pos(idxglo(idiagm), -dnm)
        end do
      end do
    end if
    !
    ! -- Return
    return

dsp_mc()

subroutine gwtdspmodule::dsp_mc	(	class(gwtdsptype)	this,
		integer(i4b), intent(in)	moffset,
		class(matrixbasetype), pointer	matrix_sln
	)

Map connections and construct iax, jax, and idxglox

Definition at line 191 of file gwt-dsp.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(GwtDspType) :: this
    integer(I4B), intent(in) :: moffset
    class(MatrixBaseType), pointer :: matrix_sln
    ! -- local
    !
    ! -- Call xt3d map connections
    if (this%ixt3d > 0) call this%xt3d%xt3d_mc(moffset, matrix_sln)
    !
    ! -- Return
    return

log_griddata()

subroutine gwtdspmodule::log_griddata	(	class(gwtdsptype)	this,
		type(gwtdspparamfoundtype), intent(in)	found
	)

Definition at line 556 of file gwt-dsp.f90.

    use gwtdspinputmodule, only: gwtdspparamfoundtype
    class(GwtDspType) :: this
    type(GwtDspParamFoundType), intent(in) :: found
 
    write (this%iout, '(1x,a)') 'Setting DSP Griddata'
 
    if (found%diffc) then
      write (this%iout, '(4x,a)') 'DIFFC set from input file'
    end if
 
    if (found%alh) then
      write (this%iout, '(4x,a)') 'ALH set from input file'
    end if
 
    if (found%alv) then
      write (this%iout, '(4x,a)') 'ALV set from input file'
    end if
 
    if (found%ath1) then
      write (this%iout, '(4x,a)') 'ATH1 set from input file'
    end if
 
    if (found%ath2) then
      write (this%iout, '(4x,a)') 'ATH2 set from input file'
    end if
 
    if (found%atv) then
      write (this%iout, '(4x,a)') 'ATV set from input file'
    end if
 
    write (this%iout, '(1x,a,/)') 'End Setting DSP Griddata'
 

log_options()

subroutine gwtdspmodule::log_options	(	class(gwtdsptype)	this,
		type(gwtdspparamfoundtype), intent(in)	found
	)

Definition at line 510 of file gwt-dsp.f90.

    use gwtdspinputmodule, only: gwtdspparamfoundtype
    class(GwTDspType) :: this
    type(GwtDspParamFoundType), intent(in) :: found
 
    write (this%iout, '(1x,a)') 'Setting DSP Options'
    write (this%iout, '(4x,a,i0)') 'XT3D formulation [0=INACTIVE, 1=ACTIVE, &
                                   &3=ACTIVE RHS] set to: ', this%ixt3d
    write (this%iout, '(1x,a,/)') 'End Setting DSP Options'
    ! -- Return
    return

source_griddata()

subroutine gwtdspmodule::source_griddata

(

class(gwtdsptype)

this

)

Definition at line 593 of file gwt-dsp.f90.

    ! -- modules
    use simmodule, only: count_errors, store_error
    use memorymanagermodule, only: mem_reallocate, mem_reassignptr
    use memorymanagerextmodule, only: mem_set_value
    use constantsmodule, only: linelength
    use gwtdspinputmodule, only: gwtdspparamfoundtype
    ! -- dummy
    class(GwtDsptype) :: this
    ! -- locals
    character(len=LINELENGTH) :: errmsg
    type(GwtDspParamFoundType) :: found
    integer(I4B), dimension(:), pointer, contiguous :: map
    ! -- formats
    !
    ! -- set map
    map => null()
    if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%diffc, 'DIFFC', this%input_mempath, map, found%diffc)
    call mem_set_value(this%alh, 'ALH', this%input_mempath, map, found%alh)
    call mem_set_value(this%alv, 'ALV', this%input_mempath, map, found%alv)
    call mem_set_value(this%ath1, 'ATH1', this%input_mempath, map, found%ath1)
    call mem_set_value(this%ath2, 'ATH2', this%input_mempath, map, found%ath2)
    call mem_set_value(this%atv, 'ATV', this%input_mempath, map, found%atv)
    !
    ! -- set active flags
    if (found%diffc) this%idiffc = 1
    if (found%alh) this%ialh = 1
    if (found%alv) this%ialv = 1
    if (found%ath1) this%iath1 = 1
    if (found%ath2) this%iath2 = 1
    if (found%atv) this%iatv = 1
    !
    ! -- reallocate diffc if not found
    if (.not. found%diffc) then
      call mem_reallocate(this%diffc, 0, 'DIFFC', trim(this%memoryPath))
    end if
    !
    ! -- set this%idisp flag
    if (found%alh) this%idisp = this%idisp + 1
    if (found%alv) this%idisp = this%idisp + 1
    if (found%ath1) this%idisp = this%idisp + 1
    if (found%ath2) this%idisp = this%idisp + 1
    !
    ! -- manage dispersion arrays
    if (this%idisp > 0) then
      if (.not. (found%alh .and. found%ath1)) then
        write (errmsg, '(a)') &
          'if dispersivities are specified then ALH and ATH1 are required.'
        call store_error(errmsg)
      end if
      ! -- If alv not specified then point it to alh
      if (.not. found%alv) &
        call mem_reassignptr(this%alv, 'ALV', trim(this%memoryPath), &
                             'ALH', trim(this%memoryPath))
      ! -- If ath2 not specified then point it to ath1
      if (.not. found%ath2) &
        call mem_reassignptr(this%ath2, 'ATH2', trim(this%memoryPath), &
                             'ATH1', trim(this%memoryPath))
      ! -- If atv not specified then point it to ath2
      if (.not. found%atv) &
        call mem_reassignptr(this%atv, 'ATV', trim(this%memoryPath), &
                             'ATH2', trim(this%memoryPath))
    else
      call mem_reallocate(this%alh, 0, 'ALH', trim(this%memoryPath))
      call mem_reallocate(this%alv, 0, 'ALV', trim(this%memoryPath))
      call mem_reallocate(this%ath1, 0, 'ATH1', trim(this%memoryPath))
      call mem_reallocate(this%ath2, 0, 'ATH2', trim(this%memoryPath))
      call mem_reallocate(this%atv, 0, 'ATV', trim(this%memoryPath))
    end if
    !
    ! -- log griddata
    if (this%iout > 0) then
      call this%log_griddata(found)
    end if
    !
    ! -- Return
    return

Here is the call graph for this function:

source_options()

subroutine gwtdspmodule::source_options

(

class(gwtdsptype)

this

)

Definition at line 525 of file gwt-dsp.f90.

    ! -- modules
    !use KindModule, only: LGP
    use memorymanagerextmodule, only: mem_set_value
    use gwtdspinputmodule, only: gwtdspparamfoundtype
    ! -- dummy
    class(GwtDspType) :: this
    ! -- locals
    type(GwtDspParamFoundType) :: found
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%ixt3doff, 'XT3D_OFF', this%input_mempath, &
                       found%xt3d_off)
    call mem_set_value(this%ixt3drhs, 'XT3D_RHS', this%input_mempath, &
                       found%xt3d_rhs)
    !
    ! -- set xt3d state flag
    if (found%xt3d_off) this%ixt3d = 0
    if (found%xt3d_rhs) this%ixt3d = 2
    !
    ! -- log options
    if (this%iout > 0) then
      call this%log_options(found)
    end if
    !
    ! -- Return
    return