search

rniswon / mfnwt

MODFLOW-NWT - Version: 1.1.0 Newton Formulation for MODFLOW-2005 For more information go to http://water.usgs.gov/ogw/modflow-nwt
4 stars 7 forks source link

CHD (and possibly FHB) called after LMT header values set #42

Open emorway-usgs opened 5 years ago

emorway-usgs commented 5 years ago

This issue was originally reported by @rbwinst-usgs:

"I have been trying to reproduce the UZT_Disp_Lamb01_TVD example in MT3D-USGS, I ran into a possible bug in lmt8_NWT.f. When I looked at the ftl file created by my version of the example, I noticed that the value printed for "MTCHD" was 0 for my version of the model even though I had defined a CHD cell via the CHD package. In

The code for setting MTCHD on lines 487-494 is as follows.

MTCHD=0    !loop through the entire grid to get
DO K=1,NLAY    !total number of constant-head cells
  DO I=1,NROW
    DO J=1,NCOL
      IF(IBOUND(J,I,K).LT.0) MTCHD=MTCHD+1
    ENDDO
  ENDDO
ENDDO

In your version of the model the CHD package is not used and the specified head cell is specified via the IBOUND array. I'm not sure whether or not this makes a difference."

emorway-usgs commented 5 years ago

Further down, in LMT8UPW1(), the following occurs on (or near) lines 1798-1805:

C--FOR EACH CELL IF IT IS CONSTANT HEAD COMPUTE FLOW ACROSS 6
C--FACES.
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
C
C--IF CELL IS NOT CONSTANT HEAD SKIP IT & GO ON TO NEXT CELL.
            IF(IBOUND(J,I,K).GE.0) CYCLE
            NCNH=NCNH+1

But at this point, IBOUND has been adjusted by CHD, which is to say that MTCHD in the previous post was wrong. These screen captures illustrate the issue: Header: header Flow terms: cnh

One solution might be to change the order of function calls in MF_NWT.f so that CHD (and possibly FHB) has had the opportunity to modify IBOUND and the header will be set correctly in MT3D-USGS. For example, Line 157 of MF_NWT.f is:

IF(IUNIT(49).GT.0) CALL LMT8BAS7AR(INUNIT,CUNIT,IGRID)

But IBOUND isn't modified by CHD until the call to GWF2CHD7RP, which doesn't happen until line 198, which is inside the stress period loop of the main code.

One solution might be to move LMT8BAS7AR() to line 541 of the main code and wrap it in an if statement. So this 

CLMT----CALL LINK-MT3DMS SUBROUTINES TO SAVE FLOW-TRANSPORT LINK FILE
CLMT----FOR USE BY MT3DMS FOR TRANSPORT SIMULATION
CLMT
          IF(IUNIT(49).GT.0) CALL LMT8BD(KKSTP,KKPER,IGRID)
CLMT

would be come

CLMT----CALL LINK-MT3DMS SUBROUTINES TO SAVE FLOW-TRANSPORT LINK FILE
CLMT----FOR USE BY MT3DMS FOR TRANSPORT SIMULATION
CLMT
          IF(IUNIT(49).GT.0) THEN
            IF(KPER.EQ.1.AND.KSTP.EQ.1) THEN 
              CALL LMT8BD(KKSTP,KKPER,IGRID)
            END IF
            CALL LMT8BD(KKSTP,KKPER,IGRID)
          END IF
CLMT

@langevin-usgs, I suspect this is happening in MF2005 as well.