exg-gwfgwe.f90

Go to the documentation of this file.

module gwfgweexchangemodule
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: lenpackagename
  use listsmodule, only: basemodellist, baseexchangelist, &
                         baseconnectionlist
  use simmodule, only: store_error
  use simvariablesmodule, only: errmsg
  use baseexchangemodule, only: baseexchangetype, addbaseexchangetolist
  use spatialmodelconnectionmodule, only: spatialmodelconnectiontype, &
                                          get_smc_from_list
  use gwegweconnectionmodule, only: gwegweconnectiontype, castasgwegweconnection
  use gwfgwfconnectionmodule, only: gwfgwfconnectiontype, castasgwfgwfconnection
  use gwfgwfexchangemodule, only: gwfexchangetype, &
                                  getgwfexchangefromlist
  use basemodelmodule, only: basemodeltype, getbasemodelfromlist
  use gwfmodule, only: gwfmodeltype
  use gwemodule, only: gwemodeltype
  use bndmodule, only: bndtype, getbndfromlist
 
  implicit none
  public :: gwfgweexchangetype
  public :: gwfgwe_cr
 
  type, extends(baseexchangetype) :: gwfgweexchangetype
 
    integer(I4B), pointer :: m1_idx => null() !< index into the list of base exchanges for model 1
    integer(I4B), pointer :: m2_idx => null() !< index into the list of base exchanges for model 2
 
  contains
 
    procedure :: exg_df
    procedure :: exg_ar
    procedure :: exg_da
    procedure, private :: set_model_pointers
    procedure, private :: allocate_scalars
    procedure, private :: gwfbnd2gwefmi
    procedure, private :: gwfconn2gweconn
    procedure, private :: link_connections
 
  end type gwfgweexchangetype
 
