mf6core.f90

Go to the documentation of this file.

!> @brief Core MODFLOW 6 module
!!
!! This module contains the core components for MODFLOW 6. This module
!! is used by the stand-alone executable and the share object versions
!! of MODFLOW 6.
!!
!<
module mf6coremodule
  use kindmodule, only: i4b, lgp
  use listsmodule, only: basesolutionlist, solutiongrouplist, &
                         basemodellist, baseexchangelist, &
                         baseconnectionlist
  use basemodelmodule, only: basemodeltype, getbasemodelfromlist
  use baseexchangemodule, only: baseexchangetype, getbaseexchangefromlist
  use spatialmodelconnectionmodule, only: spatialmodelconnectiontype, &
                                          get_smc_from_list
  use basesolutionmodule, only: basesolutiontype, getbasesolutionfromlist
  use solutiongroupmodule, only: solutiongrouptype, getsolutiongroupfromlist
  use runcontrolmodule, only: runcontroltype
  use simstagesmodule
  use profilermodule
  implicit none
 
  class(runcontroltype), pointer :: run_ctrl => null() !< the run controller for this simulation
 
contains
 
  !> @brief Main controller
    !!
    !! This subroutine is the main controller for MODFLOW 6.
    !!
  !<
  subroutine mf6run
    ! -- modules
    use commandarguments, only: getcommandlinearguments
    use tdismodule, only: endofsimulation
    ! -- local
    logical(LGP) :: hasConverged
    !
    ! -- parse any command line arguments
    call getcommandlinearguments()
    !
    ! initialize simulation
    call mf6initialize()
    !
    ! -- time loop
    do while (.not. endofsimulation)
 
      ! perform a time step
      hasconverged = mf6update()
 
      ! if not converged, break
      if (.not. hasconverged) exit
 
    end do
    !
    ! -- finalize simulation
    call mf6finalize()
 
  end subroutine mf6run
 
  !> @brief Initialize a simulation
    !!
    !! This subroutine initializes a MODFLOW 6 simulation. The subroutine:
    !!   - creates the simulation
    !!   - defines
    !!   - allocates and reads static data
    !!
  !<
  subroutine mf6initialize()
    ! -- modules
    use runcontrolfactorymodule, only: create_run_control
    use simulationcreatemodule, only: simulation_cr
    use sourceloadmodule, only: export_cr
 
    call g_prof%pre_init()
 
    ! -- get the run controller for sequential or parallel builds
    run_ctrl => create_run_control()
    call run_ctrl%start()
 
    ! -- print info and start timer
    call print_info()
 
    ! -- create mfsim.lst
    call create_lstfile()
 
    ! -- load input context
    call static_input_load()
 
    ! init timer and start
    call g_prof%initialize()
    call g_prof%start("Initialize", g_prof%tmr_init)
 
    ! -- create
    call simulation_cr()
 
    ! -- define
    call simulation_df()
 
    ! -- allocate and read
    call simulation_ar()
 
    ! -- create model exports
    call export_cr()
 
    ! -- stop the timer
    call g_prof%stop(g_prof%tmr_init)
 
  end subroutine mf6initialize
 
  !> @brief Run a time step
    !!
    !! This function runs a single time step to completion.
    !!
    !! @return  hasConverged   boolean indicating if convergence was achieved for the time step
    !!
  !<
  function mf6update() result(hasConverged)
    logical(LGP) :: hasconverged
    ! start timer
    call g_prof%start("Update", g_prof%tmr_update)
    !
    ! -- prepare timestep
    call mf6preparetimestep()
    !
    ! -- do timestep
    call mf6dotimestep()
    !
    ! -- after timestep
    hasconverged = mf6finalizetimestep()
 
    ! stop timer
    call g_prof%stop(g_prof%tmr_update)
 
  end function mf6update
 
  !> @brief Finalize the simulation
    !!
    !! This subroutine finalizes a simulation. Steps include:
    !!   - final processing
    !!   - deallocate memory
    !!
  !<
  subroutine mf6finalize()
    ! -- modules
    use, intrinsic :: iso_fortran_env, only: output_unit
    use listsmodule, only: lists_da
    use simulationcreatemodule, only: simulation_da
    use tdismodule, only: tdis_da
    use idmloadmodule, only: idm_da
    use sourceloadmodule, only: export_da
    use simvariablesmodule, only: iout
    ! -- local variables
    integer(I4B) :: im
    integer(I4B) :: ic
    integer(I4B) :: is
    integer(I4B) :: isg
    class(solutiongrouptype), pointer :: sgp => null()
    class(basesolutiontype), pointer :: sp => null()
    class(basemodeltype), pointer :: mp => null()
    class(baseexchangetype), pointer :: ep => null()
    class(spatialmodelconnectiontype), pointer :: mc => null()
    integer(I4B) :: tmr_dealloc
 
    ! start timer
    call g_prof%start("Finalize", g_prof%tmr_finalize)
 
    !
    ! -- FINAL PROCESSING (FP)
    ! -- Final processing for each model
    do im = 1, basemodellist%Count()
      mp => getbasemodelfromlist(basemodellist, im)
      call mp%model_fp()
    end do
    !
    ! -- Final processing for each exchange
    do ic = 1, baseexchangelist%Count()
      ep => getbaseexchangefromlist(baseexchangelist, ic)
      call ep%exg_fp()
    end do
    !
    ! -- Final processing for each solution
    do is = 1, basesolutionlist%Count()
      sp => getbasesolutionfromlist(basesolutionlist, is)
      call sp%sln_fp()
    end do
 
    ! start timer for deallocation
    tmr_dealloc = -1
    call g_prof%start("Deallocate", tmr_dealloc)
 
    !
    ! -- DEALLOCATE (DA)
    ! -- Deallocate tdis
    call tdis_da()
    !
    ! -- Deallocate for each model
    do im = 1, basemodellist%Count()
      mp => getbasemodelfromlist(basemodellist, im)
      call mp%model_da()
      deallocate (mp)
    end do
    !
    ! -- Deallocate for each exchange
    do ic = 1, baseexchangelist%Count()
      ep => getbaseexchangefromlist(baseexchangelist, ic)
      call ep%exg_da()
      deallocate (ep)
    end do
    !
    ! -- Deallocate for each connection
    do ic = 1, baseconnectionlist%Count()
      mc => get_smc_from_list(baseconnectionlist, ic)
      call mc%exg_da()
      deallocate (mc)
    end do
    !
    ! -- Deallocate for each solution
    do is = 1, basesolutionlist%Count()
      sp => getbasesolutionfromlist(basesolutionlist, is)
      call sp%sln_da()
      deallocate (sp)
    end do
    !
    ! -- Deallocate solution group and simulation variables
    do isg = 1, solutiongrouplist%Count()
      sgp => getsolutiongroupfromlist(solutiongrouplist, isg)
      call sgp%sgp_da()
      deallocate (sgp)
    end do
    !
    call idm_da(iout)
    call export_da()
    call simulation_da()
    call lists_da()
 
    ! stop timer
    call g_prof%stop(tmr_dealloc)
 
    ! finish gently (No calls after this)
    ! timer is stopped inside because this call does not return
    call run_ctrl%finish()
 
  end subroutine mf6finalize
 
  !> @brief print initial message
  !<
  subroutine print_info()
    use simmodule, only: initial_message
    use timermodule, only: print_start_time
 
    ! print initial message
    call initial_message()
 
    ! get start time
    call print_start_time()
 
  end subroutine print_info
 
  !> @brief Set up mfsim list file output logging
    !!
    !! This subroutine creates the mfsim list file
    !! and writes the header.
    !!
  !<
  subroutine create_lstfile()
    use constantsmodule, only: linelength
    use simvariablesmodule, only: proc_id, nr_procs, simlstfile, iout
    use inputoutputmodule, only: getunit, openfile, append_processor_id
    use messagemodule, only: write_message
    use versionmodule, only: write_listfile_header
    character(len=LINELENGTH) :: line
    !
    ! -- Open simulation list file
    iout = getunit()
    !
    if (nr_procs > 1) then
      call append_processor_id(simlstfile, proc_id)
    end if
    !
    call openfile(iout, 0, simlstfile, 'LIST', filstat_opt='REPLACE')
    !
    ! -- write simlstfile to stdout
    write (line, '(2(1x,A))') 'Writing simulation list file:', &
      trim(adjustl(simlstfile))
    !
    call write_message(line)
    call write_listfile_header(iout)
  end subroutine create_lstfile
 
  !> @brief Create simulation input context
    !!
    !! This subroutine creates the simulation input context
    !!
  !<
  subroutine static_input_load()
    ! -- modules
    use constantsmodule, only: lenmempath
    use simvariablesmodule, only: iout
    use idmloadmodule, only: simnam_load, simtdis_load, &
                             load_models, load_exchanges
    use memoryhelpermodule, only: create_mem_path
    use memorymanagermodule, only: mem_setptr, mem_allocate
    use simvariablesmodule, only: iparamlog
    !
    ! -- load simnam input context
    call simnam_load(iparamlog)
    !
    ! -- load tdis to input context
    call simtdis_load()
    !
    ! -- load in scope models
    call load_models(iout)
    !
    ! -- load in scope exchanges
    call load_exchanges(iout)
  end subroutine static_input_load
 
  !> @brief Define the simulation
    !!
    !! This subroutine defined the simulation. Steps include:
    !!   - define each model
    !!   - define each solution
    !!
  !<
  subroutine simulation_df()
    ! -- modules
    use idmloadmodule, only: idm_df
    ! -- local variables
    integer(I4B) :: im
    integer(I4B) :: ic
    integer(I4B) :: is
    class(basesolutiontype), pointer :: sp => null()
    class(basemodeltype), pointer :: mp => null()
    class(baseexchangetype), pointer :: ep => null()
    class(spatialmodelconnectiontype), pointer :: mc => null()
 
    ! -- init virtual data environment
    call run_ctrl%at_stage(stg_bfr_mdl_df)
 
    ! -- Define each model
    do im = 1, basemodellist%Count()
      mp => getbasemodelfromlist(basemodellist, im)
      call mp%model_df()
    end do
    !
    ! -- synchronize
    call run_ctrl%at_stage(stg_aft_mdl_df)
    !
    ! -- Define each exchange
    do ic = 1, baseexchangelist%Count()
      ep => getbaseexchangefromlist(baseexchangelist, ic)
      call ep%exg_df()
    end do
    !
    ! -- synchronize
    call run_ctrl%at_stage(stg_aft_exg_df)
    !
    ! -- when needed, this is where the interface models are
    !     created and added to the numerical solutions
    call connections_cr()
    !
    ! -- synchronize
    call run_ctrl%at_stage(stg_aft_con_cr)
    !
    ! -- synchronize
    call run_ctrl%at_stage(stg_bfr_con_df)
    !
    ! -- Define each connection
    do ic = 1, baseconnectionlist%Count()
      mc => get_smc_from_list(baseconnectionlist, ic)
      call mc%exg_df()
    end do
    !
    ! -- synchronize
    call run_ctrl%at_stage(stg_aft_con_df)
    !
    ! -- Define each solution
    do is = 1, basesolutionlist%Count()
      sp => getbasesolutionfromlist(basesolutionlist, is)
      call sp%sln_df()
    end do
 
    ! idm df
    call idm_df()
 
  end subroutine simulation_df
 
  !> @brief Simulation allocate and read
    !!
    !! This subroutine allocates and reads static data for the simulation.
    !! Steps include:
    !!   - allocate and read for each model
    !!   - allocate and read for each exchange
    !!   - allocate and read for each solution
    !!
  !<
  subroutine simulation_ar()
    use distvariablemodule
    ! -- local variables
    integer(I4B) :: im
    integer(I4B) :: ic
    integer(I4B) :: is
    class(basesolutiontype), pointer :: sp => null()
    class(basemodeltype), pointer :: mp => null()
    class(baseexchangetype), pointer :: ep => null()
    class(spatialmodelconnectiontype), pointer :: mc => null()
 
    ! -- Allocate and read each model
    do im = 1, basemodellist%Count()
      mp => getbasemodelfromlist(basemodellist, im)
      call mp%model_ar()
    end do
    !
    ! -- Allocate and read each exchange
    do ic = 1, baseexchangelist%Count()
      ep => getbaseexchangefromlist(baseexchangelist, ic)
      call ep%exg_ar()
    end do
    !
    ! -- Synchronize
    call run_ctrl%at_stage(stg_bfr_con_ar)
    !
    ! -- Allocate and read all model connections
    do ic = 1, baseconnectionlist%Count()
      mc => get_smc_from_list(baseconnectionlist, ic)
      call mc%exg_ar()
    end do
    !
    ! -- Synchronize
    call run_ctrl%at_stage(stg_aft_con_ar)
    !
    ! -- Allocate and read each solution
    do is = 1, basesolutionlist%Count()
      sp => getbasesolutionfromlist(basesolutionlist, is)
      call sp%sln_ar()
    end do
    !
  end subroutine simulation_ar
 
  !> @brief Create the model connections from the exchanges
    !!
    !! This will upgrade the numerical exchanges in the solution,
    !! whenever the configuration requires this, to Connection
    !! objects. Currently we anticipate:
    !!
    !!   GWF-GWF => GwfGwfConnection
    !!   GWT-GWT => GwtGwtConecction
  !<
  subroutine connections_cr()
    use connectionbuildermodule
    use simvariablesmodule, only: iout
    use versionmodule, only: idevelopmode
    integer(I4B) :: isol
    type(connectionbuildertype) :: connectionBuilder
    class(basesolutiontype), pointer :: sol => null()
    integer(I4B) :: status
    character(len=16) :: envvar
 
    write (iout, '(/a)') 'PROCESSING MODEL CONNECTIONS'
 
    if (baseexchangelist%Count() == 0) then
      ! if this is not a coupled simulation in any way,
      ! then we will not need model connections
      return
    end if
 
    if (idevelopmode == 1) then
      call get_environment_variable('DEV_ALWAYS_USE_IFMOD', &
                                    value=envvar, status=status)
      if (status == 0 .and. envvar == '1') then
        connectionbuilder%dev_always_ifmod = .true.
        write (iout, '(/a)') "Development option: forcing interface model"
      end if
    end if
 
    do isol = 1, basesolutionlist%Count()
      sol => getbasesolutionfromlist(basesolutionlist, isol)
      call connectionbuilder%processSolution(sol)
    end do
 
    write (iout, '(a)') 'END OF MODEL CONNECTIONS'
  end subroutine connections_cr
 
  !> @brief Read and prepare time step
    !!
    !! This subroutine reads and prepares period data for the simulation.
    !! Steps include:
    !!   - read and prepare for each model
    !!   - read and prepare for each exchange
    !!   - reset convergence flag
    !!   - calculate maximum time step for each model
    !!   - calculate maximum time step for each exchange
    !!   - calculate maximum time step for each solution
    !!   - set time discretization timestep using smallest maximum timestep
    !!
  !<
  subroutine mf6preparetimestep()
    ! -- modules
    use kindmodule, only: i4b
    use constantsmodule, only: linelength, mnormal, mvalidate
    use tdismodule, only: tdis_set_counters, tdis_set_timestep, &
                          kstp, kper
    use listsmodule, only: basemodellist, baseexchangelist
    use basemodelmodule, only: basemodeltype, getbasemodelfromlist
    use baseexchangemodule, only: baseexchangetype, getbaseexchangefromlist
    use basesolutionmodule, only: basesolutiontype, getbasesolutionfromlist
    use simmodule, only: converge_reset
    use simvariablesmodule, only: isim_mode
    use idmloadmodule, only: idm_rp
    use sourceloadmodule, only: export_post_prepare
    ! -- local variables
    class(basemodeltype), pointer :: mp => null()
    class(baseexchangetype), pointer :: ep => null()
    class(spatialmodelconnectiontype), pointer :: mc => null()
    class(basesolutiontype), pointer :: sp => null()
    character(len=LINELENGTH) :: line
    character(len=LINELENGTH) :: fmt
    integer(I4B) :: im
    integer(I4B) :: ie
    integer(I4B) :: ic
    integer(I4B) :: is
 
    ! start timer
    call g_prof%start("Prepare time step", g_prof%tmr_prep_tstp)
 
    !
    ! -- initialize fmt
    fmt = "(/,a,/)"
    !
    ! -- period update
    call tdis_set_counters()
    !
    ! -- set base line
    write (line, '(a,i0,a,i0,a)') &
      'start timestep kper="', kper, '" kstp="', kstp, '" mode="'
    !
    ! -- evaluate simulation mode
    select case (isim_mode)
    case (mvalidate)
      line = trim(line)//'validate"'
    case (mnormal)
      line = trim(line)//'normal"'
    end select
 
    ! -- load dynamic input
    call idm_rp()
 
    ! -- Read and prepare each model
    do im = 1, basemodellist%Count()
      mp => getbasemodelfromlist(basemodellist, im)
      call mp%model_message(line, fmt=fmt)
      call mp%model_rp()
    end do
    !
    ! -- Synchronize
    call run_ctrl%at_stage(stg_bfr_exg_rp)
    !
    ! -- Read and prepare each exchange
    do ie = 1, baseexchangelist%Count()
      ep => getbaseexchangefromlist(baseexchangelist, ie)
      call ep%exg_rp()
    end do
    !
    ! -- Read and prepare each connection
    do ic = 1, baseconnectionlist%Count()
      mc => get_smc_from_list(baseconnectionlist, ic)
      call mc%exg_rp()
    end do
    !
    ! -- Synchronize
    call run_ctrl%at_stage(stg_aft_con_rp)
    !
    ! -- reset simulation convergence flag
    call converge_reset()
    !
    ! -- time update for each model
    do im = 1, basemodellist%Count()
      mp => getbasemodelfromlist(basemodellist, im)
      call mp%model_dt()
    end do
    !
    ! -- time update for each exchange
    do ie = 1, baseexchangelist%Count()
      ep => getbaseexchangefromlist(baseexchangelist, ie)
      call ep%exg_dt()
    end do
    !
    ! -- time update for each connection
    do ic = 1, baseconnectionlist%Count()
      mc => get_smc_from_list(baseconnectionlist, ic)
      call mc%exg_dt()
    end do
    !
    ! -- time update for each solution
    do is = 1, basesolutionlist%Count()
      sp => getbasesolutionfromlist(basesolutionlist, is)
      call sp%sln_dt()
    end do
    !
    ! -- update exports
    call export_post_prepare()
    !
    ! -- set time step
    call tdis_set_timestep()
 
    ! stop timer
    call g_prof%stop(g_prof%tmr_prep_tstp)
 
  end subroutine mf6preparetimestep
 
  !> @brief Run time step
    !!
    !! This subroutine runs a single time step for the simulation.
    !! Steps include:
    !!   - formulate the system of equations for each model and exchange
    !!   - solve each solution
    !!
  !<
  subroutine mf6dotimestep()
    ! --  modules
    use kindmodule, only: i4b
    use listsmodule, only: solutiongrouplist
    use simvariablesmodule, only: ifailedstepretry
    use solutiongroupmodule, only: solutiongrouptype, getsolutiongroupfromlist
    use idmloadmodule, only: idm_ad
    ! -- local variables
    class(solutiongrouptype), pointer :: sgp => null()
    integer(I4B) :: isg
    logical :: finishedTrying
 
    ! start timer
    call g_prof%start("Do time step", g_prof%tmr_do_tstp)
 
    ! -- By default, the solution groups will be solved once, and
    !    may fail.  But if adaptive stepping is active, then
    !    the solution groups may be solved over and over with
    !    progressively smaller time steps to see if convergence
    !    can be obtained.
    ifailedstepretry = 0
    retryloop: do
 
      ! -- idm advance
      call idm_ad()
 
      do isg = 1, solutiongrouplist%Count()
        sgp => getsolutiongroupfromlist(solutiongrouplist, isg)
        call sgp%sgp_ca()
      end do
 
      call sim_step_retry(finishedtrying)
      if (finishedtrying) exit retryloop
      ifailedstepretry = ifailedstepretry + 1
 
    end do retryloop
 
    ! stop timer
    call g_prof%stop(g_prof%tmr_do_tstp)
 
  end subroutine mf6dotimestep
 
  !> @brief Rerun time step
    !!
    !! This subroutine reruns a single time step for the simulation when
    !! the adaptive time step option is used.
    !!
  !<
  subroutine sim_step_retry(finishedTrying)
    ! -- modules
    use kindmodule, only: dp
    use simvariablesmodule, only: laststepfailed
    use simmodule, only: converge_reset
    use tdismodule, only: kstp, kper, delt, tdis_delt_reset
    use adaptivetimestepmodule, only: ats_reset_delt
    ! -- dummy variables
    logical, intent(out) :: finishedTrying !< boolean that indicates if no
    ! additional reruns of the time step are required
    !
    ! -- Check with ats to reset delt and keep trying
    finishedtrying = .true.
    call ats_reset_delt(kstp, kper, laststepfailed, delt, finishedtrying)
    !
    if (.not. finishedtrying) then
      !
      ! -- Reset delt, which requires updating pertim, totim
      !    and end of period and simulation indicators
      call tdis_delt_reset(delt)
      !
      ! -- Reset state of the simulation convergence flag
      call converge_reset()
 
    end if
  end subroutine sim_step_retry
 
  !> @brief Finalize time step
    !!
    !! This function finalizes a single time step for the simulation
    !! and writes output for the time step. Steps include:
    !!   - write output for each model
    !!   - write output for each exchange
    !!   - write output for each solutions
    !!   - perform a final convergence check and whether the simulation
    !!     can continue if convergence was not achieved
    !!
    !! @return  hasConverged   boolean indicating if convergence was achieved for the time step
    !!
  !<
  function mf6finalizetimestep() result(hasConverged)
    ! -- modules
    use kindmodule, only: i4b
    use constantsmodule, only: linelength, mnormal, mvalidate
    use listsmodule, only: basesolutionlist, basemodellist, baseexchangelist
    use basemodelmodule, only: basemodeltype, getbasemodelfromlist
    use baseexchangemodule, only: baseexchangetype, getbaseexchangefromlist
    use basesolutionmodule, only: basesolutiontype, getbasesolutionfromlist
    use simmodule, only: converge_check
    use simvariablesmodule, only: isim_mode
    use sourceloadmodule, only: export_post_step
    ! -- return variable
    logical(LGP) :: hasconverged
    ! -- local variables
    class(basesolutiontype), pointer :: sp => null()
    class(basemodeltype), pointer :: mp => null()
    class(baseexchangetype), pointer :: ep => null()
    class(spatialmodelconnectiontype), pointer :: mc => null()
    character(len=LINELENGTH) :: line
    character(len=LINELENGTH) :: fmt
    integer(I4B) :: im
    integer(I4B) :: ix
    integer(I4B) :: ic
    integer(I4B) :: is
    !
    ! -- initialize format and line
    fmt = "(/,a,/)"
    line = 'end timestep'
 
    ! start timer
    call g_prof%start("Finalize time step", g_prof%tmr_final_tstp)
 
    !
    ! -- evaluate simulation mode
    select case (isim_mode)
    case (mvalidate)
      !
      ! -- Write final message for timestep for each model
      do im = 1, basemodellist%Count()
        mp => getbasemodelfromlist(basemodellist, im)
        call mp%model_message(line, fmt=fmt)
      end do
    case (mnormal)
 
      call g_prof%start("Write output", g_prof%tmr_output)
      !
      ! -- Write output and final message for timestep for each model
      do im = 1, basemodellist%Count()
        mp => getbasemodelfromlist(basemodellist, im)
        call mp%model_ot()
        call mp%model_message(line, fmt=fmt)
      end do
      !
      ! -- Write output for each exchange
      do ix = 1, baseexchangelist%Count()
        ep => getbaseexchangefromlist(baseexchangelist, ix)
        call ep%exg_ot()
      end do
      !
      ! -- Write output for each connection
      do ic = 1, baseconnectionlist%Count()
        mc => get_smc_from_list(baseconnectionlist, ic)
        call mc%exg_ot()
      end do
      !
      ! -- Write output for each solution
      do is = 1, basesolutionlist%Count()
        sp => getbasesolutionfromlist(basesolutionlist, is)
        call sp%sln_ot()
      end do
      !
      ! -- update exports
      call g_prof%start("NetCDF export", g_prof%tmr_nc_export)
      call export_post_step()
      call g_prof%stop(g_prof%tmr_nc_export)
 
      call g_prof%stop(g_prof%tmr_output)
    end select
    !
    ! -- Check if we're done
    call converge_check(hasconverged)
 
    ! stop timer
    call g_prof%stop(g_prof%tmr_final_tstp)
 
  end function mf6finalizetimestep
 
end module mf6coremodule