DisConnExchange.f90

Go to the documentation of this file.

module disconnexchangemodule
  use kindmodule, only: i4b, dp, lgp
  use simvariablesmodule, only: errmsg
  use constantsmodule, only: lenauxname, lenboundname, linelength
  use simmodule, only: store_error, count_errors, store_error_filename
  use characterstringmodule
  use listmodule, only: listtype
  use memorymanagermodule, only: mem_allocate, mem_reallocate
  use numericalmodelmodule, only: numericalmodeltype
  use virtualmodelmodule, only: virtualmodeltype
  use numericalexchangemodule, only: numericalexchangetype
  implicit none
 
  private
  public :: disconnexchangetype
  public :: castasdisconnexchangeclass, adddisconnexchangetolist, &
            getdisconnexchangefromlist
 
  !> Exchange based on connection between discretizations of DisBaseType.
  !! The data specifies the connections, similar to the information stored
  !! in the connections object: DisBaseType%con
  !<
  type, extends(numericalexchangetype) :: disconnexchangetype
    character(len=LINELENGTH), pointer :: filename => null() !< name of the input file
 
    class(numericalmodeltype), pointer :: model1 => null() !< model 1
    class(numericalmodeltype), pointer :: model2 => null() !< model 2
    class(virtualmodeltype), pointer :: v_model1 => null() !< virtual model 1
    class(virtualmodeltype), pointer :: v_model2 => null() !< virtual model 2
    logical(LGP) :: is_datacopy !< when true, this exchange is just a data copy on another process and
                                !! not responsible for controlling movers, observations, ...
 
    integer(I4B), pointer :: nexg => null() !< number of exchanges
    integer(I4B), dimension(:), pointer, contiguous :: nodem1 => null() !< node numbers in model 1
    integer(I4B), dimension(:), pointer, contiguous :: nodem2 => null() !< node numbers in model 2
    integer(I4B), dimension(:), pointer, contiguous :: ihc => null() !< horizontal connection indicator array, size: nexg
    real(dp), dimension(:), pointer, contiguous :: cl1 => null() !< connection length 1, size: nexg
    real(dp), dimension(:), pointer, contiguous :: cl2 => null() !< connection length 2, size: nexg
    real(dp), dimension(:), pointer, contiguous :: hwva => null() !< horizontal widths, vertical flow areas, size: nexg
    integer(I4B), pointer :: naux => null() !< number of auxiliary variables
    character(len=LENBOUNDNAME), dimension(:), &
      pointer, contiguous :: boundname => null() !< boundnames
 
    character(len=LENAUXNAME), dimension(:), &
      pointer, contiguous :: auxname => null() !< vector of auxname
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: auxname_cst => null() !< copy of vector auxname that can be stored in memory manager
    real(dp), dimension(:, :), pointer, contiguous :: auxvar => null() !< array of auxiliary variable values
    integer(I4B), pointer :: ianglex => null() !< flag indicating anglex was read, if read, ianglex is index in auxvar
    integer(I4B), pointer :: icdist => null() !< flag indicating cdist was read, if read, icdist is index in auxvar
    integer(I4B), pointer :: iprpak => null() !< print input flag
    integer(I4B), pointer :: iprflow => null() !< print flag for cell by cell flows
    integer(I4B), pointer :: ipakcb => null() !< save flag for cell by cell flows
    integer(I4B), pointer :: inamedbound => null() !< flag to read boundnames
 
    integer(I4B), pointer :: ixt3d => null() !< flag indicating if XT3D should be applied on the interface: 0 = off, 1 = lhs, 2 = rhs
    logical(LGP), pointer :: dev_ifmod_on !< development option, forces interface model for this exchange
 
  contains
 
    procedure :: allocate_scalars
    procedure :: allocate_arrays
    procedure :: disconnex_da
    procedure :: use_interface_model
 
    ! protected
    procedure, pass(this) :: source_options
    procedure, pass(this) :: source_dimensions
    procedure, pass(this) :: source_data
    procedure, pass(this) :: noder
    procedure, pass(this) :: cellstr
 
  end type disconnexchangetype
 
  !> @ brief DisConnExchangeFoundType
  !!
  !!  This type is used to simplify the tracking of common parameters
  !!  that are sourced from the input context.
  !<
  type disconnexchangefoundtype
    logical :: naux = .false.
    logical :: ipakcb = .false.
    logical :: iprpak = .false.
    logical :: iprflow = .false.
    logical :: boundnames = .false.
    logical :: auxiliary = .false.
    logical :: dev_ifmod_on = .false.
    logical :: nexg = .false.
  end type disconnexchangefoundtype
 
