!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !> mfmsdfm_mapping: Mapping module for !> i*M*OD*F*LOW / *M*eta*Sw*ap / *D*Flow-*FM* !> interaction !!! !!! (c) Deltares, may 2018 !!! !!! stef.Hummel@deltares.nl !!! !!! File organization: !!! !!! DECLARATIONS !!! - Data Type Definitions !!! - Module variables !!! !!! IMPLEMENTATION !!! !!! - Public Functions: !!! !!! . Initialization / general functions !!! - initialize mapping configs !!! - last error, finalize !!! . Map functions !!! - set source and target ids !!! - do mapping !!! !!! - Private Functions: !!! !!! . Parse mapping config file !!! . (TODO: dump configuration) !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module mfmsdfm_mapping implicit none !! !! Public constants/enumerations !! ! Mapper types: integer, parameter :: dmm_FM2MSW_WL = 1 ! water level D-Flow FM to MetaSWAP integer, parameter :: dmm_MSW2FM_DPV = 2 ! Delta Ponding Volume MetaSWAP to D-Flow FM integer, parameter :: dmm_MSW2FM_SPR = 3 ! Sprinkling Demand MetaSWAP to D-Flow FM integer, parameter :: dmm_FM2MF_WL = 4 ! water level D-Flow FM to iMODLOW integer, parameter :: dmm_MF2FM_Q = 5 ! fluxes iMODFLOW D-Flow FM integer, parameter :: dmmSrc = 1 ! index for source size in data structures integer, parameter :: dmmTgt = 2 ! index for target size in data structures integer, parameter :: DMMSRCTGT = 2 ! #source/target indices integer, parameter :: DMM_ID_LEN = 50 ! id-len for laterals !! !! Private constants/enumerations !! integer, parameter, private :: dmm_operat_sum = 1 ! operations for mapping functions integer, parameter, private :: dmm_operat_average = 2 integer, parameter, private :: dmm_operat_devide = 3 integer, parameter, private :: DMM_LK_BLOCK_SIZE = 10000! size for reallocing links integer, parameter, private :: DMM_SL_BLOCK_SIZE = 2 ! size for reallocing source_items in links integer, parameter, private :: dmmIndexInvalid = -99 ! item in source or target array not mapped ! ! link info target side ! type t_dmm_target integer :: id ! integer identifier (default) integer :: index ! index of target point in target array double precision :: totalDemand ! total demand towards target double precision :: x ! x-coord at target side double precision :: y ! y-coord at target side end type t_dmm_target ! ! link info source side ! type t_dmm_source integer :: id ! integer identifiers (default) integer :: index ! index of source point in source array double precision :: weight ! weight of source point (read from config file) double precision :: demand ! demand at source side double precision :: x ! x-coord at source side double precision :: y ! y-coord at source side end type t_dmm_source ! ! link info (source points -> target point) ! type t_dmm_src_tgt type(t_dmm_target) :: tgt ! point at target side (id, index) type(t_dmm_source), dimension(:), pointer :: srcs ! points at source side (id, index, weight, demand) integer :: numSrcs ! #points at source side integer :: maxNumSrcs ! max #points at source side that the link can store (needed for reallocing) logical :: logMapping ! should the exchange at this point be logged? type(t_dmm_src_tgt) , pointer :: subLink ! additional link, for combining mappers towards points ! (for combining modflow fluxes with metaswap sprinkling demands) end type t_dmm_src_tgt ! TODO: rename to t_dmm_link ! ! mapping configuration ! type t_dmm_mapping character(Len=DMM_ID_LEN), & dimension(dmmSrcTgt) :: componentID ! missing value at source and target side double precision, & dimension(dmmSrcTgt) :: missingValue ! missing value at source and target side integer :: mapperType character(Len=DMM_ID_LEN) :: varName logical :: hasXY type(t_dmm_src_tgt), pointer, & dimension(:) :: links ! #links in the mapper integer :: numLinks ! #links in the mapper integer :: maxNumLinks ! max #links that the mapper can store (needed for reallocing) character(Len=256) :: configFile = '' ! Config File that has been read logical :: dumpMapping ! log values that are exchanged by mapper? character(Len=256) :: dumpFileName ! file name for logging values integer :: lundumpFile ! handle to file for logging values integer :: tStart ! start time step for logging (modflow time steps) integer :: tEnd ! end time step for logging (modflow time steps) logical :: useArea ! mapping area specified? double precision :: xLow ! x-coord of lower left corner of mapping area double precision :: yLow ! y-coord of lower left corner of mapping area double precision :: xUpp ! x-coord of upper right corner of mapping area double precision :: yUpp ! y-coord of upper right corner of mapping area integer, pointer, dimension(:) :: index_lookup ! mapping from id => cell index (TODO: rename to link_idx_lookup) end type t_dmm_mapping ! v integer, parameter :: DMM_MAX_MAPPINGS = 20 type(t_dmm_mapping), dimension(DMM_MAX_MAPPINGS), target :: mappings integer :: numMappings character(Len=1024) :: dmmLastError ! TODO: panic shortcut, make this smart/efficient integer, parameter :: LOOKUP_BLOCKSIZE = 200000 contains subroutine dmmGetIdent(retVal) ! arguments character(Len=*) :: retval ! locals ! body retVal = 'Copy of dmm (ds Rev.@@@)' ! TODO call getfullversionstring_dmm(retval) end subroutine dmmGetIdent !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! PUBLIC functions, INITIALIZATION / GENERAL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine dmmInitialize() ! locals integer :: m, st ! loop counters for mappings and source/target-side ! body numMappings = 0 do m = 1, DMM_MAX_MAPPINGS mappings(m) % configFile = ' ' mappings(m) % mapperType = -1 mappings(m) % varName = ' ' mappings(m) % hasXY = .false. mappings(m) % links => NULL() mappings(m) % numLinks = 0 mappings(m) % dumpMapping = .false. mappings(m) % dumpFileName = ' ' mappings(m) % lundumpFile = -1 mappings(m) % tStart = 0 mappings(m) % tEnd = 1E+6 mappings(m) % links => NULL() do st = 1, DMMSRCTGT mappings(m) % missingValue(st) = -99999.0D0 mappings(m) % componentID(st) = ' ' enddo allocate(mappings(m) % index_lookup(LOOKUP_BLOCKSIZE)) mappings(m) % index_lookup = -1 enddo dmmLastError = '' end subroutine dmmInitialize function dmmCreateMapping(mappingFile, mapperType, srcCompName, tgtCompName) result(mh) ! result integer :: mh ! mapping handle. >0: OK, <=0: error ! arguments character(len=*) , intent(in) :: mappingFile ! mapping configuration file integer , intent(in) :: mapperType ! type of mapper (see enum at the start of this file) character(Len=*) , intent(in) :: srcCompName ! source component name character(Len=*) , intent(in) :: tgtCompName ! target component name ! locals integer :: fileNameLen ! length of mapping file name integer :: retVal ! result of functions for reading settings and logical :: hasXY ! does mapping configuration file contain x/y's? ! (is the case when the file name ends with .dmm) hasXY = .false. fileNameLen = len_trim(mappingFile) if (mappingFile(fileNameLen-3:fileNameLen) == '.dmm' .or. mappingFile(fileNameLen-3:fileNameLen) == '.DMM') then hasXY = .true. endif mh = -1 if(numMappings -1) then close(mappings(m) % lundumpFile) endif if (associated(mappings(m) % links)) then do lk = 1, mappings(m) % numLinks if (associated(mappings(m) % links(lk) % srcs)) then deallocate(mappings(m) % links(lk) % srcs) endif enddo endif enddo call dmmInitialize() end subroutine dmmFinish !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! PUBLIC functions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Mapping: set ids at source side (and build index mapping) ! function dmmSetSourceIds(mh, ids) result(success) ! return value logical :: success ! .false.: error ! arguments integer, intent(in) :: mh ! handle to mapping integer, dimension(:), intent(in) :: ids ! ids at source side ! locals integer :: i, lk, sl ! loop counters (ids, links, srsIds in link) integer, dimension(:), allocatable :: idxlookup ! lookup table for indices integer :: maxidx integer :: tid ! id for which to lookup index integer :: allocstat ! status for alloc ! build lookup table id => index maxidx = maxval(ids) allocate(idxlookup(maxidx), STAT=allocstat) if (allocstat /= 0) then write(*,*) 'ERROR allocating lookup table with size: ', maxidx stop endif idxlookup = -1 do i = 1, size(ids) idxlookup(ids(i)) = i enddo do lk = 1, mappings(mh) % numLinks do sl = 1, mappings(mh) % links(lk) % numSrcs tid = mappings(mh) % links(lk) % srcs(sl) % id ! check if id is mapped to index: if ( tid == dmmIndexInvalid .or. tid == -1 ) then write(*,*) 'ERROR in mh:', mh, ' lk:', lk, 'sl', sl, ' id:', tid else if (idxlookup(tid) > 0) then mappings(mh) % links(lk) % srcs(sl) % index = idxlookup(tid) endif endif enddo enddo do lk = 1, mappings(mh) % numLinks do sl = 1, mappings(mh) % links(lk) % numSrcs if ( mappings(mh) % links(lk) % srcs(sl) % index == dmmIndexInvalid ) then write(*,*) 'ERROR in mh:', mh, ' lk:', lk, 'sl', sl, ' id:', mappings(mh) % links(lk) % srcs(sl) % id endif enddo enddo deallocate(idxlookup) success = .true. end function dmmSetSourceIds ! ! Mapping: set ids at target side (and build index mapping) ! function dmmSetTargetIds(mh, ids) result(success) ! return value logical :: success ! .false.: error ! arguments integer, intent(in) :: mh ! handle to mapping integer, dimension(:), intent(in) :: ids ! ids at target side ! locals integer :: i, lk ! loop counters (ids, links, srsIds in link) integer, dimension(:), allocatable :: idxlookup ! lookup table for indices integer :: maxidx integer :: tid ! id for which to lookup index integer :: allocstat ! status for alloc ! build lookup table id => index maxidx = maxval(ids) allocate(idxlookup(maxidx), STAT=allocstat) if (allocstat /= 0) then write(*,*) 'ERROR allocating lookup table with size: ', maxidx stop endif idxlookup = -1 do i = 1, size(ids) idxlookup(ids(i)) = i enddo do lk = 1, mappings(mh) % numLinks tid = mappings(mh) % links(lk) % tgt % id ! check if id is mapped to index: if (idxlookup(tid) > 0) then mappings(mh) % links(lk) % tgt % index = idxlookup(tid) endif enddo do lk = 1, mappings(mh) % numLinks if ( mappings(mh) % links(lk) % tgt % index == dmmIndexInvalid ) then write(*,*) 'ERROR in mh:', mh, ' lk:', lk, ' id:', mappings(mh) % links(lk) % tgt % id endif enddo success = .true. end function dmmSetTargetIds function dmmCombineMappers(mhFlux1dMapper, mhSprinklingMapper) result(success) ! return value logical :: success ! .false.: error ! arguments integer, intent(in) :: mhFlux1dMapper ! handle to mapper for fluxes between iMod and D-Flow FM integer, intent(in) :: mhSprinklingMapper ! handle to mapper for sprinkling demands MetaSwap to D-Flow FM ! locals integer :: fLink, sLink ! loop counter for links and source items integer :: tgtId ! target id type(t_dmm_mapping), pointer :: flux1DMapper type(t_dmm_mapping), pointer :: sprinkMapper ! body flux1DMapper => mappings(mhFlux1dMapper) sprinkMapper => mappings(mhSprinklingMapper) do fLink = 1, flux1DMapper % numLinks tgtId = flux1DMapper % links(fLink) % tgt % id do sLink = 1, sprinkMapper % numLinks if (sprinkMapper % links(sLink) % tgt % id == tgtId) then flux1DMapper % links(fLink) % subLink => sprinkMapper % links(sLink) endif enddo end do success = .true. end function dmmCombineMappers ! ! Mapping: perform mapping ! subroutine dmmMapValues(mh, sourceValues, targetValues) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(in) :: sourceValues ! values at source side double precision, dimension(:), intent(out) :: targetValues ! computed values for target side ! locals integer :: link, s ! loop counter for links and source items integer :: indxTgt ! index at target side, index double precision :: srcVal ! value at source side character(len=32):: srcXyString, tgtXyString ! body if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) 'MAPPING ', trim(mappings(mh) % varName), ' from ', trim(mappings(mh) % componentID(dmmSrc)), ' to ', trim(mappings(mh) % componentID(dmmTgt)) targetValues = mappings(mh) % missingValue(dmmTgt) do link = 1, mappings(mh) % numLinks indxTgt = mappings(mh) % links(link) % tgt % index targetValues(indxTgt) = 0 do s = 1, mappings(mh) % links(link) % numSrcs srcVal = sourceValues(mappings(mh) % links(link) % srcs(s) % index) if(.not. dabs(srcVal - mappings(mh) % missingValue(dmmSrc))<1.0d-9) then if (dabs(targetValues(indxTgt) - mappings(mh) % missingValue(dmmTgt))<1.0d-9) then ! first time source delivers a value targetValues(indxTgt) = 0 endif targetValues(indxTgt) = targetValues(indxTgt) + & srcVal * mappings(mh) % links(link) % srcs(s) % weight mappings(mh) % links(link) % srcs(s) % demand = srcVal * mappings(mh) % links(link) % srcs(s) % weight if(mappings(mh) % mapperType == dmm_MF2FM_Q) then ! convert iMODFLOW volume to m3/s, swap sign mappings(mh) % links(link) % srcs(s) % demand = -mappings(mh) % links(link) % srcs(s) % demand / 86400.0D0 ! TODO: make modflow time step dynamic endif else ! No action (skip missing value) endif end do if (.not. dabs(targetValues(indxTgt) - mappings(mh) % missingValue(dmmTgt))<1.0d-9) then if(mappings(mh) % mapperType == dmm_MF2FM_Q) then ! convert iMODFLOW volume to m3/s, swap sign targetValues(indxTgt) = - targetValues(indxTgt) / 86400.0D0 ! TODO: make modflow time step dynamic endif if(mappings(mh) % mapperType == dmm_MSW2FM_DPV) then ! convert MetaSWAP volume to m3/s targetValues(indxTgt) = targetValues(indxTgt) / 86400.0D0 ! TODO: make metaswap time step dynamic endif endif tgtXyString = ' -,-' srcXyString = ' -,-' mappings(mh) % links(link) % tgt % totalDemand = targetValues(indxTgt) if (mappings(mh) % dumpMapping .and. mappings(mh) % links(link) % logMapping) then if (mappings(mh) % hasXY .and. (mappings(mh) % mapperType == dmm_MSW2FM_DPV .or. dmm_MF2FM_Q)) then write(tgtXyString, '(F10.3,A,F10.3)') mappings(mh) % links(link) % tgt % x, & ',', mappings(mh) % links(link) % tgt % y endif do s = 1, mappings(mh) % links(link) % numSrcs if (mappings(mh) % hasXY .and. (mappings(mh) % mapperType == dmm_FM2MSW_WL .or. dmm_FM2MF_WL)) then write(srcXyString, '(F10.3,A,F10.3)') mappings(mh) % links(link) % srcs(s) % x, & ',', mappings(mh) % links(link) % srcs(s) % y endif if (s < mappings(mh) % links(link) % numSrcs) then write(mappings(mh) % lunDumpFile,'(A,A,A,I6,A,F16.8)') & 'x/y: ', trim(srcXyString), 'src(', mappings(mh) % links(link) % srcs(s) % id, ') =', sourceValues(mappings(mh) % links(link) % srcs(s) % index) else write(mappings(mh) % lunDumpFile,'(A,A,A,I6,A,F16.8,A,I6,A,F16.8,A,A)') & 'x/y: ', trim(srcXyString), ' src(', mappings(mh) % links(link) % srcs(s) % id, ') =', sourceValues(mappings(mh) % links(link) % srcs(s) % index), & ' tgt(', mappings(mh) % links(link) % tgt % id, ') =', targetValues(indxTgt), ' x/y: ', trim(tgtXyString) endif enddo call flush(mappings(mh) % lunDumpFile) endif end do if (mappings(mh) % dumpMapping) then write(mappings(mh) % lunDumpFile,*) 'END MAPPING from ', trim(mappings(mh) % componentID(dmmSrc)), ' to ', trim(mappings(mh) % componentID(dmmTgt)) write(mappings(mh) % lunDumpFile,*) '' call flush(mappings(mh) % lunDumpFile) endif end subroutine dmmMapValues subroutine dmmCombinedMapperSetFluxes(mh, sourceValues) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(in) :: sourceValues ! values at source side ! locals integer :: link, s ! loop counter for links and source items integer :: indxTgt ! index at target side, index double precision :: srcVal ! value at source side ! body if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) 'MAPPING ', trim(mappings(mh) % varName), ' from ', trim(mappings(mh) % componentID(dmmSrc)), ' to ', trim(mappings(mh) % componentID(dmmTgt)) ! targetValues = mappings(mh) % missingValue(dmmTgt) do link = 1, mappings(mh) % numLinks indxTgt = mappings(mh) % links(link) % tgt % index do s = 1, mappings(mh) % links(link) % numSrcs srcVal = sourceValues(mappings(mh) % links(link) % srcs(s) % index) if(.not. dabs(srcVal - mappings(mh) % missingValue(dmmSrc))<1.0d-9) then mappings(mh) % links(link) % srcs(s) % demand = srcVal * mappings(mh) % links(link) % srcs(s) % weight if(mappings(mh) % mapperType == dmm_MF2FM_Q) then ! convert iMODFLOW volume to m3/s, swap sign mappings(mh) % links(link) % srcs(s) % demand = -mappings(mh) % links(link) % srcs(s) % demand / 86400.0D0 ! TODO: make modflow time step dynamic endif mappings(mh) % links(link) % tgt % totalDemand = mappings(mh) % links(link) % tgt % totalDemand + & mappings(mh) % links(link) % srcs(s) % demand else ! No action (skip missing value) endif end do if (mappings(mh) % dumpMapping .and. mappings(mh) % links(link) % logMapping) & write(mappings(mh) % lunDumpFile,'(A,I6,A,F16.8,A,I6,A,F16.8)') & 'src(', mappings(mh) % links(link) % srcs(1) % id, ') =', sourceValues(mappings(mh) % links(link) % srcs(1) % index), & ' tgt(', mappings(mh) % links(link) % tgt % id, ') =', mappings(mh) % links(link) % tgt % totalDemand if (mappings(mh) % dumpMapping) call flush(mappings(mh) % lunDumpFile) end do if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) 'END MAPPING from ', trim(mappings(mh) % componentID(dmmSrc)), ' to ', trim(mappings(mh) % componentID(dmmTgt)) if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) '' if (mappings(mh) % dumpMapping) call flush(mappings(mh) % lunDumpFile) end subroutine dmmCombinedMapperSetFluxes subroutine dmmCombinedMapperGetFluxes(mh, targetValues) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(out) :: targetValues ! computed values for target side ! locals integer :: link ! loop counter for links integer :: indxTgt ! index at target side, index ! body if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) 'COMBINING ', trim(mappings(mh) % varName), ' to ', trim(mappings(mh) % componentID(dmmTgt)) targetValues = mappings(mh) % missingValue(dmmTgt) do link = 1, mappings(mh) % numLinks indxTgt = mappings(mh) % links(link) % tgt % index targetValues(indxTgt) = mappings(mh) % links(link) % tgt % totalDemand if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,'(A,I6,A,F16.8,A,I6,A,F16.8)') 'ID(', mappings(mh) % links(link) % tgt % id, & '), Flux1: ', mappings(mh) % links(link) % tgt % totalDemand, & ', Flux2: ' , mappings(mh) % links(link) % subLink % tgt % totalDemand, & ' total: ' , mappings(mh) % links(link) % tgt % id, ') =', targetValues(indxTgt) if (mappings(mh) % dumpMapping) call flush(mappings(mh) % lunDumpFile) end do if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) 'END COMBINING to ', trim(mappings(mh) % componentID(dmmTgt)) if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) '' if (mappings(mh) % dumpMapping) call flush(mappings(mh) % lunDumpFile) end subroutine dmmCombinedMapperGetFluxes subroutine dmmCombinedMapperSetRealizedFluxes(mh, realizedValues) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(out) :: realizedValues ! computed values for target side ! locals integer :: link ! loop counter for links integer :: indxTgt ! index at target side, index ! body if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) 'COMBINING ', trim(mappings(mh) % varName), ' to ', trim(mappings(mh) % componentID(dmmTgt)) realizedValues = mappings(mh) % missingValue(dmmTgt) do link = 1, mappings(mh) % numLinks indxTgt = mappings(mh) % links(link) % tgt % index realizedValues(indxTgt) = mappings(mh) % links(link) % tgt % totalDemand if (associated(mappings(mh) % links(link) % subLink)) then realizedValues(indxTgt) = realizedValues(indxTgt) + & mappings(mh) % links(link) % subLink % tgt % totalDemand endif if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,'(A,I6,A,F16.8,A,I6,A,F16.8)') 'ID(', mappings(mh) % links(link) % tgt % id, & '), Flux1: ', mappings(mh) % links(link) % tgt % totalDemand, & ', Flux2: ' , mappings(mh) % links(link) % subLink % tgt % totalDemand, & ' total: ' , mappings(mh) % links(link) % tgt % id, ') =', realizedValues(indxTgt) if (mappings(mh) % dumpMapping) call flush(mappings(mh) % lunDumpFile) end do if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) 'END COMBINING to ', trim(mappings(mh) % componentID(dmmTgt)) if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) '' if (mappings(mh) % dumpMapping) call flush(mappings(mh) % lunDumpFile) end subroutine dmmCombinedMapperSetRealizedFluxes ! ! Mapping: perform mapping back (realized values) ! subroutine dmmMapValuesBack(mh, valuesFromTargetComponent, realizedSourceValues) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(in) :: valuesFromTargetComponent ! values as realized at the target side double precision, dimension(:), intent(out) :: realizedSourceValues ! mapped realized values at the source side ! locals double precision, dimension(:), & allocatable :: realizedTargetValues ! copy of valuesFromTargetComponent, for swapping sign integer :: link, s ! loop counter for links and source items integer :: indxSrc, indxTgt ! index at source and target side double precision :: totalRealizedValue ! total realized value at target side double precision :: totalDemandedValue ! total demanded value towards target double precision :: allOverCorrection ! all over correction flux double precision :: sourceDemandFraction ! fraction per source at the link side double precision :: epsilon = 1.0D-10 ! epsilon for comparing demanded and realized flux ! body if (mappings(mh) % dumpMapping) write(mappings(mh) % lunDumpFile,*) 'MAPPING realized ', trim(mappings(mh) % varName), ' from ', trim(mappings(mh) % componentID(dmmTgt)), ' back to ', trim(mappings(mh) % componentID(dmmSrc)) allocate(realizedTargetValues(size(valuesFromTargetComponent))) realizedTargetValues = valuesFromTargetComponent if (mappings(mh) % mapperType == dmm_MF2FM_Q) then ! swap sign realizedTargetValues = -realizedTargetValues endif realizedSourceValues = mappings(mh) % missingValue(dmmSrc) allOverCorrection = 0.0D+0 do link = 1, mappings(mh) % numLinks indxTgt = mappings(mh) % links(link) % tgt % index totalRealizedValue = valuesFromTargetComponent(indxTgt) totalDemandedValue = mappings(mh) % links(link) % tgt % totalDemand if (totalDemandedValue < 0 .and. totalDemandedValue < totalRealizedValue - epsilon) then ! distribute shortage do s = 1, mappings(mh) % links(link) % numSrcs indxSrc = mappings(mh) % links(link) % srcs(s) % index sourceDemandFraction = mappings(mh) % links(link) % srcs(s) % demand / totalDemandedValue if (mappings(mh) % mapperType == dmm_MF2FM_Q) then ! receiving party (modflow) expects correction instead of realized value realizedSourceValues(indxSrc) = & sourceDemandFraction * (totalRealizedValue - totalDemandedValue) allOverCorrection = allOverCorrection + realizedSourceValues(indxSrc) else realizedSourceValues(indxSrc) = & sourceDemandFraction * totalRealizedValue endif if (mappings(mh) % mapperType == dmm_MF2FM_Q) then ! convert m3/s to iMODFLOW volume, swap sign realizedSourceValues(indxSrc) = - realizedSourceValues(indxSrc) * 86400.0D0 ! TODO: make modflow time step dynamic else if (mappings(mh) % mapperType == dmm_MSW2FM_DPV) then ! convert m3/s to MetaSWAP volume realizedSourceValues(indxSrc) = realizedSourceValues(indxSrc) * 86400.0D0 ! TODO: make dtsw dynamic endif if (mappings(mh) % dumpMapping .and. mappings(mh) % links(link) % logMapping) then if (s == 1) then write(mappings(mh) % lunDumpFile,'(A,I6,A,F16.8,A,I6,A,F16.8)') & 'tgt(', mappings(mh) % links(link) % tgt % id, ')C=', valuesFromTargetComponent(indxTgt), & ' src(', mappings(mh) % links(link) % srcs(s) % id, ')C=', realizedSourceValues(mappings(mh) % links(link) % srcs(s) % index) else write(mappings(mh) % lunDumpFile,'(A,I6,A,F16.8)') & ' src(', mappings(mh) % links(link) % srcs(s) % id, ')C=', & realizedSourceValues(mappings(mh) % links(link) % srcs(s) % index) endif endif end do else ! no shortage, all 1D-fluxes and sprinkling demands realized do s = 1, mappings(mh) % links(link) % numSrcs indxSrc = mappings(mh) % links(link) % srcs(s) % index if (mappings(mh) % mapperType == dmm_MF2FM_Q) then ! receiving party (modflow) expects correction instead of realized value realizedSourceValues(indxSrc) = 0.0D+0 else realizedSourceValues(indxSrc) = & mappings(mh) % links(link) % srcs(s) % demand endif if (mappings(mh) % mapperType == dmm_MSW2FM_DPV) then realizedSourceValues(indxSrc) = realizedSourceValues(indxSrc) * 86400.0D0 ! TODO: make dtsw time step dynamic endif if (mappings(mh) % dumpMapping .and. mappings(mh) % links(link) % logMapping) then if (s == 1) then write(mappings(mh) % lunDumpFile,'(A,I6,A,F16.8,A,I6,A,F16.8)') & 'tgt(', mappings(mh) % links(link) % tgt % id, ') =', valuesFromTargetComponent(indxTgt), & ' src(', mappings(mh) % links(link) % srcs(s) % id, ') =', realizedSourceValues(mappings(mh) % links(link) % srcs(s) % index) else write(mappings(mh) % lunDumpFile,'(A,I6,A,F16.8)') & ' src(', mappings(mh) % links(link) % srcs(s) % id, ') =', & realizedSourceValues(mappings(mh) % links(link) % srcs(s) % index) endif endif end do if (mappings(mh) % dumpMapping .and. mappings(mh) % links(link) % logMapping) call flush(mappings(mh) % lunDumpFile) endif end do if (mappings(mh) % dumpMapping) then write(mappings(mh) % lunDumpFile,*) 'END MAPPING realized ', trim(mappings(mh) % varName), & ' from ', trim(mappings(mh) % componentID(dmmTgt)), ' back to ', trim(mappings(mh) % componentID(dmmSrc)) write(mappings(mh) % lunDumpFile,*) '' call flush(mappings(mh) % lunDumpFile) endif if (mappings(mh) % mapperType == dmm_MF2FM_Q) then write(2323,*) 'Total correction flux: ', allOverCorrection call flush(2323) endif deallocate(realizedTargetValues) end subroutine dmmMapValuesBack !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! private functions, read config !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Read configuration file ! function readMappingConfigFile(mh, confFileName, mapperType, hasXY) result(retVal) use dflowfmmod ! for determining node and cell indices from x/y coordinates ! return value integer :: retVal ! 0: everything OK ! arguments integer , intent(in) :: mh ! handle to mapping character(Len=*), intent(in) :: confFileName ! mapping config file integer , intent(in) :: mapperType ! type of mapper (see enum at the start of this file) logical , intent(in) :: hasXY ! does mapping configuration file contain x/y's (or indices only)? ! locals integer :: lunConfFile ! Handle to config file integer :: read_res ! result of read call integer :: srcId, tgtId ! id's at source and target size double precision :: x, y ! x/y coordinates source or target size character(Len=256) :: line ! line in config file logical :: eof ! yet at end of file? integer :: lk, s ! link counter, source index double precision :: weight ! weight at source side type(t_dmm_src_tgt), pointer :: link ! pointer to current link type(t_dmm_source), & dimension (:), pointer :: oldSrcs ! pointer to old sources (needed for reallocing) integer :: numOldSrcs ! #old sources (needed for reallocing) integer :: status ! allocation status flag logical :: success ! result of x/y->id mapping ! body write(*, '(2A)') 'dmm: Reading Config file ', trim(confFileName) ! open file open(newunit=lunConfFile,file=confFileName,status='old', readonly, iostat=read_res) ! read first line (TARGET, SOURCE, etc.) read(lunConfFile, '(a)', iostat=read_res) line ! read all configuration line eof = .false. do while (.not. eof) x = -777.7 y = -777.7 if (hasXY) then success = .false. if (mapperType == dmm_FM2MSW_WL) then ! file contains (FM)X, (FM)Y, SVAT-ID ! FM-cell is source, svat is target read(lunConfFile, *, iostat=read_res) x, y, tgtId if (read_res /= 0) then eof = .true. exit else success = dflowfm_MapCoordinateTo2DCellId(x, y, srcId) if (.not. success) then write(1920,*) 'ERROR getting source cell-id from FM, x:', x, ', y:', y, ', tgtId:', tgtId cycle endif endif ! currently no weights for MSW/DFM (1:1) weight = 1.0D+0 else if (mapperType == dmm_MSW2FM_DPV) then ! file contains SVAT-ID (FM)X, (FM)Y, ! svat is source, FM-cell is target read(lunConfFile, *, iostat=read_res) srcId, x, y if (read_res /= 0) then eof = .true. exit else success = dflowfm_MapCoordinateTo2DCellId(x, y, tgtId) if (.not. success) then write(1920,*) 'ERROR getting target cell-id from FM, x:', x, ', y:', y, ', srcId:', srcId cycle endif endif ! currently no weights for MSW/DFM (1:1) weight = 1.0D+0 else if (mapperType == dmm_FM2MF_WL) then ! file contains RIV-ID, (FM)X, (FM)Y, weight read(lunConfFile, *, iostat=read_res) tgtId, x, y, weight if (read_res /= 0) then eof = .true. exit else success = dflowfm_MapCoordinateTo1DCellId(x, y, srcId) if (.not. success) then write(1920,*) 'ERROR getting source node-id from FM, x:', x, ', y:', y, ', tgtId:', tgtId cycle endif endif else if (mapperType == dmm_MF2FM_Q) then ! file contains (FM)X, (FM)Y, RIV-ID read(lunConfFile, *, iostat=read_res) x, y, srcId if (read_res /= 0) then eof = .true. exit else success = dflowfm_MapCoordinateTo1DCellId(x, y, tgtId) if (.not. success) then write(1920,*) 'ERROR getting target node-id from FM, x:', x, ', y:', y, ', srcId:', srcId cycle endif ! no weights for discharges (will be summed) weight = 1.0D+0 endif endif else read(lunConfFile, *, iostat=read_res) tgtId, srcId, weight if (read_res /= 0) then eof = .true. exit endif endif ! this part will only be reached if (.not. eof .and. success) ! find or add link lk = dmmFindOrAddLink(mh, tgtId, .true.) ! TODO: handle many one to if (lk>0) then link => mappings(mh) % links(lk) link % tgt % id = tgtId if (mapperType == dmm_MF2FM_Q .or. mapperType == dmm_MSW2FM_DPV) then link % tgt % x = x link % tgt % y = y endif ! (re)alloc sources array if needed s = dmmIndexInvalid if (link % numSrcs == 0) then ! first source allocate(link % srcs(DMM_SL_BLOCK_SIZE), stat=status) if ( status /= 0 ) then write(dmmLastError, '(A)') 'dmm: readMappingConfigFile: Could not allocate initial block of source items' else link % maxNumSrcs = DMM_SL_BLOCK_SIZE s = 1 endif else if (link % numSrcs == link % maxNumSrcs) then ! realloc oldSrcs => link % srcs numOldSrcs = link % numSrcs link % maxNumSrcs = link % maxNumSrcs + DMM_SL_BLOCK_SIZE allocate(link % srcs(link % maxNumSrcs), stat=status) if ( status /= 0 ) then write(dmmLastError, '(2A,I4,A)') 'dmm: readMappingConfigFile: ', & 'Could not allocate ', link % maxNumSrcs, ' source items' else link % srcs(1:numOldSrcs) = oldSrcs(1:numOldSrcs) s = link % numSrcs + 1 deallocate(oldSrcs) endif else s = link % numSrcs + 1 endif ! add source if ( s /= dmmIndexInvalid ) then link % srcs(s) % id = srcId link % srcs(s) % weight = weight link % srcs(s) % index = dmmIndexInvalid link % srcs(s) % demand = mappings(mh) % missingValue(dmmSrc) if (mapperType == dmm_MF2FM_Q .or. mapperType == dmm_MSW2FM_DPV) then link % srcs(s) % x = x link % srcs(s) % y = y endif link % numSrcs = s endif if (hasXY .and. mappings(mh) % useArea) then if (x < mappings(mh) % xLow .or. x > mappings(mh) % xUpp .or. & y < mappings(mh) % yLow .or. y > mappings(mh) % yUpp) then link % logMapping = .false. endif endif endif end do close(lunConfFile) retVal = 0 end function readMappingConfigFile ! ! Read settings file (for specifying logging) ! function readMappingSettingsFile(mh, confFileName) result(retVal) ! return value integer :: retVal ! 0: everything OK ! arguments integer , intent(in) :: mh ! handle to mapping character(Len=*), intent(in) :: confFileName ! mapping config file ! locals character(Len=256) :: settFileName ! mapping settings file integer :: lunSettFile ! Handle to settings file integer :: read_res ! result of read call character(Len=256) :: line ! line in config file logical :: exists ! is there a settings file? logical :: eof ! yet at end of file? integer :: stringPos ! Position of a charachter in a string character(Len=128) :: key, val ! 'key = value'-pair in settings file logical :: lower, upper ! have lower left and upper right corner been set? ! body retVal = -99 ! open the settings file if there is one stringPos = index(confFileName, '.', back=.true.) if (stringPos == 0) then stringPos = len_trim(confFileName) endif settFileName = confFileName(1:stringPos-1)//'.settings' inquire(file=settFileName, exist=exists) if (.not. exists) then retVal = 0 return endif open(newunit=lunSettFile,file=settFileName,status='old', readonly, iostat=read_res) if (read_res /= 0) then retVal = -1 return endif mappings(mh) % dumpmapping = .true. mappings(mh) % dumpFileName = confFileName(1:stringPos-1)//'.log' ! read all settinfs line write(*, '(2A)') 'dmm: Reading settings file ', trim(settFileName) eof = .false. lower = .false. upper = .false. do while (.not. eof) read(lunSettFile, '(a)', iostat=read_res) line if (read_res == 0) then stringPos = index(line, '=') if (stringPos /= 0) then key = trim(line(1:stringPos-1)) val = trim(line(stringPos+1:)) if (key == 'dumpFile') then mappings(mh) % dumpFileName = val endif if (key == 'tStart') then read(val, *) mappings(mh) % tStart endif if (key == 'tEnd') then read(val, *) mappings(mh) % tEnd endif if (key == 'llCorner') then stringPos = index(val, ',') key = trim(val(1:stringPos-1)) val = trim(val(stringPos+1:)) read(key, *) mappings(mh) % xLow read(val, *) mappings(mh) % yLow lower = .true. endif if (key == 'urCorner') then stringPos = index(val, ',') key = trim(val(1:stringPos-1)) val = trim(val(stringPos+1:)) read(key, *) mappings(mh) % xUpp read(val, *) mappings(mh) % yUpp upper = .true. endif endif else eof = .true. endif end do close(lunSettFile) mappings(mh) % useArea = lower .and. upper retVal = 0 end function readMappingSettingsFile ! ! Add a link to a mapping ! function dmmFindOrAddLink(mh, tgtId, nToOne) result(link) ! return value integer :: link ! handle to new link ! arguments integer, intent(in) :: mh ! connection group integer, intent(in) :: tgtId ! id at target side logical :: nToOne ! many to one mapping, 'find' needed ! locals type(t_dmm_src_tgt) , & pointer, dimension(:) :: oldLinks ! old links (before re-alloc) integer :: numOldLinks ! #old links (before re-alloc) integer :: status ! allocation status flag ! body status = 0 link = dmmIndexInvalid if ( nToOne) then ! find link with same target id if (mappings(mh) % index_lookup(tgtId) > 0) then link = mappings(mh) % index_lookup(tgtId) return endif endif if ( link /= dmmIndexInvalid ) then ! write(*,*) 'ERROR' endif if ( link == dmmIndexInvalid ) then ! not found, add link ! first re-allocate if necessary if ( mappings(mh) % numLinks == 0 ) then allocate(mappings(mh) % links(DMM_LK_BLOCK_SIZE), stat=status) if ( status /= 0 ) then write(dmmLastError, '(A)') 'dmm: dmmFindOrAddLink: Could not allocate initial block of links' else mappings(mh) % maxNumLinks = DMM_LK_BLOCK_SIZE link = 1 endif else if ( mappings(mh) % maxNumLinks == mappings(mh) % numLinks ) then oldLinks => mappings(mh) % links numOldLinks = mappings(mh) % numLinks mappings(mh) % maxNumLinks = mappings(mh) % maxNumLinks + DMM_LK_BLOCK_SIZE allocate(mappings(mh) % links(mappings(mh) % maxNumLinks), stat=status) if ( status /= 0 ) then write(dmmLastError, '(2A,I4,A)') 'dmm: dmmFindOrAddLink: ', & 'Could not allocate ', mappings(mh) % maxNumLinks, 'items' else mappings(mh) % links(1:numOldLinks) = oldLinks(1:numOldLinks) link = numOldLinks+1 deallocate(oldLinks) endif else link = mappings(mh) % numLinks + 1 endif endif if (link /= dmmIndexInvalid) then mappings(mh) % numLinks = link mappings(mh) % links(mappings(mh) % numLinks) % tgt % id = tgtId mappings(mh) % links(mappings(mh) % numLinks) % subLink => NULL() mappings(mh) % links(mappings(mh) % numLinks) % tgt % index = dmmIndexInvalid mappings(mh) % links(mappings(mh) % numLinks) % tgt % x = -777.7 mappings(mh) % links(mappings(mh) % numLinks) % tgt % y = -777.7 mappings(mh) % links(mappings(mh) % numLinks) % numSrcs = 0 mappings(mh) % links(mappings(mh) % numLinks) % logMapping = .true. mappings(mh) % index_lookup(tgtId) = link endif end function dmmFindOrAddLink function dmmDetermineVarName(mapperType) result(varName) ! result character(len=DMM_ID_LEN) :: varName ! 'UNKNOWN VAR' => error ! arguments integer , intent(in) :: mapperType ! type of mapper select case (mapperType) case (dmm_FM2MSW_WL) varName = 'Water level 2D' case (dmm_MSW2FM_DPV) varName = 'Delta Ponding Volume' case (dmm_MSW2FM_SPR) varName = 'Sprinkling Demand' case (dmm_FM2MF_WL) varName = 'Water level 1D' case (dmm_MF2FM_Q) varName = 'Flux' case default varName = 'UNKNOWN VAR' end select end function dmmDetermineVarName end module mfmsdfm_mapping