gwecndmodule Module Reference

Data Types
type	gwecndtype

Functions/Subroutines
subroutine, public	cnd_cr (cndobj, name_model, input_mempath, inunit, iout, fmi, eqnsclfac, gwecommon)
	Create a new CND object. More...
subroutine	cnd_df (this, dis, cndOptions)
	Define CND object. More...
subroutine	cnd_ac (this, moffset, sparse)
	Add connections to CND. More...
subroutine	cnd_mc (this, moffset, matrix_sln)
	Map CND connections. More...
subroutine	cnd_ar (this, ibound, porosity)
	Allocate and read method for package. More...
subroutine	cnd_ad (this)
	Advance method for the package. More...
subroutine	cnd_fc (this, kiter, nodes, nja, matrix_sln, idxglo, rhs, cnew)
	Fill coefficient method for package. More...
subroutine	cnd_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	cnd_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 CND 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 gwecndmodule::allocate_arrays	(	class(gwecndtype)	this,
		integer(i4b), intent(in)	nodes
	)

Method to allocate arrays for the package.

Definition at line 420 of file gwe-cnd.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    use constantsmodule, only: dzero
    ! -- dummy
    class(GweCndType) :: this
    integer(I4B), intent(in) :: nodes
    !
    ! -- 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%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))
    call mem_allocate(this%ktw, nodes, 'KTW', trim(this%memoryPath))
    call mem_allocate(this%kts, nodes, 'KTS', 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 gwecndmodule::allocate_scalars ( class(gwecndtype)  this )

private

Method to allocate scalar variables for the package.

Definition at line 366 of file gwe-cnd.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    use constantsmodule, only: dzero
    ! -- dummy
    class(GweCndType) :: this
    !
    ! -- allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- Allocate
    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)
    call mem_allocate(this%iktw, 'IKTW', this%memoryPath)
    call mem_allocate(this%ikts, 'IKTS', this%memoryPath)
    !
    ! -- Initialize
    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
    this%iktw = 1
    this%ikts = 1
    !
    ! -- Return
    return

calcdispcoef()

subroutine gwecndmodule::calcdispcoef ( class(gwecndtype)  this )

private

Definition at line 810 of file gwe-cnd.f90.

    ! -- modules
    use hgeoutilmodule, only: hyeff
    use gwfconductanceutilsmodule, only: staggered_thkfrac
    ! -- dummy
    class(GweCndType) :: 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%ibdgwfsat0(n)
        satm = this%fmi%ibdgwfsat0(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 gwecndmodule::calcdispellipse

(

class(gwecndtype)

this

)

Definition at line 690 of file gwe-cnd.f90.

    ! -- modules
    ! -- dummy
    class(GweCndType) :: 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) :: ktbulk ! TODO: Implement additional options for characterizing ktbulk (see Markle refs)
    real(DP) :: dstar
    real(DP) :: qsw
    !
    ! -- 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
      !
      ! -- calculate
      ktbulk = dzero
      if (this%iktw > 0) ktbulk = ktbulk + this%porosity(n) * this%ktw(n) * &
                                  this%fmi%gwfsat(n)
      if (this%ikts > 0) ktbulk = ktbulk + (done - this%porosity(n)) * this%kts(n)
      !
      ! -- The division by rhow*cpw below is undertaken to render dstar in the
      !    form of a thermal diffusivity, and not because the governing
      !    equation is scaled by rhow*cpw. Because of this conceptual
      !    distinction, ktbulk is divided by the explicitly calculated product
      !    rhow*cpw, and not by the equivalent scale factor eqnsclfac, even
      !    though it should make no practical difference in the result.
      dstar = ktbulk / (this%gwecommon%gwecpw * this%gwecommon%gwerhow)
      !
      ! -- 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
      qsw = q * this%fmi%gwfsat(n) * this%eqnsclfac
      this%d11(n) = al * qsw + ktbulk
      this%d22(n) = at1 * qsw + ktbulk
      this%d33(n) = at2 * qsw + ktbulk
      !
      ! -- Angles of rotation if velocity based dispersion tensor
      if (this%idisp > 0) then
        !
        ! -- 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

cnd_ac()

subroutine gwecndmodule::cnd_ac ( class(gwecndtype)  this, integer(i4b), intent(in)  moffset, type(sparsematrix), intent(inout)  sparse  )

private

Add connections for extended neighbors to the sparse matrix

Definition at line 182 of file gwe-cnd.f90.

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

cnd_ad()

