swfevpmodule Module Reference

This module contains the evaporation (EVP) package methods. More...

Data Types
type	swfevptype

Functions/Subroutines
subroutine, public	evp_create (packobj, id, ibcnum, inunit, iout, namemodel, pakname, mempath, dis, dfw, cxs)
	Create a Evaporation Package. More...
subroutine	evp_allocate_scalars (this)
	Allocate scalar members. More...
subroutine	evp_allocate_arrays (this, nodelist, auxvar)
	Allocate package arrays. More...
subroutine	evp_source_options (this)
	Source options specific to EVPType. More...
subroutine	log_evp_options (this, found_readasarrays)
	Log options specific to SwfEvpType. More...
subroutine	evp_source_dimensions (this)
	Source the dimensions for this package. More...
subroutine	evp_read_initial_attr (this)
	Part of allocate and read. More...
subroutine	evp_rp (this)
	Read and Prepare. More...
subroutine	evp_ck (this)
	Ensure evaporation is positive. More...
subroutine	evp_cf (this)
	Formulate the HCOF and RHS terms. More...
real(dp) function	get_qevp (this, node, rlen, snew, sold, evaporation)
	Calculate qevp. More...
real(dp) function	get_evap_reduce_mult (this, stage, bottom)
	Calculate multiplier to reduce evap as depth goes to zero. More...
subroutine	evp_fc (this, rhs, ia, idxglo, matrix_sln)
	Copy rhs and hcof into solution rhs and amat. More...
subroutine	evp_da (this)
	Deallocate memory. More...
subroutine	evp_define_listlabel (this)
	Define the list heading that is written to iout when PRINT_INPUT option is used. More...
subroutine	default_nodelist (this)
	Assign default nodelist when READASARRAYS is specified. More...
logical function	evp_obs_supported (this)
	Overrides BndTypebnd_obs_supported() More...
subroutine	evp_df_obs (this)
	Implements bnd_df_obs. More...
real(dp) function	evp_bound_value (this, col, row)
	Return requested boundary value. More...
real(dp) function, dimension(:), pointer	reach_length_pointer (this)

Variables
character(len=lenftype)	ftype = 'EVP'
character(len=lenpackagename)	text = ' EVP'

Detailed Description

This module can be used to represent evaporation onto streams and overland flow cells.

Function/Subroutine Documentation

default_nodelist()

subroutine swfevpmodule::default_nodelist ( class(swfevptype)  this )

private

Definition at line 554 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType) :: this
    ! local
    integer(I4B) :: nodeu, noder
 
    ! This is only called for readasarrays, so nodelist will be the size of
    ! the user grid, and will have a value of 0 for any entries where idomain
    ! is not 1
    do nodeu = 1, this%maxbound
      noder = this%dis%get_nodenumber(nodeu, 0)
      this%nodelist(nodeu) = noder
    end do
 
    ! Assign nbound
    this%nbound = this%maxbound
 

evp_allocate_arrays()

subroutine swfevpmodule::evp_allocate_arrays ( class(swfevptype)  this, integer(i4b), dimension(:), optional, pointer, contiguous  nodelist, real(dp), dimension(:, :), optional, pointer, contiguous  auxvar  )

private

Definition at line 147 of file swf-evp.f90.

    ! modules
    use memorymanagermodule, only: mem_setptr, mem_checkin
    ! dummy
    class(SwfEvpType) :: this
    integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
    real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
 
    ! allocate base arrays
    call this%BndExtType%allocate_arrays(nodelist, auxvar)
 
    ! set input context pointers
    call mem_setptr(this%evaporation, 'EVAPORATION', this%input_mempath)
 
    ! checkin input context pointers
    call mem_checkin(this%evaporation, 'EVAPORATION', this%memoryPath, &
                     'EVAPORATION', this%input_mempath)

evp_allocate_scalars()

subroutine swfevpmodule::evp_allocate_scalars ( class(swfevptype), intent(inout)  this )

private

Definition at line 127 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType), intent(inout) :: this
 
    ! allocate base scalars
    call this%BndExtType%allocate_scalars()
 
    ! allocate internal members
    call mem_allocate(this%iflowred, 'IFLOWRED', this%memoryPath)
    call mem_allocate(this%reduction_depth, 'REDUCTION_DEPTH', this%memoryPath)
    allocate (this%read_as_arrays)
 
    ! Set values
    this%iflowred = 1
    this%reduction_depth = dem6
    this%read_as_arrays = .false.