contains
 
  !> @brief Create a new GWF to GWE exchange object
  !<
  subroutine gwfgwe_cr(filename, id, m1_id, m2_id)
    ! -- modules
    use simvariablesmodule, only: model_loc_idx
    ! -- dummy
    character(len=*), intent(in) :: filename
    integer(I4B), intent(in) :: id
    integer(I4B), intent(in) :: m1_id
    integer(I4B), intent(in) :: m2_id
    ! -- local
    class(baseexchangetype), pointer :: baseexchange => null()
    type(gwfgweexchangetype), pointer :: exchange => null()
    character(len=20) :: cint
    !
    ! -- Create a new exchange and add it to the baseexchangelist container
    allocate (exchange)
    baseexchange => exchange
    call addbaseexchangetolist(baseexchangelist, baseexchange)
    !
    ! -- Assign id and name
    exchange%id = id
    write (cint, '(i0)') id
    exchange%name = 'GWF-GWE_'//trim(adjustl(cint))
    exchange%memoryPath = exchange%name
    !
    ! -- allocate scalars
    call exchange%allocate_scalars()
    !
    ! -- NB: convert from id to local model index in base model list
    exchange%m1_idx = model_loc_idx(m1_id)
    exchange%m2_idx = model_loc_idx(m2_id)
    !
    ! -- set model pointers
    call exchange%set_model_pointers()
  end subroutine gwfgwe_cr
 
  !> @brief Allocate and read
  !<
  subroutine set_model_pointers(this)
    ! -- dummy
    class(gwfgweexchangetype) :: this
    ! -- local
    class(basemodeltype), pointer :: mb => null()
    type(gwfmodeltype), pointer :: gwfmodel => null()
    type(gwemodeltype), pointer :: gwemodel => null()
    !
    ! -- set gwfmodel
    gwfmodel => null()
    mb => getbasemodelfromlist(basemodellist, this%m1_idx)
    select type (mb)
    type is (gwfmodeltype)
      gwfmodel => mb
    end select
    !
    ! -- set gwemodel
    gwemodel => null()
    mb => getbasemodelfromlist(basemodellist, this%m2_idx)
    select type (mb)
    type is (gwemodeltype)
      gwemodel => mb
    end select
    !
    ! -- Verify that gwf model is of the correct type
    if (.not. associated(gwfmodel)) then
      write (errmsg, '(3a)') 'Problem with GWF-GWE exchange ', trim(this%name), &
        '.  Specified GWF Model does not appear to be of the correct type.'
      call store_error(errmsg, terminate=.true.)
    end if
    !
    ! -- Verify that gwe model is of the correct type
    if (.not. associated(gwemodel)) then
      write (errmsg, '(3a)') 'Problem with GWF-GWE exchange ', trim(this%name), &
        '.  Specified GWF Model does not appear to be of the correct type.'
      call store_error(errmsg, terminate=.true.)
    end if
    !
    ! -- Tell transport model fmi flows are not read from file
    gwemodel%fmi%flows_from_file = .false.
    !
    ! -- Set a pointer to the GWF bndlist.  This will allow the transport model
    !    to look through the flow packages and establish a link to GWF flows
    gwemodel%fmi%gwfbndlist => gwfmodel%bndlist
  end subroutine set_model_pointers
 
  !> @brief Define the GwfGwe Exchange object
  !<
  subroutine exg_df(this)
    ! -- modules
    use memorymanagermodule, only: mem_checkin
    ! -- dummy
    class(gwfgweexchangetype) :: this
    ! -- local
    class(basemodeltype), pointer :: mb => null()
    type(gwfmodeltype), pointer :: gwfmodel => null()
    type(gwemodeltype), pointer :: gwemodel => null()
    !
    ! -- set gwfmodel
    mb => getbasemodelfromlist(basemodellist, this%m1_idx)
    select type (mb)
    type is (gwfmodeltype)
      gwfmodel => mb
    end select
    !
    ! -- set gwemodel
    mb => getbasemodelfromlist(basemodellist, this%m2_idx)
    select type (mb)
    type is (gwemodeltype)
      gwemodel => mb
    end select
    !
    ! -- Check to make sure that flow is solved before transport and in a
    !    different IMS solution
    if (gwfmodel%idsoln >= gwemodel%idsoln) then
      write (errmsg, '(3a)') 'Problem with GWF-GWE exchange ', trim(this%name), &
        '.  The GWF model must be solved by a different IMS than the GWE model. &
        &Furthermore, the IMS specified for GWF must be listed in mfsim.nam &
        &before the IMS for GWE.'
      call store_error(errmsg, terminate=.true.)
    end if
    !
    ! -- Set pointer to flowja
    gwemodel%fmi%gwfflowja => gwfmodel%flowja
    call mem_checkin(gwemodel%fmi%gwfflowja, &
                     'GWFFLOWJA', gwemodel%fmi%memoryPath, &
                     'FLOWJA', gwfmodel%memoryPath)
 
    !
    ! -- Set the npf flag so that specific discharge is available for
    !    transport calculations if dispersion is active
    if (gwemodel%incnd > 0) then
      gwfmodel%npf%icalcspdis = 1
    end if
  end subroutine exg_df
 
  !> @brief Allocate and read
  !<
  subroutine exg_ar(this)
    ! -- modules
    use memorymanagermodule, only: mem_checkin
    use dismodule, only: distype
    use disvmodule, only: disvtype
    use disumodule, only: disutype
    ! -- dummy
    class(gwfgweexchangetype) :: this
    ! -- local
    class(basemodeltype), pointer :: mb => null()
    type(gwfmodeltype), pointer :: gwfmodel => null()
    type(gwemodeltype), pointer :: gwemodel => null()
    ! -- formats
    character(len=*), parameter :: fmtdiserr = &
      "('GWF and GWE Models do not have the same discretization for exchange&
      & ',a,'.&
      &  GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
      &  GWE Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
      &  Ensure discretization packages, including IDOMAIN, are identical.')"
    character(len=*), parameter :: fmtidomerr = &
      "('GWF and GWE Models do not have the same discretization for exchange&
      & ',a,'.&
      &  GWF Model and GWE Model have different IDOMAIN arrays.&
      &  Ensure discretization packages, including IDOMAIN, are identical.')"
    !
    ! -- set gwfmodel
    mb => getbasemodelfromlist(basemodellist, this%m1_idx)
    select type (mb)
    type is (gwfmodeltype)
      gwfmodel => mb
    end select
    !
    ! -- set gwemodel
    mb => getbasemodelfromlist(basemodellist, this%m2_idx)
    select type (mb)
    type is (gwemodeltype)
      gwemodel => mb
    end select
    !
    ! -- Check to make sure sizes are identical
    if (gwemodel%dis%nodes /= gwfmodel%dis%nodes .or. &
        gwemodel%dis%nodesuser /= gwfmodel%dis%nodesuser) then
      write (errmsg, fmtdiserr) trim(this%name), &
        gwfmodel%dis%nodesuser, &
        gwfmodel%dis%nodes, &
        gwemodel%dis%nodesuser, &
        gwemodel%dis%nodes
      call store_error(errmsg, terminate=.true.)
    end if
    !
    ! -- Make sure idomains are identical
    select type (gwfdis => gwfmodel%dis)
    type is (distype)
      select type (gwedis => gwemodel%dis)
      type is (distype)
        if (.not. all(gwfdis%idomain == gwedis%idomain)) then
          write (errmsg, fmtidomerr) trim(this%name)
          call store_error(errmsg, terminate=.true.)
        end if
      end select
    type is (disvtype)
      select type (gwedis => gwemodel%dis)
      type is (disvtype)
        if (.not. all(gwfdis%idomain == gwedis%idomain)) then
          write (errmsg, fmtidomerr) trim(this%name)
          call store_error(errmsg, terminate=.true.)
        end if
      end select
    type is (disutype)
      select type (gwedis => gwemodel%dis)
      type is (disutype)
        if (.not. all(gwfdis%idomain == gwedis%idomain)) then
          write (errmsg, fmtidomerr) trim(this%name)
          call store_error(errmsg, terminate=.true.)
        end if
      end select
    end select
    !
    ! -- setup pointers to gwf variables allocated in gwf_ar
    gwemodel%fmi%gwfhead => gwfmodel%x
    call mem_checkin(gwemodel%fmi%gwfhead, &
                     'GWFHEAD', gwemodel%fmi%memoryPath, &
                     'X', gwfmodel%memoryPath)
    gwemodel%fmi%gwfsat => gwfmodel%npf%sat
    call mem_checkin(gwemodel%fmi%gwfsat, &
                     'GWFSAT', gwemodel%fmi%memoryPath, &
                     'SAT', gwfmodel%npf%memoryPath)
    gwemodel%fmi%gwfspdis => gwfmodel%npf%spdis
    call mem_checkin(gwemodel%fmi%gwfspdis, &
                     'GWFSPDIS', gwemodel%fmi%memoryPath, &
                     'SPDIS', gwfmodel%npf%memoryPath)
    gwemodel%fmi%igwfspdis = gwfmodel%npf%icalcspdis
    !
    ! -- setup pointers to the flow storage rates. GWF strg arrays are
    !    available after the gwf_ar routine is called.
    if (gwemodel%inest > 0) then
      if (gwfmodel%insto > 0) then
        gwemodel%fmi%gwfstrgss => gwfmodel%sto%strgss
        gwemodel%fmi%igwfstrgss = 1
        if (gwfmodel%sto%iusesy == 1) then
          gwemodel%fmi%gwfstrgsy => gwfmodel%sto%strgsy
          gwemodel%fmi%igwfstrgsy = 1
        end if
      end if
    end if
    !
    ! -- Set a pointer to conc in buy
    if (gwfmodel%inbuy > 0) then
      call gwfmodel%buy%set_concentration_pointer(gwemodel%name, gwemodel%x, &
                                                  gwemodel%ibound)
    end if
    !
    ! -- Set a pointer to conc (which could be a temperature) in vsc
    if (gwfmodel%invsc > 0) then
      call gwfmodel%vsc%set_concentration_pointer(gwemodel%name, gwemodel%x, &
                                                  gwemodel%ibound, 1)
    end if
    !
    ! -- transfer the boundary package information from gwf to gwe
    call this%gwfbnd2gwefmi()
    !
    ! -- if mover package is active, then set a pointer to it's budget object
    if (gwfmodel%inmvr /= 0) then
      gwemodel%fmi%mvrbudobj => gwfmodel%mvr%budobj
    end if
    !
    ! -- connect Connections
    call this%gwfconn2gweconn(gwfmodel, gwemodel)
  end subroutine exg_ar
 
  !> @brief Link GWE connections to GWF connections or exchanges
  !<
  subroutine gwfconn2gweconn(this, gwfModel, gweModel)
    ! -- modules
    use simmodule, only: store_error
    use simvariablesmodule, only: iout
    use memorymanagermodule, only: mem_checkin
    ! -- dummy
    class(gwfgweexchangetype) :: this !< this exchange
    type(gwfmodeltype), pointer :: gwfModel !< the flow model
    type(gwemodeltype), pointer :: gweModel !< the energy transport model
    ! -- local
    class(spatialmodelconnectiontype), pointer :: conn => null()
    class(*), pointer :: objPtr => null()
    class(gwegweconnectiontype), pointer :: gweConn => null()
    class(gwfgwfconnectiontype), pointer :: gwfConn => null()
    class(gwfexchangetype), pointer :: gwfEx => null()
    integer(I4B) :: ic1, ic2, iex
    integer(I4B) :: gwfConnIdx, gwfExIdx
    logical(LGP) :: areEqual
    !
    ! loop over all connections
    gweloop: do ic1 = 1, baseconnectionlist%Count()
      !
      conn => get_smc_from_list(baseconnectionlist, ic1)
      if (.not. associated(conn%owner, gwemodel)) cycle gweloop
      !
      ! start with a GWE conn.
      objptr => conn
      gweconn => castasgwegweconnection(objptr)
      gwfconnidx = -1
      gwfexidx = -1
      !
      ! find matching GWF conn. in same list
      gwfloop: do ic2 = 1, baseconnectionlist%Count()
        conn => get_smc_from_list(baseconnectionlist, ic2)
        !
        if (associated(conn%owner, gwfmodel)) then
          objptr => conn
          gwfconn => castasgwfgwfconnection(objptr)
          !
          ! A model can have multiple exchanges, even connecting the same two
          ! models. Require that the GWF models match those referenced by the
          ! GWE exchange before comparing node lists, otherwise a GWE exchange
          ! could be linked to the wrong GWF-GWF exchange.
          if (gwfconn%gwfExchange%v_model1%name /= &
              gweconn%gweExchange%gwfmodelname1) cycle
          if (gwfconn%gwfExchange%v_model2%name /= &
              gweconn%gweExchange%gwfmodelname2) cycle
          !
          ! then, connecting the same nodes nrs will be
          ! sufficient evidence of equality
          areequal = (gwfconn%gwfExchange%nexg == gweconn%gweExchange%nexg)
          if (areequal) then
            areequal = all(gwfconn%prim_exchange%nodem1 == &
                           gweconn%prim_exchange%nodem1)
            areequal = areequal .and. all(gwfconn%prim_exchange%nodem2 == &
                                          gweconn%prim_exchange%nodem2)
          end if
          if (areequal) then
            ! same DIS, same exchange: link and go to next GWE conn.
            write (iout, '(/6a)') 'Linking exchange ', &
              trim(gweconn%prim_exchange%name), &
              ' to ', trim(gwfconn%prim_exchange%name), &
              ' (using interface model) for GWE model ', &
              trim(gwemodel%name)
            gwfconnidx = ic2
            call this%link_connections(gweconn, gwfconn)
            exit gwfloop
          end if
        end if
      end do gwfloop
      !
      ! fallback option: coupling to old gwfgwf exchange,
      ! (this will go obsolete at some point)
      if (gwfconnidx == -1) then
        gwfloopexg: do iex = 1, baseexchangelist%Count()
          gwfex => getgwfexchangefromlist(baseexchangelist, iex)
          !
          ! -- There is no guarantee that iex is a gwfExg, in which case
          !    it will return as null.  cycle if so.
          if (.not. associated(gwfex)) cycle gwfloopexg
          !
          if (associated(gwfex%model1, gwfmodel) .or. &
              associated(gwfex%model2, gwfmodel)) then
            !
            ! require that the GWF models match those referenced by the GWE
            ! exchange, otherwise a GWE exchange could be linked to the wrong
            ! GWF-GWF exchange
            if (gwfex%v_model1%name /= gweconn%gweExchange%gwfmodelname1) cycle
            if (gwfex%v_model2%name /= gweconn%gweExchange%gwfmodelname2) cycle
 
            ! check exchanges have same node counts
            areequal = size(gwfex%nodem1) == size(gweconn%prim_exchange%nodem1)
            ! then, connecting the same nodes nrs will be
            ! sufficient evidence of equality
            if (areequal) &
              areequal = all(gwfex%nodem1 == gweconn%prim_exchange%nodem1)
            if (areequal) &
              areequal = all(gwfex%nodem2 == gweconn%prim_exchange%nodem2)
            if (areequal) then
              ! link exchange to connection
              write (iout, '(/6a)') 'Linking exchange ', &
                trim(gweconn%prim_exchange%name), &
                ' to ', trim(gwfex%name), ' for GWE model ', &
                trim(gwemodel%name)
              gwfexidx = iex
              if (gweconn%owns_exchange) then
                gweconn%gweExchange%gwfsimvals => gwfex%simvals
                call mem_checkin(gweconn%gweExchange%gwfsimvals, &
                                 'GWFSIMVALS', gweconn%gweExchange%memoryPath, &
                                 'SIMVALS', gwfex%memoryPath)
              end if
              !
              !cdl link up mvt to mvr
              if (gwfex%inmvr > 0) then
                if (gweconn%owns_exchange) then
                  !cdl todo: check and make sure gweEx has mvt active
                  call gweconn%gweExchange%mvt%set_pointer_mvrbudobj( &
                    gwfex%mvr%budobj)
                end if
              end if
              !
              if (associated(gwfex%model2, gwfmodel)) gweconn%exgflowSign = -1
              gweconn%gweInterfaceModel%fmi%flows_from_file = .false.
              !
              exit gwfloopexg
            end if
          end if
          !
        end do gwfloopexg
      end if
      !
      if (gwfconnidx == -1 .and. gwfexidx == -1) then
        ! none found, report
        write (errmsg, '(/6a)') 'Missing GWF-GWF exchange when connecting GWE'// &
          ' model ', trim(gwemodel%name), ' with exchange ', &
          trim(gweconn%prim_exchange%name), ' to GWF model ', &
          trim(gwfmodel%name)
        call store_error(errmsg, terminate=.true.)
      end if
      !
    end do gweloop
  end subroutine gwfconn2gweconn
 
  !> @brief Links a GWE connection to its GWF counterpart
  !<
  subroutine link_connections(this, gweConn, gwfConn)
    ! -- modules
    use memorymanagermodule, only: mem_checkin
    ! -- dummy
    class(gwfgweexchangetype) :: this !< this exchange
    class(gwegweconnectiontype), pointer :: gweConn !< GWE connection
    class(gwfgwfconnectiontype), pointer :: gwfConn !< GWF connection
    !
    !gweConn%exgflowja => gwfConn%exgflowja
    if (gweconn%owns_exchange) then
      gweconn%gweExchange%gwfsimvals => gwfconn%gwfExchange%simvals
      call mem_checkin(gweconn%gweExchange%gwfsimvals, &
                       'GWFSIMVALS', gweconn%gweExchange%memoryPath, &
                       'SIMVALS', gwfconn%gwfExchange%memoryPath)
    end if
    !
    !cdl link up mvt to mvr
    if (gwfconn%gwfExchange%inmvr > 0) then
      if (gweconn%owns_exchange) then
        !cdl todo: check and make sure gweEx has mvt active
        call gweconn%gweExchange%mvt%set_pointer_mvrbudobj( &
          gwfconn%gwfExchange%mvr%budobj)
      end if
    end if
    !
    if (associated(gwfconn%gwfExchange%model2, gwfconn%owner)) then
      gweconn%exgflowSign = -1
    end if
    !
    ! fmi flows are not read from file
    gweconn%gweInterfaceModel%fmi%flows_from_file = .false.
    !
    ! set concentration pointer for buoyancy
    !call gwfConn%gwfInterfaceModel%buy%set_concentration_pointer( &
    !  gweConn%gweModel%name, &
    !  gweConn%conc, &
    !  gweConn%icbound)
  end subroutine link_connections
 
  !> @brief Deallocate memory
  !<
  subroutine exg_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(gwfgweexchangetype) :: this
    !
    call mem_deallocate(this%m1_idx)
    call mem_deallocate(this%m2_idx)
  end subroutine exg_da
 
  !> @brief Allocate GwfGwe exchange scalars
  !<
  subroutine allocate_scalars(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(gwfgweexchangetype) :: this
    !
    call mem_allocate(this%m1_idx, 'M1ID', this%memoryPath)
    call mem_allocate(this%m2_idx, 'M2ID', this%memoryPath)
    this%m1_idx = 0
    this%m2_idx = 0
  end subroutine allocate_scalars
 
  !> @brief Call routines in FMI that will set pointers to the necessary flow
  !! data (SIMVALS and SIMTOMVR) stored within each GWF flow package
  !<
  subroutine gwfbnd2gwefmi(this)
    ! -- dummy
    class(gwfgweexchangetype) :: this
    ! -- local
    integer(I4B) :: ngwfpack, ip, iterm, imover
    class(basemodeltype), pointer :: mb => null()
    type(gwfmodeltype), pointer :: gwfmodel => null()
    type(gwemodeltype), pointer :: gwemodel => null()
    class(bndtype), pointer :: packobj => null()
    !
    ! -- set gwfmodel
    mb => getbasemodelfromlist(basemodellist, this%m1_idx)
    select type (mb)
    type is (gwfmodeltype)
      gwfmodel => mb
    end select
    !
    ! -- set gwemodel
    mb => getbasemodelfromlist(basemodellist, this%m2_idx)
    select type (mb)
    type is (gwemodeltype)
      gwemodel => mb
    end select
    !
    ! -- Call routines in FMI that will set pointers to the necessary flow
    !    data (SIMVALS and SIMTOMVR) stored within each GWF flow package
    ngwfpack = gwfmodel%bndlist%Count()
    iterm = 1
    do ip = 1, ngwfpack
      packobj => getbndfromlist(gwfmodel%bndlist, ip)
      call gwemodel%fmi%gwfpackages(iterm)%set_pointers( &
        'SIMVALS', &
        packobj%memoryPath, packobj%input_mempath)
      iterm = iterm + 1
      !
      ! -- If a mover is active for this package, then establish a separate
      !    pointer link for the mover flows stored in SIMTOMVR
      imover = packobj%imover
      if (packobj%isadvpak /= 0) imover = 0
      if (imover /= 0) then
        call gwemodel%fmi%gwfpackages(iterm)%set_pointers( &
          'SIMTOMVR', &
          packobj%memoryPath, packobj%input_mempath)
        iterm = iterm + 1
      end if
    end do
  end subroutine gwfbnd2gwefmi
 
end module gwfgweexchangemodule