subroutine gwecndmodule::cnd_ad ( class(gwecndtype)  this )

private

Definition at line 243 of file gwe-cnd.f90.

    ! -- modules
    use tdismodule, only: kstp, kper
    ! -- dummy
    class(GweCndType) :: this
    !
    ! -- 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, .true.)
      end if
    end if
    !
    ! -- Return
    return

cnd_ar()

subroutine gwecndmodule::cnd_ar	(	class(gwecndtype)	this,
		integer(i4b), dimension(:), pointer, contiguous	ibound,
		real(dp), dimension(:), pointer, contiguous	porosity
	)

Method to allocate and read static data for the package.

Definition at line 221 of file gwe-cnd.f90.

    ! -- modules
    ! -- dummy
    class(GweCndType) :: this
    integer(I4B), dimension(:), pointer, contiguous :: ibound
    real(DP), dimension(:), pointer, contiguous :: porosity
    ! -- local
    ! -- formats
    character(len=*), parameter :: fmtcnd = &
      "(1x,/1x,'CND-- THERMAL CONDUCTION AND DISPERSION PACKAGE, VERSION 1, ', &
      &'5/01/2023, INPUT READ FROM UNIT ', i0, //)"
    !
    ! -- cnd pointers to arguments that were passed in
    this%ibound => ibound
    this%porosity => porosity
    !
    ! -- Return
    return

cnd_cq()

subroutine gwecndmodule::cnd_cq ( class(gwecndtype)  this, real(dp), dimension(:), intent(inout)  cnew, real(dp), dimension(:), intent(inout)  flowja  )

private

Method to calculate dispersion contribution to flowja

Definition at line 330 of file gwe-cnd.f90.

    ! -- modules
    ! -- dummy
    class(GweCndType) :: this
    real(DP), intent(inout), dimension(:) :: cnew
    real(DP), intent(inout), dimension(:) :: flowja
    ! -- local
    integer(I4B) :: n, m, ipos, isympos
    real(DP) :: dnm, qnm
    !
    ! -- 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)
!!          qnm = dnm * (cnew(m) - cnew(n)) * this%eqnsclfac
          qnm = dnm * (cnew(m) - cnew(n))
          flowja(ipos) = flowja(ipos) + qnm
        end do
      end do
    end if
    !
    ! -- Return
    return

cnd_cr()

subroutine, public gwecndmodule::cnd_cr	(	type(gwecndtype), pointer	cndobj,
		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,
		real(dp), intent(in), pointer	eqnsclfac,
		type(gweinputdatatype), intent(in), optional, target	gwecommon
	)

Create a new CND package

Parameters

[in]	eqnsclfac	governing equation scale factor
[in]	gwecommon	shared data container for use by multiple GWE packages

Definition at line 84 of file gwe-cnd.f90.

    ! -- modules
    use kindmodule, only: lgp
    use memorymanagerextmodule, only: mem_set_value
    ! -- dummy
    type(GweCndType), pointer :: cndobj
    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
    real(DP), intent(in), pointer :: eqnsclfac !< governing equation scale factor
    type(GweInputDataType), intent(in), target, optional :: gwecommon !< shared data container for use by multiple GWE packages
    ! -- formats
    character(len=*), parameter :: fmtcnd = &
      "(1x,/1x,'CND-- THERMAL CONDUCTION AND DISPERSION PACKAGE, VERSION 1, ', &
      &'01/01/2024, INPUT READ FROM MEMPATH ', A, //)"
    !
    ! -- Create the object
    allocate (cndobj)
    !
    ! -- create name and memory path
    call cndobj%set_names(1, name_model, 'CND', 'CND', input_mempath)
    !
    ! -- Allocate scalars
    call cndobj%allocate_scalars()
    !
    ! -- Set variables
    cndobj%inunit = inunit
    cndobj%iout = iout
    cndobj%fmi => fmi
    cndobj%eqnsclfac => eqnsclfac
    if (present(gwecommon)) then
      cndobj%gwecommon => gwecommon
    end if
    !
    if (cndobj%inunit > 0) then
      !
      ! -- Print a message identifying the dispersion package.
      if (cndobj%iout > 0) then
        write (cndobj%iout, fmtcnd) input_mempath
      end if
    end if
    !
    ! -- Return
    return

Here is the caller graph for this function:

cnd_da()

subroutine gwecndmodule::cnd_da

(

class(gwecndtype)

this

)

Method to deallocate memory for the package.

