!! 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. !> Module for the coupling MODFLOW-MetaSWAP with MOZART. module dflowfmmod use iso_c_binding use kdtree2_module implicit none real(kdkind) :: search_radius type(kdtree2), pointer :: coordinate1D_tree ! kdtree for node (1D) id lookup by coordinate type(kdtree2), pointer :: coordinate2D_tree ! kdtree for grid cell (2D) id lookup by coordinate double precision, allocatable :: & ! cumulative volumes from FM. From this and the values cumulative_volumes_in(:) ! at the previous time step the realized flux is computed double precision, allocatable :: & cumulative_volumes_in_prev_MF_timestep(:) double precision, allocatable :: & cumulative_volumes_in_prev_MSW_timestep(:) double precision :: lastTimestep ! TODO: differentiate between MF and MSW private :: getCumulativeVolumes ! private method for getting the volumes from FM ! integer(c_size_t) :: dllHandle = -1 ! handle to dflowfm.dll or dimr_dll.dll character(len=16) :: modelPrefix = '' ! model prefix in BMI call ! (empty in case of FM only, 'dflowfm/' in case of dimr-run) ! todo: make configurable logical :: doDebugLogging = .false. ! todo: make configurable integer :: logTimestep = 1 ! counter for debu-logging contains !----------------------------- General ------------------------------ function dflowfm_LoadLibrary(dllName) result(retVal) implicit none integer :: retVal ! arguments character(len=*) :: dllName ! externals integer, external :: open_shared_library retVal = open_shared_library(dllHandle, dllName) if (dllName == 'dimr') then modelPrefix = 'FlowFM/' endif end function dflowfm_LoadLibrary !> Initialize a model, based on main input file (*.mdu). !! After this call, model is ready for timesteps. function dflowfm_InitSimulation(mduOrDimrConfigFile) result(retVal) implicit none integer :: retVal character(len=*) :: mduOrDimrConfigFile ! externals integer, external :: bmi_initialize integer :: numnod, ndx2d, ndxi retVal = bmi_initialize(dllHandle, mduOrDimrConfigFile) if (retVal /= 0) then ! print *, 'Error in bmi_initialize: ', retVal, ', config file: ', mduOrDimrConfigFile ! stop' FATAL' endif numnod = dflowfm_PutNumberOf1DNodes() ndx2d = dflowfm_PutNumberOf2DGridCells() ndxi = dflowfm_PutNumberOf1D2DGridCells() if (allocated(cumulative_volumes_in_prev_MSW_timestep)) deallocate(cumulative_volumes_in_prev_MSW_timestep) allocate(cumulative_volumes_in_prev_MSW_timestep(ndxi)) cumulative_volumes_in_prev_MSW_timestep = 0.0d0 if (allocated(cumulative_volumes_in_prev_MF_timestep)) deallocate(cumulative_volumes_in_prev_MF_timestep) allocate(cumulative_volumes_in_prev_MF_timestep(ndxi)) cumulative_volumes_in_prev_MF_timestep = 0.0d0 if (allocated(cumulative_volumes_in)) deallocate(cumulative_volumes_in) allocate(cumulative_volumes_in(ndxi)) cumulative_volumes_in = 0.0d0 call init_kdtree() end function dflowfm_InitSimulation subroutine init_kdtree() integer :: ierr, i integer, external :: bmi_get_var integer :: ndx1d, ndx2d double precision, dimension(:), allocatable :: xz, yz double precision, dimension(:,:), allocatable :: coordinates ! get nr. of nodes 1D + 2D: ndx1d = dflowfm_PutNumberOf1DNodes() ndx2d = dflowfm_PutNumberOf2DGridCells() ! get (x,y) allocate(xz(ndx2d + ndx1d), yz(ndx2d + ndx1d)) allocate(coordinates(2,ndx2d + ndx1d)) ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'xz', coordinates(1,:), ndx2d + ndx1d) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(xz): ', ierr ! stop' FATAL' endif ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'yz', coordinates(2,:), ndx2d + ndx1d) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(yz): ', ierr ! stop' FATAL' endif coordinate1D_tree => kdtree2_create(coordinates(:,ndx2d+1:),sort=.false.,rearrange=.true.) ! creating the tree for 1D nodes coordinate2D_tree => kdtree2_create(coordinates(:,1:ndx2d),sort=.false.,rearrange=.true.) ! creating the tree for 2D grid cells ! set radius for matching: search_radius = 0.0001d+0 end subroutine init_kdtree !> Performs (one or more) timesteps to reach the requested targetTime. subroutine dflowfm_performSurfacewaterTimestep(currentTime, targetTime, idtsw) implicit none double precision, intent(in) :: currentTime, targetTime !< period to be computed, times in MJD integer , intent(in) :: idtsw !< msw top system step within modflow step ! externals integer, external :: bmi_update double precision :: timestep logical :: success integer :: ierr ! NOTE: assume that given currentTime is correct (responsibility of caller) timestep = (86400.0D+0)*(targetTime - currentTime) ! get volumes from prev timestep success = getCumulativeVolumes(cumulative_volumes_in_prev_MSW_timestep) if (idtsw == 1) then success = getCumulativeVolumes(cumulative_volumes_in_prev_MF_timestep) endif ! perform time step ierr = bmi_update(dllHandle, timestep) if (ierr /= 0) then ! print *, 'Error in bmi_update: ', ierr, ', currentTime/targetTime: ', & ! currentTime, ' / ', targetTime ! stop' FATAL' endif lastTimeStep = timestep logTimestep = logTimestep + 1 end subroutine dflowfm_performSurfacewaterTimestep !> Returns the user time step, specific for dflowfm. double precision function dflowfm_PutUserTimestep() result(res) ! externals integer, external :: bmi_get_time_step integer :: ierr ierr = bmi_get_time_step(dllHandle, res) if (ierr /= 0) then ! print *, 'Error in bmi_get_time_step: ', ierr ! stop' FATAL' endif end function dflowfm_PutUserTimestep double precision function dflowfm_PutEndTime() result(res) ! externals integer, external :: bmi_get_end_time integer :: ierr ierr = bmi_get_end_time(dllHandle, res) if (ierr /= 0) then ! print *, 'Error in bmi_get_end_time: ', ierr ! stop' FATAL' endif end function dflowfm_PutEndTime !> Finalizes the current dflowfm surface water timestep(s). subroutine dflowfm_finishSurfacewaterTimestep() implicit none ! NOTE: if performSurfacewaterTimestep calls update(), then no finalize is needed. return end subroutine dflowfm_finishSurfacewaterTimestep !> Finishes the current model simulation. subroutine dflowfm_finishSimulation() implicit none ! externals integer, external :: bmi_finalize integer :: ierr if (allocated(cumulative_volumes_in_prev_MSW_timestep)) deallocate(cumulative_volumes_in_prev_MSW_timestep) if (allocated(cumulative_volumes_in_prev_MF_timestep)) deallocate(cumulative_volumes_in_prev_MF_timestep) if (allocated(cumulative_volumes_in)) deallocate(cumulative_volumes_in) ierr = bmi_finalize(dllHandle) if (ierr /= 0) then ! print *, 'Error in bmi_finalize: ', ierr ! stop' FATAL' endif end subroutine dflowfm_finishSimulation !----------------------------- 1D nodes------------------------------ ! Intended use: infiltration to and drainage from groundwater, combined with ! sprinkling demand/allocation from unsaturated zone !> Returns the number of 1D nodes. function dflowfm_PutNumberOf1DNodes() result(nFlowCells) implicit none ! externals integer, external :: bmi_get_var integer :: nFlowCells integer :: ierr integer :: ndxi ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ndxi', ndxi, 1) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ndxi): ', ierr ! stop' FATAL' endif nFlowCells = ndxi - dflowfm_PutNumberOf2DGridCells() !nFlowCells = 21 end function dflowfm_PutNumberOf1DNodes !> Returns the integer node ids for all 1D nodes. !! NOTE: This is effectively the increasing range of 1D indices. function dflowfm_Put1DNodeIDs(ids) result(success) implicit none logical :: success integer, dimension(:), intent(out) :: ids ! len(ids) == number of 2D cells integer :: n, ndx2d, ndx1d ndx1d = dflowfm_PutNumberOf1DNodes() ndx2d = dflowfm_PutNumberOf2DGridCells() do n=1,ndx1d ids(n) = ndx2d + n end do !ids = (/ 119, 118, 117, 116, 115, 114, & ! 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, & ! 113, 112, 111, 110, 109 /) success = .true. end function dflowfm_Put1DNodeIDs !> Returns the surface water levels for all 1D nodes. function dflowfm_Put1DNodeSurfacewaterLevels(levels) result(success) implicit none logical :: success double precision, dimension(:), & intent(out) :: levels !< len(ids) == number of 1D Nodes !< units: m AD NAP !< values ordered according to ids1Dnodes array ! externals integer, external :: bmi_get_var integer :: ierr integer :: n, ndxi, ndx2d, ndx1d double precision, allocatable :: bmivalues(:) ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ndxi', ndxi, 1) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ndxi): ', ierr ! stop' FATAL' endif allocate(bmivalues(ndxi)) ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'s1', bmivalues, ndxi) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(s1): ', ierr ! stop' FATAL' endif ndx1d = dflowfm_PutNumberOf1DNodes() ndx2d = dflowfm_PutNumberOf2DGridCells() do n=1,ndx1d levels(n) = bmivalues(ndx2d+n) end do !do c = 1, dflowfm_PutNumberOf2DGridCells() ! levels = 10.D0 + 0.001D0 * c !enddo deallocate(bmivalues) success = .true. end function dflowfm_Put1DNodeSurfacewaterLevels !> Sets in dflowfm the requested fluxes 1D function dflowfm_Get1DNodeFlux(fluxes, missingValue) result(success) logical :: success double precision, dimension(:), & intent(in) :: fluxes !< len(fluxes) == number of 1D Nodes !< units: m3/s ! values ordered according to ids1Dnodes array real :: missingValue ! externals integer, external :: bmi_set_var_slice integer :: ierr, i, qinextFileHandle integer :: ndx2d, ndx1d character(len=128) :: qinextFileName ndx1d = dflowfm_PutNumberOf1DNodes() ndx2d = dflowfm_PutNumberOf2DGridCells() call filterMissingValues(fluxes, dble(missingValue), .false.) ierr = bmi_set_var_slice(dllHandle, trim(modelPrefix)//'qext', fluxes, ndx2d, ndx1d) if (ierr /= 0) then ! print *, 'Error in bmi_set_var_slice(qext - 1D): ', ierr ! stop' FATAL' endif if (doDebugLogging) then write(qinextFileName, '(A,I4.4,A)') '1d-node-fluxes-step-', logTimestep, '.log' open(newunit=qinextFileHandle, file=qinextFileName) write(qinextFileHandle, '(A)') 'bmi_set_var_slice(dllHandle, ''qext'', fluxes, ndx2d, ndx1d)' write(qinextFileHandle, '(A,I6,A,I6)') ' ndx2d =' , ndx2d, ' ndx1d =' , ndx1d do i = 1, ndx1d write(qinextFileHandle, '(F16.10)') fluxes(i) enddo close(qinextFileHandle) endif success = .true. end function dflowfm_Get1DNodeFlux !> Returns the realised fluxes for all 1D grid nodes. function dflowfm_PutRealised1DNodeFlux(flux) result(success) implicit none logical :: success double precision, dimension(:), intent(inout) :: flux !< m3/s, length: number of 1D grid points integer :: n, ndx1d, ndx2d integer :: ierr ierr = getCumulativeVolumes(cumulative_volumes_in) ndx1d = dflowfm_PutNumberOf1DNodes() ndx2d = dflowfm_PutNumberOf2DGridCells() do n=1,ndx1d flux(n) = (cumulative_volumes_in(ndx2d+n) - cumulative_volumes_in_prev_MF_timestep(ndx2d+n)) / lastTimestep enddo success = .true. end function dflowfm_PutRealised1DNodeFlux !> Returns the realised fluxes for all 2D grid cells. function dflowfm_PutRealised2DCellFlux(flux) result(success) implicit none logical :: success double precision, dimension(:), intent(inout) :: flux !< m3/s, length: number of 2D grid cells integer :: n, ndx2d integer :: ierr ierr = getCumulativeVolumes(cumulative_volumes_in) ndx2d = dflowfm_PutNumberOf2DGridCells() do n=1,ndx2d flux(n) = (cumulative_volumes_in(n) - cumulative_volumes_in_prev_MSW_timestep(n)) / lastTimestep enddo success = .true. end function dflowfm_PutRealised2DCellFlux !> Gets the cumulative volumes from external discharges for all grid nodes ('vincum') function getCumulativeVolumes(volumes) result(success) logical :: success double precision, dimension(:), intent(inout) :: volumes !< m3, length: number of 1D+2D grid cells ('ndxi') integer, external :: bmi_get_var integer :: ndxi, ierr ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ndxi', ndxi, 1) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ndxi): ', ierr ! stop' FATAL' endif ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'vextcum ', volumes, ndxi) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(vextcum): ', ierr ! stop' FATAL' endif success = .true. end function getCumulativeVolumes function dflowfm_PutNumberOf1D2DGridCells() result(n1D2DCells) implicit none ! externals integer, external :: bmi_get_var integer :: n1D2DCells integer :: ierr integer :: ndxi ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ndxi', ndxi, 1) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ndxi): ', ierr ! stop' FATAL' endif n1D2DCells = ndxi end function dflowfm_PutNumberOf1D2DGridCells !----------------------------- 2D ------------------------------ ! Intended use: ponding/infiltration on 2D surface level. function dflowfm_PutNumberOf2DGridCells() result(nFlowCells) implicit none ! externals integer, external :: bmi_get_var integer :: nFlowCells integer :: ierr integer :: ndx2d ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ndx2d', ndx2d, 1) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ndx2d): ', ierr ! stop' FATAL' endif nFlowCells = ndx2d end function dflowfm_PutNumberOf2DGridCells !> Returns the integer ids for all 2D nodes. !> Returns the integer node ids for all 1D nodes. !! NOTE: This is effectively the increasing range of 1D indices. function dflowfm_Put2DGridCellIDs(ids) result(success) implicit none logical :: success integer, dimension(:), intent(out) :: ids ! len(ids) == number of 2D cells ! externals integer, external :: bmi_get_var integer :: ierr integer :: n, ndxi ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ndxi', ndxi, 1) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ndxi): ', ierr ! stop' FATAL' endif do n = 1, dflowfm_PutNumberOf2DGridCells() ids(n) = n enddo success = .true. end function dflowfm_Put2DGridCellIDs !> returns internal id for 1D grid cell closest to (x,y), and -1 when fails function dflowfm_MapCoordinateTo1DCellId(x, y, id) result(success) implicit none logical :: success double precision, intent(in) :: x double precision, intent(in) :: y integer, intent(out) :: id integer :: ndx2d call find_nearest_node(coordinate1D_tree, x, y, id) if (id == -1) then success = .false. return end if ! offset: 1D coordinates are stacked after 2D coordinates ndx2d = dflowfm_PutNumberOf2DGridCells() id = id + ndx2d success = .true. end function !> returns internal id for 2D grid cell closest to (x,y), and -1 when fails function dflowfm_MapCoordinateTo2DCellId(x, y, id) result(success) implicit none logical :: success double precision, intent(in) :: x double precision, intent(in) :: y integer, intent(out) :: id call find_nearest_node(coordinate2D_tree, x, y, id) if (id == -1) then success = .false. return end if success = .true. end function !> Find nearest node in kdtree within tolerance. !> Lookup failed will return id = -1 subroutine find_nearest_node(tree, x, y, id) implicit none type(kdtree2), pointer, intent(in) :: tree ! the kdtree with grid coordinates double precision, intent(in) :: x double precision, intent(in) :: y integer, intent(out) :: id integer :: nf ! number of hits, should be 1 in our case real(kdkind) :: coordinate(2) ! coordinate to match type(kdtree2_result) :: results(1) ! search results ! search crashes on empty tree... if (tree.n == 0) then id = -1 return end if coordinate = (/ x, y /) call kdtree2_r_nearest(tp=tree, qv=coordinate, r2=search_radius, nfound=nf, nalloc=1, results=results) if (nf == 0) then id = -1 return end if id = results(1)%idx end subroutine !> Returns the cell areas for all 2D grid cells. function dflowfm_Put2DGridCellAreas(areas) result(success) implicit none logical :: success double precision, dimension(:), intent(out) :: areas ! len(areas) == number of 2D grid cells ! externals integer, external :: bmi_get_var integer :: ierr integer :: n, ndxi double precision, dimension(:), allocatable :: bmivalues ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ndxi', ndxi, 1) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ndxi): ', ierr ! stop' FATAL' endif ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ba', bmivalues, ndxi) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ba): ', ierr ! stop' FATAL' endif do n=1,dflowfm_PutNumberOf2DGridCells() areas(n) = bmivalues(n) end do !areas = -1D0 success = .true. end function dflowfm_Put2DGridCellAreas !> Returns the current water levels for all 2D grid cells. function dflowfm_Put2DWaterlevels(SurfacewaterLevels) result(success) implicit none logical :: success double precision, dimension(:), intent(out) :: SurfacewaterLevels ! externals integer, external :: bmi_get_var integer :: ierr integer :: n, ndxi double precision, dimension(:), allocatable :: bedLevels, waterLevels ! TODO: get bedLevels once ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'ndxi', ndxi, 1) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(ndxi): ', ierr ! stop' FATAL' endif allocate(bedLevels(ndxi)) allocate(waterLevels(ndxi)) ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'bl', bedLevels, ndxi) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(bl): ', ierr ! stop' FATAL' endif ierr = bmi_get_var(dllHandle, trim(modelPrefix)//'s1', waterLevels, ndxi) if (ierr /= 0) then ! print *, 'Error in bmi_get_var(s1): ', ierr ! stop' FATAL' endif do n=1,dflowfm_PutNumberOf2DGridCells() ! TODO: decide on epsilon if (waterLevels(n) > (bedLevels(n) + 0.001D+0)) then SurfacewaterLevels(n) = waterLevels(n) else SurfacewaterLevels(n) = bedLevels(n) ! TODO missing value? -99999.0D+0 endif end do !do c = 1, dflowfm_PutNumberOf2DGridCells() ! SurfacewaterLevels(c) = 1.5d0+0.0003d0*c !enddo deallocate(bedLevels) deallocate(waterLevels) success = .true. end function dflowfm_Put2DWaterlevels !> Sets into dflowfm the requested Delta volumes for ponding on all 2D grid cells. logical function DflowFM_GetFluxes2D(delvol, missVal) result(success) double precision, dimension(:), intent(in) :: delvol !< len(ids) == number of 2D grid cells !< units: m3 real, intent(in) :: missVal ! missingValue, TODO: handle in a generic way ! values ordered according to grid cell ID array ! externals integer, external :: bmi_set_var_slice integer :: ierr, i, qinextFileHandle integer :: ndx2d character(len=128) :: qinextFileName call filterMissingValues(delvol, dble(missVal), .false.) ndx2d = dflowfm_PutNumberOf2DGridCells() ierr = bmi_set_var_slice(dllHandle, trim(modelPrefix)//'qext', delvol, 0, ndx2d) if (ierr /= 0) then ! print *, 'Error in bmi_set_var_slice(qext - 2D): ', ierr ! stop' FATAL' endif if (doDebugLogging) then write(qinextFileName, '(A,I4.4,A)') '2d-cell-fluxes-step-', logTimestep, '.log' open(newunit=qinextFileHandle, file=qinextFileName) write(qinextFileHandle, '(A)') 'bmi_set_var_slice(dllHandle, ''qext'', delvol, 0, ndx2d)' write(qinextFileHandle, '(A,I6)') ' ndx2d =' , ndx2d do i = 1, ndx2d write(qinextFileHandle, '(F16.10)') delvol(i) enddo close(qinextFileHandle) endif success = .true. end function DflowFM_GetFluxes2D !> replace missing values by 0.0 subroutine filterMissingValues(arr, mv, setZero) double precision, dimension(:) :: arr double precision :: mv logical :: setZero integer :: i do i = 1, size(arr) if (setZero) then arr(i) = 0.0D+0 else if (dabs(arr(i) - dble(mv)) < 0.001D+0) then arr(i) = 0.0D+0 endif endif enddo end subroutine end module dflowfmmod