!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !> 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: ! Note 2: dmm_MF2FM_DRN, dmm_MF2FM_RIV2 and dmm_MSW2FM_RUN have to be ! numbered as 1 2 3, because they are used as index for sub links, and ! because the first two are use in a small loop integer, parameter :: dmm_FM2MSW_WL = 5 ! water level D-Flow FM cells to MetaSWAP integer, parameter :: dmm_MSW2FM_DPV = 6 ! Delta Ponding Volume MetaSWAP to D-Flow FM cells integer, parameter :: dmm_FM2MF_WL = 7 ! water level D-Flow FM 1D nodes to iMODLOW integer, parameter :: dmm_MF2FM_Q = 8 ! iMODFLOW RIV to D-Flow FM 1D nodes integer, parameter :: dmm_MF2FM_DRN = 1 ! iMODFLOW-DRN drainage to D-Flow FM nodes (sub mapper) integer, parameter :: dmm_MF2FM_RIV2 = 2 ! iMODFLOW -> D-Flow FM 1D, drainage from non-coupled rivers (sub mapper) integer, parameter :: dmm_MSW2FM_RUN = 3 ! MetaSwap runnoff to D-Flow FM nodes (sub mapper) integer, parameter :: dmm_MSW2FM_SPR = 4 ! Sprinkling Demand MetaSWAP to D-Flow FM 1D nodes (sub mapper) integer, parameter :: dmm_NUM_SUB_TYPES = 4 ! number of sub mapper types (and thus sub links) 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 at target side ! type t_dmm_target integer :: id ! integer identifier (default) integer :: index ! index of target point in target array double precision :: targetDemand ! summed demand towards target (summed source demands) double precision :: totalRealized ! total realized at target side double precision :: x ! x-coord at target side double precision :: y ! y-coord at target side end type t_dmm_target ! ! link info at 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 (can be positive, i.e. providing water, ! or negative (i.e. asking for water) double precision :: realized ! demand at source side (can be positive, i.e. providing water, ! or negative (i.e. asking for water) double precision :: demandFullDay ! total demand at source side per iMODFLOW time step ! (needed for providing correction flux) double precision :: x ! x-coord at source side double precision :: y ! y-coord at source side end type t_dmm_source type t_sublink_ptr type(t_dmm_link), pointer :: link ! pointer to additional link for combining mapping of fluxes towards on FM node integer :: mh ! handle to the mapping that the sublink is part of end type t_sublink_ptr ! ! Link info (source points -> target point). Use for both the top level (MF-RIV -> FM) link and the sublinks ! type t_dmm_link integer :: mapType ! in case of sublink: the mapper type that the link belongs to 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? double precision :: demandPos ! total positive 'demand' per Dtsw-step at source side (i.e. providing water) double precision :: demandNeg ! total negative demand per Dtsw-step at source side (i.e. asking for water) ! Only used in 'main link' (iMODFLOW RIV-element fluxes): type(t_sublink_ptr), & dimension(dmm_NUM_SUB_TYPES) :: subLinks ! additional links, for combining mapping of fluxes towards one FM node ! (for combining modflow fluxes with metaswap sprinkling demands and modflow drainage) double precision :: totalDemand ! total demand (sum of target demands of main links and sub links) double precision :: totalRealized ! total realized (sum of demands of main links and sub links) end type 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 ! type of mapper character(Len=DMM_ID_LEN) :: varName logical :: hasXY ! TODO: remove (in new version, all configs have X,Y) type(t_dmm_link), 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 ! global variables integer, parameter :: DMM_MAX_MAPPINGS = 20 type(t_dmm_mapping), & dimension(DMM_MAX_MAPPINGS), & target :: mappings integer :: numMappings character(Len=1024) :: dmmLastError double precision :: dmmDtsw = 1.0D0 ! dT (time step size) MetaSWAP (can set) double precision :: dmmDtMf = 1.0D0 ! dT (time step size) iMODFLOW (todo: make settable) integer :: dmmMfStep = 0 ! current iMODFLOW time step (for logging) integer :: dmmMswStep = 0 ! current iMODFLOW time step (for logging) ! TODO: panic shortcut (MR), make this smart/efficient integer, parameter :: LOOKUP_BLOCKSIZE = 200000 double precision, parameter :: numSecPerDay = 86400.0D+0 ! # seconds in a Day integer :: lunFluxes1D = -1 ! log file handle for combined 1D fluxes integer :: lunTotCorr = -1 ! log file handle for total correcton towards iMODFLOW-rivs 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 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! PUBLIC functions, perform mapping for one exchange type !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Mapping: perform mapping. ! This method is used for mapping water levels (1D and 1D) ! and for mapping delta ponding values ! ! 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 type(t_dmm_link), pointer & :: pLink ! pointer to link ! body targetValues = mappings(mh) % missingValue(dmmTgt) do link = 1, mappings(mh) % numLinks pLink => mappings(mh) % links(link) indxTgt =pLink % tgt % index targetValues(indxTgt) = 0 pLink % demandNeg = 0.0D+0 pLink % demandPos = 0.0D+0 do s = 1,pLink % numSrcs srcVal = sourceValues(pLink % srcs(s) % index) if(.not. dabs(srcVal - mappings(mh) % missingValue(dmmSrc))<1.0d-9) then srcVal = srcVal * pLink % srcs(s) % weight targetValues(indxTgt) = targetValues(indxTgt) + srcVal plink % srcs(s) % demand = srcVal if (srcVal < 0.0D+0) then pLink % demandNeg = pLink % demandNeg + srcVal else pLink % demandPos = pLink % demandPos + srcVal endif else ! No action (skip missing value) endif end do pLink % tgt % targetDemand = targetValues(indxTgt) if (.not. dabs(targetValues(indxTgt) - mappings(mh) % missingValue(dmmTgt))<1.0d-9) then if(mappings(mh) % mapperType == dmm_MSW2FM_DPV) then ! convert MetaSWAP volume to m3/s targetValues(indxTgt) = targetValues(indxTgt) / (dmmDtsw * numSecPerDay) else ! water level mapping. write logging. if (mappings(mh) % dumpMapping .and. plink % logMapping) then do s = 1,pLink % numSrcs if(mappings(mh) % mapperType == dmm_FM2MF_WL) then write(mappings(mh) % lunDumpFile,'(I4,A,F16.5,A,F16.5,A,F8.4,A,F6.3,A,I6,A,F8.4)') & dmmMfStep, ',', plink % srcs(s) % x, ',', plink % srcs(s) % y, ',', & plink % srcs(s) % demand, ',', plink % srcs(s) % weight, & ',', plink % tgt % id, ',', plink % tgt % targetDemand else if(mappings(mh) % mapperType == dmm_FM2MSW_WL) then write(mappings(mh) % lunDumpFile,'(I4,A,I2,A,F16.5,A,F16.5,A,F8.4,A,F6.3,A,I6,A,F8.4)') & dmmMfStep, ',', dmmMswStep, ',', plink % srcs(s) % x, ',', plink % srcs(s) % y, ',', & plink % srcs(s) % demand, ',', plink % srcs(s) % weight, & ',', plink % tgt % id, ',', plink % tgt % targetDemand else stop 'unknown mapper type in dmmMapValues' endif enddo call flush(mappings(mh) % lunDumpFile) endif endif endif end do end subroutine dmmMapValues ! ! Mapping: perform mapping back (realized values) ! This method is used for mapping realized delta ponding values back. ! It distributes a deficit according to the reatio of the requesters. ! 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 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 :: targetDemandValue ! total demanded value towards target double precision :: epsilon = 1.0D-10 ! epsilon for comparing demanded and realized flux type(t_dmm_link), pointer :: pLink ! pointer to link double precision :: shortage ! non realized part of the total demand on a link double precision :: notRealizedFraction ! non realized fraction of negative demand on a link ! body realizedSourceValues = mappings(mh) % missingValue(dmmSrc) do link = 1, mappings(mh) % numLinks plink => mappings(mh) % links(link) indxTgt = plink % tgt % index totalRealizedValue = valuesFromTargetComponent(indxTgt) ! convert m3/s to volume totalRealizedValue = totalRealizedValue * dmmDtsw * numSecPerDay targetDemandValue = plink % tgt % targetDemand ! TODO: combine if's if (targetDemandValue < 0 .and. totalRealizedValue > targetDemandValue + epsilon) then ! distribute shortage shortage = targetDemandValue - totalRealizedValue notRealizedFraction = shortage / pLink % demandNeg do s = 1, plink % numSrcs indxSrc = plink % srcs(s) % index if (pLink % srcs(s) % demand < 0) then realizedSourceValues(indxSrc) = pLink % srcs(s) % demand * (1.0D+0 - notRealizedFraction) else realizedSourceValues(indxSrc) = pLink % srcs(s) % demand endif if (mappings(mh) % dumpMapping .and. plink % logMapping) then write(mappings(mh) % lunDumpFile,'(I4,A,I2,A,I6,A,F16.5,A,F16.5,A,F12.4,A,F12.4,A,F16.10,A,F16.10)') & dmmMfStep, ',', dmmMswStep, ',', plink % srcs(s) % id, ',', plink % tgt % x, ',', plink % tgt % y, & ',', plink % srcs(s) % demand, ',', realizedSourceValues(plink % srcs(s) % index), & ',', plink % tgt % targetDemand / (dmmDtsw * numSecPerDay), ',', valuesFromTargetComponent(indxTgt) endif end do else ! no shortage, all delta ponding volumes realized do s = 1, plink % numSrcs indxSrc = plink % srcs(s) % index realizedSourceValues(indxSrc) = plink % srcs(s) % demand if (mappings(mh) % dumpMapping .and. plink % logMapping) then write(mappings(mh) % lunDumpFile,'(I4,A,I2,A,I6,A,F16.5,A,F16.5,A,F12.4,A,F12.4,A,F16.10,A,F16.10)') & dmmMfStep, ',', dmmMswStep, ',', plink % srcs(s) % id, ',', plink % tgt % x, ',', plink % tgt % y, & ',', plink % srcs(s) % demand, ',', realizedSourceValues(plink % srcs(s) % index), & ',', plink % tgt % targetDemand / (dmmDtsw * numSecPerDay), ',', valuesFromTargetComponent(indxTgt) endif end do endif if (mappings(mh) % dumpMapping .and. plink % logMapping) call flush(mappings(mh) % lunDumpFile) end do end subroutine dmmMapValuesBack !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! PUBLIC functions, perform mapping for combined fluxes !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! This function combines the four different flux mappers. ! The mhFlux1dMapper (iModflow -> D-Flow FM, infiltration/seepage) is see as the main mapper ! For each of its links, the links of the other mappers are added as sub links ! function dmmCombineMappers(mhFlux1dMapper , mapperSprinkling , & mapperMfDrnToDfm1D , mapperMswRunoffToDfm1D, & mapperMfRivDrnToDfm1D ) result(success) ! return value logical :: success ! .false.: error ! arguments integer, intent(in) :: mhFlux1dMapper ! handle to mapper iModflow -> D-Flow FM (infiltration/seepage) integer, intent(in) :: mapperSprinkling ! handle to mapper MetaSWAP -> D-Flow FM (sprinkling demand) integer, intent(in) :: mapperMfDrnToDfm1D ! handle to mapper iModflow -> D-Flow FM (DRN) integer, intent(in) :: mapperMswRunoffToDfm1D ! handle to mapper MetaSWAP -> D-Flow FM (Runoff) integer, intent(in) :: mapperMfRivDrnToDfm1D ! handle to mapper MetaSWAP -> D-Flow FM (Drainage from small rivers) ! locals integer :: m, link, sLink ! loop counter for links and source items integer :: tgtId ! target id type(t_dmm_mapping), pointer :: flux1DMapper ! mapper for fluxes between iModflow riv elements and FM 1D nodes type(t_dmm_mapping), pointer :: addMapper ! mapper for additional fluxes ! body flux1DMapper => mappings(mhFlux1dMapper) do m = 1,4 addMapper => NULL() select case(m) case(1) if(mapperSprinkling > 0) then addMapper => mappings(mapperSprinkling) endif case(2) if(mapperMfDrnToDfm1D > 0) then addMapper => mappings(mapperMfDrnToDfm1D) endif case(3) if(mapperMswRunoffToDfm1D > 0) then addMapper => mappings(mapperMswRunoffToDfm1D) endif case(4) if(mapperMfRivDrnToDfm1D > 0) then addMapper => mappings(mapperMfRivDrnToDfm1D) endif end select if (associated(addMapper)) then do link = 1, flux1DMapper % numLinks tgtId = flux1DMapper % links(link) % tgt % id do sLink = 1, addMapper % numLinks if (addMapper % links(sLink) % tgt % id == tgtId) then flux1DMapper % links(link) % subLinks(addMapper % mapperType) % link => addMapper % links(sLink) endif enddo end do endif end do success = .true. end function dmmCombineMappers ! ! This function sets the requested volumes. Called by the driver for each flux (= volume per time step) mapper ! subroutine dmmCombMapperSetVolumes(mh, sourceValues) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(in) :: sourceValues ! values at source side type(t_dmm_link), pointer :: pLink ! pointer to link ! locals integer :: link, s ! loop counter for links and source items double precision :: srcVal ! value at source side if (mappings(mh) % mapperType == dmm_MF2FM_Q) then ! method was called before dtsw-loop. reset msw time step counter dmmMswStep = 0 else if (mappings(mh) % mapperType == dmm_MSW2FM_SPR) then ! First 1D D-Flow FM <-> MetaSWAP call dtsw-loop. increment counter dmmMswStep = dmmMswStep + 1 endif do link = 1, mappings(mh) % numLinks pLink => mappings(mh) % links(link) pLink % tgt % targetDemand = 0.D+0 pLink % demandPos = 0.D+0 pLink % demandNeg = 0.D+0 do s = 1, pLink % numSrcs srcVal = sourceValues(pLink % srcs(s) % index) if(.not. dabs(srcVal - mappings(mh) % missingValue(dmmSrc))<1.0d-9) then pLink % srcs(s) % demand = srcVal * pLink % srcs(s) % weight ! MetaSWAP provides outgoing delta ponding volume and runoff (TODO: check) as positive ! iMODFLOW provides outgoing volumes as negative (TODO: check) ! So for iMODFLOW out towards D-FM in, swap sign. ! Also, for those mappers, devide the iMODFLOW full day demand ! by the number Dtsw steps per day pLink % srcs(s) % demandFullDay = pLink % srcs(s) % demand if (mappings(mh) % mapperType == dmm_MF2FM_Q .or. & mappings(mh) % mapperType == dmm_MF2FM_DRN .or. & mappings(mh) % mapperType == dmm_MF2FM_RIV2 ) then pLink % srcs(s) % demand = -pLink % srcs(s) % demand / nint(dmmDtMf/dmmDtsw) endif ! In case of runoff from MSW directly to FM, no negative MetaSWAP delta ponding allowed. if (mappings(mh) % mapperType == dmm_MSW2FM_RUN .and. & pLink % srcs(s) % demand < -1.0D-7 ) then stop 'Negative runoff from MSW' endif if (pLink % srcs(s) % demand < 0.0D+0) then pLink % demandNeg = pLink % demandNeg + pLink % srcs(s) % demand else pLink % demandPos = pLink % demandPos + pLink % srcs(s) % demand endif pLink % tgt % targetDemand = & pLink % tgt % targetDemand + pLink % srcs(s) % demand else ! No action (skip missing value) endif pLink % srcs(s) % realized = 0.0D+0 end do end do end subroutine dmmCombMapperSetVolumes ! ! This method is called by the driver to provided all summed volumes to FM, as a flux (volume over the dtsw step) per 1D grid node ! subroutine dmmCombMapperGetFluxes(mh, targetValues) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(out) :: targetValues ! computed values for target side ! locals integer :: indxTgt ! index at target side, index integer :: link, sLink ! loop counter for links and subLinks type(t_dmm_link), pointer :: pLink, pSubLink ! pointer to link and sub link ! body if (mappings(mh) % mapperType /= dmm_MF2FM_Q) then stop 'wrong caller for dmmCombMapperGetFluxes' endif targetValues = mappings(mh) % missingValue(dmmTgt) do link = 1, mappings(mh) % numLinks pLink => mappings(mh) % links(link) ! First add the requested or delivered volume of the main link (iMODFLOW RIV -> FM 1D) pLink % totalDemand = pLink % tgt % targetDemand ! Then add the the requested or delivered volumes of the sub links ! Devide modflow demand by number of metaswap steps in a modflow step do sLink = 1, dmm_NUM_SUB_TYPES if (associated(pLink % subLinks(sLink) % link)) then pSubLink => pLink % subLinks(sLink) % link pLink % totalDemand = pLink % totalDemand + pSubLink % tgt % targetDemand endif end do indxTgt = pLink % tgt % index targetValues(indxTgt) = pLink % totalDemand / (dmmDtsw * numSecPerDay) end do end subroutine dmmCombMapperGetFluxes ! ! This function is called by the driver to set the fluxees as they were realized by D-Flow FM. ! It distributes the realized values over the various requesters ! subroutine dmmCombMapperSetRealizedFluxes(mh, realizedFluxes) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(in) :: realizedFluxes ! computed realized fluxes ! locals integer :: link, sLink, s ! loop counter for links, subLinks and sources integer :: indxTgt ! index at target side, index double precision :: shortage ! non realized part of the total demand on a link double precision :: notRealizedFraction ! non realized fraction of negative demand on a link double precision :: linkDemand ! demand per link double precision :: srcDemand ! demand per source double precision :: mfDemand ! demand from iMODFLOW main link double precision :: mswDemand ! demand from msw sub link double precision :: mfRealized ! realized for iMODFLOW main link (for debugging only) double precision :: mswRealized ! realized for msw sub link (for debugging only) type(t_dmm_link), pointer :: pLink ! pointer to link character(len=1024) :: logLine ! composed line to write to log file ! body do link = 1, mappings(mh) % numLinks pLink => mappings(mh) % links(link) indxTgt = pLink % tgt % index pLink % totalRealized = realizedFluxes(indxTgt) * dmmDtsw * numSecPerDay if(mappings(mh) % dumpMapping .and. plink % logMapping) then write(logLine,'(I4,A,I2,A,F16.5,A,F16.5,A,F16.10,A,F16.10)') dmmMfStep, ',', dmmMswStep, ',', & pLink % tgt % x, ',', pLink % tgt % y, ',', & pLink % totalDemand / (dmmDtsw * numSecPerDay) , ',', realizedFluxes(indxTgt) endif ! First substract the always realized drain- or runoff- fluxes do sLink = dmm_MF2FM_DRN, dmm_MSW2FM_RUN if (associated(pLink % subLinks(sLink) % link )) then call dmmLogFixedFluxLink(pLink, sLink) pLink % totalRealized = pLink % totalRealized - pLink % subLinks(sLink) % link % demandPos if(mappings(mh) % dumpMapping .and. plink % logMapping) then write(logLine,'(A,A,F12.4)') trim(logLine), ',', pLink % subLinks(sLink) % link % demandPos endif else if(mappings(mh) % dumpMapping .and. plink % logMapping) then write(logLine,'(A,A,F12.4)') trim(logLine), ',', 0.D+0 endif endif enddo ! Now check if both iMODFLOW and MetaSWAP demand could be realized, which is the case if: ! - the resulting realization for these fluxes is positive ! or ! - the resulting realization is negative, and the negative value is not greater than ! (or in absolute sense not smaller than) the summed demand of iMODFLOW and MetaSwap mfDemand = pLink % demandPos + pLink % demandNeg mswDemand = 0.D+0 linkDemand = mfDemand if (associated(pLink % subLinks(dmm_MSW2FM_SPR) % link)) then mswDemand = pLink % subLinks(dmm_MSW2FM_SPR) % link % demandPos + pLink % subLinks(dmm_MSW2FM_SPR) % link % demandNeg linkDemand = linkDemand + mswDemand endif mswRealized = 0.0D+0 ! for debugging mfRealized = 0.0D+0 ! " " if (pLink % totalRealized > 0 .or. (pLink % totalRealized <= 0 .and. pLink % totalRealized <= linkDemand) ) then ! All has been realized ! MetaSWAP sprinkling -> D-FM sub link if (associated(pLink % subLinks(dmm_MSW2FM_SPR) % link)) then do s = 1, pLink % subLinks(dmm_MSW2FM_SPR) % link % numSrcs pLink % subLinks(dmm_MSW2FM_SPR) % link % srcs(s) % realized = pLink % subLinks(dmm_MSW2FM_SPR) % link % srcs(s) % demand end do endif mswRealized = mswDemand ! iMODFLOW RIV -> D-FM link (i.e. the present 'main' link). ! Add the realized demand for this Dtsw step to the full day realized value do s = 1, pLink % numSrcs pLink % srcs(s) % realized = pLink % srcs(s) % realized + pLink % srcs(s) % demand end do mfRealized = mfRealized + mfDemand else if (pLink % totalRealized < mfDemand) then ! iMODFLOW demand could be realized. ! Add the realized demand for this Dtsw step to the full day realized value do s = 1, pLink % numSrcs pLink % srcs(s) % realized = pLink % srcs(s) % realized + pLink % srcs(s) % demand end do pLink % totalRealized = pLink % totalRealized - mfDemand mfRealized = mfRealized + mfDemand ! Provide remainder to MetaSWAP ! Distribute negative shortage over the negative demand ! TODO: simplify (no positive sprinkling!) if (associated(pLink % subLinks(dmm_MSW2FM_SPR) % link)) then shortage = mswDemand - pLink % totalRealized notRealizedFraction = shortage / pLink % subLinks(dmm_MSW2FM_SPR) % link % demandNeg do s = 1, pLink % subLinks(dmm_MSW2FM_SPR) % link % numSrcs srcDemand = pLink % subLinks(dmm_MSW2FM_SPR) % link % srcs(s) % demand if (srcDemand < 0) then pLink % subLinks(dmm_MSW2FM_SPR) % link % srcs(s) % realized = srcDemand * (1.0D+0 - notRealizedFraction) else pLink % subLinks(dmm_MSW2FM_SPR) % link % srcs(s) % realized = srcDemand endif mswRealized = mswRealized + pLink % subLinks(dmm_MSW2FM_SPR) % link % srcs(s) % realized end do endif else ! iMODFLOW demand could not be met. Distribute shortages. ! Add the realized fraction for this Dtsw step to the full day realized value. ! Set MetaSwap allocation to zere. shortage = mfDemand - pLink % totalRealized notRealizedFraction = shortage / pLink % demandNeg do s = 1, pLink % numSrcs srcDemand = pLink % srcs(s) % demand if (srcDemand < 0) then pLink % srcs(s) % realized = pLink % srcs(s) % realized + pLink % srcs(s) % demand * (1.0D+0 - notRealizedFraction) else pLink % srcs(s) % realized = pLink % srcs(s) % realized + pLink % srcs(s) % demand endif mfRealized = mfRealized + pLink % srcs(s) % realized end do if (associated(pLink % subLinks(dmm_MSW2FM_SPR) % link)) then do s = 1, pLink % subLinks(dmm_MSW2FM_SPR) % link % numSrcs pLink % subLinks(dmm_MSW2FM_SPR) % link % srcs(s) % realized = 0.0D+0 end do endif endif if(mappings(mh) % dumpMapping .and. plink % logMapping) then write(logLine,'(A,A,F12.4,A,F12.4,A,F12.4,A,F12.4)') trim(logLine), & ',', mfDemand, ',', mfRealized, ',', mswDemand, ',', mswRealized write(lunFluxes1D,'(A)') trim(logLine) call flush(lunFluxes1D) endif end do end subroutine dmmCombMapperSetRealizedFluxes subroutine dmmCombMapperGetRealizedFluxes(mh, realizedVolumes) ! arguments integer , intent(in) :: mh ! handle to mapping double precision, dimension(:), intent(out) :: realizedVolumes ! computed and distributed realized volumes ! locals integer :: link, s ! loop counter for links and sources integer :: indxSrc ! Index at source site type(t_dmm_link), pointer :: pLink ! pointer to link and sub link double precision :: totalCorrection ! all over correction flux to iMODFLOW RIV totalCorrection = 0.0D+0 ! body realizedVolumes = mappings(mh) % missingValue(dmmTgt) do link = 1, mappings(mh) % numLinks pLink => mappings(mh) % links(link) do s = 1, pLink % numSrcs indxSrc = pLink % srcs(s) % index realizedVolumes(indxSrc) = pLink % srcs(s) % realized if (mappings(mh) % mapperType == dmm_MF2FM_Q) then ! Swap sign, and provide correction instead of realized volumes realizedVolumes(indxSrc) = -realizedVolumes(indxSrc) realizedVolumes(indxSrc) = realizedVolumes(indxSrc) - pLink % srcs(s) % demandFullDay totalCorrection = totalCorrection + realizedVolumes(indxSrc) endif if (mappings(mh) % dumpMapping .and. pLink % logMapping) then if(mappings(mh) % mapperType == dmm_MSW2FM_SPR) then write(mappings(mh) % lunDumpFile,'(I4,A,I2,A,F16.5,A,F16.5,A,I6,A,F12.4,A,F12.4)') & dmmMfStep, ',', dmmMswStep, ',', plink % tgt % x, ',', plink % tgt % y, ',', & plink % srcs(s) % id, ',', plink % srcs(s) % demand, ',', realizedVolumes(indxSrc) else if(mappings(mh) % mapperType == dmm_MF2FM_Q) then write(mappings(mh) % lunDumpFile,'(I4,A,F16.5,A,F16.5,A,I6,A,F12.4,A,F12.4)') & dmmMfStep, ',', plink % tgt % x, ',', plink % tgt % y, ',', & plink % srcs(s) % id, ',', plink % srcs(s) % demand, ',', realizedVolumes(indxSrc) else stop 'unknown mapper type in dmmCombMapperGetRealizedFluxes' endif endif end do if (mappings(mh) % dumpMapping .and. pLink % logMapping) call flush(mappings(mh) % lunDumpFile) end do if (mappings(mh) % mapperType == dmm_MF2FM_Q) then write(lunTotCorr,'(I4,A,F18.8)') dmmMfStep, ',', totalCorrection call flush(lunTotCorr) endif end subroutine dmmCombMapperGetRealizedFluxes !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! 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_link), 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 retVal = -9 write(*, '(2A)') 'dmm: Reading Config file ', trim(confFileName) ! open file open(newunit=lunConfFile,file=confFileName,form='formatted', status='old', readonly, iostat=read_res) if (read_res/=0) then write(1920,*) 'readMappingConfigFile: file ', trim(confFileName), ' could not be opened' return endif ! read first line (TARGET, SOURCE, etc.) read(lunConfFile, '(a)', iostat=read_res) line if (read_res/=0) then write(1920,*) 'readMappingConfigFile: file ', trim(confFileName), ' could not be read' return endif ! read all configuration line eof = .false. do while (.not. eof) x = -777.7 y = -777.7 srcId = -1 tgtId = -1 weight= 1.0D+0 if (hasXY) then success = .false. read(lunConfFile,'(A)', iostat=read_res) line if (read_res /= 0) then eof = .true. exit endif if (mapperType == dmm_FM2MSW_WL .or. mapperType == dmm_FM2MF_WL) then ! file contains SVAT/MF-RIV ID, (FM)X, (FM)Y, WEIGHT ! FM-cell is source, svat is target read(line, *, iostat=read_res) tgtId, x, y, weight ! read(line, '(F,F,I,F)', iostat=read_res) x, y, tgtId, weight if (read_res /= 0) then write(1920,*) 'ERROR reading SVAT/MF-RIV ID, (FM)X, (FM)Y, WEIGHT from ', trim(confFileName) else if (mapperType == dmm_FM2MSW_WL) then 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 else success = dflowfm_MapCoordinateTo1DCellId(x, y, srcId) if (.not. success) then write(1920,*) 'ERROR getting source gridpoint-id from FM, x:', x, ', y:', y, ', tgtId:', tgtId cycle endif endif endif else if (mapperType == dmm_MSW2FM_DPV .or. mapperType == dmm_MF2FM_Q .or. & mapperType == dmm_MF2FM_DRN .or. mapperType == dmm_MF2FM_RIV2 .or. & mapperType == dmm_MSW2FM_RUN .or. mapperType == dmm_MSW2FM_SPR ) then ! file contains (FM)X, (FM)Y, SVAT/MF-(RIV/DRN) ID ! svat is source, FM-cell is target read(line, *, iostat=read_res) x, y, srcId if (read_res /= 0) then write(1920,*) 'ERROR reading(FM)X, (FM)Y, SVAT/MF-(RIV/DRN) ID from ', trim(confFileName) else if (mapperType == dmm_MSW2FM_DPV) then 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 else success = dflowfm_MapCoordinateTo1DCellId(x, y, tgtId) if (.not. success) then write(1920,*) 'ERROR getting target gridpoint-id from FM, x:', x, ', y:', y, ', srcId:', srcId cycle endif endif endif else stop 'unknown mapper type in readMappingConfigFile' 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) if (lk>0) then link => mappings(mh) % links(lk) link % tgt % id = tgtId if (mapperType == dmm_MSW2FM_DPV .or. mapperType == dmm_MF2FM_Q .or. & mapperType == dmm_MF2FM_DRN .or. mapperType == dmm_MF2FM_RIV2 .or. & mapperType == dmm_MSW2FM_RUN .or. mapperType == dmm_MSW2FM_SPR ) 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_FM2MSW_WL .or. mapperType == dmm_FM2MF_WL) 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)//'.dms' 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)//'.csv' ! 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) 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 ! locals type(t_dmm_link) , & pointer, dimension(:) :: oldLinks ! old links (before re-alloc) integer :: numOldLinks ! #old links (before re-alloc) integer :: status ! allocation status flag integer :: sLink ! loop over sublinks type(t_dmm_mapping), pointer :: mapper ! pointer to mapper ! body status = 0 link = dmmIndexInvalid mapper => mappings(mh) ! find link with same target id if (mapper % index_lookup(tgtId) > 0) then link = mapper % index_lookup(tgtId) return endif if ( link /= dmmIndexInvalid ) then ! write(*,*) 'ERROR' endif if ( link == dmmIndexInvalid ) then ! not found, add link ! first re-allocate if necessary if ( mapper % numLinks == 0 ) then allocate(mapper % links(DMM_LK_BLOCK_SIZE), stat=status) if ( status /= 0 ) then write(dmmLastError, '(A)') 'dmm: dmmFindOrAddLink: Could not allocate initial block of links' else mapper % maxNumLinks = DMM_LK_BLOCK_SIZE link = 1 endif else if ( mapper % maxNumLinks == mapper % numLinks ) then oldLinks => mapper % links numOldLinks = mapper % numLinks mapper % maxNumLinks = mapper % maxNumLinks + DMM_LK_BLOCK_SIZE allocate(mapper % links(mapper % maxNumLinks), stat=status) if ( status /= 0 ) then write(dmmLastError, '(2A,I4,A)') 'dmm: dmmFindOrAddLink: ', & 'Could not allocate ', mapper % maxNumLinks, 'items' else mapper % links(1:numOldLinks) = oldLinks(1:numOldLinks) link = numOldLinks+1 deallocate(oldLinks) endif else link = mapper % numLinks + 1 endif endif if (link /= dmmIndexInvalid) then mapper % numLinks = link mapper % links(mapper % numLinks) % tgt % id = tgtId do sLink = 1, dmm_NUM_SUB_TYPES mapper % links(mapper % numLinks) % subLinks(sLink) % link => NULL() enddo mapper % links(mapper % numLinks) % tgt % index = dmmIndexInvalid mapper % links(mapper % numLinks) % tgt % x = -777.7 mapper % links(mapper % numLinks) % tgt % y = -777.7 mapper % links(mapper % numLinks) % numSrcs = 0 mapper % links(mapper % numLinks) % logMapping = .true. mapper % 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, dmm_MF2FM_RIV2, dmm_MF2FM_DRN, dmm_MSW2FM_RUN) varName = 'Volumes' case default varName = 'UNKNOWN VAR' end select end function dmmDetermineVarName subroutine dmmLogFixedFluxLink(pLink, subLinkIndex) ! arguments type(t_dmm_link), pointer :: pLink ! pointer to link in main (MF->FM Q) mapper integer, intent(in) :: subLinkIndex ! index of sublink ! (dmm_MF2FM_DRN, dmm_MF2FM_RIV2 or dmm_MSW2FM_RUN) ! locals integer :: mh ! handle to sub mapper integer :: s ! loop counter for sources type(t_dmm_link), pointer :: pSubLink ! pointer to sub link mh = dmmGetMapperHandle(subLinkIndex) if (.not. (mappings(mh) % dumpMapping .and. plink % logMapping)) then return endif pSubLink => pLink % subLinks(subLinkIndex) % link do s = 1, pSubLink % numSrcs if(subLinkIndex == dmm_MF2FM_DRN .or. subLinkIndex == dmm_MF2FM_RIV2) then ! log iMODFLOW drainage only at last MSW step if (dmmMswStep == nint(dmmDtMf/dmmDtsw)) then write(mappings(mh) % lunDumpFile,'(I4,A,F16.5,A,F16.5,A,I6,A,F12.4)') & dmmMfStep, ',', pSubLink % tgt % x, ',', pSubLink % tgt % y, ',', & pSubLink % srcs(s) % id, ',', pSubLink % srcs(s) % demand endif else if(subLinkIndex == dmm_MSW2FM_RUN) then write(mappings(mh) % lunDumpFile,'(I4,A,I2,A,F16.5,A,F16.5,A,I6,A,F12.4)') & dmmMfStep, ',', dmmMswStep, ',', pSubLink % tgt % x, ',', pSubLink % tgt % y, ',', & pSubLink % srcs(s) % id, ',', pSubLink % srcs(s) % demand endif enddo call flush(mappings(mh) % lunDumpFile) end subroutine dmmLogFixedFluxLink function dmmGetMapperHandle(mapperType) result(mh) ! result integer :: mh ! handle to mapper (index in mappings array) ! arguments integer, intent(in) :: mapperType ! type of mapper ! body integer :: m ! loop counter mh = -1 do m = 1, numMappings if (mappings(m) % mapperType == mapperType) then mh = m endif enddo ! TODO: if == -1 => error end function dmmGetMapperHandle end module mfmsdfm_mapping