program LumbricusUnitTests use mfmsdfm_mapping use LumbricusUnitTestsData implicit none ! variables ! 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 ! size for for exchange integer, parameter :: nModRiv = 3 ! #1D-rivers in modflow schematization integer, parameter :: nModRiv2 = 5 ! #not connected 1D-rivers in modflow schematization integer, parameter :: nModDrn = 4 ! #DRN+systeem3 items in modflow schematization integer, parameter :: nMswSprnk = 2 ! #sprinkling demands in MetaSWAP schematization integer, parameter :: nMswSvats = 2 ! #svats in MetaSWAP schematization integer, parameter :: nFlowNodes= 1 ! #1D-nodes (= 1D-gridpoints) in D-Flow FM schematization integer, parameter :: nFlowCells= 1 ! #2D-cells in D-Flow FM schematization ! id-arrays and vars for setting up exchange integer, dimension(nModRiv) :: rivLocIds = (/ 1, 2, 3 /) integer, dimension(nModRiv2) :: riv2LocIds = (/ 1, 2, 3, 4, 5 /) integer, dimension(nModDrn) :: drnLocIds = (/ 1, 2, 3, 4 /) integer, dimension(nMswSprnk) :: sprlocIds = (/ 1, 2 /) integer, dimension(nMswSvats) :: svtlocIds = (/ 1, 2 /) integer, dimension(nFlowNodes) :: fmNodeIds = (/ 1 /) integer, dimension(nFlowCells) :: fmCellIds = (/ 1 /) logical :: success, ok integer :: retVal ! time stepping and loop counters integer :: timestep ! time step counter double precision :: dtMf ! iMODFLOW time step size integer :: numDtswTimeSteps ! #MetaSWAP steps per day integer :: idtsw ! MetaSWAP top system time step counter double precision :: dtsw ! MetaSWAP top system time step size ! ! value-arrays for the (mapped) exchanged values ! double precision, parameter :: mv = -99999.0D+0 ! ! ! 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, DRN package: double precision, dimension(:), allocatable :: XchFromMfFluxesInDrains ! iMODFLOW -> D-Flow FM 1D, drainage from non-coupled rivers: double precision, dimension(:), allocatable :: XchFromMfFluxesInRivers2 ! MetaSwap -> D-Flow FM 1D, runoff from svats: double precision, dimension(:), allocatable :: XchFromMswRunoffOnSvats ! ! END variables for coupling with dflow-fm ! call dmmInitialize() ! Read mapper configs for FM <-> Msw exchange (2D FM nodes <-> MetaSWAP svats) mapperWatLev2D = dmmCreateMapping('Dfm2DWatLevToMsw_H.dmm' , dmm_FM2MSW_WL , 'DFLOWFM', 'METASWAP') mapperDeltaPonding = dmmCreateMapping('MswPondingToDfm2D_DV.dmm' , dmm_MSW2FM_DPV, 'METASWAP', 'DFLOWFM') ! Read mapper configs for 1D exchange (waterlevel from FM, combined fluxes to FM) mapperWatLev1D = dmmCreateMapping('Dfm1DWatLevToMfRiv_H.dmm', dmm_FM2MF_WL , 'DFLOWFM' , 'MODFLOW') mapperMfRivToDfm1D = dmmCreateMapping('MfRivToDfm1D_Q.dmm' , dmm_MF2FM_Q , 'MODFLOW' , 'DFLOWFM') mapperMfRiv2ToDfm1D = dmmCreateMapping('MfRiv2ToDfm1D_Q.dmm' , dmm_MF2FM_RIV2, 'MODFLOW' , 'DFLOWFM') mapperSprinkling = dmmCreateMapping('MswSprinkToDfm1D_Q.dmm' , dmm_MSW2FM_SPR, 'METASWAP', 'DFLOWFM') mapperMfDrnToDfm1D = dmmCreateMapping('MfDrnToDfm1D_Q.dmm' , dmm_MF2FM_DRN , 'MODFLOW' , 'DFLOWFM') mapperMswRunoffToDfm1D = dmmCreateMapping('MswRunoffToDfm1D_Q.dmm' , dmm_MSW2FM_RUN, 'METASWAP', 'DFLOWFM') if (mapperWatLev2D < 0 .or. & mapperDeltaPonding < 0 .or. & mapperWatLev1D < 0 .or. & mapperMfRivToDfm1D < 0 .or. & mapperSprinkling < 0 .or. & mapperMfDrnToDfm1D < 0 .or. & mapperMswRunoffToDfm1D < 0 ) then stop 'Error reading on or more of the mapping config files' endif ! Feed mappers with source and target ids: ! Water level FM 2D-cells -> MetaSWAP svats retVal = dmmSetSourceIds(mapperWatLev2D, fmCellIds) retVal = dmmSetTargetIds(mapperWatLev2D, svtlocIds) ! Delta Ponding Volume MetaSWAP svats -> FM 2D-cells retVal = dmmSetSourceIds(mapperDeltaPonding, svtlocIds) retVal = dmmSetTargetIds(mapperDeltaPonding, fmCellIds) ! Water level D-Flow FM 1D-nodes -> Modflow river elements retVal = dmmSetSourceIds(mapperWatLev1D, fmNodeIds) retVal = dmmSetTargetIds(mapperWatLev1D, rivLocIds) ! Infiltration/seepage Modflow river elements -> D-Flow FM 1D-nodes retVal = dmmSetSourceIds(mapperMfRivToDfm1D, rivLocIds) retVal = dmmSetTargetIds(mapperMfRivToDfm1D, fmNodeIds) if (mapperSprinkling > 0) then ! Sprinkling demand MetaSWAP svats -> D-Flow FM 1D-nodes retVal = dmmSetSourceIds(mapperSprinkling, sprlocIds) retVal = dmmSetTargetIds(mapperSprinkling, fmNodeIds) endif if (mapperMfRiv2ToDfm1D > 0) then ! Drainage from Modflow non-coupled rivers retVal = dmmSetSourceIds(mapperMfRiv2ToDfm1D, riv2LocIds) retVal = dmmSetTargetIds(mapperMfRiv2ToDfm1D, fmNodeIds) endif if (mapperMfDrnToDfm1D > 0) then ! Drainage Modflow DRN/system3 -> D-Flow FM 1D-nodes retVal = dmmSetSourceIds(mapperMfDrnToDfm1D, drnLocIds) retVal = dmmSetTargetIds(mapperMfDrnToDfm1D, fmNodeIds) endif if (mapperMswRunoffToDfm1D > 0) then ! Drainage MetaSWAP svats -> D-Flow FM 1D-nodes retVal = dmmSetSourceIds(mapperMswRunoffToDfm1D, svtlocIds) retVal = dmmSetTargetIds(mapperMswRunoffToDfm1D, fmNodeIds) 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: allocate(XchFromDfmWaterlevels1D(nFlowNodes)) allocate(XchToModflowWaterlevelsOnRivers(nModRiv)) ! iMODFLOW/MetaSWAP -> D-Flow FM, fluxes and realized fluxes in 1D nodes: allocate(XchToDfmFluxes1D(nFlowNodes)) allocate(XchFromDfmRealizedFluxes1D(nFlowNodes)) ! iMODFLOW -> D-Flow FM, infiltration/seepage: allocate(XchFromMfFluxesInRivers(nModRiv)) allocate(XchToMfCorrectionFluxesInRivers(nModRiv)) ! iMODFLOW -> D-Flow FM, drains + system3: if (mapperMfDrnToDfm1D > 0) then allocate(XchFromMfFluxesInDrains(nModDrn)) endif ! iMODFLOW -> D-Flow FM, small non-coupled rivers: if (mapperMfRiv2ToDfm1D > 0) then allocate(XchFromMfFluxesInRivers2(nModRiv2)) endif ! MetaSWAP -> D-Flow FM, sprinkling: if (mapperSprinkling > 0) then allocate(XchFromMswSprinklingDemand(nMswSvats)) allocate(XchToMswRealizedSprinklingDemand(nMswSvats)) endif ! MetaSWAP -> D-Flow FM, runoff on svats: if (mapperMswRunoffToDfm1D > 0) then allocate(XchFromMswRunoffOnSvats(nMswSvats)) endif ! combine the iMod fluxes (RIV, DRN) and MetaSWAP fluxes (sprinkling, runoff) retVal = dmmCombineMappers(mapperMfRivToDfm1D , mapperSprinkling , & mapperMfDrnToDfm1D , mapperMswRunoffToDfm1D, & mapperMfRiv2ToDfm1D) ! ... timestep timestep = 0 if (.not. dmmDoLogTimeStepCheck(timestep)) then stop endif dtMf = 1.0d+0 ! 1 julian day ! initalize MetaSWAP surface water ok = DflowFM_Put2DWaterlevels(XchFromDfmWaterlevelsOnCells) call dmmMapValues(mapperWatLev2D, XchFromDfmWaterlevelsOnCells, XchToMswWaterlevelsOnSvats) ok = MetaSWAP_Get2DWaterlevels(XchToMswWaterlevelsOnSvats,mv) !!D!! call MetaSWAP_sw_initcomponent() !!!! Remark: getDTSW has been moved down for testing purposes !!!! (needed to be ableto vary the dtsw per modflow day) !!D!! ! get the MetaSWAP dtsw !!D!! ok = MetaSWAP_PutSurfacewaterTimestep(dtsw) !!D!! ! compute the river fluxes for intial heads !!D!! call mf2005_RivBudgetShead() timestep = 0 do while (timestep < dmmTestGetNumTimeSteps()) ! Main time loop (= test number loop) loop (steps of 1 day) timestep = timestep + 1 call dmmTestSetModflowStep(timestep) if (.not. dmmDoLogTimeStepCheck(timestep)) then stop endif !!!! Remark: getDTSW has been moved to here for testing purposes !!!! (needed to be ableto vary the dtsw per modflow day) ok = MetaSWAP_PutSurfacewaterTimestep(dtsw) call dmmSetDtsw(dtsw) !!D!! iMODFLOW/MetaSWAP exchange at BeginOfTimestep ! FM's 1D water levels at start of MF-step; units: m AD NAP ok = dflowfm_Put1DNodeSurfacewaterLevels(XchFromDfmWaterlevels1D); call dmmMapValues(mapperWatLev1D, XchFromDfmWaterlevels1D, XchToModflowWaterlevelsOnRivers) ok = mf2005_Get1DWaterlevels(XchToModflowWaterlevelsOnRivers,mv); !!D!! ! ... init timestep for MODFLOW !!D!! ! call mf2005_prepareTimestep(currentTime,.false.,retValMF2005) !!D!! ! ! compute MODFLOW river fluxes based on !BM: calculating V1D using FM-levels !!D!! ! ! current groundwater heads and current FM-levels !BM: calculating V1D using FM-levels !!D!! ! if (rt.eq.rtmodsimdflowfm) then !BM: calculating V1D using FM-levels !!D!! ! call mf2005_RivBudgetShead() !BM: calculating V1D using FM-levels !!D!! ! endif !BM: calculating V1D using FM-levels ! provide MODFLOW fluxes to the combined 1D fluxes mapper ok = mf2005_Put1DVolumes(XchFromMfFluxesInRivers); call dmmCombMapperSetVolumes(mapperMfRivToDfm1D, XchFromMfFluxesInRivers) if (mapperMfRiv2ToDfm1D > 0) then ok = mf2005_Put1DFlux2(XchFromMfFluxesInRivers2) call dmmCombMapperSetVolumes(mapperMfRiv2ToDfm1D, XchFromMfFluxesInRivers2) endif if (mapperMfDrnToDfm1D > 0) then ok = mf2005_PutDrnFlux(XchFromMfFluxesInDrains) call dmmCombMapperSetVolumes(mapperMfDrnToDfm1d, XchFromMfFluxesInDrains) endif !!D!! call mf2005_initTimeStep(currentTime,retValMF2005) !!D!! call mf2005_getEndOfCurrentTimeStep(endOfCurrentTimeStepMF2005,retValMF2005) ! compute the number of MetaSWAP surface water iterations numDtswTimeSteps = nint(dtMf/dtsw) ! let MetaSWAP compute sprinkling demand !!D!! MetaSWAP_initTimestep(dtMetaSWAP) ! MetaSWAP-TOP + DflowFM loop do idtsw = 1, numDtswTimeSteps call dmmTestSetDtswStep(idtsw) !!D!! call MetaSWAP_performSurfacewaterTimestep(idtsw) ok = MetaSWAP_PutPondingDeltaVolumes(XchFromMswDeltaVolumesOnSvats) call dmmMapValues(mapperDeltaPonding, XchFromMswDeltaVolumesOnSvats, XchToDfmFluxes2D) ok = DflowFM_GetFluxes2D(XchToDfmFluxes2D, mv) if (mapperSprinkling > 0) then ! Add MetaSWAP sprinkling demands to the MODFLOW fluxes (in the combined 1D fluxes mapper) ok = MetaSWAP_PutSprinklingDemand(XchFromMswSprinklingDemand) call dmmCombMapperSetVolumes(mapperSprinkling, XchFromMswSprinklingDemand) endif if ( mapperMswRunoffToDfm1D > 0) then ! Add MetaSWAP SVAT runnoff to the MODFLOW fluxes (in the combined 1D fluxes mapper) ok = MetaSWAP_PutRunoffVolumes(XchFromMswRunoffOnSvats) call dmmCombMapperSetVolumes(mapperMswRunoffToDfm1D, XchFromMswRunoffOnSvats) endif ! Provide combined 1D fluxes to FM call dmmCombMapperGetFluxes(mapperMfRivToDfm1D, XchToDfmFluxes1D) ok = dflowfm_Get1DNodeFlux(XchToDfmFluxes1D, mv); ! Let FM compute over a dtsw time step !!D!! call dflowfm_performSurfacewaterTimestep(tStartTop, tEndTop) ! Provide realized FM 1D fluxes to combined mapper ok = dflowfm_PutRealised1DNodeFlux(XchFromDfmRealizedFluxes1D) call dmmCombMapperSetRealizedFluxes(mapperMfRivToDfm1D, XchFromDfmRealizedFluxes1D) ! Provide realized FM 2D fluxes to MetaSWAP ok = dflowfm_PutRealised2DCellFlux(XchFromDfmRealizedFluxes2D) call dmmMapValuesBack(mapperDeltaPonding, XchFromDfmRealizedFluxes2D, XchToMswRealizedDeltaVolumesOnSvats) ok = MetaSWAP_GetDeltaPondingAllocation(XchToMswRealizedDeltaVolumesOnSvats, mv) ! Provide realized sprinkling demand to MetaSWAP if (mapperSprinkling > 0) then call dmmCombMapperGetRealizedFluxes(mapperSprinkling, XchToMswRealizedSprinklingDemand) ok = MetaSWAP_GetSprinklingAllocation(XchToMswRealizedSprinklingDemand, mv) endif ! Put the water levels for the next step to MetaSWAP time step ! (MetaSWAP needs them for finishing its current time step) ok = DflowFM_Put2DWaterlevels(XchFromDfmWaterlevelsOnCells) call dmmMapValues(mapperWatLev2D, XchFromDfmWaterlevelsOnCells, XchToMswWaterlevelsOnSvats) ok = MetaSWAP_Get2DWaterlevels(XchToMswWaterlevelsOnSvats,mv) !!D!! call MetaSWAP_finishSurfacewaterTimestep(idtsw) enddo call dmmCombMapperGetRealizedFluxes(mapperMfRivToDfm1D, XchToMfCorrectionFluxesInRivers) ok = mf2005_Get1DCorrectionFlux(XchToMfCorrectionFluxesInRivers, mv) !!D!! MetaSWAP / iMODFLOW iteration loop enddo ! Main time loop end program LumbricusUnitTests