!! Copyright (C) Stichting Deltares, 2005-2014.
!!
!! This file is part of iMOD.
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see .
!!
!! Contact: imod.support@deltares.nl
!! Stichting Deltares
!! P.O. Box 177
!! 2600 MH Delft, The Netherlands.
!> Module for the coupling MODFLOW-MetaSWAP with dgflow.
module dgflowmod
implicit none
character(len=100), parameter :: mflpref = 'mf2dgf'
logical, save :: first = .true.
character(len=100), parameter :: dgratesignal = 'mf2dgf_flux_done.txt' !< DFLOW signal file for rates done
character(len=100), parameter :: dgexitsignal = 'mf2dgf_exit.txt' !< DFLOW signal file for exiting
character(len=100), parameter :: mfinitsignal = 'mf2dgf_init_done.txt' !< MODFLOW signal file for initialization done
character(len=100), parameter :: mfbndsignal = 'mf2dgf_bnd_done.txt' !< MODFLOW signal file for boundary done
character(len=100), parameter :: mfheadsignal = 'mf2dgf_head_done.txt' !< MODFLOW signal file for heads done
!... Parameters
integer, parameter :: tsleep = 1 !< time in seconds to check coupling signal file
double precision :: beginOfCurrentTimeStep, endOfCurrentTimeStep
integer, parameter :: mv = -1234.0 !< general missing value
end module
!! =============================================================================
! !> Initialize DGFLOW: read command line
! subroutine dgflow_initComponent(cla_dgflow)
!! =============================================================================
!!... modules
! use dgflowmod
!
! implicit none
!
!!... arguments
! character (len=*), intent(in) :: cla_dgflow !< command line for dgflow
!
!!... functions
! integer :: cfn_n_elem
!
!!... local variables
! character(len=50), dimension(:), allocatable :: str
! integer :: iarg, nargs
!!.......................................................................
!
! write(*,*) 'dgf: *** initialization ***'
! write(*,*) 'dgf: reading command line arguments'
!
! nargs = cfn_n_elem(' ',1,cla_dgflow)
!
! if (nargs.eq.0) then
! write(*,*) 'dgf: error, dgflow commandline arguments not found!'
! call exit(1)
! end if
! allocate(str(nargs))
! read(cla_dgflow,*)(str(iarg),iarg=1,nargs)
!
! ! read(str(1),*) dgfsignal ! read DGFLOW signal file name
!
! deallocate(str)
!
! end subroutine
!
!! =============================================================================
! !> Initialize DGFLOW interface: allocate arrays
! subroutine dgflow_initSimulation()
!! =============================================================================
!!... modules
! use dgflowmod
!
! implicit none
!!.......................................................................
!
! end subroutine
!
!! =============================================================================
! subroutine dgflow_prepareTimeStep(endOfSimulation,convergeddgflow,currentTime)
!! =============================================================================
!!... modules
! use dgflowmod
!
! implicit none
!
!!... arguments
! logical, intent(inout) :: endOfSimulation
! logical, intent(out) :: convergeddgflow
! double precision, intent(inout) :: currentTime
!
!!... functions
! integer :: cfn_getlun
!
!!... locals
! logical :: lcont, lexit
! character(len=1) :: dora
! character(len=15) :: time
! integer :: lun, ios, date, year, month, day, hour, minute, second
! double precision :: dt
!!.......................................................................
!
!!... wait until the dgflow signal file is found
! write(*,*) 'dgf: waiting for signal from dgflow ',dgfsignal(1:len_trim(dgfsignal)), ' or ',&
! dgfexitsignal(1:len_trim(dgfexitsignal)),' to be found ...'
! do while(.true.)
! inquire(file=dgfsignal,exist=lcont)
! inquire(file=dgfexitsignal,exist=lexit)
! if (lcont.or.lexit) exit
! call sleep(tsleep)
! end do
!
! endOfSimulation = .false.
! if (lexit) then
! write(*,*) 'dgf: found ',dgfexitsignal(1:len_trim(dgfexitsignal))
! write(*,*) 'dgf: dgflow finished simulation'
! endOfSimulation = .true.
! return
! end if
!
! !... open the dgflow signal file, read, and delete file
! write(*,*) 'dgf: found ',dgfsignal(1:len_trim(dgfsignal)),' and reading'
! lun = cfn_getlun(10,99)
! open(unit=lun,file=dgfsignal,status='old',form='formatted',iostat=ios)
! if (ios.ne.0) write(*,*) 'Error opening file...'
!
!!... read dgflow signal file
! read(lun,*) dora, time, dt
! write(*,*) 'dgf: read from signal file: ', dora, time, dt
!
!!... delete signal file
! close(lun,status='delete')
!
! end subroutine
!
!! =============================================================================
! subroutine dgflow_initTimeStep(iter)
!! =============================================================================
! use dgflowmod
!
! implicit none
!
!!... arguments
! integer, intent(in) :: iter
!!.......................................................................
!
! end subroutine
!
!! =============================================================================
! subroutine dgflow_finishTimeStep(iter,dt,usetransol)
!! =============================================================================
!!... modules
!
! implicit none
!
!!... arguments
! integer, intent(in) :: iter
! double precision, intent(in) :: dt
! logical, intent(in) :: usetransol
!!.......................................................................
!
! ! write the signal file
! call dgflow_write_sig()
!
! end subroutine
!
! !> Write dgflow signal file.
! subroutine dgflow_write_sig()
!
! use dgflowmod
!
! implicit none
!
!!... functions
! integer :: cfn_getlun
!
!!... local variables
! integer :: lun, ios
! character(len=100) :: tmpfile
!!.......................................................................
!
!!... write (empty) MODFLOW-MetaSWAP file
! write(*,*) 'dgf: writing signal file for dgflow (',trim(mfmssignal),')'
!
! lun = cfn_getlun(10,99)
! write(tmpfile,'(2a)') mfmssignal(1:len_trim(mfmssignal)), 'tmp'
! open(unit=lun,file=tmpfile,status='new',form='formatted',share='denyrd')
! close(lun)
! call osd_rename(tmpfile,mfmssignal,ios) ! rename file
! if (ios.ne.0) write(*,*) 'Error renaming file...'
!
! end subroutine
!
!> write geometry
logical function dgflow_GetModDgfGeom(cells,geom,ncells,nlay,nrow,ncol,nper,perlen,dt,ts,clay)
!... modules
use dgflowmod
use rf2mf_module, only: dis
implicit none
!... arguments
integer, intent(in) :: ncells
integer, dimension(3,ncells), intent(in) :: cells
real, dimension(13,ncells), intent(in) :: geom
integer, intent(in) :: nlay
integer, intent(in) :: nrow
integer, intent(in) :: ncol
integer, intent(in) :: nper
real, dimension(nper), intent(in) :: perlen
real, intent(in) :: dt
double precision, intent(in) :: ts
integer, intent(in) :: clay
!... locals
integer :: lun, cfn_getlun, iper, i, j, ilay, irow, icol, n, mcol, mrow, mlay
character(len=256), dimension(10) :: sa
character(len=256) :: s, mflfile
integer, dimension(:,:), allocatable :: idx
double precision :: t, dt2
integer :: date, hour, minute, second, nt
!.......................................................................
mlay = maxval(cells(1,:))-minval(cells(1,:)) + 1
mrow = maxval(cells(2,:))-minval(cells(2,:)) + 1
mcol = maxval(cells(3,:))-minval(cells(3,:)) + 1
lun = cfn_getlun(10,99)
write(mflfile,'(2a)') trim(mflpref)//'_init.mfl'
open(unit=lun,file=mflfile)
! nx 150 ny 100 nz 19
write(sa(1),*) mcol
write(sa(2),*) mrow
write(sa(3),*) mlay*2-1
! if (mlay.ne.nlay)then
! write(sa(3),*) mlay*2
! else
! write(sa(3),*) mlay*2-1
! end if
write(s,'(6(a,1x),a)') 'nx',trim(adjustl(sa(1))),'ny',trim(adjustl(sa(2))),'nz',trim(adjustl(sa(3))),''
write(lun,'(a)') trim(s)
allocate(idx(2,nlay))
idx(1,:) = ncells; idx(2,:) = 0
do i = 1, ncells
ilay = cells(1,i)
idx(1,ilay) = min(i,idx(1,ilay))
idx(2,ilay) = max(i,idx(2,ilay))
end do
do ilay = 1, nlay
if (idx(2,ilay).eq.0) cycle
! aquifer
do i = idx(1,ilay), idx(2,ilay)
irow = cells(2,i)
icol = cells(3,i)
n = icol + (irow-1)*ncol + (ilay-1)*nrow*ncol
write(sa(1),*) n
write(sa(2),*) geom(1,i) ! x
write(sa(3),*) geom(2,i) ! y
write(sa(4),*) geom(3,i) ! z1
write(sa(5),*) geom(5,i) ! dx
write(sa(6),*) geom(6,i) ! dy
write(sa(7),*) geom(7,i) ! dz1
write(sa(8),*) geom(9,i) ! kh1
write(sa(9),*) geom(10,i)! kv1
write(sa(10),*) geom(13,i) ! head
write(lun,'(9(a,1x),a)')( trim(adjustl(sa(j))),j=1,10)
end do
! aquitard
if (ilay.ne.min(clay,nlay)) then
do i = idx(1,ilay), idx(2,ilay)
irow = cells(2,i)
icol = cells(3,i)
n = icol + (irow-1)*ncol + (ilay-1)*nrow*ncol
write(sa(1),*) -n
write(sa(2),*) geom(1,i) ! x
write(sa(3),*) geom(2,i) ! y
write(sa(4),*) geom(4,i) ! z2
write(sa(5),*) geom(5,i) ! dx
write(sa(6),*) geom(6,i) ! dy
write(sa(7),*) geom(8,i) ! dz2
write(sa(8),*) geom(11,i) ! kh2
write(sa(9),*) geom(12,i) ! kv2
write(sa(10),*) geom(13,i) ! head
write(lun,'(9(a,1x),a)')( trim(adjustl(sa(j))),j=1,10)
end do
end if
end do
! time step information
write(sa(1),*) nper
write(sa(2),*) dt
write(lun,'(a,1x,a)')( trim(adjustl(sa(i))),i=1,2)
t = ts
do iper = 1, nper
nt = nint(perlen(iper)/dt)
dt2 = dble(perlen(iper))
call cfn_mjd2datehms(t,date,hour,minute,second)
write(sa(1),'(i8,3(a,i2.2))') date,' ',abs(hour),':',minute,':',second
t = t + dt2
call cfn_mjd2datehms(t,date,hour,minute,second)
write(sa(2),'(i8,3(a,i2.2))') date,' ',abs(hour),':',minute,':',second
write(sa(3),*) nt
if (dis%sp(iper)%writeoc) then
write(sa(4),*) 1
else
write(sa(4),*) 0
end if
write(lun,'(3(a,1x),a)')( trim(adjustl(sa(i))),i=1,4)
end do
close(lun)
! write the signalfile
lun = cfn_getlun(10,99)
open(unit=lun,file=mfinitsignal)
write(lun,'(a)') trim(mflfile)
close(lun)
deallocate(idx)
dgflow_GetModDgfGeom = .true.
end function
logical function dgflow_GetModDgfRiversDrains(ts,t,geom,ncells,nlay,nrow,ncol,rivr,nriver,irivsys,nrivsys,drai,ndrain,idrnsys,ndrnsys)
!... modules
use dgflowmod
implicit none
!... arguments
integer, intent(in) :: ts,ncells,nlay,nrow,ncol,nriver,irivsys,nrivsys,ndrain,idrnsys,ndrnsys
double precision, intent(in) :: t
integer, dimension(3,ncells) :: geom
real, dimension(8,*) :: rivr
real, dimension(6,*) :: drai
!... locals
integer :: lun, cfn_getlun, ig, i, isys, iact, ilay, irow, icol, n, m, ilist
integer :: date, hour, minute, second
character(len=256) :: bndfile, s
character(len=256), dimension(10) :: sa
integer, dimension(:,:,:), allocatable :: iwrk
!.......................................................................
dgflow_GetModDgfRiversDrains = .true.
call cfn_mjd2datehms(t,date,hour,minute,second)
write(s,*) ts-1
write(bndfile,'(3a)') trim(mflpref),'_bnd.cp',trim(adjustl(s))
lun = cfn_getlun(10,99)
open(unit=lun,file=bndfile)
if (nriver.le.0.and.ndrain.le.0)then
write(lun,'(a)') '0'
else
write(sa(1),*) nrivsys
write(sa(2),*), ndrnsys
write(sa(3),'(i8,3(a,i2.2))') date,' ',abs(hour),':',minute,':',second
write(lun,'(a,1x,a,1x,a)')(trim(adjustl(sa(i))),i=1,3)
allocate(iwrk(ncol,nrow,nlay))
! rivers
do isys = 1, nrivsys
iwrk = 0
do ilist = 1, nriver
if (int(rivr(irivsys,ilist)).eq.isys) then
ilay = rivr(1,ilist)
irow = rivr(2,ilist)
icol = rivr(3,ilist)
iwrk(icol,irow,ilay) = ilist
end if
end do
do iact = 1, 2
m = 0
do ig = 1, ncells
ilay = geom(1,i)
irow = geom(2,i)
icol = geom(3,i)
ilist = iwrk(icol,irow,ilay)
if (ilist.gt.0) then
m = m + 1
if (iact.eq.2) then
n = icol + (irow-1)*ncol + (ilay-1)*nrow*ncol
write(sa(1),*) n
write(sa(2),*) rivr(4,ilist) ! stage
write(sa(3),*) rivr(5,ilist) ! cond
write(sa(4),*) rivr(6,ilist) ! rbot
write(sa(5),*) rivr(7,ilist) ! inf
write(lun,'(4(a,1x),a)')(trim(adjustl(sa(i))),i=1,5)
end if
end if
end do
if (iact.eq.1) then
write(sa(1),*) m
write(lun,'(a)') trim(adjustl(sa(1)))
end if
end do ! iact
end do ! rivsys
! drains
do isys = 1, ndrnsys
iwrk = 0
do ilist = 1, ndrain
if (int(drai(idrnsys,ilist)).eq.isys) then
ilay = drai(1,ilist)
irow = drai(2,ilist)
icol = drai(3,ilist)
iwrk(icol,irow,ilay) = ilist
end if
end do
do iact = 1, 2
m = 0
do ig = 1, ncells
ilay = geom(1,ig)
irow = geom(2,ig)
icol = geom(3,ig)
ilist = iwrk(icol,irow,ilay)
if (ilist.gt.0) then
m = m + 1
if (iact.eq.2) then
n = icol + (irow-1)*ncol + (ilay-1)*nrow*ncol
write(sa(1),*) n
write(sa(2),*) drai(4,ilist) ! stage
write(sa(3),*) drai(5,ilist) ! cond
write(lun,'(2(a,1x),a)')(trim(adjustl(sa(i))),i=1,3)
end if
end if
end do
if (iact.eq.1) then
write(sa(1),*) m
write(lun,'(a)') trim(adjustl(sa(1)))
end if
end do ! iact
end do ! drnsys
deallocate(iwrk)
end if
close(lun)
! write the signalfile
lun = cfn_getlun(10,99)
open(unit=lun,file=mfbndsignal)
write(sa(1),'(i8,3(a,i2.2))') date,' ',abs(hour),':',minute,':',second
write(sa(2),'(a)') bndfile
write(lun,'(a)') trim(adjustl(sa(1)))
write(lun,'(a)') trim(adjustl(sa(2)))
close(lun)
end function
logical function dgflow_GetHeads(ts, t, bnd, nbnd, heads, nlay, nrow, ncol)
use dgflowmod
implicit none
!... arguments
integer, intent(in) :: ts
double precision, intent(in) :: t
integer, intent(in) :: nbnd
integer, dimension(4,nbnd), intent(in) :: bnd
double precision, dimension(nbnd), intent(in) :: heads
integer, intent(in) :: nlay
integer, intent(in) :: nrow
integer, intent(in) :: ncol
!... locals
integer :: lun, cfn_getlun, ibnd, i, n, icol, irow, ilay, flg
integer :: date, hour, minute, second
character(len=256) :: headfile, s
character(len=256), dimension(10) :: sa
!.......................................................................
dgflow_GetHeads = .true.
call cfn_mjd2datehms(t,date,hour,minute,second)
write(s,*) ts
write(headfile,'(3a)') trim(mflpref),'_head.cp',trim(adjustl(s))
lun = cfn_getlun(10,99)
open(unit=lun,file=headfile)
write(sa(1),*) nbnd
write(sa(2),'(i8,3(a,i2.2))') date,' ',abs(hour),':',minute,':',second
write(lun,'(a,1x,a)')(trim(adjustl(sa(i))),i=1,2)
do ibnd = 1, nbnd
ilay = bnd(1,ibnd)
irow = bnd(2,ibnd)
icol = bnd(3,ibnd)
flg = bnd(4,ibnd)
n = icol + (irow-1)*ncol + (ilay-1)*nrow*ncol
if (flg.eq.-1) n = -n
write(sa(1),*) n
write(sa(2),*) heads(ibnd)
! write(sa(2),*) 0.D0
write(lun,'(a,1x,a)')(trim(adjustl(sa(i))),i=1,2)
end do
close(lun)
! write the signalfile
lun = cfn_getlun(10,99)
open(unit=lun,file=mfheadsignal)
write(sa(1),'(i8,3(a,i2.2))') date,' ',abs(hour),':',minute,':',second
write(sa(2),'(a)') headfile
write(lun,'(a)') trim(adjustl(sa(1)))
write(lun,'(a)') trim(adjustl(sa(2)))
close(lun)
end function
logical function dgflow_PutRatesMf(rates)
use dgflowmod
implicit none
!... arguments
real, dimension(*), intent(out) :: rates
!... locals
logical :: lexit, lcont
character(len=256) :: ratefile
integer :: lun, cfn_getlun, i, n, idum
real :: rate
!.......................................................................
dgflow_PutRatesMf = .true.
! wait for signal file
100 continue ! check (TEMPORARY)
do while(.true.)
write(*,'( 4(a,1x),a)') 'Waiting for',trim(dgratesignal),'or',trim(dgexitsignal),'...'
inquire(file=dgratesignal,exist=lcont)
inquire(file=dgexitsignal,exist=lexit)
if (lcont.or.lexit) exit
call sleep(tsleep)
end do
if (lexit) then
dgflow_PutRatesMf = .false.
return
end if
lun = cfn_getlun(10,99)
open(unit=lun,file=dgratesignal)
read(lun,'(a)') ratefile
close(lun,status='delete')
! check (TEMPORARY)
!if (index(ratefile,'.cp0').gt.0) goto 100
! read the rate file
if(.true.)then
lun = cfn_getlun(10,99)
write(*,*) 'dg: opening '//trim(ratefile)//'...'
open(unit=lun,file=ratefile)
read(lun,*) n
do i = 1, n
read(lun,*) idum, rate
rates(i) = rate
end do
close(lun)
else
write(*,*) '!!!! Not using fluxes from DgFlow !!!!'
end if
end function