exg-gwfgwt.f90

Go to the documentation of this file.

module gwfgwtexchangemodule
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: lenpackagename, linelength
  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 gwtgwtconnectionmodule, only: gwtgwtconnectiontype, castasgwtgwtconnection
  use gwfgwfconnectionmodule, only: gwfgwfconnectiontype, castasgwfgwfconnection
  use gwfgwfexchangemodule, only: gwfexchangetype, &
                                  getgwfexchangefromlist
  use gwtgwtexchangemodule, only: gwtexchangetype
  use basemodelmodule, only: basemodeltype, getbasemodelfromlist
  use gwfmodule, only: gwfmodeltype
  use gwtmodule, only: gwtmodeltype
  use bndmodule, only: bndtype, getbndfromlist
 
  implicit none
  public :: gwfgwtexchangetype
  public :: gwfgwt_cr
 
  type, extends(baseexchangetype) :: gwfgwtexchangetype
 
    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
    character(len=LINELENGTH) :: filename !< the input file for the GWF-GWT exchange
 
  contains
 
    procedure :: exg_df
    procedure :: exg_ar
    procedure :: exg_da
    procedure, private :: set_model_pointers
    procedure, private :: allocate_scalars
    procedure, private :: gwfbnd2gwtfmi
    procedure, private :: gwfconn2gwtconn
    procedure, private :: link_connections
 
  end type gwfgwtexchangetype
 