evp_bound_value()

real(dp) function swfevpmodule::evp_bound_value ( class(swfevptype), intent(inout)  this, integer(i4b), intent(in)  col, integer(i4b), intent(in)  row  )

private

Parameters

[in,out]

this

BndExtType object

Definition at line 604 of file swf-evp.f90.

    ! modules
    use constantsmodule, only: dzero
    ! dummy
    class(SwfEvpType), intent(inout) :: this !< BndExtType object
    integer(I4B), intent(in) :: col
    integer(I4B), intent(in) :: row
    ! result
    real(DP) :: bndval
 
    select case (col)
    case (1)
      if (this%iauxmultcol > 0) then
        bndval = this%evaporation(row) * this%auxvar(this%iauxmultcol, row)
      else
        bndval = this%evaporation(row)
      end if
    case default
      errmsg = 'Programming error. EVP bound value requested column '&
               &'outside range of ncolbnd (1).'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end select

Here is the call graph for this function:

evp_cf()

subroutine swfevpmodule::evp_cf ( class(swfevptype)  this )

private

Skip if no evaporation. Otherwise, calculate hcof and rhs

Definition at line 314 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType) :: this
    ! local
    integer(I4B) :: i
    integer(I4B) :: node
    integer(I4B) :: inwt
    real(DP) :: q
    real(DP) :: qeps
    real(DP) :: eps
    real(DP) :: derv
    real(DP) :: evap
    real(DP) :: rlen
    real(DP), dimension(:), pointer :: reach_length
 
    ! Return if no evaporation
    if (this%nbound == 0) return
 
    ! Set pointer to reach_length for 1d
    reach_length => this%reach_length_pointer()
    rlen = dzero
 
    ! Calculate hcof and rhs for each evaporation entry
    do i = 1, this%nbound
 
      ! Find the node number
      node = this%nodelist(i)
 
      ! cycle if nonexistent bound
      if (node <= 0) then
        this%hcof(i) = dzero
        this%rhs(i) = dzero
        cycle
      end if
 
      ! cycle if dry or constant head
      if (this%ibound(node) <= 0) then
        this%hcof(i) = dzero
        this%rhs(i) = dzero
        cycle
      end if
 
      ! Initialize hcof
      this%hcof(i) = dzero
 
      ! assign evap rate in length per time and multiply by auxvar
      evap = this%evaporation(i)
      if (this%iauxmultcol > 0) then
        evap = evap * this%auxvar(this%iauxmultcol, i)
      end if
 
      ! get reach length for 1d channel
      if (this%dis%is_1d()) then
        rlen = reach_length(node)
      end if
 
      ! Calculate volumetric evaporation flow in L^3/Td and add to rhs
      q = -this%get_qevp(node, rlen, this%xnew(node), this%xold(node), evap)
      this%rhs(i) = -q
 
      ! Code for adding newton terms
      inwt = 1
      if (inwt == 1) then
 
        ! calculate perturbed q
        eps = get_perturbation(this%xnew(node))
        qeps = -this%get_qevp(node, rlen, this%xnew(node) + eps, &
                              this%xold(node), evap)
 
        ! calculate derivative
        derv = (qeps - q) / eps
 
        ! add derivative to hcof and update rhs with derivate contribution
        this%hcof(i) = derv
        this%rhs(i) = this%rhs(i) + derv * this%xnew(node)
      end if
 
    end do

Here is the call graph for this function:

evp_ck()

subroutine swfevpmodule::evp_ck

(

class(swfevptype), intent(inout)

this

)

Definition at line 283 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType), intent(inout) :: this
    ! local
    character(len=30) :: nodestr
    integer(I4B) :: i, nr
    character(len=*), parameter :: fmterr = &
      &"('Specified stress ',i0, &
      &' evaporation (',g0,') is less than zero for cell', a)"
 
    ! Ensure evaporation rates are positive
    do i = 1, this%nbound
      nr = this%nodelist(i)
      if (nr <= 0) cycle
      if (this%evaporation(i) < dzero) then
        call this%dis%noder_to_string(nr, nodestr)
        write (errmsg, fmt=fmterr) i, this%evaporation(i), trim(nodestr)
        call store_error(errmsg)
      end if
    end do
 
    ! write summary of package error messages
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if

Here is the call graph for this function:

evp_create()

