search

MODFLOW-USGS / mt3d-usgs

MT3D-USGS Repository
23 stars 12 forks source link

'DRYCELL' keyword in BTN causes mass accounting error in SSM1BD() #65

Closed emorway-usgs closed 5 years ago

emorway-usgs commented 5 years ago

@vivekbedekar The first chunk of code below is from lines 1206-1250 in SSM. Note that when DRYON is .TRUE. as a result of using the DRYCELL keyword in BTN, MT3D-USGS will drop into the highest-level ELSEIF(..) statement for inactive cells and by virtue of DRYON will proceed to calculate contributions to the overall budget terms [RMASIO(8,x,icomp)] using whatever value of CTMP the user has provided. For cases when the user supplies a non-zero value for CINACT, this is a problem.

A test model from Milo Simonic demonstrates the issue. In that model, the following conditions exist in a subset of cells:

And as a result, RMASIO(8,2,icomp) is non-zero, which is wrong. The second chunk of code below, from the original MT3DMS source code, shows that historically a cycle would occur owing to ICBUND(J,I,K,ICOMP)==0 for these cells. I think the fix is to check if CEVT(J,I,ICOMP)==0 immediately after the IF(DRYON) THEN and if so, CYCLE. Please confirm.

      DO I=1,NROW
        DO J=1,NCOL
          K=IEVT(J,I)
          IF(K.EQ.0) CYCLE
          IF(ICBUND(J,I,K,ICOMP).GT.0) THEN
            CTMP=CEVT(J,I,ICOMP)
            IF(EVTR(J,I).LT.0.AND.(CTMP.LT.0 .or.
     &                             CTMP.GE.CNEW(J,I,K,ICOMP))) THEN
              CTMP=CNEW(J,I,K,ICOMP)
            ELSEIF(CTMP.LT.0) THEN        
              CTMP=0.
            ENDIF
            IF(EVTR(J,I).GT.0) THEN
              RMASIO(8,1,ICOMP)=RMASIO(8,1,ICOMP)+EVTR(J,I)*CTMP*DTRANS*
     &                          DELR(J)*DELC(I)*DH(J,I,K)
            ELSE
              RMASIO(8,2,ICOMP)=RMASIO(8,2,ICOMP)+EVTR(J,I)*CTMP*DTRANS*
     &                          DELR(J)*DELC(I)*DH(J,I,K)
            ENDIF
C
          ELSEIF(ICBUND(J,I,K,ICOMP).EQ.0) THEN
            IF(DRYON) THEN
              CTMP=CEVT(J,I,ICOMP)
              IF(EVTR(J,I).LT.0.AND.(CTMP.LT.0 .or.
     &                         CTMP.GE.CNEW(J,I,K,ICOMP))) THEN
                CTMP=CNEW(J,I,K,ICOMP)
              ELSEIF(CTMP.LT.0) THEN        
                CTMP=0.
              ENDIF
              VOLAQU=DELR(J)*DELC(I)*DH(J,I,K)
              IF(ABS(VOLAQU).LE.1.E-8) VOLAQU=1.E-8
              IF(EVTR(J,I).LT.0) THEN
                RMASIO(8,2,ICOMP)=RMASIO(8,2,ICOMP)+EVTR(J,I)*CTMP*
     &                            DTRANS*ABS(VOLAQU)
                QC7(J,I,K,9)=QC7(J,I,K,9)-EVTR(J,I)*ABS(VOLAQU)
              ELSE
                RMASIO(8,1,ICOMP)=RMASIO(8,1,ICOMP)+EVTR(J,I)*CTMP*
     &                            DTRANS*ABS(VOLAQU)
                QC7(J,I,K,7)=QC7(J,I,K,7)-EVTR(J,I)*ABS(VOLAQU)*CTMP
                QC7(J,I,K,8)=QC7(J,I,K,8)-EVTR(J,I)*ABS(VOLAQU)
              ENDIF
            ENDIF
          ENDIF
        ENDDO
      ENDDO

Original MT3DMS code for this section:

C
C--(EVAPOTRANSPIRATION)
  100 IF(.NOT.FEVT .AND. .NOT.FETS) GOTO 200
C
      DO I=1,NROW
        DO J=1,NCOL
          K=IEVT(J,I)
          IF(K.EQ.0 .OR. ICBUND(J,I,K,ICOMP).LE.0) CYCLE
          CTMP=CEVT(J,I,ICOMP)
          IF(EVTR(J,I).LT.0.AND.(CTMP.LT.0 .or.
     &                           CTMP.GE.CNEW(J,I,K,ICOMP))) THEN
            CTMP=CNEW(J,I,K,ICOMP)
          ELSEIF(CTMP.LT.0) THEN        
            CTMP=0.
          ENDIF
          IF(EVTR(J,I).GT.0) THEN
            RMASIO(8,1,ICOMP)=RMASIO(8,1,ICOMP)+EVTR(J,I)*CTMP*DTRANS*
     &       DELR(J)*DELC(I)*DH(J,I,K)
          ELSE
            RMASIO(8,2,ICOMP)=RMASIO(8,2,ICOMP)+EVTR(J,I)*CTMP*DTRANS*
     &       DELR(J)*DELC(I)*DH(J,I,K)
          ENDIF
        ENDDO
      ENDDO
C
emorway-usgs commented 5 years ago

Passed on Travis, closing. Also, Vivek confirmed suggested fix by way of email.