IVSD and IVC for NLAY=1

langevin-usgs commented 7 years ago

I thought that in mfusg if we had a model with cells that were overlying and underlying one another, that we could put all the cells into one layer, and just set IVC appropriately. I think I had done that in the past. But a colleague in the Netherlands, just pointed out that IVC is only read if IVSD = 1 AND NLAY > 1. Code is below. Is this a new requirement? Is there any reason we don't read IVC in this case? The NLAY > 1 requirement does not seem to be documented in the instructions.

C10-----FILL VERTICAL CONNECTION ARRAY IVC.
      IF(IVSD.LE.0)THEN
C10A------compute IVC for IVSD. LE. 0
        DO K=1,NLAY
          NNDLAY = NODLAY(K)
          NSTRT = NODLAY(K-1)+1
          DO N=NSTRT,NNDLAY
C10A1-----LOOP OVER CONNECTIONS OF NODE N AND FILL
          DO II = IA(N)+1,IA(N+1)-1
            JJ = JA(II)
            IIS = JAS(II)
            IF(JJ.LE.N.OR.JJ.GT.NODES) CYCLE
            IF(JJ.GT.NNDLAY)THEN
              IVC(IIS) = 1  ! LAYER IS BELOW
            ELSE
              IVC(IIS) = 0
            ENDIF
          ENDDO
          ENDDO
        ENDDO
      ELSE
C10B----read IVC, for IVSD. GT. 0.
        IF(NLAY.GT.1) 
     *  CALL U1DINTNJA(IVC,IATMP,ANAME,NJATMP,INDIS,IOUT,IDSYMRD)
      ENDIF

sorab commented 7 years ago

Sorry. Get rid of it. I am doing it here on USG-Beta too.

jdhughes-usgs commented 6 years ago

It does not seem to be as easy. Doesn't look like vertical flow is calculated if NLAY is 1. Also there are some array bounds issues in:

C
C5A-------IF THERE IS NO FLOW TO CALCULATE THROUGH THIS FACE, THEN GO ON
C5A-------TO NEXT FACE.
          IF(ICHFLG.EQ.0) THEN
            IF((IBOUND(N).LE.0) .AND. (IBOUND(JJ).LE.0)) GO TO 30
          ELSE
            IF(IBOUND(JJ).EQ.0) GO TO 30
          END IF
          HD=HNEW(JJ)
          IF(IVC(IIS).EQ.1.AND.JJ.GT.N)THEN !VERTICAL DIRECTION DOWN
            IF(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 122 !CSP 2 SHOULD BE 1
            TMP=HD
            IF(NOVFC.EQ.0)THEN
              IF(TMP.LT.TOP(JJ)) HD=TOP(JJ)
            ENDIF
          ENDIF

with LAYCON(K+1).

sorab commented 5 years ago

The LAYCON(K+1) issue is also readily solved (done) but not sure how this will affect LAKE package where NALY.eq.1 is also considered in several places.

sorab commented 5 years ago

We could document that LAKE package may want the layers separated and close the issue.

MODFLOW-USGS / mfusg

IVSD and IVC for NLAY=1 #29