Definition at line 459 of file gwe-cnd.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    use memorymanagerextmodule, only: memorylist_remove
    use simvariablesmodule, only: idm_context
    ! -- dummy
    class(GweCndType) :: this
    ! -- local
    !
    ! -- Deallocate input memory
    call memorylist_remove(this%name_model, 'CND', 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%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%ktw)
      call mem_deallocate(this%kts)
      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)
    nullify (this%gwecommon)
    !
    ! -- deallocate scalars
    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)
    call mem_deallocate(this%iktw)
    call mem_deallocate(this%ikts)
    !
    ! -- deallocate variables in NumericalPackageType
    call this%NumericalPackageType%da()
    !
    ! -- Return
    return

Here is the call graph for this function:

cnd_df()

subroutine gwecndmodule::cnd_df	(	class(gwecndtype)	this,
		class(disbasetype), pointer	dis,
		type(gwecndoptionstype), intent(in), optional	cndOptions
	)

Parameters

[in]

cndoptions

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

Definition at line 135 of file gwe-cnd.f90.

    ! -- modules
    ! -- dummy
    class(GweCndType) :: this
    class(DisBaseType), pointer :: dis
    type(GweCndOptionsType), optional, intent(in) :: cndOptions !< the optional CND options, used when not
                                                                !! creating CND from file
    !
    ! -- Store pointer to dis
    this%dis => dis
    !
    !
    ! -- set default xt3d representation to on and lhs
    this%ixt3d = 1
    !
    ! -- Read dispersion options
    if (present(cndoptions)) then
      this%ixt3d = cndoptions%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:

cnd_fc()

subroutine gwecndmodule::cnd_fc	(	class(gwecndtype)	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 281 of file gwe-cnd.f90.

    ! -- dummy
    class(GweCndType) :: 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

cnd_mc()

subroutine gwecndmodule::cnd_mc	(	class(gwecndtype)	this,
		integer(i4b), intent(in)	moffset,
		class(matrixbasetype), pointer	matrix_sln
	)

Map connections and construct iax, jax, and idxglox

Definition at line 202 of file gwe-cnd.f90.

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

log_griddata()

subroutine gwecndmodule::log_griddata	(	class(gwecndtype)	this,
		type(gwecndparamfoundtype), intent(in)	found
	)

Definition at line 566 of file gwe-cnd.f90.

    use gwecndinputmodule, only: gwecndparamfoundtype
    class(GweCndType) :: this
    type(GweCndParamFoundType), intent(in) :: found
    !
    write (this%iout, '(1x,a)') 'Setting CND Griddata'
    !
    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
    !
    if (found%ktw) then
      write (this%iout, '(4x,a)') 'KTW set from input file'
    end if
    !
    if (found%kts) then
      write (this%iout, '(4x,a)') 'KTS set from input file'
    end if
    !
    write (this%iout, '(1x,a,/)') 'End Setting CND Griddata'
    !
    ! -- Return
    return

log_options()

subroutine gwecndmodule::log_options	(	class(gwecndtype)	this,
		type(gwecndparamfoundtype), intent(in)	found
	)

Definition at line 521 of file gwe-cnd.f90.

    use gwecndinputmodule, only: gwecndparamfoundtype
    class(GweCndType) :: this
    type(GweCndParamFoundType), intent(in) :: found
    !
    write (this%iout, '(1x,a)') 'Setting CND 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 CND Options'
    ! -- Return
    return

source_griddata()

subroutine gwecndmodule::source_griddata

(

class(gwecndtype)

this

)

Definition at line 609 of file gwe-cnd.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: lenmempath, linelength
    use gwecndinputmodule, only: gwecndparamfoundtype
    ! -- dummy
    class(GweCndType) :: this
    ! -- locals
    character(len=LINELENGTH) :: errmsg
    type(GweCndParamFoundType) :: 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%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)
    call mem_set_value(this%ktw, 'KTW', this%input_mempath, map, found%ktw)
    call mem_set_value(this%kts, 'KTS', this%input_mempath, map, found%kts)
    !
    ! -- set active flags
    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
    if (found%ktw) this%iktw = 1
    if (found%kts) this%ikts = 1
    !
    ! -- 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, '(1x,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 gwecndmodule::source_options

(

class(gwecndtype)

this

)

Definition at line 536 of file gwe-cnd.f90.

    ! -- modules
    use memorymanagerextmodule, only: mem_set_value
    use gwecndinputmodule, only: gwecndparamfoundtype
    ! -- dummy
    class(GweCndType) :: this
    ! -- locals
    type(GweCndParamFoundType) :: 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