!! Copyright (C) Stichting Deltares, 2005-2017.
!!
!! This file is part of iMOD.
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see .
!!
!! Contact: imod.support@deltares.nl
!! Stichting Deltares
!! P.O. Box 177
!! 2600 MH Delft, The Netherlands.
program driver
! modules
use driver_module
use IMOD_UTL, only : imod_utl_capf
use imod_idf_par, only : idfobj
use m_main_info
use imod_utl, only: imod_utl_closeunits, imod_utl_has_ext, imod_utl_printtext,imod_utl_getunit, imod_utl_getslash
use pks_imod_utl, only: pks_imod_utl_idfmerge_init, pks_imod_utl_write_idfmergefile, pks_imod_utl_idfmerge ! PKS
use mod_pest, only: pest1_meteo_metaswap, pest1alpha_metaswap, pest1appendlogfile, pestnext, pestdumpfct, PEST1INIT, PEST1CLOSELOGFILES
use PESTVAR, only : IUPESTOUT
use pks_imod_utl, only: pks_imod_utl_iarmwp_xch_read
use, intrinsic :: ieee_exceptions
use dflowfmmod
use mfmsdfm_mapping
implicit none
! general
logical :: mf_steadystate
integer :: mf_ngrid, mf_igrid, iout, mf_minid, mf_maxid
! iPEST with MetaSWAP
integer :: mf_nrow, mf_ncol, mf_nlay
! parameters
integer, parameter :: nsubmax = 1000
! functions
integer osd_open2,cfn_length
double precision cfn_mjd_nodata
logical pks7mpimasterwrite
! local variables
integer ios,lunc,lunsts,lc,exitcode,ivcl,jvcl,iarg,jarg,narg
logical :: ok
logical endOfSimulation
character del*1,cont*1,comm*1
character (len=1024) :: root, wd, modwd1, modwd2, simwd1, simwd2, dfmwd1, infile, rfopt, pfile, rfroot
! integer :: privsys, srivsys, trivsys, bdrnsys, odrnsys ! KEEP: needed for data exchange arrays
character (len=1024) :: record, modrecord
character (len=32) :: compcls,time_string
character (len=1024) :: compfile
integer :: arglen
logical :: usemodflow,usemetaswap,usedflowfm,usests,usestsmodflow
logical :: converged,convergedMF2005,convergedMetaSwap,convergedMozart,convergedPest
logical :: doTimeLoop
integer :: retValMF2005
double precision :: currentTime,currentTimeMF2005, &
endOfCurrentTimeStep,endOfCurrentTimeStepMF2005, &
nodataTime,BeginOfNextTimeStepMF2005,&
endOfCurrentTimeStepFM
double precision :: dtMetaSwap ! variable for metaswap
double precision :: dtsw ! MetaSWAP surface water time step
integer :: iterMetaSwap ! variable for metaswap
character tekens(1)*1
integer tsc ! debug variable
integer :: date, hour, minute, second
integer :: numTopTimeStepsPerMfTimeStep
double precision :: tStartTop, tEndTop ! start and end time stamp of inner MetaSWAT/DflowFM computation
integer :: idtsw ! timestep counter inner MetaSWAT/DflowFM steps
logical :: lrunfile, lnamfile, llpf, lipest, lpwt, lss, lrf, psolved
logical :: lwstdo, lpks, lidfmerge
character(len=1024) :: idfmergefile
character(len=256) :: mduFile, dimr_config
integer :: isub, nsub, nnsub
character(len=50), dimension(nsubmax) :: submstr=''
real :: hnoflo
type(idfobj) :: idf
character(len=1024) :: str
real, dimension(4) :: imodusebox
character(len=1) :: slash
! coupling with Dflow-FM
! . parse mapping lines in config file
character(len=32 ) :: fromComp, toComp, locType
character(len=256) :: mapConfig
! . handles to mapping objects
integer :: mapperDfm2dToMsw = -1
integer :: mapperMswToDfm1d = -1
integer :: mapperMswToDfm2d = -1
integer :: mapperDfm1dToMfRiv = -1
integer :: mapperMfRivToDfm1d = -1
! arrays for exchange
integer :: nSimSvats ! #svats in metaswap schematization
integer :: nmodriv ! #1D-rivers in modflow schematization
integer :: nFlowNodes ! #1D-nodes in D-Flow FM schematization
integer :: nFlowCells ! #2D-nodes in D-Flow FM schematization
! id-arrays for setting up exchange
integer, dimension(:), allocatable :: locationIds
logical :: success
integer :: retVal
integer :: riv
!
! value-arrays for the (mapped) exchanged values
!
! D-Flow FM -> MetaSwap, waterlevels:
double precision, dimension(:), allocatable :: XchFromDfmWaterlevelsOnCells
double precision, dimension(:), allocatable :: XchToMetaswapWaterlevelsOnSvats
! MetaSwap -> D-Flow FM, ponding volumes:
double precision, dimension(:), allocatable :: XchFromMetaswapDeltaVolumesOnSvats
double precision, dimension(:), allocatable :: XchToDfmDeltaVolumesOnCells
double precision, dimension(:), allocatable :: XchFromDfmRealizedFluxes2D
double precision, dimension(:), allocatable :: XchToMswRealizedDeltaVolumesOnSvats
! D-Flow FM -> iMODFLOW, waterlevels in 1D nodes:
double precision, dimension(:), allocatable :: XchFromDfmWaterlevelsOnRivers
double precision, dimension(:), allocatable :: XchToModflowWaterlevelsOnRivers
! iMODFLOW -> D-Flow FM, fluxes and realized fluxes in 1D nodes:
double precision, dimension(:), allocatable :: XchFromModflowFluxesInRivers
double precision, dimension(:), allocatable :: XchToDfmFluxesInRivers
double precision, dimension(:), allocatable :: XchFromDfmRealizedFluxesInRivers
double precision, dimension(:), allocatable :: XchToModflowCorrectedFluxesInRivers
! MetaSwap -> D-Flow FM, sprinkling:
! (Toward DFM the sprinkling demands will be add to the fluxes on 1D nodes
! the realized sprinkling demands will be derived from the realized fluxes on 1D nodes)
double precision, dimension(:), allocatable :: XchFromMetaSwapSprinklingDemand
double precision, dimension(:), allocatable :: XchToMetaSwapRealizedSprinklingDemand
! debug
integer ncvgerr
! program section
! ------------------------------------------------------------------------------
CALL IEEE_SET_HALTING_MODE (IEEE_DIVIDE_BY_ZERO, .TRUE.)
call dmmInitialize()
call pks7mpiini1(lwstdo) ! PKS
call pks7mpiactive(lpks) ! PKS
! Evaluate iMOD license
if(pks7mpimasterwrite()) call imod_license()
call pks7mpibarrier() ! PKS
! ... init
exitcode = 0 ! when 0: ok
lunc = 9 ! unit number for components file
lunsts = 8 ! <=0: sts2init() will assign a free unit number
! >0: unit number will be used
comm = '!' ! comment sign for components input file
cont = char(92) ! continuation mark for components input file records
nodataTime = cfn_mjd_nodata()
nsub = 0
NCVGERR = 0
! ... process initialisation
! ... init components
! ios: 0: OK
! -1: commandline argument not found
! <-1: ERROR
! define virtual commandline
call cfn_vcl_set(' ',ivcl)
! get the root name
call osd_getcwd(root)
modwd1 = root
simwd1 = root
dfmwd1 = root
! get slash
call imod_utl_getslash(slash)
! run-file
lrunfile = .false.
lnamfile = .false.
usemodflow = .false.
usemetaswap = .false.
lipest = .false.
lidfmerge = .false.
call cfn_vcl_narg(ivcl,narg)
if (narg.eq.0) then
call imod_utl_printtext(' ',3)
call imod_utl_printtext('Driver program valid arguments:',3)
call imod_utl_printtext(' 1: ',3)
call imod_utl_printtext(' 2: ',3)
call imod_utl_printtext(' 3: -components ',3)
call imod_utl_printtext(' 4: -pksmergeidf ',3)
call imod_utl_printtext(' 5: -ipest ',3)
call pks7mpifinalize()! PKS
call exit(0)
end if
! option for merging PKS output IDF files
call cfn_vcl_fndc(ivcl,iarg,'-pksmergeidf',.true.,idfmergefile,1)
if (iarg.gt.0) then
call pks_imod_utl_idfmerge(idfmergefile)
stop
end if
nmodriv = 0
call cfn_vcl_fndi(ivcl,iarg,'-nriv',.true.,nmodriv,1)
call cfn_vcl_arg(ivcl,1,infile,arglen)
if (imod_utl_has_ext(infile,'nam')) lnamfile= .true.
if (imod_utl_has_ext(infile,'run')) lrunfile= .true.
if (lrunfile .or. lnamfile) then ! run-file or nam-file
usemodflow=.true.
modrecord=''
if (lnamfile) write(modrecord,'(a,1x,3a)') '-namfile','"',trim(infile),'"'
if (lrunfile) then
write(modrecord,'(a,1x,3a)') '-runfile','"',trim(infile),'"'
if (narg.eq.3) then ! runfile options
call cfn_vcl_arg(ivcl,2,rfopt,arglen)
write(modrecord,'(a,1x,a,1x,a)') trim(modrecord),'-rfopt',trim(rfopt)
call cfn_vcl_arg(ivcl,3,rfopt,arglen)
write(modrecord,'(a,1x,a)') trim(modrecord),trim(rfopt)
end if
end if
wd = ''
call rf2mf_prg(lrunfile,lipest,lidfmerge,modrecord,usemetaswap,submstr,nsub,nsubmax,wd,imodusebox,rfroot,nmodriv)
modwd1 = trim(wd)
call osd_s_filename(modwd1)
if (usemetaswap) then
simwd1 = trim(wd)
end if
end if
usests = .false.
usestsmodflow = .false.
usedflowfm = .false.
! get explicit arguments
call cfn_vcl_fndc(ivcl,iarg,'-components',.true.,compfile,1)
if (iarg.gt.0) then
! open file and initialise components
ios=osd_open2(lunc,0,compfile,'readonly,shared')
if (ios.eq.0) then
! read all component information and initialise all components
call cfn_getrec(lunc,record,comm,cont)
do while (cfn_length(record).gt.0)
! extract the first argument of variable record,
! leave the rest of the record for initialisation of the component
tekens(1) = ' '
call cfn_par_ext(compcls,record,tekens,1,del)
call cfn_token(compcls,'tu')
lc=cfn_length(compcls)
select case( compcls(1:lc) )
case( 'MODFLOW' )
modrecord = record
usemodflow=.true.
! set virtual command line
call cfn_vcl_set(modrecord,jvcl)
! set working directory
call cfn_vcl_fndc(jvcl,jarg,'-wd',.true.,wd,1)
if (jarg.gt.0) modwd1 = trim(modwd1)//'\'//wd
call osd_chdir(modwd1)
! convert iMOD run-file to MODFLOW
call rf2mf_prg(lrf,lipest,lidfmerge,modrecord,usemetaswap,submstr,nsub,nsubmax,modwd1,imodusebox,rfroot,nmodriv)
if (lrf) then
modwd1 = trim(modwd1)//slash//'mf2005_tmp'
call osd_s_filename(modwd1)
end if
#ifdef INCLUDE_METASWAP
case( 'METASWAP' )
usemetaswap=.true.
! set working directory
call cfn_vcl_set(record,jvcl)
call cfn_vcl_fndc(jvcl,jarg,'-wd',.true.,wd,1)
if (jarg.gt.0) simwd1 = trim(simwd1)//'\'//wd
#endif
#ifdef INCLUDE_DFLOWFM
case( 'DFLOWFM')
usedflowfm = .true.
! set working directory
call cfn_vcl_set(record,jvcl)
call cfn_vcl_fndc(jvcl,jarg,'-wd',.true.,wd,1)
if (jarg.gt.0) dfmwd1= trim(dfmwd1)//'\'//wd
call cfn_vcl_fndc(jvcl,jarg,'-mdu',.true.,mduFile,1)
if (jarg<=0) then
! TODO error
endif
retVal = dflowfm_LoadLibrary('dflowfm')
print *, 'LOAD dflowfm Library result: ', retVal
call osd_chdir(dfmwd1)
retVal = dflowfm_InitSimulation(mduFile)
call osd_chdir(modwd1)
case( 'DIMR')
usedflowfm = .true.
! set working directory
call cfn_vcl_set(record,jvcl)
call cfn_vcl_fndc(jvcl,jarg,'-wd',.true.,wd,1)
if (jarg.gt.0) dfmwd1= trim(dfmwd1)//'\'//wd
call cfn_vcl_fndc(jvcl,jarg,'-config',.true.,dimr_config,1)
if (jarg<=0) then
! TODO error
endif
retVal = dflowfm_LoadLibrary('dimr')
print *, 'LOAD dimr Library result: ', retVal
call osd_chdir(dfmwd1)
retVal = dflowfm_InitSimulation(dimr_config)
call osd_chdir(modwd1)
case( 'MAPPING')
! set working directory
wd = '.'
call cfn_vcl_set(record,jvcl)
call cfn_vcl_fndc(jvcl,jarg,'-wd',.true.,wd,1)
call cfn_vcl_fndc(jvcl,jarg,'from',.true.,fromComp,1)
call cfn_vcl_fndc(jvcl,jarg,'to',.true.,toComp,1)
call cfn_vcl_fndc(jvcl,jarg,'-locType',.true.,locType,1)
call cfn_vcl_fndc(jvcl,jarg,'-config',.true.,mapConfig,1)
mapConfig = trim(root)//'\'//trim(wd)//'\'//trim(mapConfig)
if (locType=='FM2MSW_WL') then
mapperDfm2dToMsw = dmmCreateMapping(mapConfig, dmm_FM2MSW_WL, 'DFLOWFM', 'METASWAP')
else if (locType=='MSW2FM_DPV') then
mapperMswToDfm2d = dmmCreateMapping(mapConfig, dmm_MSW2FM_DPV, 'METASWAP', 'DFLOWFM')
else if (locType=='sprinkling') then
mapperMswToDfm1d = dmmCreateMapping(mapConfig, dmm_MSW2FM_SPR, 'METASWAP', 'DFLOWFM')
else if (locType=='FM2MF_WL') then
mapperDfm1dToMfRiv = dmmCreateMapping(mapConfig, dmm_FM2MF_WL, 'DFLOWFM', 'MODFLOW')
else if (locType=='MF2FM_Q') then
mapperMfRivToDfm1d = dmmCreateMapping(mapConfig, dmm_MF2FM_Q, 'MODFLOW', 'DFLOWFM')
else if (.false.) then
! TODO: 1D-laterals
endif
#endif
case default
! ERROR, component unknown
if (pks7mpimasterwrite()) write(*,*) ' Warning: component ',compcls(1:cfn_length(compcls)),' unknown.'
end select
! get next component
call cfn_getrec(lunc,record,comm,cont)
enddo
! close components file
close(lunc)
else
call driverChk(.false.,'Could not open component description file '//compfile(1:cfn_length(compfile)))
call pks7mpifinalize()! PKS
call exit(1)
endif
endif
! check if user has specified -ipest
if(lnamfile)then
lipest = .false.
call cfn_vcl_fndc(ivcl,iarg,'-ipest',.true.,pfile,1)
if (iarg.gt.0) then
lipest=.true.
CALL PEST1INIT(1,pfile,IUPESTOUT,modwd1,idf,0) !'CODE OF HET UIT RUN- OF NAMFILE KOMT')
endif
endif
rt = driverGetRunType(usemodflow,usemetaswap,usedflowfm)
str = 'Running mode: '//trim(rtdesc(rt))
if(lpks) str = trim(str)//' (Parallel Krylov Solver activated)'
if (pks7mpimasterwrite()) write(*,'(79(''*''),/,(1x,a),/,79(''*''))') trim(str)
! check state-save options
if (usests.and..not.usestsmodflow) call driverChk(.false.,'MODFLOW state-save not activated.')
if (.not.usests.and.usestsmodflow) call driverChk(.false.,'MODFLOW state-save activated.')
nnsub = max(1,nsub)
! ######################################### SUBMODEL LOOP #########################################
submodelloop: do isub = 1,nnsub
! ######################################### PEST LOOP #########################################
call pks_imod_utl_idfmerge_init() ! PKS
convergedPest = .false.
pestloop: do while (.not.convergedPest)
if (lrunfile) then
if (nsub.gt.0) then
write(modwd2,'(5a)') trim(modwd1),slash,trim(submstr(isub)),slash,'mf2005_tmp'
if (usemetaswap) then
write(simwd2,'(5a)') trim(simwd1),slash,trim(submstr(isub)),slash,'metaswap'
end if
write(*,'(50(''+''),/,1x,a,1x,a,1x,a,i3.3,a,i3.3,a,/,50(''+''))') 'Computing for submodel:',trim(submstr(isub)),'(',isub,'/',nsub,')'
else
write(modwd2,'(3a)') trim(modwd1),slash,'mf2005_tmp'
if (usemetaswap) then
write(simwd2,'(3a)') trim(simwd1),slash,'metaswap'
end if
end if
else
modwd2 = modwd1
simwd2 = simwd1
end if
! PKS IARMWP option
call pks_imod_utl_iarmwp_xch_read(modwd2)
timestep = 1
! init return values
retValMF2005 = 0
! time values
currentTime = nodataTime
currentTimeMF2005 = 0.d0
! convergence
convergedMF2005 = .true.
convergedMetaSwap = .true.
convergedMozart = .true.
endOfSimulation = .false.
! init all components TODO!!!!!!!
if (usemodflow) then
call osd_chdir(modwd2)
if (usestsmodflow) call sts2init(usestsmodflow,lunsts)
call mf2005_initComponent(modrecord,retValMF2005)
! get number of grids
ok = mf2005_PutNumberOfGrids(mf_ngrid); call driverChk(ok,'mf2005_PutNumberOfGrids')
if (mf_ngrid.ne.1) call driverChk(.false.,'More than one MODFLOW grids is not supported.'); mf_igrid = 1
if (lipest) then
ok = mf2005_GetPestFlag(lipest); call driverChk(ok,'mf2005_GetPestFlag')
end if
end if
if (usemetaswap) then
call osd_chdir(simwd2)
if (lipest) then
mf_nrow = 0; mf_ncol = 0; lpwt = .false.
ok = mf2005_PutGridDimensions(mf_igrid, mf_nrow, mf_ncol, mf_nlay); call driverChk(ok,'mf2005_PutGridDimensions')
ok = mf2005_PutPWTActive(mf_igrid, lpwt); call driverChk(ok,'mf2005_PutPWTActive')
call pest1_meteo_metaswap()
call pest1alpha_metaswap(mf_nrow,mf_ncol,lpwt)
end if
call MetaSwap_initComponent()
end if
if (usedflowfm) then
! retVal = dflowfm_LoadLibrary()
! call osd_chdir(dfmwd1)
! retVal = dflowfm_InitSimulation(mduFile)
! call osd_chdir(root)
end if
! append the PEST log-file
if (lipest) then
call pest1appendlogfile(modwd1)
call pest1log()
CALL PEST1CLOSELOGFILES()
end if
! check if MODFLOW is activated
if (retValMF2005.ne.0) call driverChk(.false.,'MODFLOW initialization failed.')
! get starting date for this simulation
call mf2005_getCurrentTime(currentTimeMF2005,retValMF2005)
currentTime = currentTimeMF2005
! start sts
if (usestsmodflow) call sts2start(currentTime)
! read data for entire simulation
call osd_chdir(modwd2)
call mf2005_initSimulation(currentTime,retValMF2005)
if (retValMF2005.ne.0) exitcode = -12
! Coupling MODFLOW-MetaSWAP phase 1
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
ok = mf2005_PutModSimNumberOfIDs(mf_igrid, mf_minid, mf_maxid, XchModSimModNID); call driverChk(ok,'mf2005_PutModSimNumberOfIDs') ! MODFLOW puts the number of MOD-SIM IDs
call driverXchInitModSim1() ! allocate and init
ok = mf2005_PutModSimIDs(mf_igrid,XchModSimModIds); call driverChk(ok,'mf2005_PutModSimIDs') ! MODFLOW puts the MOD-SIM exchange IDs
ok = mf2005_PutModSimCells(mf_igrid,XchModSimModCells); call driverChk(ok,'mf2005_PutModSimCells') ! MODFLOW puts the MOD-SIM cells
call pks7mpimssetids(XchModSimModIds,XchModSimModNID,mf_minid, mf_maxid) ! PKS
call osd_chdir(simwd2)
call metaswap_initSimulation(currentTime)
end if
! #### BEGIN EXCHANGE: BeforeInitSimulationMozart #############################
ok = mf2005_PutLPFActive(mf_igrid, llpf); call driverChk(ok,'mf2005_PutLPFActive') ! get flag if BCF is used or LPF
ok = mf2005_PutSimulationType(mf_igrid, lss); call driverChk(ok,'mf2005_PutSimulationType') ! get flag if simulation is steady-state or transient
ok = mf2005_PutHeadNoFlo(mf_igrid, hnoflo); call driverChk(ok,'mf2005_PutHeadNoFlo') ! get hnoflo
! append the PEST log-file
if (lipest) then
call mf2005_returnIOUT(iout)
CALL PESTDUMPFCT(modwd1,iout)
endif
!#### TIMESERIES ####
call osd_chdir(root)
call tserie1init1(lipest,lss,hnoflo)
ok = mf2005_TimeserieInit(mf_igrid); call driverChk(ok,'mf2005_TimeserieInit')
call tserie1init2(lipest,lss,hnoflo,modwd1,submstr(isub))
if (rt.eq.rtmodsimdflowfm) then
! call osd_chdir(dfmwd1)
! todo
!convergedMozart = .false.
endOfSimulation = .false.
!call mozart_prepareTimeStep(endOfSimulation,convergedMozart,currentTime) ! read signal file
end if
! #### BEGIN EXCHANGE: BeforeTimestep #########################################
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
ok = mf2005_PutSimulationType(mf_igrid, mf_steadystate); call driverChk(ok,'mf2005_PutSimulationType') ! MODFLOW puts the simulation type (transient/steady-state)
if (mf_steadystate) call driverChk(ok,'MODFLOW may not be steady-state')
end if
! Coupling MODFLOW-MetaSWAP phase 2
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
ok = metaswap_PutModSimNumberOfIDs(XchModSimSimNID); call driverChk(ok,'metaswap_PutModSimNumberOfIDs') ! MetaSWAP puts the number of MOD-SIM IDs
call driverXchInitModSim2() ! allocate and init
ok = metaswap_PutModSimIDs(XchModSimSimIds); call driverChk(ok,'metaswap_PutModSimIDs') ! get exchange id's
end if
! Coupling iMODFLOW / MetaSWAP / D-Flow FM
if (rt.eq.rtmodsimdflowfm) then
end if
! Create mappings
!call driverXchIniMapDummy()
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) call driverXchIniMapModSim()
if (rt.eq.rtmodsimdflowfm) then
! set up D-Flow FM / Metaswap mapping for 2D-cells
! #locations
nFlowCells = dflowfm_PutNumberOf2DGridCells()
success = metaswap_PutNumberOfSvats(nSimSvats)
! get the ids of DFM 2D-cells and feed them to the mappers
allocate(locationIds(nFlowCells))
success = dflowfm_Put2DGridCellIDs(locationIds)
retVal = dmmSetSourceIds(mapperDfm2dToMsw, locationIds)
retVal = dmmSetTargetIds(mapperMswToDfm2d, locationIds)
deallocate(locationIds)
! get the ids of the metaswap SVATS and feed them to the mappers
allocate(locationIds(nSimSvats))
success = metaswap_PutSvatIDs(locationIds)
retVal = dmmSetSourceIds(mapperMswToDfm2d, locationIds)
retVal = dmmSetTargetIds(mapperDfm2dToMsw, locationIds)
if (mapperMswToDfm1d > 0) then
retVal = dmmSetSourceIds(mapperMswToDfm1d, locationIds)
endif
deallocate(locationIds)
! allocatie exchange variables for the D-Flow FM / Metaswap coupling
allocate(XchFromMetaswapDeltaVolumesOnSvats(nSimSvats))
allocate(XchToDfmDeltaVolumesOnCells(nFlowCells))
allocate(XchFromDfmRealizedFluxes2D(nFlowCells))
allocate(XchToMswRealizedDeltaVolumesOnSvats(nSimSvats))
allocate(XchFromDfmWaterlevelsOnCells(nFlowCells))
allocate(XchToMetaswapWaterlevelsOnSvats(nSimSvats))
allocate(XchFromMetaSwapSprinklingDemand(nSimSvats))
allocate(XchToMetaSwapRealizedSprinklingDemand(nSimSvats))
! set up D-Flow FM / mf2005 mapping for 1D-nodes/rivers
! #locations
nFlowNodes = dflowfm_PutNumberOf1DNodes()
! get the ids of the DFM 1D nodes and feed them to the mappers
allocate(locationIds(nFlowNodes))
success = dflowfm_Put1DNodeIDs(locationIds)
retVal = dmmSetSourceIds(mapperDfm1dToMfRiv, locationIds)
retVal = dmmSetTargetIds(mapperMfRivToDfm1d, locationIds)
if (mapperMswToDfm1d > 0) then
retVal = dmmSetTargetIds(mapperMswToDfm1d, locationIds)
endif
deallocate(locationIds)
! get (=generate) the ids of the Modflow river elements, feed them to the mappers
ok = mf2005_PutNumberOfRivers(nmodriv); call DriverChk(ok,'mf2005_PutNumberOfRivers') !
allocate(locationIds(nmodriv))
do riv= 1, nmodriv
locationIds(riv) = riv
enddo
retVal = dmmSetSourceIds(mapperMfRivToDfm1d, locationIds)
retVal = dmmSetTargetIds(mapperDfm1dToMfRiv, locationIds)
deallocate(locationIds)
! allocatie exchange variables for the D-Flow FM / Modflow coupling
allocate(XchFromDfmWaterlevelsOnRivers(nFlowNodes))
allocate(XchToModflowWaterlevelsOnRivers(nmodriv))
allocate(XchFromModflowFluxesInRivers(nmodriv))
allocate(XchToDfmFluxesInRivers(nFlowNodes))
allocate(XchFromDfmRealizedFluxesInRivers(nFlowNodes))
allocate(XchToModflowCorrectedFluxesInRivers(nmodriv))
! combine the iMod fluxes and MetaSwap sprinkling demands
if (mapperMswToDfm1d > 0) then
retVal = dmmCombineMappers(mapperMfRivToDfm1d, mapperMswToDfm1d)
endif
end if
! ###### END EXCHANGE: BeforeTimestep #########################################
! ... timestep
timestep = 0
if (.not. dmmCheckLogTimeStep(timestep)) then
stop
endif
tsc=0 ! debug: TimeStepCycle
if (rt.eq.rtmodsimdflowfm) then
! initalize MetaSWAP surface water
ok = DflowFM_Put2DWaterlevels(XchFromDfmWaterlevelsOnCells); call DriverChk(ok,'DflowFM_Put2DWaterlevels') ! initial exchange of water levels; units: m depth
call dmmMapValues(mapperDfm2dToMsw, XchFromDfmWaterlevelsOnCells, XchToMetaswapWaterlevelsOnSvats)
ok = metaswap_Get2DWaterlevels(XchToMetaswapWaterlevelsOnSvats,mv); call DriverChk(ok,'metaswap_Get2DWaterlevels') ! units: m depth
call metaswap_sw_initcomponent()
! get the metaswap dtsw
ok = metaswap_PutSurfacewaterTimestep(dtsw); call DriverChk(ok,'metaswap_PutSurfacewaterTimestep')
endif
! compute the river fluxes for intial heads
if (rt.eq.rtmodsimdflowfm) then
call mf2005_RivBudgetShead()
endif
do while (.not.endOfSimulation .and. exitcode.eq.0)
tsc=tsc+1
timestep = timestep + 1
if (.not. dmmCheckLogTimeStep(timestep)) then
stop
endif
call cfn_mjd2datehms(currentTime,date,hour,minute,second)
!##### BEGIN EXCHANGE: BeginOfTimestep ########################################
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
ok = mf2005_PutHeads(mf_igrid,XchModSimModCells,XchModSimModNID,XchModSimModHeads,mv); call DriverChk(ok,'mf2005_PutHeads') ! Put groundwater heads
ok = metaswap_GetHeads(XchModSimModHeads,XchModSimSimNID,XchMod2SimIdx,XchMod2SimOff,mv); call DriverChk(ok,'metaswap_GetHeads') ! Get groundwater heads
XchModSimModHeads = mv
end if
!####### END EXCHANGE: BeginOfTimestep ########################################
if (rt.eq.rtmodsimdflowfm) then
! provide FM's water levels to MODFLOW
ok = dflowfm_Put1DNodeSurfacewaterLevels(XchFromDfmWaterlevelsOnRivers); call DriverChk(ok,'dflowfm_Put1DNodeSurfacewaterLevels') ! 1D water levels at start of MF-step; units: m AD NAP
call dmmMapValues(mapperDfm1dToMfRiv, XchFromDfmWaterlevelsOnRivers, XchToModflowWaterlevelsOnRivers)
ok = mf2005_Get1DWaterlevels(XchToModflowWaterlevelsOnRivers,mv); call DriverChk(ok,'mf2005_Get1DWaterlevels') ! units: m AD NAP
endif
! ... init timestep for MODFLOW
call mf2005_prepareTimestep(currentTime,.false.,retValMF2005)
! compute MODFLOW river fluxes based on !BM: calculating V1D using FM-levels
! current groundwater heads and current FM-levels !BM: calculating V1D using FM-levels
if (rt.eq.rtmodsimdflowfm) then !BM: calculating V1D using FM-levels
call mf2005_RivBudgetShead() !BM: calculating V1D using FM-levels
endif !BM: calculating V1D using FM-levels
if (rt.eq.rtmodsimdflowfm) then
! provide MODFLOW fluxes to FM
ok = mf2005_Put1DFluxRealised(XchFromModflowFluxesInRivers); call DriverChk(ok,'mf2005_Put1DFluxRealised') ! 1D fluxes; units: m3 for next dtgw (from computed previous dtgw)
! TODO: vervangen door:
! hier: dmmCombinedMapperSetFluxes1(mapperMfRivToDfm1d, XchFromModflowFluxesInRivers)
! verderop: dmmCombinedMapperSetFluxes1(mapperMswToDfm1d, XchFromMetaSwapSprinklingDemand)
! (reg.718) dmmMapCombined values(mapperMfRivToDfm1d, mapperMswToDfm1d, XchToDfmFluxesInRivers)
! dflowfm_Get1DNodeFlux(XchToDfmFluxesInRivers, mv)
call dmmMapValues(mapperMfRivToDfm1d, XchFromModflowFluxesInRivers, XchToDfmFluxesInRivers)
ok = dflowfm_Get1DNodeFlux(XchToDfmFluxesInRivers, mv); call DriverChk(ok,'dflowfm_Get1DNodeFlux')
endif
call osd_chdir(modwd2)
call mf2005_initTimeStep(currentTime,retValMF2005)
if (retValMF2005.ne.0) exitcode = -15
! get end of current timestep
!!!!!!!!!!!!!!!! Only MODFLOW information is used !!!!!!!!!!!!!!!!!!!!!!!
call mf2005_getEndOfCurrentTimeStep(endOfCurrentTimeStepMF2005,retValMF2005)
if (retValMF2005.ne.0) exitcode = -17
endOfCurrentTimeStep = endOfCurrentTimeStepMF2005
! compute the number of MetaSWAP surface water iterations
if (rt.eq.rtmodsimdflowfm) then
numTopTimeStepsPerMfTimeStep = nint((endOfCurrentTimeStep-currentTime)/dtsw)
else
numTopTimeStepsPerMfTimeStep = 0
endif
! let MetaSwap compute sprinkling demand
if (exitcode.eq.0) then
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
dtMetaSwap = 1.D0 ! fixed
call osd_chdir(simwd2)
call MetaSwap_initTimestep(dtMetaSwap)
iterMetaSwap = 0
end if
end if
! MetaSWAP-TOP + DflowFM loop
do idtsw = 1, numTopTimeStepsPerMfTimeStep
tStartTop = currentTime + (idtsw-1)*dtsw
tEndTop = currentTime + idtsw*dtsw
call MetaSWAP_performSurfacewaterTimestep(idtsw)
ok = metaswap_PutPondingDeltaVolumes(XchFromMetaswapDeltaVolumesOnSvats); call DriverChk(ok,'metaswap_PutPondingDeltaVolumes') ! units: m3 for dtsw
call dmmMapValues(mapperMswToDfm2d, XchFromMetaswapDeltaVolumesOnSvats, XchToDfmDeltaVolumesOnCells)
ok = DflowFM_GetPondingDeltaVolumes(XchToDfmDeltaVolumesOnCells, mv); call DriverChk(ok,'DflowFM_GetPondingDeltaVolumes') ! units: m3 for dtsw
if (mapperMswToDfm1d > 0) then
ok = MetaSWAP_PutSprinklingDemand(XchFromMetaSwapSprinklingDemand); call DriverChk(ok,'') ! TODO; units: m3 for next dtgw
call dmmMapValues(mapperMswToDfm1d, XchFromMetaSwapSprinklingDemand, XchToDfmFluxesInRivers)
ok = dflowfm_Get1DNodeFlux(XchToDfmFluxesInRivers, mv); call DriverChk(ok,'') ! TODO; units: m3/s
endif
print *, 'FM TS START'
call dflowfm_performSurfacewaterTimestep(tStartTop, tEndTop)
endOfCurrentTimeStepFM = 0.0D+0 ! TODO: conversion
print *, 'FM TS DONE'
ok = dflowfm_PutRealised2DCellFlux(XchFromDfmRealizedFluxes2D); call DriverChk(ok,'dflowfm_PutRealised2DCellFlux')
call dmmMapValuesBack(mapperMswToDfm2d, XchFromDfmRealizedFluxes2D, XchToMswRealizedDeltaVolumesOnSvats)
ok = MetaSWAP_GetDeltaPondingAllocation(XchToMswRealizedDeltaVolumesOnSvats, mv)
if (mapperMswToDfm1d > 0) then
ok = dflowfm_PutRealised1DNodeFlux(XchFromDfmRealizedFluxesInRivers); call DriverChk(ok,'dflowfm_PutRealisedLateralFlux')
call dmmMapValues(mapperMswToDfm1d, XchFromDfmRealizedFluxesInRivers, XchToMetaSwapRealizedSprinklingDemand)
ok = MetaSWAP_GetSprinklingAllocation(XchToMetaSwapRealizedSprinklingDemand, mv); call DriverChk(ok,'MetaSWAP_GetSprinklingAllocation')
else
XchToMetaSwapRealizedSprinklingDemand = 0.0D+0
ok = MetaSWAP_GetSprinklingAllocation(XchToMetaSwapRealizedSprinklingDemand, mv)
endif
! Put the water levels for the next MetaSwap time step
! (MetaSwap needs them for finishing its current time step)
ok = DflowFM_Put2DWaterlevels(XchFromDfmWaterlevelsOnCells); call DriverChk(ok,'DflowFM_Put2DWaterlevels') ! m AD NAP
call dmmMapValues(mapperDfm2dToMsw, XchFromDfmWaterlevelsOnCells, XchToMetaswapWaterlevelsOnSvats)
ok = metaswap_Get2DWaterlevels(XchToMetaswapWaterlevelsOnSvats,mv); call DriverChk(ok,'metaswap_Get2DWaterlevels') ! m AD NAP
call MetaSWAP_finishSurfacewaterTimestep(idtsw)
enddo
if (rt.eq.rtmodsimdflowfm) then
ok = dflowfm_PutRealised1DNodeFlux(XchFromDfmRealizedFluxesInRivers); call DriverChk(ok,'dflowfm_Put1DFluxRealised') ! 1D fluxes; units: m3 for next dtgw (from computed previous dtgw)
call dmmMapValuesBack(mapperMfRivToDfm1d, XchFromDfmRealizedFluxesInRivers, XchToModflowCorrectedFluxesInRivers)
ok = mf2005_Get1DCorrectionFlux(XchToModflowCorrectedFluxesInRivers, mv); call DriverChk(ok,'mf2005_Get1DNodeFlux') ! 1D water levels at start of MF-step; units: m3/s
! ok = dflowfm_Put1DNodeSurfacewaterLevels(XchFromDfmWaterlevelsOnRivers); call DriverChk(ok,'dflowfm_Put1DNodeSurfacewaterLevels') ! 1D water levels at start of MF-step; units: m AD NAP
! call dmmMapValues(mapperDfm1dToMfRiv, XchFromDfmWaterlevelsOnRivers, XchToModflowWaterlevelsOnRivers)
! ok = mf2005_Get1DWaterlevels(XchToModflowWaterlevelsOnRivers,mv); call DriverChk(ok,'mf2005_Get1DWaterlevels') ! units: m AD NAP
endif
! extra check
if (endOfCurrentTimeStep.eq.nodataTime) call driverChk(.false.,'endOfCurrentTimeStep = nodataTime')
if (pks7mpimasterwrite()) call mf2005_writeTimeStep(tsc,date,hour,minute,second)
! one timestep for each cycle
! ... iteration
converged =.false.
do while (.not. converged .and. exitcode.eq.0)
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
call MetaSwap_prepareIter(iterMetaSwap) ! only action: iter=iter+1
call MetaSwap_performIter(iterMetaSwap)
!##### BEGIN EXCHANGE: AfterSolveMetaSwap #####################################
ok = metaswap_PutModSimUnsaturatedZoneFlux(XchModSimSimUnsZFlux,mv); call DriverChk(ok,'metaswap_PutUnsaturatedZoneFlux') ! Put unsaturated zone flux
ok = mf2005_GetUnsaturatedZoneFlux(mf_igrid,&
XchModSimModNID,XchModSimSimUnsZFlux,XchSim2ModIdx,XchSim2ModOff,mv); call DriverChk(ok,'mf2005_GetUnsaturatedZoneFlux') ! Get unsaturated zone flux
XchModSimSimUnsZFlux = mv
if (.not.mf_steadystate) then
ok = metaswap_PutStorageFactor(XchModSimSimStrFct,mv); call DriverChk(ok,'metaswap_PutStorageFactor') ! Put storage factor
if (.not.llpf) then ! BCF
ok = mf2005_GetStorageFactor(mf_igrid,XchModSimSimStrFct,&
XchModSimModNID,XchSim2ModIdx,XchSim2ModOff,mv); call DriverChk(ok,'mf2005_GetStorageFactor') ! Get storage factor
else ! LPF
ok = mf2005_GetStorageFactorLPF(mf_igrid,XchModSimSimStrFct,&
XchModSimModNID,XchSim2ModIdx,XchSim2ModOff,mv); call DriverChk(ok,'mf2005_GetStorageFactorLPF') ! Get storage factor
end if
XchModSimSimStrFct = mv
end if
end if
!####### END EXCHANGE: AfterSolveMetaSwap #####################################
call mf2005_performIter(retValMF2005,psolved)
if (retValMF2005.ne.0) exitcode = -26
!##### BEGIN EXCHANGE: AfterSolve #############################################
if (exitcode.eq.0) then
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
ok = mf2005_PutHeads(mf_igrid,XchModSimModCells,XchModSimModNID,XchModSimModHeads,mv); call DriverChk(ok,'mf2005_PutHeads') ! MODFLOW put groundwater heads
ok = metaswap_GetHeads(XchModSimModHeads,XchModSimSimNID,XchMod2SimIdx,XchMod2SimOff,mv); call DriverChk(ok,'metaswap_GetHeads') ! MetaSWAP gets the groundwater heads
XchModSimModHeads = mv
end if
end if
!####### END EXCHANGE: AfterSolve #############################################
! finish iteration
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
call MetaSwap_finishIter(iterMetaSwap,convergedMetaSwap)
endif
call mf2005_finishIter(convergedMF2005,retValMF2005)
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
if (pks7mpimasterwrite()) write(*,*) ' convergence MODFLOW - MetaSWAP: ',convergedMF2005, convergedMetaSwap
end if
! get next iteration information
converged = .true.
converged = converged .and. convergedMF2005
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm)&
converged = converged .and. convergedMetaSwap
if (retValMF2005.ne.0 ) exitcode = -261
enddo
if (.not.converged) then
if (pks7mpimasterwrite()) write(*,*) ' Model did not converge!',exitcode
end if
! ... write results
if (converged .and. exitcode.eq.0) then
call osd_chdir(modwd2)
call mf2005_finishTimestep(retValMF2005,psolved,ncvgerr)
if (retValMF2005.ne.0) exitcode = -31
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
call osd_chdir(simwd2)
call MetaSwap_finishTimestep(currentTime)
end if
!#### TIMESERIES #####
call osd_chdir(root)
ok = mf2005_TimeserieGetHead(mf_igrid); call DriverChk(ok,'mf2005_TimeserieGetHead')
call gwf2getcurrentdate(mf_igrid,time_string) !,issflg(kkper)
call tserie1write(0,lss,currentTime,hnoflo,usests,modwd1,submstr(isub),time_string)
!#### TIMESERIES #####
endif
! ... next timestep
if (exitcode.ne.0) then
! ERROR, do not continue
endOfSimulation=.true.
else
! get new time information
! Let Modflow be leading. don't ask the other components
call mf2005_getBeginOfNextTimeStep(BeginOfNextTimeStepMF2005,retValMF2005)
! find errors
if (retValMF2005.eq.-1) then
! end of simulation for Modflow
endOfSimulation=.true.
retValMF2005 =0
endif
if (retValMF2005.ne.0) exitcode = -31
if (exitcode.ne.0) endOfSimulation=.true.
! find start time of next timestep
! this may be a time in the 'future', the 'past' or the 'current' time
! nodata value of time-value means: this was the last timestep for this component
if (exitcode.eq.0) then
currentTime=nodataTime
! check for each component (ok, only modflow available now)
if (BeginOfNextTimeStepMF2005.ne.nodataTime) then
if (currentTime.eq.nodataTime) then
currentTime=BeginOfNextTimeStepMF2005
else
currentTime=min(currentTime,BeginOfNextTimeStepMF2005)
endif
endif
if (currentTime.eq.nodataTime.and.rt.ne.rtmod) endOfSimulation=.true.
endif
endif
if (rt.eq.rtmodsimdflowfm) then
doTimeLoop = currenttime.lt.endOfCurrentTimeStepFM .and. .not.endOfSimulation
else
doTimeLoop = .not.endOfSimulation
end if
enddo ! time loop
! ... end
call osd_chdir(modwd2)
if (exitcode.eq.0) call mf2005_finishSimulation(retValMF2005)
if (retValMF2005.ne.0 ) exitcode = -61
if (rt.eq.rtmodsim .or. rt.eq.rtmodsimdflowfm) then
call osd_chdir(simwd2)
call MetaSWAP_finishSimulation()
end if
!#### TIMESERIES #####
call osd_chdir(root)
call tserie1write(1,lss,currentTime,hnoflo,usests,modwd1,submstr(isub),time_string)
call tserie1close()
!#### TIMESERIES #####
! list routine timer info
call cfn_rtt_list(0)
! Deallocate coupling arrays
call driverXchDeallocate()
! exit
if (exitcode.ne.0) then
if (pks7mpimasterwrite()) write(unit=*,fmt=*) 'ERROR, exit code: ',exitcode
call pks7mpifinalize()! PKS
call exit(10)
endif
!C10-----END OF PROGRAM.
IF(NCVGERR.GT.0) THEN
if (pks7mpimasterwrite()) WRITE(*,*) 'FAILED TO MEET SOLVER CONVERGENCE CRITERIA ',NCVGERR,' TIME(S)'
END IF
! next pest iteration
call imod_utl_closeunits()
if (.not.lipest) then
convergedPest=.true.
else
convergedPest=pestnext(lss,modwd1)
end if
! call imod_utl_closeunits()
end do pestloop
! PKS IDF merge
call osd_chdir(root)
call pks_imod_utl_write_idfmergefile(modwd2,idfmergefile)
if (lidfmerge.and.lpks) call pks_imod_utl_idfmerge(idfmergefile)
end do submodelloop
! timing output
call osd_chdir(root)
call timing_stat()
call pks7mpifinalize()! PKS
end program