subroutine, public swfevpmodule::evp_create	(	class(bndtype), pointer	packobj,
		integer(i4b), intent(in)	id,
		integer(i4b), intent(in)	ibcnum,
		integer(i4b), intent(in)	inunit,
		integer(i4b), intent(in)	iout,
		character(len=*), intent(in)	namemodel,
		character(len=*), intent(in)	pakname,
		character(len=*), intent(in)	mempath,
		class(disbasetype), intent(inout), pointer	dis,
		type(swfdfwtype), intent(in), pointer	dfw,
		type(swfcxstype), intent(in), pointer	cxs
	)

Parameters

	packobj	pointer to default package type
[in]	id	package id
[in]	ibcnum	boundary condition number
[in]	inunit	unit number of CDB package input file
[in]	iout	unit number of model listing file
[in]	namemodel	model name
[in]	pakname	package name
[in]	mempath	input mempath
[in,out]	dis	the pointer to the discretization
[in]	dfw	the pointer to the dfw package
[in]	cxs	the pointer to the cxs package

Definition at line 78 of file swf-evp.f90.

    ! dummy
    class(BndType), pointer :: packobj !< pointer to default package type
    integer(I4B), intent(in) :: id !< package id
    integer(I4B), intent(in) :: ibcnum !< boundary condition number
    integer(I4B), intent(in) :: inunit !< unit number of CDB package input file
    integer(I4B), intent(in) :: iout !< unit number of model listing file
    character(len=*), intent(in) :: namemodel !< model name
    character(len=*), intent(in) :: pakname !< package name
    character(len=*), intent(in) :: mempath !< input mempath
    class(DisBaseType), pointer, intent(inout) :: dis !< the pointer to the discretization
    type(SwfDfwType), pointer, intent(in) :: dfw !< the pointer to the dfw package
    type(SwfCxsType), pointer, intent(in) :: cxs !< the pointer to the cxs package
    ! local
    type(SwfEvpType), pointer :: evpobj
 
    ! allocate evaporation object and scalar variables
    allocate (evpobj)
    packobj => evpobj
 
    ! create name and memory path
    call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath)
    packobj%text = text
 
    ! allocate scalars
    call evpobj%evp_allocate_scalars()
 
    ! initialize package
    call packobj%pack_initialize()
 
    packobj%inunit = inunit
    packobj%iout = iout
    packobj%id = id
    packobj%ibcnum = ibcnum
    packobj%ictMemPath = create_mem_path(namemodel, 'DFW')
 
    ! store pointer to dis
    evpobj%dis => dis
 
    ! store pointer to dfw
    evpobj%dfw => dfw
 
    ! store pointer to cxs
    evpobj%cxs => cxs

Here is the call graph for this function:

Here is the caller graph for this function:

evp_da()

subroutine swfevpmodule::evp_da ( class(swfevptype)  this )

private

Definition at line 503 of file swf-evp.f90.

    ! modules
    ! dummy
    class(SwfEvpType) :: this
 
    ! Deallocate parent package
    call this%BndExtType%bnd_da()
 
    ! scalars
    call mem_deallocate(this%iflowred)
    call mem_deallocate(this%reduction_depth)
    deallocate (this%read_as_arrays)
 
    ! arrays
    call mem_deallocate(this%evaporation, 'EVAPORATION', this%memoryPath)
 
    ! pointers
    nullify (this%dis)
    nullify (this%dfw)
    nullify (this%cxs)

evp_define_listlabel()

subroutine swfevpmodule::evp_define_listlabel ( class(swfevptype), intent(inout)  this )

private

Definition at line 528 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType), intent(inout) :: this
    !
    ! create the header list label
    this%listlabel = trim(this%filtyp)//' NO.'
    if (this%dis%ndim == 3) then
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
    elseif (this%dis%ndim == 2) then
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
    else
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
    end if
    write (this%listlabel, '(a, a16)') trim(this%listlabel), 'EVAPORATION'
!    if(this%multindex > 0) &
!      write(this%listlabel, '(a, a16)') trim(this%listlabel), 'MULTIPLIER'
    if (this%inamedbound == 1) then
      write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
    end if

evp_df_obs()

subroutine swfevpmodule::evp_df_obs ( class(swfevptype)  this )

private

Store observation type supported by EVP package. Overrides BndTypebnd_df_obs