contains
 
  !> @brief Source options from input context
  !<
  subroutine source_options(this, iout)
    ! -- modules
    use memorymanagerextmodule, only: mem_set_value
    use arrayhandlersmodule, only: ifind
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    integer(I4B), intent(in) :: iout !< for logging
    ! -- local
    type(disconnexchangefoundtype) :: found
    integer(I4B) :: ival, n
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%naux, 'NAUX', this%input_mempath, found%naux)
    call mem_set_value(this%ipakcb, 'IPAKCB', this%input_mempath, found%ipakcb)
    call mem_set_value(this%iprpak, 'IPRPAK', this%input_mempath, found%iprpak)
    call mem_set_value(this%iprflow, 'IPRFLOW', this%input_mempath, found%iprflow)
    call mem_set_value(this%inamedbound, 'BOUNDNAMES', this%input_mempath, &
                       found%boundnames)
    call mem_set_value(this%dev_ifmod_on, 'DEV_IFMOD_ON', this%input_mempath, &
                       found%dev_ifmod_on)
    !
    ! -- reallocate aux arrays if aux variables provided
    if (found%naux .and. this%naux > 0) then
      call mem_reallocate(this%auxname, lenauxname, this%naux, &
                          'AUXNAME', this%memoryPath)
      call mem_reallocate(this%auxname_cst, lenauxname, this%naux, &
                          'AUXNAME_CST', this%memoryPath)
      call mem_set_value(this%auxname_cst, 'AUXILIARY', this%input_mempath, &
                         found%auxiliary)
      !
      do n = 1, this%naux
        this%auxname(n) = this%auxname_cst(n)
      end do
      !
      ! -- If ANGLDEGX is an auxiliary variable, then anisotropy can be
      !    used in either model.  Store ANGLDEGX position in this%ianglex
      ival = ifind(this%auxname, 'ANGLDEGX')
      if (ival > 0) then
        this%ianglex = ival
      end if
      !
      ival = ifind(this%auxname, 'CDIST')
      if (ival > 0) then
        this%icdist = ival
      end if
    end if
    !
    if (found%ipakcb) then
      this%ipakcb = -1
      write (iout, '(4x,a)') &
        'EXCHANGE FLOWS WILL BE SAVED TO BINARY BUDGET FILES.'
    end if
    !
    if (found%iprpak) then
      write (iout, '(4x,a)') &
        'THE LIST OF EXCHANGES WILL BE PRINTED.'
    end if
    !
    if (found%iprflow) then
      write (iout, '(4x,a)') &
        'EXCHANGE FLOWS WILL BE PRINTED TO LIST FILES.'
    end if
    !
    if (found%boundnames) then
      write (iout, '(4x,a)') 'EXCHANGE BOUNDARIES HAVE NAMES IN LAST COLUMN'
    end if
    !
    if (found%dev_ifmod_on) then
      write (iout, '(4x,2a)') 'Interface model coupling approach manually &
        &activated for ', trim(this%name)
    end if
    !
    ! -- Return
    return
  end subroutine source_options
 
  !> @brief Source dimension from input context
  !<
  subroutine source_dimensions(this, iout)
    ! -- modules
    use memorymanagerextmodule, only: mem_set_value
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    integer(I4B), intent(in) :: iout !< for logging
    ! -- local
    type(disconnexchangefoundtype) :: found
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%nexg, 'NEXG', this%input_mempath, found%nexg)
    !
    write (iout, '(1x,a)') 'PROCESSING EXCHANGE DIMENSIONS'
    !
    if (found%nexg) then
      write (iout, '(4x,a,i0)') 'NEXG = ', this%nexg
    end if
    !
    write (iout, '(1x,a)') 'END OF EXCHANGE DIMENSIONS'
    !
    ! -- return
    return
  end subroutine source_dimensions
 
  !> @brief Returns reduced node number from user
  !< specified cell id.
  function noder(this, model, cellid, iout)
    ! -- modules
    use geomutilmodule, only: get_node
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    class(numericalmodeltype), pointer, intent(in) :: model
    integer(I4B), dimension(:), pointer, intent(in) :: cellid
    integer(I4B), intent(in) :: iout !< the output file unit
    integer(I4B) :: noder, node
    !
    if (model%dis%ndim == 1) then
      node = cellid(1)
    elseif (model%dis%ndim == 2) then
      node = get_node(cellid(1), 1, cellid(2), &
                      model%dis%mshape(1), 1, &
                      model%dis%mshape(2))
    else
      node = get_node(cellid(1), cellid(2), cellid(3), &
                      model%dis%mshape(1), &
                      model%dis%mshape(2), &
                      model%dis%mshape(3))
    end if
    noder = model%dis%get_nodenumber(node, 0)
    !
    ! -- return
    return
  end function noder
 
  !> @brief
  !<
  function cellstr(this, ndim, cellid, iout)
    ! -- modules
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    integer(I4B) :: ndim !< model DIS dimension
    integer(I4B), dimension(:), pointer, intent(in) :: cellid
    integer(I4B), intent(in) :: iout !< the output file unit
    character(len=20) :: cellstr
    character(len=*), parameter :: fmtndim1 = &
                                   "('(',i0,')')"
    character(len=*), parameter :: fmtndim2 = &
                                   "('(',i0,',',i0,')')"
    character(len=*), parameter :: fmtndim3 = &
                                   "('(',i0,',',i0,',',i0,')')"
    !
    cellstr = ''
    !
    select case (ndim)
    case (1)
      write (cellstr, fmtndim1) cellid(1)
    case (2)
      write (cellstr, fmtndim2) cellid(1), cellid(2)
    case (3)
      write (cellstr, fmtndim3) cellid(1), cellid(2), cellid(3)
    case default
    end select
    !
    ! -- return
    return
  end function cellstr
 
  !> @brief Source exchange data from input context
  !<
  subroutine source_data(this, iout)
    ! -- modules
    use memorymanagermodule, only: mem_setptr
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    integer(I4B), intent(in) :: iout !< the output file unit
    ! -- local
    integer(I4B), dimension(:, :), contiguous, pointer :: cellidm1
    integer(I4B), dimension(:, :), contiguous, pointer :: cellidm2
    integer(I4B), dimension(:), contiguous, pointer :: ihc
    real(DP), dimension(:), contiguous, pointer :: cl1
    real(DP), dimension(:), contiguous, pointer :: cl2
    real(DP), dimension(:), contiguous, pointer :: hwva
    real(DP), dimension(:, :), contiguous, pointer :: auxvar
    type(characterstringtype), dimension(:), contiguous, pointer :: boundname
    integer(I4B) :: ndim1, ndim2
    character(len=20) :: cellstr1, cellstr2
    character(len=2) :: cnfloat
    integer(I4B) :: nerr, iaux
    integer(I4B) :: iexg, nodem1, nodem2
    ! -- format
    character(len=*), parameter :: fmtexglabel = "(1x, 3a10, 50(a16))"
    character(len=*), parameter :: fmtexgdata = &
                                   "(5x, a, 1x, a ,I10, 50(1pg16.6))"
    character(len=40) :: fmtexgdata2
    !
    call mem_setptr(cellidm1, 'CELLIDM1', this%input_mempath)
    call mem_setptr(cellidm2, 'CELLIDM2', this%input_mempath)
    call mem_setptr(ihc, 'IHC', this%input_mempath)
    call mem_setptr(cl1, 'CL1', this%input_mempath)
    call mem_setptr(cl2, 'CL2', this%input_mempath)
    call mem_setptr(hwva, 'HWVA', this%input_mempath)
    call mem_setptr(auxvar, 'AUXVAR', this%input_mempath)
    call mem_setptr(boundname, 'BOUNDNAME', this%input_mempath)
    ndim1 = size(cellidm1, dim=1)
    ndim2 = size(cellidm2, dim=1)
    !
    write (iout, '(1x,a)') 'PROCESSING EXCHANGEDATA'
    !
    if (this%iprpak /= 0) then
      if (this%inamedbound == 0) then
        write (iout, fmtexglabel) 'NODEM1', 'NODEM2', 'IHC', &
          'CL1', 'CL2', 'HWVA', (adjustr(this%auxname(iaux)), &
                                 iaux=1, this%naux)
      else
        write (iout, fmtexglabel) 'NODEM1', 'NODEM2', 'IHC', 'CL1', 'CL2', &
          'HWVA', (adjustr(this%auxname(iaux)), iaux=1, this%naux), &
          ' BOUNDNAME      '
        ! Define format suitable for writing input data,
        ! any auxiliary variables, and boundname.
        write (cnfloat, '(i0)') 3 + this%naux
        fmtexgdata2 = '(5x, a, 1x, a, i10, '//trim(cnfloat)// &
                      '(1pg16.6), 1x, a)'
      end if
    end if
    !
    do iexg = 1, this%nexg
      !
      if (associated(this%model1)) then
        !
        ! -- Determine reduced node number
        nodem1 = this%noder(this%model1, cellidm1(:, iexg), iout)
        this%nodem1(iexg) = nodem1
        !
      else
        this%nodem1(iexg) = -1
      end if
      !
      if (associated(this%model2)) then
        !
        ! -- Determine reduced node number
        nodem2 = this%noder(this%model2, cellidm2(:, iexg), iout)
        this%nodem2(iexg) = nodem2
        !
      else
        this%nodem2(iexg) = -1
      end if
      !
      ! -- Read rest of input line
      this%ihc(iexg) = ihc(iexg)
      this%cl1(iexg) = cl1(iexg)
      this%cl2(iexg) = cl2(iexg)
      this%hwva(iexg) = hwva(iexg)
      do iaux = 1, this%naux
        this%auxvar(iaux, iexg) = auxvar(iaux, iexg)
      end do
      if (this%inamedbound == 1) then
        this%boundname(iexg) = boundname(iexg)
      end if
      !
      ! -- Write the data to listing file if requested
      if (this%iprpak /= 0) then
        cellstr1 = this%cellstr(ndim1, cellidm1(:, iexg), iout)
        cellstr2 = this%cellstr(ndim2, cellidm2(:, iexg), iout)
        if (this%inamedbound == 0) then
          write (iout, fmtexgdata) trim(cellstr1), trim(cellstr2), &
            this%ihc(iexg), this%cl1(iexg), this%cl2(iexg), &
            this%hwva(iexg), &
            (this%auxvar(iaux, iexg), iaux=1, this%naux)
        else
          write (iout, fmtexgdata2) trim(cellstr1), trim(cellstr2), &
            this%ihc(iexg), this%cl1(iexg), this%cl2(iexg), &
            this%hwva(iexg), &
            (this%auxvar(iaux, iexg), iaux=1, this%naux), &
            trim(this%boundname(iexg))
        end if
      end if
      !
      ! -- Check to see if nodem1 is outside of active domain
      if (associated(this%model1)) then
        if (nodem1 <= 0) then
          cellstr1 = this%cellstr(ndim1, cellidm1(:, iexg), iout)
          write (errmsg, *) &
            trim(adjustl(this%model1%name))// &
            ' Cell is outside active grid domain ('// &
            trim(adjustl(cellstr1))//').'
          call store_error(errmsg)
        end if
      end if
      !
      ! -- Check to see if nodem2 is outside of active domain
      if (associated(this%model2)) then
        if (nodem2 <= 0) then
          cellstr2 = this%cellstr(ndim2, cellidm2(:, iexg), iout)
          write (errmsg, *) &
            trim(adjustl(this%model2%name))// &
            ' Cell is outside active grid domain ('// &
            trim(adjustl(cellstr2))//').'
          call store_error(errmsg)
        end if
      end if
    end do
    !
    write (iout, '(1x,a)') 'END OF EXCHANGEDATA'
    !
    ! -- Stop if errors
    nerr = count_errors()
    if (nerr > 0) then
      call store_error('Errors encountered in exchange input file.')
      call store_error_filename(this%filename)
    end if
    !
    ! -- Return
    return
  end subroutine source_data
 
  !> @brief Allocate scalars and initialize to defaults
  !<
  subroutine allocate_scalars(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    !
    allocate (this%filename)
    this%filename = ''
    !
    call mem_allocate(this%nexg, 'NEXG', this%memoryPath)
    call mem_allocate(this%naux, 'NAUX', this%memoryPath)
    call mem_allocate(this%ianglex, 'IANGLEX', this%memoryPath)
    call mem_allocate(this%icdist, 'ICDIST', this%memoryPath)
    call mem_allocate(this%ixt3d, 'IXT3D', this%memoryPath)
    call mem_allocate(this%iprpak, 'IPRPAK', this%memoryPath)
    call mem_allocate(this%iprflow, 'IPRFLOW', this%memoryPath)
    call mem_allocate(this%ipakcb, 'IPAKCB', this%memoryPath)
    call mem_allocate(this%inamedbound, 'INAMEDBOUND', this%memoryPath)
    call mem_allocate(this%dev_ifmod_on, 'DEV_IFMOD_ON', this%memoryPath)
 
    call mem_allocate(this%auxname, lenauxname, 0, &
                      'AUXNAME', this%memoryPath)
    call mem_allocate(this%auxname_cst, lenauxname, 0, &
                      'AUXNAME_CST', this%memoryPath)
    !
    this%nexg = 0
    this%naux = 0
    this%ianglex = 0
    this%icdist = 0
    this%ixt3d = 0
    this%iprpak = 0
    this%iprflow = 0
    this%ipakcb = 0
    this%inamedbound = 0
    !
    this%dev_ifmod_on = .false.
    !
    ! -- Return
    return
  end subroutine allocate_scalars
 
  !> @brief Allocate array data, using the number of
  !! connected nodes @param nexg
  !<
  subroutine allocate_arrays(this)
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    !
    call mem_allocate(this%nodem1, this%nexg, 'NODEM1', this%memoryPath)
    call mem_allocate(this%nodem2, this%nexg, 'NODEM2', this%memoryPath)
    call mem_allocate(this%ihc, this%nexg, 'IHC', this%memoryPath)
    call mem_allocate(this%cl1, this%nexg, 'CL1', this%memoryPath)
    call mem_allocate(this%cl2, this%nexg, 'CL2', this%memoryPath)
    call mem_allocate(this%hwva, this%nexg, 'HWVA', this%memoryPath)
    ! NB: auxname array is allocated while parsing
    call mem_allocate(this%auxvar, this%naux, this%nexg, &
                      'AUXVAR', this%memoryPath)
    !
    ! allocate boundname
    if (this%inamedbound == 1) then
      allocate (this%boundname(this%nexg))
    else
      allocate (this%boundname(1))
    end if
    this%boundname(:) = ''
    !
    ! -- Return
    return
  end subroutine allocate_arrays
 
  !> @brief Should interface model be used to handle these
  !! exchanges, to be overridden for inheriting types
  !<
  function use_interface_model(this) result(use_im)
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    ! -- return
    logical(LGP) :: use_im !< flag whether interface model should be used
                          !! for this exchange instead
    !
    ! use im when one of the models is not local
    use_im = .not. (this%v_model1%is_local .and. this%v_model2%is_local)
    !
    ! -- Return
    return
  end function use_interface_model
 
  !> @brief Clean up all scalars and arrays
  !<
  subroutine disconnex_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(disconnexchangetype) :: this !< instance of exchange object
    !
    ! arrays
    call mem_deallocate(this%nodem1)
    call mem_deallocate(this%nodem2)
    call mem_deallocate(this%ihc)
    call mem_deallocate(this%cl1)
    call mem_deallocate(this%cl2)
    call mem_deallocate(this%hwva)
    call mem_deallocate(this%auxvar)
    !
    deallocate (this%boundname)
    !
    ! scalars
    call mem_deallocate(this%nexg)
    call mem_deallocate(this%naux)
    call mem_deallocate(this%auxname, 'AUXNAME', this%memoryPath)
    call mem_deallocate(this%auxname_cst, 'AUXNAME_CST', this%memoryPath)
    call mem_deallocate(this%ianglex)
    call mem_deallocate(this%icdist)
    call mem_deallocate(this%ixt3d)
    call mem_deallocate(this%iprpak)
    call mem_deallocate(this%iprflow)
    call mem_deallocate(this%ipakcb)
    call mem_deallocate(this%inamedbound)
    call mem_deallocate(this%dev_ifmod_on)
    !
    ! -- Return
    return
  end subroutine disconnex_da
 
  function castasdisconnexchangeclass(obj) result(res)
    implicit none
    ! -- dummy
    class(*), pointer, intent(inout) :: obj
    ! -- return
    class(disconnexchangetype), pointer :: res
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    class is (disconnexchangetype)
      res => obj
    end select
    !
    ! -- Return
    return
  end function castasdisconnexchangeclass
 
  subroutine adddisconnexchangetolist(list, exchange)
    implicit none
    ! -- dummy
    type(listtype), intent(inout) :: list
    class(disconnexchangetype), pointer, intent(in) :: exchange
    ! -- local
    class(*), pointer :: obj
    !
    obj => exchange
    call list%Add(obj)
    !
    ! -- Return
    return
  end subroutine adddisconnexchangetolist
 
  function getdisconnexchangefromlist(list, idx) result(res)
    implicit none
    ! -- dummy
    type(listtype), intent(inout) :: list
    integer(I4B), intent(in) :: idx
    ! -- return
    class(disconnexchangetype), pointer :: res
    ! -- local
    class(*), pointer :: obj
    !
    obj => list%GetItem(idx)
    res => castasdisconnexchangeclass(obj)
    !
    ! -- Return
    return
  end function getdisconnexchangefromlist
 
end module disconnexchangemodule