!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !> 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 1: The first 5 mappers must be number 1 to 5 ! they are used as index for the sub links in a cluster. ! Note 2: dmm_MF2FM_DRN/dmm_MF2FM_RIV2/dmm_MSW2FM_RUN must have ! this order, because htat are used in a small loop integer, parameter :: dmm_MF2FM_Q = 1 ! iMODFLOW RIV to D-Flow FM 1D nodes integer, parameter :: dmm_MF2FM_DRN = 2 ! iMODFLOW-DRN drainage to D-Flow FM nodes (sub mapper) integer, parameter :: dmm_MF2FM_RIV2 = 3 ! iMODFLOW -> D-Flow FM 1D, drainage from non-coupled rivers (sub mapper) integer, parameter :: dmm_MSW2FM_RUN = 4 ! MetaSwap runnoff to D-Flow FM nodes (sub mapper) integer, parameter :: dmm_MSW2FM_SPR = 5 ! Sprinkling Demand MetaSWAP to D-Flow FM 1D nodes (sub mapper) integer, parameter :: dmm_NUM_SUB_TYPES = 5 ! number of sub mapper types (and thus sub links) integer, parameter :: dmm_FM2MSW_WL = 6 ! water level D-Flow FM cells to MetaSWAP integer, parameter :: dmm_MSW2FM_DPV = 7 ! Delta Ponding Volume MetaSWAP to D-Flow FM cells integer, parameter :: dmm_FM2MF_WL = 8 ! water level D-Flow FM 1D nodes to iMODLOW 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 = 1000 ! size for reallocing links integer, parameter, private :: DMM_SL_BLOCK_SIZE = 5 ! size for reallocing source_items in links integer, parameter, private :: DMM_CL_BLOCK_SIZE = 1000 ! size for reallocing clusters 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 ! realization 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) double precision :: totalDemand ! total demand (sum of target demands the clustered links) double precision :: totalRealized ! total realized (sum of demands of main links and sub links) logical :: connected ! has the link in the sub mapper indeed been connected? end type t_dmm_link ! ! mapper (N sources to 1 target) ! 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 ! name of the variable that is exchanged by this mapper 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) integer :: numConnected ! #links in a mapper that have been connected to a cluster (should be all) 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(:) :: lk_lookup ! mapping from inter id to index, to avoid inefficient search integer :: lk_lookup_length ! current length of the lookup table end type t_dmm_mapping ! ! Cluster of mapper links that all have one FM 1D-node as target. ! The cluster contains the total of 1D-demands towards that FM 1D-node for all 1D-flux types: ! iMODFLOW: . seepage/drainage in river package elements that are interactively connected with an FM 1D-node ! . drainages from drainage package (one way flux) ! . river fluxes from river package elements that are not interactively connected (one way flux) ! MetaSWAP: . sprinkling demand ! . 2D runoff (in case FM's 2D-part is de-activated) ! type t_dmm_cluster integer :: tgtIndex ! (FM-) index of target point in target array double precision :: tgtX, tgtY ! (FM-) coordinates of target point in target array (for logging) type(t_sublink_ptr), & dimension(dmm_NUM_SUB_TYPES) :: subLinks ! clustered links, for combining the mapping of fluxes towards one FM 1D-node ! (modflow river fluxes, metaswap sprinkling demands, modflow drainage, ! modflow river2 fluxes, metaswap runoff) double precision :: totalDemand ! total demand towards FM (sum of target demands of the sub links) double precision :: totalRealized ! total realized by FM 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) logical :: logMapping ! do log the 1D-fluxes for this cluster? end type t_dmm_cluster ! 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) double precision, parameter :: numSecPerDay = 86400.0D+0 ! # seconds in a Day integer :: lunFluxes1D = -1 ! log file handle for combined 1D fluxes integer :: lunSummed = -1 ! log file handle for sum of 1D fluxes per timestep integer :: lunMswRunoff = -1 ! log file handle for detailed logging of MSW's runoff ! ! Clusters of mapper links that all have one FM 1D-node as target ! type(t_dmm_cluster), & pointer, dimension(:) :: clusters ! the clusters integer :: maxNumClusters = DMM_CL_BLOCK_SIZE ! maximum number of clusters (realloc if needed) integer :: numClusters = 0 ! actual number of clusters integer, & pointer, dimension(:) :: cl_lookup ! lookup table: (integer) identifier to index in clusters-array integer :: cl_lookup_length = DMM_CL_BLOCK_SIZE ! maximum number of clusters (realloc if needed) double precision :: clusterMissingValue = -99999.0D0 ! missing value identifier at target side (i.e. FM) ! (will be set to FM's miss.val. in MF->FM-1D mapper) logical :: dumpClusterFluxes = .false. ! do dump 1D-fluxes? ! (will be set to: do dump MF->FM-1D mapper?) ! FOR LOGGING: character(len=128) :: demandLineFm ! D-Flow Fm demand part of log line character(len=128) :: demandLine(dmm_NUM_SUB_TYPES) ! Mf/Drn/Riv2/Run/Spr demand parts of log line character(len=128) :: realizedLineFm ! D-Flow Fm realized part of log line character(len=128) :: realizedLine(dmm_NUM_SUB_TYPES) ! Mf/Drn/Riv2/Run/Spr realized part of log line character(len=128) :: connectedLine ! MSW connected / not connected volumes line logical, & dimension(dmm_NUM_SUB_TYPES) :: subTypeUsed ! sub link type used in one of the clusters? 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) % maxNumLinks = 0 mappings(m) % numConnected = 0 mappings(m) % dumpMapping = .false. mappings(m) % dumpFileName = ' ' mappings(m) % lundumpFile = -1 mappings(m) % tStart = 0 mappings(m) % tEnd = 1E+6 mappings(m) % useArea = .false. mappings(m) % links => NULL() do st = 1, DMMSRCTGT mappings(m) % missingValue(st) = -99999.0D0 mappings(m) % componentID(st) = ' ' enddo allocate(mappings(m) % lk_lookup(DMM_LK_BLOCK_SIZE)) mappings(m) % lk_lookup = -1 enddo allocate(clusters(DMM_CL_BLOCK_SIZE), cl_lookup(DMM_CL_BLOCK_SIZE)) cl_lookup = -1 dmmLastError = '' demandLineFm = '' demandLine = ' 0.00000, -0.00000, 0.00000, 0.00000,' realizedLineFm = '' realizedLine = ' 0.00000, -0.00000, 0.00000, 0.00000,' connectedLine = ' 0.00000, -0.00000,' 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(mapperFlux1D , mapperSprinkling , & mapperMfDrnToDfm1D , mapperMswRunoffToDfm1D, & mapperMfRivDrnToDfm1D ) result(success) ! return value logical :: success ! .false.: error ! arguments integer, intent(in) :: mapperFlux1D ! 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, sLink ! loop counter for mappers and the sub links type(t_dmm_mapping), pointer :: addMapper ! mapper for additional fluxes integer :: clusterIndex ! cluster index for exc type(t_dmm_cluster), pointer :: pCluster ! pointer to a cluster ! body do m = 1,dmm_NUM_SUB_TYPES addMapper => NULL() select case(m) case(dmm_MF2FM_Q) if(mapperFlux1D > 0) then addMapper => mappings(mapperFlux1D) endif case(dmm_MF2FM_DRN) if(mapperMfDrnToDfm1D > 0) then addMapper => mappings(mapperMfDrnToDfm1D) endif case(dmm_MF2FM_RIV2) if(mapperMfRivDrnToDfm1D > 0) then addMapper => mappings(mapperMfRivDrnToDfm1D) endif case(dmm_MSW2FM_RUN) if(mapperMswRunoffToDfm1D > 0) then addMapper => mappings(mapperMswRunoffToDfm1D) endif case(dmm_MSW2FM_SPR) if(mapperSprinkling > 0) then addMapper => mappings(mapperSprinkling) endif end select if (associated(addMapper)) then addMapper % numConnected = 0 do sLink = 1, addMapper % numLinks addMapper % links(sLink) % connected = .false. clusterIndex = dmmFindOrAddCluster(addMapper % links(sLink) % tgt % id) if (clusterIndex == -1) then ! todo error else pCluster => clusters(clusterIndex) pCluster % subLinks(addMapper % mapperType) % link => addMapper % links(sLink) addMapper % links(sLink) % connected = .true. addMapper % numConnected = addMapper % numConnected + 1 pCluster % tgtIndex = addMapper % links(sLink) % tgt % index pCluster % tgtX = addMapper % links(sLink) % tgt % x pCluster % tgtY = addMapper % links(sLink) % tgt % y pCluster % logMapping = pCluster % logMapping .or.addMapper % links(sLink) % logMapping if (addMapper % mapperType == dmm_MF2FM_Q) then dumpClusterFluxes = addMapper % dumpMapping endif endif enddo endif end do success = .true. end function dmmCombineMappers ! ! This function sets the requested volumes for a mapper that handles 1D fluxes. ! Called by the driver for each flux (= volume per time step) mapper. ! subroutine dmm1DMapperSetVolumes(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 ! For logging: double precision :: totSrcIn ! total volume of all sourceValues when coming in double precision :: totDemand ! total demand over all links and sources double precision :: totPos ! total positive target demand over all links and sources double precision :: totNeg ! total negative target demand over all links and sources 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 ! log total demand for the 1D-fluxes mapper type, en demands per link totSrcIn = 0.0D+0 do s = 1, size(sourceValues) totSrcIn = totSrcIn + sourceValues(s) enddo totDemand = 0.0D+0 totPos = 0.0D+0 totNeg = 0.0D+0 do link = 1, mappings(mh) % numLinks pLink => mappings(mh) % links(link) totPos = totPos + pLink % demandPos totNeg = totNeg + pLink % demandNeg do s = 1, pLink % numSrcs totDemand = totDemand + pLink % srcs(s) % demand end do end do write(demandLine(mappings(mh) % mapperType), '(F12.5,A,F12.5,A,F12.5,A,F12.5,A)') totSrcIn, ',', totDemand, ',', totPos, ',', totNeg, ',' end subroutine dmm1DMapperSetVolumes ! ! 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 dmmClustersGetFluxes(targetValues) ! arguments double precision, dimension(:), intent(out) :: targetValues ! computed values for target side ! locals integer :: indxTgt ! index at target side, index integer :: cl, sLink ! loop counter for cluster and subLinks type(t_dmm_cluster), pointer :: pCluster ! pointer a cluster type(t_dmm_link) , pointer :: pSubLink ! pointer to the sub links of a cluster ! For logging: integer :: t ! loop counter for values in target array double precision :: totDemand ! total demand over all links and sources double precision :: totPos ! total positive demand over all links and sources double precision :: totNeg ! total negative demand over all links and sources double precision :: totTarget ! total volume at all FM target nodes double precision, parameter :: epsilon = 1.0D-10 ! epsilon for comparing demanded and realized flux ! body totDemand = 0.0D+0 totPos = 0.0D+0 totNeg = 0.0D+0 targetValues = clusterMissingValue do cl = 1, numClusters pCluster => clusters(cl) pCluster % totalDemand = 0.0D+0 pCluster % demandPos = 0.0D+0 pCluster % demandNeg = 0.0D+0 ! Add the requested or delivered volumes of the sub links do sLink = 1, dmm_NUM_SUB_TYPES if (associated(pCluster % subLinks(sLink) % link)) then pSubLink => pCluster % subLinks(sLink) % link pCluster % totalDemand = pCluster % totalDemand + pSubLink % tgt % targetDemand pCluster % demandPos = pCluster % demandPos + pSubLink % demandPos pCluster % demandNeg = pCluster % demandNeg + pSubLink % demandNeg endif end do indxTgt = pCluster % tgtIndex targetValues(indxTgt) = pCluster % totalDemand / (dmmDtsw * numSecPerDay) ! Vol per MSW-step => m3/s ! for logging totDemand = totDemand + pCluster % totalDemand totPos = totPos + pCluster % demandPos totNeg = totNeg + pCluster % demandNeg end do ! log total demand for FM totTarget = 0.0D+0 do t = 1, size(targetValues) if ( dabs(targetValues(t) - clusterMissingValue) > epsilon ) then totTarget = totTarget + targetValues(t) endif enddo write(demandLineFm, '(I4,A,I2,A,F15.9,A,F12.5,A,F12.5,A,F12.5,A,F12.5,A)') & dmmMfStep, ',', dmmMswStep, ',', & totTarget, ',', totTarget * dmmDtsw * numSecPerDay, ',', & totDemand, ',', totPos, ',', totNeg, ',' end subroutine dmmClustersGetFluxes ! ! This function is called by the driver to set the fluxes as they were realized by D-Flow FM. ! It distributes the realized values over the various requesters ! subroutine dmmClustersSetRealizedFluxes(realizedFluxes) ! arguments double precision, dimension(:), intent(in) :: realizedFluxes ! computed realized fluxes ! locals integer :: cl, sLink, s ! loop counter for clusters, subLinks and sources integer :: indxTgt ! index at target side, index double precision :: stillAvailable ! Available volume, initially and after substracting sublink volumes 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 :: 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_cluster), pointer :: pCluster ! pointer to a cluster type(t_dmm_link) , pointer :: pMfLink ! pointer to a link in the iMODFLOW->FM 1D-fluxes mapper type(t_dmm_link) , pointer :: pMswLink ! pointer to a link in the MetaSWAP->FM sprinkling mapper ! For logging: integer :: link, rf ! link counter, realized flux counter type(t_dmm_link), pointer :: pLink ! pointer to link integer :: mhMswRunoff ! handle to mapper for MSW-runoff to FM-1D-nodes double precision :: totFmIn ! sum of realizedFluxes (the fluxes as realized by FM) double precision :: totFmR ! totals for the sublink realizedFluxes character(Len=1024) :: logLine ! composed log line for fluxes 1D double precision, & dimension(dmm_NUM_SUB_TYPES-1) :: totSubs ! totals per sub link in a cluster double precision :: sumNotConnected, & ! summed volume provide by MetaSWAP on links that are not connected to FM sumConnected ! body logLine = '' totFmIn = 0.0D+0 totFmR = 0.0D+0 totSubs = 0.0D+0 do rf = 1, size(realizedFluxes) totFmIn = totFmIn + realizedFluxes(rf) enddo do cl = 1, numClusters pCluster => clusters(cl) pMfLink => pCluster % subLinks(dmm_MF2FM_Q) % link pMswLink => pCluster % subLinks(dmm_MSW2FM_SPR) % link indxTgt = pCluster % tgtIndex stillAvailable = realizedFluxes(indxTgt) * dmmDtsw * numSecPerDay ! m3/s => Vol per MSW-step totFmR = totFmR + stillAvailable pCluster % totalRealized = stillAvailable if(dumpClusterFluxes .and. pCluster % logMapping) then write(logLine,'(I4,A,I2,A,F16.5,A,F16.5,A,F16.10,A,F16.10)') dmmMfStep, ',', dmmMswStep, ',', & pCluster % tgtX, ',', pCluster % tgtY, ',', & pCluster % totalDemand / (dmmDtsw * numSecPerDay) , ',', realizedFluxes(indxTgt) endif ! First substract the always realized drain / one way river / runoff fluxes do sLink = dmm_MF2FM_DRN, dmm_MSW2FM_RUN if (associated(pCluster % subLinks(sLink) % link )) then call dmmLogFixedFluxLink(pCluster, sLink) stillAvailable = stillAvailable - pCluster % subLinks(sLink) % link % demandPos if(dumpClusterFluxes .and. pCluster % logMapping) then write(logLine,'(A,A,F12.4)') trim(logLine), ',', pCluster % subLinks(sLink) % link % demandPos endif pLink => pCluster % subLinks(sLink) % link totSubs(sLink) = totSubs(sLink) + pLink % demandPos subTypeUsed(sLink) = .true. pLink % totalRealized = pLink % demandPos do s = 1, pLink % numSrcs pLink % srcs(s) % realized = pLink % srcs(s) % demand end do else if(dumpClusterFluxes .and. pCluster % logMapping) then write(logLine,'(A,A,F12.4)') trim(logLine), ',', 0.D+0 endif endif enddo ! indicate that MF->FM is always true subTypeUsed(dmm_MF2FM_Q) = .true. ! 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 = 0.D+0 mswDemand = 0.D+0 if (associated(pMfLink)) then mfDemand = pMfLink % demandPos + pMfLink % demandNeg endif if (associated(pMswLink)) then mswDemand = pMswLink % demandPos + pMswLink % demandNeg endif mswRealized = 0.0D+0 ! for debugging mfRealized = 0.0D+0 ! " " if (stillAvailable > 0 .or. (stillAvailable <= 0 .and. stillAvailable <= (mfDemand + mswDemand)) ) then ! All has been realized ! iMODFLOW RIV -> D-FM link ! Add the realized demand for this Dtsw step to the full day realized value if (associated(pMfLink)) then do s = 1, pMfLink % numSrcs pMfLink % srcs(s) % realized = pMfLink % srcs(s) % realized + pMfLink % srcs(s) % demand end do endif mfRealized = mfDemand ! MetaSWAP sprinkling -> D-FM sub link if (associated(pMswLink)) then do s = 1, pMswLink % numSrcs pMswLink % srcs(s) % realized = pMswLink % srcs(s) % demand end do endif mswRealized = mswDemand else if (stillAvailable < mfDemand) then ! iMODFLOW demand could be realized. ! Add the realized demand for this Dtsw step to the full day realized value if (associated(pMfLink)) then do s = 1, pMfLink % numSrcs pMfLink % srcs(s) % realized = pMfLink % srcs(s) % realized + pMfLink % srcs(s) % demand end do pMfLink % totalRealized = mfDemand endif mfRealized = mfDemand ! Provide remainder to MetaSWAP ! Distribute negative shortage over the negative demand ! TODO: simplify (no positive sprinkling!) stillAvailable = stillAvailable - mfDemand if (associated(pMswLink)) then pMswLink % totalRealized = stillAvailable shortage = mswDemand - stillAvailable notRealizedFraction = shortage / pMswLink % demandNeg do s = 1, pMswLink % numSrcs srcDemand = pMswLink % srcs(s) % demand if (srcDemand < 0) then pMswLink % srcs(s) % realized = srcDemand * (1.0D+0 - notRealizedFraction) else pMswLink % srcs(s) % realized = srcDemand endif mswRealized = mswRealized + pMswLink % 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 - stillAvailable if (associated(pMfLink)) then notRealizedFraction = shortage / pMfLink % demandNeg do s = 1, pMfLink % numSrcs srcDemand = pMfLink % srcs(s) % demand if (srcDemand < 0) then pMfLink % srcs(s) % realized = pMfLink % srcs(s) % realized + pMfLink % srcs(s) % demand * (1.0D+0 - notRealizedFraction) else pMfLink % srcs(s) % realized = pMfLink % srcs(s) % realized + pMfLink % srcs(s) % demand endif mfRealized = mfRealized + pMfLink % srcs(s) % realized end do endif if (associated(pMswLink)) then do s = 1, pMswLink % numSrcs pMswLink % srcs(s) % realized = 0.0D+0 end do endif endif if(pCluster % logMapping) then write(logLine,'(A,A,F12.4,A,F12.4,A,F12.4,A,F12.4)') trim(logLine), & ',', mfDemand, ',', mfRealized, ',', mswDemand, ',', mswRealized if (lunFluxes1D /= -1) then write(lunFluxes1D,'(A)') trim(logLine) call flush(lunFluxes1D) endif endif end do ! Log the always realized drain / one way river / runoff fluxes do sLink = 1, dmm_NUM_SUB_TYPES if (subTypeUsed(sLink)) then write(realizedLine(sLink),'(F12.5,A)') totSubs(sLink), ',' else ! not used write zero's (including the one for the for this mapper type ! not called dmm1DMapperGetRealizedVolumes) write(realizedLine(sLink),'(F12.5,A,F12.5,A)') 0.0D+0, ',', 0.0D+0, ',' endif enddo ! log (just to be sure) the MSW runoff values for not connected nodes mhMswRunoff = dmmGetMapperHandle(dmm_MSW2FM_RUN) if (mhMswRunoff /= -1) then sumConnected = 0.D+0 sumNotConnected = 0.D+0 do link = 1, mappings(mhMswRunoff) % numLinks pLink => mappings(mhMswRunoff) % links(link) do s = 1, pLink % numSrcs if (pLink % connected) then sumConnected = sumConnected + pLink % srcs(s) % demand else sumNotConnected = sumNotConnected + pLink % srcs(s) % demand endif end do enddo write(connectedLine,'(F12.5,A,F12.5,A)') sumConnected, ',', sumNotConnected, ',' else write(connectedLine,'(F12.5,A,F12.5,A)') 0.0D+0, ',', 0.0D+0, ',' endif ! log the realized FM volumes write(realizedLineFm, '(F15.8,A,F12.5,A,F12.5,A)') totFmIn, ',', totFmIn * (dmmDtsw * numSecPerDay), ',', totFmR, ',' end subroutine dmmClustersSetRealizedFluxes ! ! This function returns the requested volumes for a mapper that handles 1D fluxes. ! Called by the driver for each flux (= volume per time step) mapper. ! subroutine dmm1DMapperGetRealizedVolumes(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, v ! loop counter for links, sources and values integer :: indxSrc ! Index at source site type(t_dmm_link), pointer :: pLink ! pointer to link and sub link double precision :: totalRealized ! total realized values double precision :: totalCorrection ! all over correction flux to iMODFLOW RIV double precision, parameter :: epsilon = 1.0D-10 ! epsilon for comparing demanded and realized flux totalCorrection = 0.0D+0 ! body realizedVolumes = mappings(mh) % missingValue(dmmSrc) 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 dmm1DMapperGetRealizedVolumes' endif endif end do if (mappings(mh) % dumpMapping .and. pLink % logMapping) call flush(mappings(mh) % lunDumpFile) end do ! Log the realized volumes totalRealized = 0.0D+0 do v = 1, size(realizedVolumes) if ( dabs(realizedVolumes(v) - mappings(mh) % missingValue(dmmTgt)) > epsilon ) then totalRealized = totalRealized + realizedVolumes(v) endif enddo write(realizedLine(mappings(mh) % mapperType),'(A,F12.5,A)') & trim(realizedLine(mappings(mh) % mapperType)), totalRealized, ',' end subroutine dmm1DMapperGetRealizedVolumes !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! 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 ! ! Find an existing link in a mapper, or create a new link ! function dmmFindOrAddLink(mh, tgtId) result(link) ! return value integer :: link ! handle to found or new link ! arguments integer, intent(in) :: mh ! mapper handle 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 type(t_dmm_mapping), pointer :: mapper ! pointer to mapper ! body status = 0 link = dmmIndexInvalid mapper => mappings(mh) ! find link with same target id if (tgtId <= mapper % lk_lookup_length) then if (mapper % lk_lookup(tgtId) > 0) then link = mapper % lk_lookup(tgtId) return endif endif ! 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 if (link /= dmmIndexInvalid) then mapper % numLinks = link mapper % links(mapper % numLinks) % tgt % id = tgtId 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. call dmmAddLinkToLookupArray(mapper, link, tgtId) endif end function dmmFindOrAddLink ! ! Find an existing cluster, or create a new cluster ! function dmmFindOrAddCluster(tgtId) result(cluster) ! return value integer :: cluster ! handle to new cluster ! arguments integer, intent(in) :: tgtId ! id at target side (i.e. the FM 1D-node index) ! locals type(t_dmm_cluster) , & pointer, dimension(:) :: oldClusters ! old clusters (before re-alloc) integer :: numOldClusters ! #old clusters (before re-alloc) integer :: status ! allocation status flag integer :: sLink ! loop counter for sub links type(t_dmm_cluster), pointer :: pNewCluster ! pointer to newly created cluster ! body status = 0 cluster = dmmIndexInvalid ! find cluster with same target id if (tgtId <= cl_lookup_length) then if (cl_lookup(tgtId) > 0) then cluster = cl_lookup(tgtId) return endif endif ! not found, add cluster and initialize it ! first re-allocate if necessary if ( maxNumClusters == numClusters ) then oldClusters => clusters numOldClusters = numClusters maxNumClusters = maxNumClusters + DMM_LK_BLOCK_SIZE allocate(clusters(maxNumClusters), stat=status) if ( status /= 0 ) then write(dmmLastError, '(2A,I4,A)') 'dmm: dmmFindOrAddCluster: ', & 'Could not allocate ', maxNumClusters, ' items' else clusters(1:numOldClusters) = oldClusters(1:numOldClusters) cluster = numOldClusters+1 deallocate(oldClusters) endif else cluster = numClusters + 1 endif if (cluster /= dmmIndexInvalid) then numClusters = cluster pNewCluster => clusters(cluster) pNewCluster % tgtIndex = tgtId pNewCluster % tgtX = -777.7 pNewCluster % tgtY = -777.7 pNewCluster % totalDemand = 0.0D+0 pNewCluster % totalRealized = 0.0D+0 pNewCluster % demandPos = 0.0D+0 pNewCluster % demandNeg = 0.0D+0 pNewCluster % logMapping = .false. do sLink = 1, dmm_NUM_SUB_TYPES pNewCluster % subLinks(sLink) % link => NULL() enddo call dmmAddClusterToLookupArray(tgtId, cluster) endif end function dmmFindOrAddCluster subroutine dmmAddClusterToLookupArray(tgtId, clusterIndex) ! arguments integer, intent(in) :: tgtId ! target id integer, intent(in) :: clusterIndex ! index in the cluster array ! locals integer :: newLookupLength ! required lenght for lookup table integer, & pointer, dimension(:) :: oldLookups ! pointer to old lookup table (needed for copying) integer :: oldLookupLength ! old length of the lookup table integer :: status ! allocation status flag ! body ! check if the cluster index lookup array can hold the tgtId newLookupLength = cl_lookup_length do while ( tgtId > newLookupLength ) newLookupLength = newLookupLength + DMM_CL_BLOCK_SIZE enddo if (newLookupLength > cl_lookup_length) then oldLookups => cl_lookup oldLookupLength = cl_lookup_length cl_lookup_length = newLookupLength allocate(cl_lookup(cl_lookup_length), stat=status) if ( status /= 0 ) then write(dmmLastError, '(2A,I4,A)') 'dmm: dmmAddClusterToLookupArray: ', & 'Could not allocate ', cl_lookup_length, ' items' else cl_lookup = -1 cl_lookup(1:oldLookupLength) = oldLookups(1:oldLookupLength) deallocate(oldLookups) endif endif ! set target lookup index cl_lookup(tgtId) = clusterIndex end subroutine dmmAddClusterToLookupArray subroutine dmmAddLinkToLookupArray(mapper, link, tgtId) ! arguments type(t_dmm_mapping), pointer :: mapper ! pointer to mapper integer , intent(in) :: link ! link index integer , intent(in) :: tgtId ! target id ! locals integer :: newLookupLength ! required lenght for lookup table integer, & pointer, dimension(:) :: oldLookups ! pointer to old lookup table (needed for copying) integer :: oldLookupLength ! old length of the lookup table integer :: status ! allocation status flag ! body ! check if the link index lookup array can hold the tgtId newLookupLength = mapper % lk_lookup_length do while ( tgtId > newLookupLength ) newLookupLength = newLookupLength + DMM_LK_BLOCK_SIZE enddo if (newLookupLength > mapper % lk_lookup_length) then oldLookups => mapper % lk_lookup oldLookupLength = mapper % lk_lookup_length mapper % lk_lookup_length = newLookupLength allocate(mapper % lk_lookup(mapper % lk_lookup_length), stat=status) if ( status /= 0 ) then write(dmmLastError, '(2A,I4,A)') 'dmm: dmmAddLinkToLookupArray: ', & 'Could not allocate ', mapper % lk_lookup_length, ' items' else mapper % lk_lookup = -1 mapper % lk_lookup(1:oldLookupLength) = oldLookups(1:oldLookupLength) deallocate(oldLookups) endif endif ! set target lookup index mapper % lk_lookup(tgtId) = link end subroutine dmmAddLinkToLookupArray 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(pCluster, subLinkIndex) ! arguments type(t_dmm_cluster), pointer :: pCluster ! pointer to a cluster integer, intent(in) :: subLinkIndex ! index of sublink for one way flux ! (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. pCluster % logMapping)) then return endif pSubLink => pCluster % 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