Definition at line 591 of file swf-evp.f90.

    implicit none
    ! dummy
    class(SwfEvpType) :: this
    ! local
    integer(I4B) :: indx
 
    call this%obs%StoreObsType('evp', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor

Here is the call graph for this function:

evp_fc()

subroutine swfevpmodule::evp_fc ( class(swfevptype)  this, real(dp), dimension(:), intent(inout)  rhs, integer(i4b), dimension(:), intent(in)  ia, integer(i4b), dimension(:), intent(in)  idxglo, class(matrixbasetype), pointer  matrix_sln  )

private

Definition at line 481 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType) :: this
    real(DP), dimension(:), intent(inout) :: rhs
    integer(I4B), dimension(:), intent(in) :: ia
    integer(I4B), dimension(:), intent(in) :: idxglo
    class(MatrixBaseType), pointer :: matrix_sln
    ! local
    integer(I4B) :: i, n, ipos
 
    ! Copy package rhs and hcof into solution rhs and amat
    do i = 1, this%nbound
      n = this%nodelist(i)
      if (n <= 0) cycle
      rhs(n) = rhs(n) + this%rhs(i)
      ipos = ia(n)
      call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i))
    end do

evp_obs_supported()

logical function swfevpmodule::evp_obs_supported ( class(swfevptype)  this )

private

Definition at line 579 of file swf-evp.f90.

    implicit none
    ! dummy
    class(SwfEvpType) :: this
    evp_obs_supported = .true.

evp_read_initial_attr()

subroutine swfevpmodule::evp_read_initial_attr ( class(swfevptype), intent(inout)  this )

private

Definition at line 246 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType), intent(inout) :: this
 
    if (this%read_as_arrays) then
      call this%default_nodelist()
    end if

evp_rp()

subroutine swfevpmodule::evp_rp ( class(swfevptype), intent(inout)  this )

private

Read itmp and read new boundaries if itmp > 0

Definition at line 259 of file swf-evp.f90.

    ! modules
    use tdismodule, only: kper
    implicit none
    ! dummy
    class(SwfEvpType), intent(inout) :: this
 
    if (this%iper /= kper) return
 
    if (this%read_as_arrays) then
      ! no need to do anything because this%evaporation points directly to
      ! the input context evaporation, which is automatically updated by idm
    else
      call this%BndExtType%bnd_rp()
    end if
 
    ! Write the list to iout if requested
    if (this%iprpak /= 0) then
      call this%write_list()
    end if

evp_source_dimensions()

subroutine swfevpmodule::evp_source_dimensions ( class(swfevptype), intent(inout)  this )

private

Definition at line 214 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType), intent(inout) :: this
 
    if (this%read_as_arrays) then
 
      ! Set maxbound to the number of cells per layer, which is simply
      ! nrow * ncol for a dis2d grid, and nodesuser for disv2d and disv1d
      this%maxbound = this%dis%get_ncpl()
 
      ! verify dimensions were set
      if (this%maxbound <= 0) then
        write (errmsg, '(a)') &
          'MAXBOUND must be an integer greater than zero.'
        call store_error(errmsg)
        call store_error_filename(this%input_fname)
      end if
 
    else
 
      ! source maxbound
      call this%BndExtType%source_dimensions()
 
    end if
 
    ! Call define_listlabel to construct the list label that is written
    !    when PRINT_INPUT option is used.
    call this%define_listlabel()

Here is the call graph for this function:

evp_source_options()

subroutine swfevpmodule::evp_source_options

(

class(swfevptype), intent(inout)

this

)

Definition at line 168 of file swf-evp.f90.

    ! modules
    use memorymanagerextmodule, only: mem_set_value
    implicit none
    ! dummy
    class(SwfEvpType), intent(inout) :: this
    ! local
    logical(LGP) :: found_readasarrays = .false.
 
    ! source common bound options
    call this%BndExtType%source_options()
 
    ! update defaults with idm sourced values
    call mem_set_value(this%read_as_arrays, 'READASARRAYS', this%input_mempath, &
                       found_readasarrays)
 
    ! log evp params
    call this%log_evp_options(found_readasarrays)

get_evap_reduce_mult()

real(dp) function swfevpmodule::get_evap_reduce_mult ( class(swfevptype)  this, real(dp), intent(in)  stage, real(dp), intent(in)  bottom  )

private

Definition at line 461 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType) :: this
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: bottom
    ! return
    real(DP) :: qmult
    ! local
    real(DP) :: tp
 
    qmult = done
    if (this%iflowred == 1) then
      tp = bottom + this%reduction_depth
      qmult = sqsaturation(tp, bottom, stage)
    end if
 