contains
 
  !> @brief Create a new GWF to GWT exchange object
  !<
  subroutine gwfgwt_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(gwfgwtexchangetype), 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-GWT_'//trim(adjustl(cint))
    exchange%memoryPath = exchange%name
    exchange%filename = filename
    !
    ! -- 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 gwfgwt_cr
 
  !> @brief Allocate and read
  !<
  subroutine set_model_pointers(this)
    ! -- dummy
    class(gwfgwtexchangetype) :: this
    ! -- local
    class(basemodeltype), pointer :: mb => null()
    type(gwfmodeltype), pointer :: gwfmodel => null()
    type(gwtmodeltype), pointer :: gwtmodel => null()
    !
    ! -- set gwfmodel
    gwfmodel => null()
    mb => getbasemodelfromlist(basemodellist, this%m1_idx)
    select type (mb)
    type is (gwfmodeltype)
      gwfmodel => mb
    end select
    !
    ! -- set gwtmodel
    gwtmodel => null()
    mb => getbasemodelfromlist(basemodellist, this%m2_idx)
    select type (mb)
    type is (gwtmodeltype)
      gwtmodel => mb
    end select
    !
    ! -- Verify that gwf model is of the correct type
    if (.not. associated(gwfmodel)) then
      write (errmsg, '(3a)') 'Problem with GWF-GWT 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 gwt model is of the correct type
    if (.not. associated(gwtmodel)) then
      write (errmsg, '(3a)') 'Problem with GWF-GWT exchange ', trim(this%name), &
        '.  Specified GWT 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
    gwtmodel%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
    gwtmodel%fmi%gwfbndlist => gwfmodel%bndlist
  end subroutine set_model_pointers
 
  !> @brief Define the GwfGwt Exchange object
  !<
  subroutine exg_df(this)
    ! -- modules
    use memorymanagermodule, only: mem_checkin
    ! -- dummy
    class(gwfgwtexchangetype) :: this
    ! -- local
    class(basemodeltype), pointer :: mb => null()
    type(gwfmodeltype), pointer :: gwfmodel => null()
    type(gwtmodeltype), pointer :: gwtmodel => null()
    !
    ! -- set gwfmodel
    mb => getbasemodelfromlist(basemodellist, this%m1_idx)
    select type (mb)
    type is (gwfmodeltype)
      gwfmodel => mb
    end select
    !
    ! -- set gwtmodel
    mb => getbasemodelfromlist(basemodellist, this%m2_idx)
    select type (mb)
    type is (gwtmodeltype)
      gwtmodel => mb
    end select
    !
    ! -- Check to make sure that flow is solved before transport and in a
    !    different IMS solution
    if (gwfmodel%idsoln >= gwtmodel%idsoln) then
      write (errmsg, '(3a)') 'Problem with GWF-GWT exchange ', trim(this%name), &
        '.  The GWF model must be solved by a different IMS than the GWT model. &
        &Furthermore, the IMS specified for GWF must be listed in mfsim.nam &
        &before the IMS for GWT.'
      call store_error(errmsg, terminate=.true.)
    end if
    !
    ! -- Set pointer to flowja
    gwtmodel%fmi%gwfflowja => gwfmodel%flowja
    call mem_checkin(gwtmodel%fmi%gwfflowja, &
                     'GWFFLOWJA', gwtmodel%fmi%memoryPath, &
                     'FLOWJA', gwfmodel%memoryPath)
 
    !
    ! -- Set the npf flag so that specific discharge is available for
    !    transport calculations if dispersion is active
    if (gwtmodel%indsp > 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(gwfgwtexchangetype) :: this
    ! -- local
    class(basemodeltype), pointer :: mb => null()
    type(gwfmodeltype), pointer :: gwfmodel => null()
    type(gwtmodeltype), pointer :: gwtmodel => null()
    ! -- formats
    character(len=*), parameter :: fmtdiserr = &
      "('GWF and GWT Models do not have the same discretization for exchange&
      & ',a,'.&
      &  GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
      &  GWT Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
      &  Ensure discretization packages, including IDOMAIN, are identical.')"
    character(len=*), parameter :: fmtidomerr = &
      "('GWF and GWT Models do not have the same discretization for exchange&
      & ',a,'.&
      &  GWF Model and GWT 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 gwtmodel
    mb => getbasemodelfromlist(basemodellist, this%m2_idx)
    select type (mb)
    type is (gwtmodeltype)
      gwtmodel => mb
    end select
    !
    ! -- Check to make sure sizes are identical
    if (gwtmodel%dis%nodes /= gwfmodel%dis%nodes .or. &
        gwtmodel%dis%nodesuser /= gwfmodel%dis%nodesuser) then
      write (errmsg, fmtdiserr) trim(this%name), &
        gwfmodel%dis%nodesuser, &
        gwfmodel%dis%nodes, &
        gwtmodel%dis%nodesuser, &
        gwtmodel%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 (gwtdis => gwtmodel%dis)
      type is (distype)
        if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
          write (errmsg, fmtidomerr) trim(this%name)
          call store_error(errmsg, terminate=.true.)
        end if
      end select
    type is (disvtype)
      select type (gwtdis => gwtmodel%dis)
      type is (disvtype)
        if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
          write (errmsg, fmtidomerr) trim(this%name)
          call store_error(errmsg, terminate=.true.)
        end if
      end select
    type is (disutype)
      select type (gwtdis => gwtmodel%dis)
      type is (disutype)
        if (.not. all(gwfdis%idomain == gwtdis%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
    gwtmodel%fmi%gwfhead => gwfmodel%x
    call mem_checkin(gwtmodel%fmi%gwfhead, &
                     'GWFHEAD', gwtmodel%fmi%memoryPath, &
                     'X', gwfmodel%memoryPath)
    gwtmodel%fmi%gwfsat => gwfmodel%npf%sat
    call mem_checkin(gwtmodel%fmi%gwfsat, &
                     'GWFSAT', gwtmodel%fmi%memoryPath, &
                     'SAT', gwfmodel%npf%memoryPath)
    gwtmodel%fmi%gwfspdis => gwfmodel%npf%spdis
    call mem_checkin(gwtmodel%fmi%gwfspdis, &
                     'GWFSPDIS', gwtmodel%fmi%memoryPath, &
                     'SPDIS', gwfmodel%npf%memoryPath)
    gwtmodel%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 (gwtmodel%inmst > 0) then
      if (gwfmodel%insto > 0) then
        gwtmodel%fmi%gwfstrgss => gwfmodel%sto%strgss
        gwtmodel%fmi%igwfstrgss = 1
        if (gwfmodel%sto%iusesy == 1) then
          gwtmodel%fmi%gwfstrgsy => gwfmodel%sto%strgsy
          gwtmodel%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(gwtmodel%name, gwtmodel%x, &
                                                  gwtmodel%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(gwtmodel%name, gwtmodel%x, &
                                                  gwtmodel%ibound)
    end if
    !
    ! -- transfer the boundary package information from gwf to gwt
    call this%gwfbnd2gwtfmi()
    !
    ! -- if mover package is active, then set a pointer to it's budget object
    if (gwfmodel%inmvr /= 0) then
      gwtmodel%fmi%mvrbudobj => gwfmodel%mvr%budobj
    end if
    !
    ! -- connect Connections
    call this%gwfconn2gwtconn(gwfmodel, gwtmodel)
  end subroutine exg_ar
 
  !> @brief Link GWT connections to GWF connections or exchanges
  !<
  subroutine gwfconn2gwtconn(this, gwfModel, gwtModel)
    ! -- modules
    use simmodule, only: store_error, store_error_filename, count_errors
    use simvariablesmodule, only: iout
    use memorymanagermodule, only: mem_checkin
    ! -- dummy
    class(gwfgwtexchangetype) :: this !< this exchange
    type(gwfmodeltype), pointer :: gwfModel !< the flow model
    type(gwtmodeltype), pointer :: gwtModel !< the transport model
    ! -- local
    class(spatialmodelconnectiontype), pointer :: conn => null()
    class(*), pointer :: objPtr => null()
    class(gwtgwtconnectiontype), pointer :: gwtConn => null()
    class(gwfgwfconnectiontype), pointer :: gwfConn => null()
    class(gwfexchangetype), pointer :: gwfExg => null()
    class(gwtexchangetype), pointer :: gwtExg => null()
    integer(I4B) :: ic1, ic2, iex
    integer(I4B) :: gwfConnIdx, gwfExIdx
    logical(LGP) :: areEqual
    !
    ! loop over all connections
    gwtloop: do ic1 = 1, baseconnectionlist%Count()
      !
      conn => get_smc_from_list(baseconnectionlist, ic1)
      if (.not. associated(conn%owner, gwtmodel)) cycle gwtloop
      !
      ! start with a GWT conn.
      objptr => conn
      gwtconn => castasgwtgwtconnection(objptr)
      gwtexg => gwtconn%gwtExchange
      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)
          gwfexg => gwfconn%gwfExchange
          !
          ! A model can have multiple exchanges, even connecting the same two
          ! models. We have a match if
          !  1. gwtgwt%model1 is connected to gwfgwf%model1
          !  2. gwtgwt%model2 is connected to gwfgwf%model2
          !  3. the list of connected nodes (nodem1, nodem2) is equivalent, such
          !     that it contains the same nodes, appearing in the same order in the
          !     exchange data block
          !
          if (gwfexg%v_model1%name /= gwtexg%gwfmodelname1) cycle
          if (gwfexg%v_model2%name /= gwtexg%gwfmodelname2) cycle
          !
          areequal = (gwfexg%nexg == gwtexg%nexg)
          if (areequal) then
            areequal = all(gwfexg%nodem1 == gwtexg%nodem1)
            areequal = areequal .and. all(gwfexg%nodem2 == gwtexg%nodem2)
          end if
          if (areequal) then
            ! same DIS, same exchange: link and go to next GWT conn.
            write (iout, '(/6a)') 'Linking exchange ', &
              trim(gwtexg%name), ' to ', trim(gwfexg%name), &
              ' (using interface model) for GWT model ', &
              trim(gwtmodel%name)
            gwfconnidx = ic2
            call this%link_connections(gwtconn, gwfconn)
            exit gwfloop
          end if
        end if
      end do gwfloop
      !
      ! fallback option: coupling to old gwfgwf exchange,
      ! the conditions are equal to what is used above
      ! (this will go obsolete at some point)
      if (gwfconnidx == -1) then
        gwfloopexg: do iex = 1, baseexchangelist%Count()
          gwfexg => getgwfexchangefromlist(baseexchangelist, iex)
          !
          if (.not. associated(gwfexg)) cycle gwfloopexg
          !
          if (associated(gwfexg%model1, gwfmodel) .or. &
              associated(gwfexg%model2, gwfmodel)) then
            !
            if (gwfexg%v_model1%name /= gwtexg%gwfmodelname1) cycle
            if (gwfexg%v_model2%name /= gwtexg%gwfmodelname2) cycle
            !
            areequal = (gwfexg%nexg == gwtexg%nexg)
            !
            if (areequal) then
              areequal = all(gwfexg%nodem1 == gwtexg%nodem1)
              areequal = areequal .and. all(gwfexg%nodem2 == gwtexg%nodem2)
            end if
            if (areequal) then
              ! link exchange to connection
              write (iout, '(/6a)') 'Linking exchange ', &
                trim(gwtexg%name), ' to ', trim(gwfexg%name), ' for GWT model ', &
                trim(gwtmodel%name)
              gwfexidx = iex
              if (gwtconn%owns_exchange) then
                gwtexg%gwfsimvals => gwfexg%simvals
                call mem_checkin(gwtexg%gwfsimvals, &
                                 'GWFSIMVALS', gwtexg%memoryPath, &
                                 'SIMVALS', gwfexg%memoryPath)
              end if
              !
              !cdl link up mvt to mvr
              if (gwfexg%inmvr > 0) then
                if (gwtconn%owns_exchange) then
                  !cdl todo: check and make sure gwtEx has mvt active
                  call gwtexg%mvt%set_pointer_mvrbudobj(gwfexg%mvr%budobj)
                end if
              end if
              !
              if (associated(gwfexg%model2, gwfmodel)) gwtconn%exgflowSign = -1
              gwtconn%gwtInterfaceModel%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, *) 'Cannot find GWF-GWF exchange when connecting'// &
          ' GWT model ', trim(gwtmodel%name), ' with exchange ', &
          trim(gwtexg%name), ' to GWF model ', trim(gwfmodel%name), &
          '. Note: GWF-GWF and GWT-GWT need identical exchange data '// &
          '(both in value and order) for the match to succeed.'
        call store_error(errmsg)
      end if
      !
    end do gwtloop
    !
    ! -- report errors
    if (count_errors() > 0) then
      call store_error_filename(this%filename)
    end if
  end subroutine gwfconn2gwtconn
 
  !> @brief Links a GWT connection to its GWF counterpart
  !<
  subroutine link_connections(this, gwtConn, gwfConn)
    ! -- modules
    use memorymanagermodule, only: mem_checkin
    ! -- dummy
    class(gwfgwtexchangetype) :: this !< this exchange
    class(gwtgwtconnectiontype), pointer :: gwtConn !< GWT connection
    class(gwfgwfconnectiontype), pointer :: gwfConn !< GWF connection
    !
    !gwtConn%exgflowja => gwfConn%exgflowja
    if (gwtconn%owns_exchange) then
      gwtconn%gwtExchange%gwfsimvals => gwfconn%gwfExchange%simvals
      call mem_checkin(gwtconn%gwtExchange%gwfsimvals, &
                       'GWFSIMVALS', gwtconn%gwtExchange%memoryPath, &
                       'SIMVALS', gwfconn%gwfExchange%memoryPath)
    end if
    !
    !cdl link up mvt to mvr
    if (gwfconn%gwfExchange%inmvr > 0) then
      if (gwtconn%owns_exchange) then
        !cdl todo: check and make sure gwtEx has mvt active
        call gwtconn%gwtExchange%mvt%set_pointer_mvrbudobj( &
          gwfconn%gwfExchange%mvr%budobj)
      end if
    end if
    !
    if (associated(gwfconn%gwfExchange%model2, gwfconn%owner)) then
      gwtconn%exgflowSign = -1
    end if
    !
    ! fmi flows are not read from file
    gwtconn%gwtInterfaceModel%fmi%flows_from_file = .false.
    !
    ! set concentration pointer for buoyancy
    ! call gwfConn%gwfInterfaceModel%buy%set_concentration_pointer( &
    !   gwtConn%gwtModel%name, &
    !   gwtConn%conc, &
    !   gwtConn%icbound)
  end subroutine link_connections
 
  !> @brief Deallocate memory
  !<
  subroutine exg_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(gwfgwtexchangetype) :: this
    !
    call mem_deallocate(this%m1_idx)
    call mem_deallocate(this%m2_idx)
  end subroutine exg_da
 
  !> @brief Allocate package scalars
  !<
  subroutine allocate_scalars(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(gwfgwtexchangetype) :: 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
  !<
  subroutine gwfbnd2gwtfmi(this)
    ! -- dummy
    class(gwfgwtexchangetype) :: this
    ! -- local
    integer(I4B) :: ngwfpack, ip, iterm, imover
    class(basemodeltype), pointer :: mb => null()
    type(gwfmodeltype), pointer :: gwfmodel => null()
    type(gwtmodeltype), pointer :: gwtmodel => 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 gwtmodel
    mb => getbasemodelfromlist(basemodellist, this%m2_idx)
    select type (mb)
    type is (gwtmodeltype)
      gwtmodel => 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 gwtmodel%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 gwtmodel%fmi%gwfpackages(iterm)%set_pointers( &
          'SIMTOMVR', &
          packobj%memoryPath, packobj%input_mempath)
        iterm = iterm + 1
      end if
    end do
  end subroutine gwfbnd2gwtfmi
 
end module gwfgwtexchangemodule