!! 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
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
! debug
integer ncvgerr
! VARIABLES FOR COUPLING WITH DFLOW-FM
!
! needed for reading config files
!
character(len=32 ) :: exchType ! soure/target components, exchange type
character(len=256) :: mapConfig ! mapper config files
logical :: compHasErrors ! component file or one or more of the sub files have errors
!
! handles to mappers
!
integer :: mapperWatLev2D = 0
integer :: mapperDeltaPonding = 0
integer :: mapperWatLev1D = 0
integer :: mapperMfRivToDfm1D = 0
integer :: mapperMfRiv2ToDfm1D = 0
integer :: mapperSprinkling = 0
integer :: mapperMfDrnToDfm1D = 0
integer :: mapperMswRunoffToDfm1D = 0
!
! number of elements per exchanged variable
!
integer :: nModRiv ! #1D-rivers in modflow schematization
integer :: nModRiv2 ! #not connected 1D-rivers in modflow schematization
integer :: nModDrn ! #DRN+systeem3 items in modflow schematization
integer :: nMswSvats ! #svats in MetaSWAP schematization
integer :: nFlowNodes ! #1D-nodes (= 1D-gridpoints) in D-Flow FM schematization
integer :: nFlowCells ! #2D-cells in D-Flow FM schematization
!
! id-arrays for setting up mappers
!
integer, dimension(:), allocatable :: locationIds
logical :: success
integer :: retVal
integer :: riv
!
! value-arrays for the (mapped) exchanged values
!
! D-Flow FM 2D -> MetaSwap, waterlevels:
double precision, dimension(:), allocatable :: XchFromDfmWaterlevelsOnCells
double precision, dimension(:), allocatable :: XchToMswWaterlevelsOnSvats
! MetaSwap -> D-Flow FM 2D, ponding volumes:
double precision, dimension(:), allocatable :: XchFromMswDeltaVolumesOnSvats
double precision, dimension(:), allocatable :: XchToDfmFluxes2D
double precision, dimension(:), allocatable :: XchFromDfmRealizedFluxes2D
double precision, dimension(:), allocatable :: XchToMswRealizedDeltaVolumesOnSvats
! D-Flow FM 1D -> iMODFLOW, waterlevels in 1D nodes:
double precision, dimension(:), allocatable :: XchFromDfmWaterlevels1D
double precision, dimension(:), allocatable :: XchToModflowWaterlevelsOnRivers
! iMODFLOW/MetaSwap -> D-Flow FM 1D, summed fluxes and realized summed fluxes in 1D nodes:
double precision, dimension(:), allocatable :: XchToDfmFluxes1D
double precision, dimension(:), allocatable :: XchFromDfmRealizedFluxes1D
! iMODFLOW -> D-Flow FM 1D, infiltration/seepage:
double precision, dimension(:), allocatable :: XchFromMfFluxesInRivers
double precision, dimension(:), allocatable :: XchToMfCorrectionFluxesInRivers
! MetaSwap -> D-Flow FM 1D, sprinkling:
double precision, dimension(:), allocatable :: XchFromMswSprinklingDemand
double precision, dimension(:), allocatable :: XchToMswRealizedSprinklingDemand
! iMODFLOW -> D-Flow FM 1D, drainage from non-coupled rivers:
double precision, dimension(:), allocatable :: XchFromMfFluxesInRivers2
! iMODFLOW -> D-Flow FM 1D, drainage from DRN-package:
double precision, dimension(:), allocatable :: XchFromMfFluxesInDrains
! MetaSwap -> D-Flow FM 1D, runoff from svats:
double precision, dimension(:), allocatable :: XchFromMswRunoffOnSvats
! END variables for coupling with dflow-fm
! 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.
compHasErrors = .false.
! get explicit arguments: check for components file
call cfn_vcl_fndc(ivcl,iarg,'-components',.true.,compfile,1)
if (iarg.gt.0) then
! open components 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 getcwd(modwd1)
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)
else
call getcwd(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
if (retVal /= 0) stop 'Could not load dflowfm library. Check dependencies'
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,'-exchType',.true.,exchType,1)
call cfn_vcl_fndc(jvcl,jarg,'-config',.true.,mapConfig,1)
mapConfig = trim(root)//'\'//trim(wd)//'\'//trim(mapConfig)
! create mappers, one for each exchang type
if (exchType=='waterlevel2D') then
mapperWatLev2D = dmmCreateMapping(mapConfig, dmm_FM2MSW_WL, 'DFLOWFM', 'METASWAP')
if (mapperWatLev2D < 0) then
print *, 'ERRORs in waterlevel2D mapping file '//trim(mapConfig)
compHasErrors = .true.
endif
else if (exchType=='deltaPonding') then
mapperDeltaPonding = dmmCreateMapping(mapConfig, dmm_MSW2FM_DPV, 'METASWAP', 'DFLOWFM')
if (mapperDeltaPonding < 0) then
print *, 'ERRORs in deltaPonding mapping file '//trim(mapConfig)
compHasErrors = .true.
endif
else if (exchType=='waterlevel1D') then
mapperWatLev1D = dmmCreateMapping(mapConfig, dmm_FM2MF_WL, 'DFLOWFM', 'MODFLOW')
if (mapperWatLev1D < 0) then
print *, 'ERRORs in waterlevel1D mapping file '//trim(mapConfig)
compHasErrors = .true.
endif
else if (exchType=='riverFlux') then
mapperMfRivToDfm1D = dmmCreateMapping(mapConfig, dmm_MF2FM_Q, 'MODFLOW', 'DFLOWFM')
if (mapperMfRivToDfm1D < 0) then
print *, 'ERRORs in riverFlux mapping file '//trim(mapConfig)
compHasErrors = .true.
endif
else if (exchType=='sprinkling') then
mapperSprinkling = dmmCreateMapping(mapConfig, dmm_MSW2FM_SPR, 'METASWAP', 'DFLOWFM')
if (mapperSprinkling < 0) then
print *, 'ERRORs in sprinkling mapping file '//trim(mapConfig)
compHasErrors = .true.
endif
else if (exchType=='mfDrains') then
mapperMfDrnToDfm1D = dmmCreateMapping(mapConfig, dmm_MF2FM_DRN, 'MODFLOW', 'DFLOWFM')
if (mapperMfDrnToDfm1D < 0) then
print *, 'ERRORs in mfDrains mapping file '//trim(mapConfig)
compHasErrors = .true.
endif
else if (exchType=='riverFlux2') then
mapperMfRiv2ToDfm1D = dmmCreateMapping(mapConfig, dmm_MF2FM_RIV2, 'MODFLOW', 'DFLOWFM')
if (mapperMfRiv2ToDfm1D < 0) then
print *, 'ERRORs in mfDrains mapping file '//trim(mapConfig)
compHasErrors = .true.
endif
else if (exchType=='mswRunnoff') then
mapperMswRunoffToDfm1D = dmmCreateMapping(mapConfig, dmm_MSW2FM_RUN, 'METASWAP', 'DFLOWFM')
if (mapperMswRunoffToDfm1D < 0) then
print *, 'ERRORs in mswRunnoff mapping file '//trim(mapConfig)
compHasErrors = .true.
endif
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
if (mapperWatLev2D == 0) then
print *, 'Error: no waterlevel2D mapper found in config file '//trim(compfile)
compHasErrors = .true.
endif
if (mapperDeltaPonding == 0) then
print *, 'Error: no deltaPonding mapper found in config file '//trim(compfile)
compHasErrors = .true.
endif
if (compHasErrors) stop
! 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 per location type
nFlowNodes = dflowfm_PutNumberOf1DNodes()
nFlowCells = dflowfm_PutNumberOf2DGridCells()
success = metaswap_PutNumberOfSvats(nMswSvats); call DriverChk(success,'mf2005_PutNumberOfRivers')
ok = mf2005_PutNumberOfRivers(nModRiv); call DriverChk(ok,'mf2005_PutNumberOfRivers')
!## remember nModRiv2=nModRiv omdat er gewerkt wordt met een inverse-index matrix
ok = mf2005_PutNumberOfRivers2(nModRiv2); call DriverChk(ok,'mf2005_PutNumberOfRivers2')
!## get maximum number of drain-elements
ok = mf2005_PutNumberOfDrains(nModDrn); call DriverChk(ok,'mf2005_PutNumberOfDrains')
! get the ids of DFM 2D-cells and feed them to the mappers
allocate(locationIds(nFlowCells))
success = dflowfm_Put2DGridCellIDs(locationIds)
retVal = dmmSetSourceIds(mapperWatLev2D, locationIds)
retVal = dmmSetTargetIds(mapperDeltaPonding, locationIds)
if (mapperMswRunoffToDfm1D > 0) then
retVal = dmmSetTargetIds(mapperMswRunoffToDfm1D, locationIds)
endif
deallocate(locationIds)
! get the ids of the metaswap SVATS and feed them to the mappers
allocate(locationIds(nMswSvats))
success = metaswap_PutSvatIDs(locationIds)
retVal = dmmSetSourceIds(mapperDeltaPonding, locationIds)
retVal = dmmSetTargetIds(mapperWatLev2D, locationIds)
if (mapperSprinkling > 0) then
! TODO: is dit een andere set SVATS?
retVal = dmmSetSourceIds(mapperSprinkling, locationIds)
endif
if (mapperMswRunoffToDfm1D > 0) then
retVal = dmmSetSourceIds(mapperMswRunoffToDfm1D, locationIds)
endif
deallocate(locationIds)
! get the ids of the DFM 1D nodes and feed them to the mappers
allocate(locationIds(nFlowNodes))
success = dflowfm_Put1DNodeIDs(locationIds)
if (mapperMfRivToDfm1D > 0) then
retVal = dmmSetTargetIds(mapperMfRivToDfm1D, locationIds)
endif
if (mapperWatLev1D > 0) then
retVal = dmmSetSourceIds(mapperWatLev1D, locationIds)
endif
! iMODFLOW -> D-Flow FM, small non-coupled rivers:
if (mapperMfRiv2ToDfm1D > 0) then
retVal = dmmSetTargetIds(mapperMfRiv2ToDfm1D, locationIds)
endif
! iMODFLOW -> D-Flow FM, DRN package
if (mapperMfDrnToDfm1D > 0) then
retVal = dmmSetTargetIds(mapperMfDrnToDfm1D, locationIds)
endif
! MetaSwap -> D-Flow FM, sprinkling
if (mapperSprinkling > 0) then
retVal = dmmSetTargetIds(mapperSprinkling, locationIds)
endif
! MetaSwap -> D-Flow FM, 2D runnoff
if (mapperMswRunoffToDfm1D > 0) then
retVal = dmmSetTargetIds(mapperMswRunoffToDfm1D, locationIds)
endif
deallocate(locationIds)
! get (=generate) the ids of the Modflow river elements, feed them to the mappers
if (mapperMfRivToDfm1D > 0) then
allocate(locationIds(nModRiv))
do riv= 1, nmodriv
locationIds(riv) = riv
enddo
retVal = dmmSetSourceIds(mapperMfRivToDfm1d, locationIds)
retVal = dmmSetTargetIds(mapperWatLev1D, locationIds)
deallocate(locationIds)
endif
if (mapperMfRiv2ToDfm1d > 0) then
! get (=generate) the ids of the Modflow river 2 elements (small non coupled rivers)
allocate(locationIds(nModriv2))
do riv= 1, nmodriv2
locationIds(riv) = riv
enddo
retVal = dmmSetSourceIds(mapperMfRiv2ToDfm1d, locationIds)
deallocate(locationIds)
endif
if (mapperMfDrnToDfm1d > 0) then
! get (=generate) the ids of the Modflow drain elements
allocate(locationIds(nModDrn))
do riv= 1, nModDrn
locationIds(riv) = riv
enddo
retVal = dmmSetSourceIds(mapperMfDrnToDfm1d, locationIds)
deallocate(locationIds)
endif
! allocatie exchange variables for the D-Flow FM / Modflow coupling
! D-Flow FM -> MetaSWAP, waterlevels:
allocate(XchFromDfmWaterlevelsOnCells(nFlowCells))
allocate(XchToMswWaterlevelsOnSvats(nMswSvats))
! MetaSWAP -> D-Flow FM, demanded ponding volumes:
allocate(XchFromMswDeltaVolumesOnSvats(nMswSvats))
allocate(XchToDfmFluxes2D(nFlowCells))
! MetaSWAP -> D-Flow FM, realized ponding volumes:
allocate(XchFromDfmRealizedFluxes2D(nFlowCells))
allocate(XchToMswRealizedDeltaVolumesOnSvats(nMswSvats))
! D-Flow FM -> iMODFLOW, waterlevels in 1D nodes:
if (mapperWatLev1D > 0) then
allocate(XchFromDfmWaterlevels1D(nFlowNodes))
allocate(XchToModflowWaterlevelsOnRivers(nModRiv))
endif
! iMODFLOW/MetaSWAP -> D-Flow FM, fluxes and realized fluxes in 1D nodes:
! iMODFLOW -> D-Flow FM, infiltration/seepage:
if (mapperMfRivToDfm1D > 0) then
allocate(XchToDfmFluxes1D(nFlowNodes))
allocate(XchFromDfmRealizedFluxes1D(nFlowNodes))
allocate(XchFromMfFluxesInRivers(nModRiv))
allocate(XchToMfCorrectionFluxesInRivers(nModRiv))
endif
! iMODFLOW -> D-Flow FM, DRN-package:
! iMODFLOW -> D-Flow FM, small non-coupled rivers:
if (mapperMfRiv2ToDfm1D > 0) then
allocate(XchFromMfFluxesInRivers2(nModRiv2))
endif
if (mapperMfDrnToDfm1D > 0) then
allocate(XchFromMfFluxesInDrains(nModDrn))
endif
! MetaSWAP -> D-Flow FM, sprinkling:
allocate(XchFromMswSprinklingDemand(nMswSvats))
allocate(XchToMswRealizedSprinklingDemand(nMswSvats))
! MetaSWAP -> D-Flow FM, runoff on svats:
if (mapperMswRunoffToDfm1D > 0) then
allocate(XchFromMswRunoffOnSvats(nMswSvats))
endif
if(mapperMfRivToDfm1D > 0) then
! combine the iMod fluxes (RIV, DRN) and MetaSWAP fluxes (sprinkling, runoff)
retVal = dmmCombineMappers(mapperMfRivToDfm1D , mapperSprinkling , &
mapperMfDrnToDfm1D , mapperMswRunoffToDfm1D, &
mapperMfRiv2ToDfm1D)
endif
end if
! ###### END EXCHANGE: BeforeTimestep #########################################
! ... timestep
timestep = 0
if (.not. dmmDoLogTimeStepCheck(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(mapperWatLev2D, XchFromDfmWaterlevelsOnCells, XchToMswWaterlevelsOnSvats)
ok = metaswap_Get2DWaterlevels(XchToMswWaterlevelsOnSvats,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')
call dmmSetDtsw(dtsw)
endif
!BM ! compute the river fluxes for intial heads
!BM if (rt.eq.rtmodsimdflowfm) then
!BM call mf2005_RivBudgetShead()
!BM endif
do while (.not.endOfSimulation .and. exitcode.eq.0)
tsc=tsc+1
timestep = timestep + 1
if (.not. dmmDoLogTimeStepCheck(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 .and. mapperWatLev1D > 0) then
! provide FM's water levels to MODFLOW
ok = dflowfm_Put1DNodeSurfacewaterLevels(XchFromDfmWaterlevels1D); call DriverChk(ok,'dflowfm_Put1DNodeSurfacewaterLevels') ! 1D water levels at start of MF-step; units: m AD NAP
call dmmMapValues(mapperWatLev1D, XchFromDfmWaterlevels1D, 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)
call osd_chdir(modwd2)
call mf2005_initTimeStep(currentTime,retValMF2005) !BM: mf2005_initTimeStep() contains all MODFLOW-GWF2xxxRP-calls (in mf2005_comp.f), the last one is gwf2dxc1rp(igrid)
if (retValMF2005.ne.0) exitcode = -15
! 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 .and. mapperWatLev1D > 0) then
! provide MODFLOW fluxes to FM
ok = mf2005_Put1DFlux(XchFromMfFluxesInRivers); call DriverChk(ok,'mf2005_Put1DFlux')
call dmmCombMapperSetVolumes(mapperMfRivToDfm1d, XchFromMfFluxesInRivers)
if (mapperMfRiv2ToDfm1D > 0) then
ok = mf2005_Put1DFlux2(XchFromMfFluxesInRivers2); call DriverChk(ok,'mf2005_Put1DFlux2')
call dmmCombMapperSetVolumes(mapperMfRiv2ToDfm1D, XchFromMfFluxesInRivers2)
endif
if (mapperMfDrnToDfm1D > 0) then
ok = mf2005_PutDrnFlux(XchFromMfFluxesInDrains); call DriverChk(ok,'mf2005_PutDrnFlux')
call dmmCombMapperSetVolumes(mapperMfDrnToDfm1D, XchFromMfFluxesInDrains)
endif
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)
if (mapperSprinkling > 0) then
ok = MetaSWAP_PutSprinklingDemand(XchFromMswSprinklingDemand); call DriverChk(ok,'')
call dmmCombMapperSetVolumes(mapperSprinkling, XchFromMswSprinklingDemand)
endif
if ( mapperMswRunoffToDfm1D > 0) then
! MSW's Svat runoffs directly go to FM-1D gridpoints. No 2D cells in FM model
! Add MetaSWAP SVAT runnoff to the MODFLOW fluxes (in the combined 1D fluxes mapper)
ok = metaswap_PutPondingDeltaVolumes(XchFromMswRunoffOnSvats)
call dmmCombMapperSetVolumes(mapperMswRunoffToDfm1D, XchFromMswRunoffOnSvats)
else
! MSW's Svat runoffs go to 2D cells in FM model
ok = metaswap_PutPondingDeltaVolumes(XchFromMswDeltaVolumesOnSvats); call DriverChk(ok,'metaswap_PutPondingDeltaVolumes')
call dmmMapValues(mapperDeltaPonding, XchFromMswDeltaVolumesOnSvats, XchToDfmFluxes2D)
ok = DflowFM_GetFluxes2D(XchToDfmFluxes2D, mv); call DriverChk(ok,'DflowFM_GetFluxes2D')
endif
if (mapperMfRivToDfm1D > 0) then
call dmmCombMapperGetFluxes(mapperMfRivToDfm1D, XchToDfmFluxes1D)
ok = dflowfm_Get1DNodeFlux(XchToDfmFluxes1D, mv); call DriverChk(ok,'dflowfm_Get1DNodeFlux')
endif
print *, 'FM TS START'
call dflowfm_performSurfacewaterTimestep(tStartTop, tEndTop, idtsw)
endOfCurrentTimeStepFM = 0.0D+0 ! TODO: conversion
print *, 'FM TS DONE'
if (mapperMfRivToDfm1D > 0) then
! provide realized 1D node fluxes to combined mapper
ok = dflowfm_PutRealised1DNodeFlux(XchFromDfmRealizedFluxes1D); call DriverChk(ok,'dflowfm_PutRealised1DNodeFlux')
call dmmCombMapperSetRealizedFluxes(mapperMfRivToDfm1D, XchFromDfmRealizedFluxes1D)
endif
if ( mapperMswRunoffToDfm1D > 0) then
! MSW's Svat runoffs directly go to FM-1D gridpoints. No 2D cells in FM model
! Add MetaSWAP SVAT runnoff to the MODFLOW fluxes (in the combined 1D fluxes mapper)
XchToMswRealizedDeltaVolumesOnSvats = mv
ok = MetaSWAP_GetDeltaPondingAllocation(XchToMswRealizedDeltaVolumesOnSvats, mv); call DriverChk(ok,'MetaSWAP_GetDeltaPondingAllocation')
else
! MSW's Svat runoffs go to 2D cells in FM model
ok = dflowfm_PutRealised2DCellFlux(XchFromDfmRealizedFluxes2D); call DriverChk(ok,'dflowfm_PutRealised2DCellFlux')
call dmmMapValuesBack(mapperDeltaPonding, XchFromDfmRealizedFluxes2D, XchToMswRealizedDeltaVolumesOnSvats)
ok = MetaSWAP_GetDeltaPondingAllocation(XchToMswRealizedDeltaVolumesOnSvats, mv); call DriverChk(ok,'MetaSWAP_GetDeltaPondingAllocation')
endif
if (mapperSprinkling > 0) then
call dmmCombMapperGetRealizedFluxes(mapperSprinkling, XchToMswRealizedSprinklingDemand)
ok = MetaSWAP_GetSprinklingAllocation(XchToMswRealizedSprinklingDemand, mv); call DriverChk(ok,'MetaSWAP_GetSprinklingAllocation')
else
XchToMswRealizedSprinklingDemand = 0.0D+0
ok = MetaSWAP_GetSprinklingAllocation(XchToMswRealizedSprinklingDemand, mv); call DriverChk(ok,'MetaSWAP_GetSprinklingAllocation')
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(mapperWatLev2D, XchFromDfmWaterlevelsOnCells, XchToMswWaterlevelsOnSvats)
ok = metaswap_Get2DWaterlevels(XchToMswWaterlevelsOnSvats,mv); call DriverChk(ok,'metaswap_Get2DWaterlevels') ! m AD NAP
call MetaSWAP_finishSurfacewaterTimestep(idtsw)
enddo
if (rt.eq.rtmodsimdflowfm .and. mapperMfRivToDfm1D > 0) then
call dmmCombMapperGetRealizedFluxes(mapperMfRivToDfm1D, XchToMfCorrectionFluxesInRivers)
ok = mf2005_Get1DCorrectionFlux(XchToMfCorrectionFluxesInRivers, mv); call DriverChk(ok,'mf2005_Get1DCorrectionFlux')
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