Here is the call graph for this function:

get_qevp()

real(dp) function swfevpmodule::get_qevp ( class(swfevptype)  this, integer(i4b), intent(in)  node, real(dp), intent(in)  rlen, real(dp), intent(in)  snew, real(dp), intent(in)  sold, real(dp), intent(in)  evaporation  )

private

Calculate qevp for both channel and overland flow grids. Approximate the average water surface width of the channel as wavg = delta A over delta h, and then multiply wavg by reach length to come up with surface water area for the channel. Reduce evaporation when depths are small and shut it off when there is no water in the cell.

Parameters

	this	this instance
[in]	node	reduced node number
[in]	rlen	length of reach
[in]	snew	current stage in reach
[in]	sold	previous stage in reach
[in]	evaporation	evaporation rate in length per time

Definition at line 403 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType) :: this !< this instance
    integer(I4B), intent(in) :: node !< reduced node number
    real(DP), intent(in) :: rlen !< length of reach
    real(DP), intent(in) :: snew !< current stage in reach
    real(DP), intent(in) :: sold !< previous stage in reach
    real(DP), intent(in) :: evaporation !< evaporation rate in length per time
    ! return
    real(DP) :: qevp
    ! local
    integer(I4B) :: idcxs
    real(DP) :: depth
    real(DP) :: bt
    real(DP) :: area
    real(DP) :: anew
    real(DP) :: aold
    real(DP) :: denom
    real(DP) :: wavg
    real(DP) :: width_channel
    real(DP) :: dummy
    real(DP) :: qmult
 
    ! calculate depth of water
    bt = this%dis%bot(node)
 
    ! Determine the water surface area
    if (this%dis%is_2d()) then
      ! overland flow case
      area = this%dis%get_area(node)
    else if (this%dis%is_1d()) then
      ! channel case
      idcxs = this%dfw%idcxs(node)
      call this%dis%get_flow_width(node, node, 0, width_channel, dummy)
 
      depth = snew - bt
      anew = this%cxs%get_area(idcxs, width_channel, depth)
      depth = sold - bt
      aold = this%cxs%get_area(idcxs, width_channel, depth)
      wavg = this%cxs%get_wetted_top_width(idcxs, width_channel, depth)
      denom = snew - sold
      if (abs(denom) > dprec) then
        wavg = (anew - aold) / (snew - sold)
      end if
      area = rlen * wavg
 
    end if
 
    ! Reduce the evap rate as cell goes dry
    qmult = this%get_evap_reduce_mult(snew, bt)
 
    ! calculate volumetric evaporation flow in L^3/T
    qevp = evaporation * area * qmult
 

log_evp_options()

subroutine swfevpmodule::log_evp_options	(	class(swfevptype), intent(inout)	this,
		logical(lgp), intent(in)	found_readasarrays
	)

Definition at line 190 of file swf-evp.f90.

    implicit none
    ! dummy
    class(SwfEvpType), intent(inout) :: this
    logical(LGP), intent(in) :: found_readasarrays
    ! formats
    character(len=*), parameter :: fmtreadasarrays = &
      &"(4x, 'EVAPORATION INPUT WILL BE READ AS ARRAY(S).')"
 
    ! log found options
    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text)) &
      //' OPTIONS'
 
    if (found_readasarrays) then
      write (this%iout, fmtreadasarrays)
    end if
 
    ! close logging block
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' OPTIONS'

reach_length_pointer()

real(dp) function, dimension(:), pointer swfevpmodule::reach_length_pointer

(

class(swfevptype)

this

)

Parameters

this	this instance

Definition at line 629 of file swf-evp.f90.

    ! dummy
    class(SwfEvpType) :: this !< this instance
    ! return
    real(DP), dimension(:), pointer :: ptr
    ! local
    class(DisBaseType), pointer :: dis
 
    ptr => null()
    dis => this%dis
    select type (dis)
    type is (disv1dtype)
      ptr => dis%length
    end select
 

Variable Documentation

ftype

character(len=lenftype) swfevpmodule::ftype = 'EVP'

private

Definition at line 35 of file swf-evp.f90.

  character(len=LENFTYPE) :: ftype = 'EVP'

text

character(len=lenpackagename) swfevpmodule::text = ' EVP'

private

Definition at line 36 of file swf-evp.f90.

  character(len=LENPACKAGENAME) :: text = '             EVP'