c Copyright (C) Stichting Deltares, 2005-2024. 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 VDF2RCH7FM(IGRID) C ****************************************************************** C SUBTRACT RECHARGE FROM RHS C--SEAWAT: REFORMULATED FOR VD FLOW C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY:NCOL,NROW,NLAY,IBOUND,RHS USE GWFRCHMODULE,ONLY:NRCHOP,RECH,IRCH USE VDFMODULE, ONLY:DENSEREF,PS C ------------------------------------------------------------------ C C1------SET POINTERS FOR THE CURRENT GRID. CALL SGWF2RCH7PNT(IGRID) C C2------DETERMINE WHICH RECHARGE OPTION. IF(NRCHOP.EQ.1) THEN C C3------NRCHOP IS 1, SO RECHARGE IS IN TOP LAYER. LAYER INDEX IS 1. DO 10 IR=1,NROW DO 10 IC=1,NCOL C C3A-----IF CELL IS VARIABLE HEAD, SUBTRACT RECHARGE RATE FROM C3A-----RIGHT-HAND-SIDE. IF(IBOUND(IC,IR,1).GT.0) THEN DENSE=DENSEREF IF(RECH(IC,IR).LT.0) DENSE=PS(IC,IR,1) C C1B-----SUBTRACT RECHARGE RATE FROM RIGHT-HAND-SIDE. C--SEAWAT: CONSERVE MASS C RHS(IC,IR,1)=RHS(IC,IR,1)-RECH(IC,IR) RHS(IC,IR,1)=RHS(IC,IR,1)-RECH(IC,IR)*DENSE ENDIF 10 CONTINUE ELSE IF(NRCHOP.EQ.2) THEN C C4------NRCHOP IS 2, SO RECHARGE IS INTO LAYER IN INDICATOR ARRAY DO 20 IR=1,NROW DO 20 IC=1,NCOL IL=IRCH(IC,IR) C C4A-----IF THE CELL IS VARIABLE HEAD, SUBTRACT RECHARGE FROM C4A-----RIGHT-HAND-SIDE. IF(IL.EQ.0) GO TO 20 IF(IBOUND(IC,IR,IL).GT.0)THEN DENSE=DENSEREF IF(RECH(IC,IR).LT.0) DENSE=PS(IC,IR,IL) C C2C-----SUBTRACT RECHARGE FROM RIGHT-HAND-SIDE. C--SEAWAT:CONSERVE MASS C RHS(IC,IR,IL)=RHS(IC,IR,IL)-RECH(IC,IR) RHS(IC,IR,IL)=RHS(IC,IR,IL)-RECH(IC,IR)*DENSE ENDIF 20 CONTINUE ELSE C C5------NRCHOP IS 3, RECHARGE IS INTO HIGHEST VARIABLE-HEAD CELL, EXCEPT C5------CANNOT PASS THROUGH CONSTANT HEAD NODE DO 30 IR=1,NROW DO 30 IC=1,NCOL DO 28 IL=1,NLAY C C5A-----IF CELL IS CONSTANT HEAD MOVE ON TO NEXT HORIZONTAL LOCATION. IF(IBOUND(IC,IR,IL).LT.0) GO TO 30 C C5B-----IF THE CELL IS VARIABLE HEAD, SUBTRACT RECHARGE FROM C5B-----RIGHT-HAND-SIDE AND MOVE TO NEXT HORIZONTAL LOCATION. IF(IBOUND(IC,IR,IL).GT.0) THEN DENSE=DENSEREF IF(RECH(IC,IR).LT.0) DENSE=PS(IC,IR,IL) C C3C-----SUBTRACT RECHARGE FROM RIGHT-HAND-SIDE. C--SEAWAT:CONSERVE MASS C RHS(IC,IR,IL)=RHS(IC,IR,IL)-RECH(IC,IR) RHS(IC,IR,IL)=RHS(IC,IR,IL)-RECH(IC,IR)*DENSE GO TO 30 END IF 28 CONTINUE 30 CONTINUE END IF C C6------RETURN RETURN END SUBROUTINE VDF2RCH7BD(KSTP,KPER,IGRID) C ****************************************************************** C CALCULATE VOLUMETRIC BUDGET FOR RECHARGE C--SEAWAT: USE VARIABLE-DENSITY EQUATIONS C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY:IOUT,NCOL,NROW,NLAY,IBOUND,BUFF USE GWFBASMODULE,ONLY:MSUM,VBVL,VBNM,ICBCFL,DELT,PERTIM,TOTIM USE GWFRCHMODULE,ONLY:NRCHOP,IRCHCB,RECH,IRCH USE VDFMODULE, ONLY: DENSEREF,PS C DOUBLE PRECISION RATIN,RATOUT,QQ CHARACTER*16 TEXT DATA TEXT /' RECHARGE'/ C ------------------------------------------------------------------ C C1------SET POINTERS FOR THE CURRENT GRID. CALL SGWF2RCH7PNT(IGRID) C C2------CLEAR THE RATE ACCUMULATORS. ZERO=0. RATIN=ZERO RATOUT=ZERO C C3------CLEAR THE BUFFER & SET FLAG FOR SAVING CELL-BY-CELL FLOW TERMS. DO 2 IL=1,NLAY DO 2 IR=1,NROW DO 2 IC=1,NCOL BUFF(IC,IR,IL)=ZERO 2 CONTINUE IBD=0 IF(IRCHCB.GT.0) IBD=ICBCFL C C4------DETERMINE THE RECHARGE OPTION. IF(NRCHOP.EQ.1) THEN C C5------NRCHOP=1, SO RECH GOES INTO LAYER 1. PROCESS EACH HORIZONTAL C5------CELL LOCATION. DO 10 IR=1,NROW DO 10 IC=1,NCOL C C5A-----IF CELL IS VARIABLE HEAD, THEN DO BUDGET FOR IT. IF(IBOUND(IC,IR,1).GT.0) THEN DENSE=DENSEREF IF(RECH(IC,IR).LT.0) DENSE=PS(IC,IR,1) C--SEAWAT:CONSERVE MASS C Q=RECH(IC,IR) Q=RECH(IC,IR)*DENSE QQ=Q C C5B-----ADD RECH TO BUFF. C--SEAWAT:CONVERT TO VOLUMETRIC FLUX FOR OUTPUT C BUFF(IC,IR,1)=Q BUFF(IC,IR,1)=Q/DENSE C C5C-----IF RECH POSITIVE ADD IT TO RATIN, ELSE ADD IT TO RATOUT. IF(Q.GE.ZERO) THEN RATIN=RATIN+QQ ELSE RATOUT=RATOUT-QQ END IF END IF 10 CONTINUE ELSE IF(NRCHOP.EQ.2) THEN C C6------NRCHOP=2, RECH IS IN LAYER SPECIFIED IN INDICATOR ARRAY(IRCH). C6------PROCESS EACH HORIZONTAL CELL LOCATION. DO 20 IR=1,NROW DO 20 IC=1,NCOL C C6A-----GET LAYER INDEX FROM INDICATOR ARRAY(IRCH). IL=IRCH(IC,IR) C C6B-----IF CELL IS VARIABLE HEAD, THEN DO BUDGET FOR IT. IF(IL.EQ.0) GO TO 20 IF(IBOUND(IC,IR,IL).GT.0) THEN DENSE=DENSEREF IF(RECH(IC,IR).LT.0) DENSE=PS(IC,IR,IL) C--SEAWAT:CONSERVE MASS C Q=RECH(IC,IR) Q=RECH(IC,IR)*DENSE QQ=Q C C6C-----ADD RECHARGE TO BUFF. C--SEAWAT: CONVERT TO VOLUMETRIC FLUX C BUFF(IC,IR,IL)=Q BUFF(IC,IR,IL)=Q/DENSE C C6D-----IF RECHARGE IS POSITIVE ADD TO RATIN, ELSE ADD IT TO RATOUT. IF(Q.GE.ZERO) THEN RATIN=RATIN+QQ ELSE RATOUT=RATOUT-QQ END IF END IF 20 CONTINUE ELSE C C7------NRCHOP=3; RECHARGE IS INTO HIGHEST CELL IN A VERTICAL COLUMN C7------THAT IS NOT NO FLOW. PROCESS EACH HORIZONTAL CELL LOCATION. DO 30 IR=1,NROW DO 29 IC=1,NCOL C C7A-----INITIALIZE IRCH TO 1, AND LOOP THROUGH CELLS IN A VERTICAL C7A-----COLUMN TO FIND WHERE TO PLACE RECHARGE. IRCH(IC,IR)=1 DO 28 IL=1,NLAY C C7B-----IF CELL IS CONSTANT HEAD, MOVE ON TO NEXT HORIZONTAL LOCATION. IF(IBOUND(IC,IR,IL).LT.0) GO TO 29 C C7C-----IF CELL IS VARIABLE HEAD, THEN DO BUDGET FOR IT. IF (IBOUND(IC,IR,IL).GT.0) THEN C--SEAWAT: SET DENSE DENSE=DENSEREF IF(RECH(IC,IR).LT.0) DENSE=PS(IC,IR,IL) C--SEAWAT:CONSERVE MASS C Q=RECH(IC,IR) Q=RECH(IC,IR)*DENSE QQ=Q C--SEAWAT:CONVERT TO VOLUMETRIC FLUX C BUFF(IC,IR,IL)=Q BUFF(IC,IR,IL)=Q/DENSE IRCH(IC,IR)=IL C C7E-----IF RECH IS POSITIVE ADD IT TO RATIN, ELSE ADD IT TO RATOUT. IF(Q.GE.ZERO) THEN RATIN=RATIN+QQ ELSE RATOUT=RATOUT-QQ END IF GO TO 29 END IF 28 CONTINUE 29 CONTINUE 30 CONTINUE C END IF C C8------IF CELL-BY-CELL FLOW TERMS SHOULD BE SAVED, CALL APPROPRIATE C8------UTILITY MODULE TO WRITE THEM. 100 IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IRCHCB,BUFF,NCOL,NROW, 1 NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV3(KSTP,KPER,TEXT,IRCHCB,BUFF,IRCH,NRCHOP, 1 NCOL,NROW,NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND) C C9------MOVE TOTAL RECHARGE RATE INTO VBVL FOR PRINTING BY BAS1OT. ROUT=RATOUT RIN=RATIN VBVL(4,MSUM)=ROUT VBVL(3,MSUM)=RIN C C10-----ADD RECHARGE FOR TIME STEP TO RECHARGE ACCUMULATOR IN VBVL. VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT VBVL(1,MSUM)=VBVL(1,MSUM)+RIN*DELT C C11-----MOVE BUDGET TERM LABELS TO VBNM FOR PRINT BY MODULE BAS_OT. VBNM(MSUM)=TEXT C C12-----INCREMENT BUDGET TERM COUNTER. MSUM=MSUM+1 C C13-----RETURN RETURN END SUBROUTINE VDF2RCH7DA(IGRID) C Deallocate RCH DATA USE GWFRCHMODULE C DEALLOCATE(GWFRCHDAT(IGRID)%NRCHOP) DEALLOCATE(GWFRCHDAT(IGRID)%IRCHCB) DEALLOCATE(GWFRCHDAT(IGRID)%NPRCH) DEALLOCATE(GWFRCHDAT(IGRID)%IRCHPF) DEALLOCATE(GWFRCHDAT(IGRID)%RECH) if(associated(gwfrchdat(igrid)%rechbuff)) ! DLT 1 deallocate(gwfrchdat(igrid)%rechbuff) ! DLT DEALLOCATE(GWFRCHDAT(IGRID)%IRCH) DEALLOCATE(GWFRCHDAT(IGRID)%IADDRECH) ! DLT C RETURN END SUBROUTINE SVDF2RCH7PNT(IGRID) C Set RCH pointers for grid. USE GWFRCHMODULE C NRCHOP=>GWFRCHDAT(IGRID)%NRCHOP IRCHCB=>GWFRCHDAT(IGRID)%IRCHCB NPRCH=>GWFRCHDAT(IGRID)%NPRCH IRCHPF=>GWFRCHDAT(IGRID)%IRCHPF RECH=>GWFRCHDAT(IGRID)%RECH RECHBUFF=>GWFRCHDAT(IGRID)%RECHBUFF ! DLT IRCH=>GWFRCHDAT(IGRID)%IRCH IADDRECH=>GWFRCHDAT(IGRID)%IADDRECH ! DLT C RETURN END SUBROUTINE SVDF2RCH7PSV(IGRID) C Save RCH pointers for grid. USE GWFRCHMODULE C GWFRCHDAT(IGRID)%NRCHOP=>NRCHOP GWFRCHDAT(IGRID)%IRCHCB=>IRCHCB GWFRCHDAT(IGRID)%NPRCH=>NPRCH GWFRCHDAT(IGRID)%IRCHPF=>IRCHPF GWFRCHDAT(IGRID)%RECH=>RECH GWFRCHDAT(IGRID)%RECHBUFF=>RECHBUFF ! DLT GWFRCHDAT(IGRID)%IRCH=>IRCH GWFRCHDAT(IGRID)%IADDRECH=>IADDRECH ! DLT C RETURN END