c Copyright (C) Stichting Deltares, 2005-2017.
c
c This file is part of iMOD.
c
c This program is free software: you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation, either version 3 of the License, or
c (at your option) any later version.
c
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c GNU General Public License for more details.
c
c You should have received a copy of the GNU General Public License
c along with this program. If not, see .
c
c Contact: imod.support@deltares.nl
c Stichting Deltares
c P.O. Box 177
c 2600 MH Delft, The Netherlands.
c
c iMod is partly based on the USGS MODFLOW2005 source code;
c for iMOD the USGS MODFLOW2005 source code has been expanded
c and extensively modified by Stichting Deltares.
c The original USGS MODFLOW2005 source code can be downloaded from the USGS
c website http://www.usgs.gov/. The original MODFLOW2005 code incorporated
c in this file is covered by the USGS Software User Rights Notice;
c you should have received a copy of this notice along with this program.
c If not, see .
SUBROUTINE VDF2RIV7FM(IGRID)
C ******************************************************************
C ADD RIVER TERMS TO RHS AND HCOF
C--C--SEAWAT: MODIFIED TO USE A VARIABLE DENSITY FORM OF DARCY'S LAW
C ******************************************************************
C
C SPECIFICATIONS:
C ------------------------------------------------------------------
USE GLOBAL, ONLY:IBOUND,HNEW,RHS,HCOF
USE GWFRIVMODULE, ONLY:NRIVER,RIVR
1 ,IRIVRFACT,RIVAUX ! RFACT
USE VDFMODULE, ONLY: DENSEREF,PS,ELEV
C
C--SEAWAT: ADD HRIV
DOUBLE PRECISION RRBOT,HRIV
C--SEAWAT: NEED AUXILIARY VARIABLES
! COMMON /RIVCOM/RIVAUX(5)
! CHARACTER*16 RIVAUX
C ------------------------------------------------------------------
CALL SGWF2RIV7PNT(IGRID)
C
C1------IF NRIVER<=0 THERE ARE NO RIVERS. RETURN.
IF(NRIVER.LE.0)RETURN
C
C--SEAWAT: FIND RBDTHK AND RIVDEN IF EXIST AS AUX VARIABLES
LOCRIVDEN=0
DO I=1,5
IF(RIVAUX(I).EQ.'RIVDEN') LOCRIVDEN=I+6
ENDDO
C2------PROCESS EACH CELL IN THE RIVER LIST.
DO 100 L=1,NRIVER
C
C3------GET COLUMN, ROW, AND LAYER OF CELL CONTAINING REACH.
IL=RIVR(1,L)
IR=RIVR(2,L)
IC=RIVR(3,L)
C
C4------IF THE CELL IS EXTERNAL SKIP IT.
IF(IBOUND(IC,IR,IL).LE.0)GO TO 100
C
C5------SINCE THE CELL IS INTERNAL GET THE RIVER DATA.
HRIV=RIVR(4,L)
CRIV=RIVR(5,L)
RBOT=RIVR(6,L)
RRBOT=RBOT
RBDTHK=ABS(RBOT-ELEV(IC,IR,IL))
RIVDENS=PS(IC,IR,IL)
IF(LOCRIVDEN.GT.0) RIVDENS=RIVR(LOCRIVDEN,L)
HRIV=FEHEAD(HRIV,RIVDENS,RBOT+RBDTHK)
HHNEW=HNEW(IC,IR,IL)
HTEMP=SALTHEAD(HNEW(IC,IR,IL),PS(IC,IR,IL),ELEV(IC,IR,IL))
c ------APPLY INFILTRATION FACTOR ! RFACT
if (IRIVRFACT.gt.0) then ! RFACT
if (HTEMP.LE.HRIV) then ! RFACT
! situation with infiltration, apply infiltration factor ! RFACT
CRIV=CRIV*RIVR(IRIVRFACT,L) ! RFACT
endif ! RFACT
endif ! RFACT
C--SEAWAT: CALCULATE FRESHWATER HEAD AT RBOT USING HEAD IN MODEL CELL
HFRBOT=HHNEW+(PS(IC,IR,IL)-DENSEREF)/DENSEREF*
& (ELEV(IC,IR,IL)-RBOT)
C
C6------COMPARE AQUIFER HEAD TO BOTTOM OF STREAM BED.
C--SEAWAT: USE SALTHEAD FOR THIS COMPARISON
IF(HTEMP.LE.RRBOT)GO TO 96 !SALTHEAD(HNEW(IC,IR,IL),PS(IC,IR,IL),ELEV(IC,IR,IL))
! + .LE.RRBOT) GO TO 96
RHOAVG=(RIVDENS+PS(IC,IR,IL))/2
DIRECT=-((HRIV-HFRBOT)+(PS(IC,IR,IL)-DENSEREF)/
+ DENSEREF*RBDTHK)
C--SEAWAT: DIRECT IS POSITIVE, FLOW IS UP INTO RIVER
IF (DIRECT.GT.0.) RIVDENS=PS(IC,IR,IL)
C7------SINCE HEAD>BOTTOM ADD TERMS TO RHS AND HCOF.
C--SEAWAT: USE VARIABLE DENSITY FORM OF DARCY'S LAW, CONSERVE MASS
RHS(IC,IR,IL)=RHS(IC,IR,IL)-CRIV*RIVDENS*
& (HRIV-(PS(IC,IR,IL)-DENSEREF)/DENSEREF*
& (ELEV(IC,IR,IL)-RBOT)+(RHOAVG-DENSEREF)/
+ DENSEREF*RBDTHK)
HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-CRIV*RIVDENS
GO TO 100
C8------SINCE HEADGWFRIVDAT(IGRID)%NRIVER
MXRIVR=>GWFRIVDAT(IGRID)%MXRIVR
NRIVVL=>GWFRIVDAT(IGRID)%NRIVVL
IRIVCB=>GWFRIVDAT(IGRID)%IRIVCB
IPRRIV=>GWFRIVDAT(IGRID)%IPRRIV
NPRIV=>GWFRIVDAT(IGRID)%NPRIV
IRIVPB=>GWFRIVDAT(IGRID)%IRIVPB
NNPRIV=>GWFRIVDAT(IGRID)%NNPRIV
RIVAUX=>GWFRIVDAT(IGRID)%RIVAUX
RIVR=>GWFRIVDAT(IGRID)%RIVR
IRIVRFACT=>GWFRIVDAT(IGRID)%IRIVRFACT ! RFACT
irivsubsys=>gwfrivdat(igrid)%irivsubsys ! rsubsys
nrivsubsys=>gwfrivdat(igrid)%nrivsubsys ! rsubsys
rivsubsidx=>gwfrivdat(igrid)%rivsubsidx ! rsubsys
irivrconc=>gwfrivdat(igrid)%irivrconc ! rconc
ifvdl=>gwfrivdat(igrid)%ifvdl ! ifvdl
isft=>gwfrivdat(igrid)%isft ! ifvdl
sft=>gwfrivdat(igrid)%sft ! ifvdl
C
RETURN
END
SUBROUTINE SVDF2RIV7PSV(IGRID)
C Save river data for a grid.
USE GWFRIVMODULE
C
GWFRIVDAT(IGRID)%NRIVER=>NRIVER
GWFRIVDAT(IGRID)%MXRIVR=>MXRIVR
GWFRIVDAT(IGRID)%NRIVVL=>NRIVVL
GWFRIVDAT(IGRID)%IRIVCB=>IRIVCB
GWFRIVDAT(IGRID)%IPRRIV=>IPRRIV
GWFRIVDAT(IGRID)%NPRIV=>NPRIV
GWFRIVDAT(IGRID)%IRIVPB=>IRIVPB
GWFRIVDAT(IGRID)%NNPRIV=>NNPRIV
GWFRIVDAT(IGRID)%RIVAUX=>RIVAUX
GWFRIVDAT(IGRID)%RIVR=>RIVR
GWFRIVDAT(IGRID)%IRIVRFACT=>IRIVRFACT ! RFACT
gwfrivdat(igrid)%irivsubsys=>irivsubsys ! rsubsys
gwfrivdat(igrid)%nrivsubsys=>nrivsubsys ! rsubsys
gwfrivdat(igrid)%rivsubsidx=>rivsubsidx ! rsubsys
gwfrivdat(igrid)%irivrconc=>irivrconc ! rconc
gwfrivdat(igrid)%ifvdl=>ifvdl ! ifvdl
gwfrivdat(igrid)%isft=>isft ! ifvdl
gwfrivdat(igrid)%sft=>sft ! ifvdl
C
RETURN
END