modflow-usgs / mt3d-usgs Goto Github PK

When LMT writes the linker file, it calculates the cross-sectional area of the flow. When the user specifies a negative value for FLOW (a direct withdraw from the channel), it was crashing MF-NWT.

Weave the MF6 zoneBudget.f code into MT3D-USGS for supporting structured grid transport with MF6. At Chris's suggestion, also need to consider reading the binary grid file generated by MF6. Need to find out more about this from Chris, but for now, need to start with simple example and read that into existing MT3D-USGS arrays.

The code originally introduced in commit 141f400 negates the constant concentration boundary condition.

C--RESET ICBUND TO ABS(ICBUND) FOR CONSTANT CONCENTRATION CELLS
      DO INDEX=1,NCOMP
        DO KK=1,NLAY
          DO II=1,NROW
            DO JJ=1,NCOL
              IF(ICBUND(JJ,II,KK,INDEX).LT.0) THEN
                ICBUND(JJ,II,KK,INDEX)=ABS(ICBUND(JJ,II,KK,INDEX))
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ENDDO  

None of our tests used this boundary condition and so it went undetected until now. In trying to determine why this code was added, did it have something to do with transport through dry cells that can happen with the use of MF-NWT?

In cell 1, 211, 149, temperature does the following:

Mt3D-USGS encounters an error when reading SFR data from an unformated ftl file. The error occurs in the following lines of READPS.
C--MARK CELLS NEAR THE SINK/SOURCE
IF(QSTEMP.LT.0 .AND. ICBUND(J,I,K).GT.0) THEN
ICBUND(J,I,K)=1000+IQ
ELSEIF(ICBUND(J,I,K).GT.0 ) THEN
ICBUND(J,I,K)=1020+IQ
ENDIF
The error occurs because when reading SFR data, there are two real values to be read for each cell but only one of them is actually read. For example, here is the beginning of the relevant portions of the ftl file when it is saved as a formatted file.
SFR 36
1 1 1 1.028048 4500.000
1 2 2 2.706948 7000.000

"IF()" statements starting with k.ne.0.and... need to be broken apart

An email from Paul Hsieh:
"The problem seems to occur when sorption is involved in a multi-layer model. I created some runs to illustrate this. The set up consists of 4 layers, 1 row, 10 columns, with flow occurring along the row from left to right. The set ups were created with ModelMuse.

When there is no sorption, both mt3d-usgs and mt3dms give the same result.

However, when sorption (linear isotherm) is turned on, mt3d-usgs does not appear to correctly read the sorption parameters (bulk density and Kd), so the retardation factor is not correctly computed. (See the listing file example.mls in the section "SORPTION AND 1ST/0TH ORDER REACTION PARAMETERS". The retardation factor should be 1.8 in all layers.)

Running the same sorption problem using mt3dms appears to give the correct result."

Email from Mike Toews on 7/25/17:
"Hi Eric and Vivek,

I've been looking to improve the speed of MT3D-USGS, and I've attached a patch of suggested changes, which seem to compile fine with GNU, PGI and Intel compilers. Here are some notes:

Changes to adv1.f: these are preliminary, and I'm not sure if they have any improvement. I was focusing on tidying up the variable declarations as I've found that the function CFACE and subroutine SADV1U have a high number of calls when using the TVD scheme, and are a focus of future improvement.

Changes to gcg1.f: there are several improvements. The most effective was to move LU and LL to mt_module.f90 as an integer parameter, as these do not change, and don't need to continuously be set. Other changes are to reverse the order of DO loops to work with Fortran's row-major array ordering. For instance, around Line 783. I've changed L to be a fixed array with 19 items, as it does not need to be a dynamic array with allocation/deallocation. Other changes are to use modern variable declarations for MVPRD, MTVPRD, QSOLVE and MIC, which generally help them read better. These subroutines also have a high number of calls while running any type of model.
I may look at other improvements in the near future, and will pass these along. Please review and take any changes you feel are improvements. Also, if you need the patch in a different format, or whatever, let me know."
mt_optim1.txt

A recent test model with LKT active but not SSM was attempting to run but should have exited with error message that already exists in MT3D-USGS ERROR: THE SSM PACKAGE MUST BE USED IN THE CURRENT SIMULATION BECAUSE THE FLOW MODEL INCLUDES A SINK/SOURCE PACKAGE."

The documentation of NTMP in LKT states that NTMP must be greater than zero for the first stress period. However, in the example model for LKT NTMP is zero.

James Rumbaugh [email protected]
10:55 AM (15 minutes ago)

to me
Hi Chris,

The data input guide for UZT says that IUZFBND array is to be read for each layer. The code seems to only read one layer of data. I assume the guide is incorrect?

Thanks

Jim

Also in UZT, the data input guide says that records 3, 4, and 5 are read free format. However the code uses U2DINT to read #3 (IUZFBND) and RARRAY for the other 2.

In Dec. 2018, we looked at SSM FMI and thought that multiple BCs in the same model cell may not be handled separately.

Noticed that the From-To node listing wasn't as expected in the linker file. Note the 3rd line (in context) in some example lines from linker file below:

36 37 1 55692.22 1.208313
37 38 1 59427.54 1.283066
38 39 1 65349.79 1.370677
-999 40 0 -4.108000 6.5532312E-02
40 41 1 502.3974 9.3378410E-02
41 42 1 2249.491 0.2373254

As currently constructed, the line looks like it might work, but need to verify. Based on the node order, it is suggesting flow is from the sink (-999) to the first reach of iseg 4 (the 40th item in the REACHINPUT list), but because the flow is negative the order is actually reversed. Just want to make sure this is OK.

The following error occurs in the attached model.
Transport Step: 4 Step Size: 2.0688E+05 Total Elapsed Time: 2.62800E+06
Outer Iter. 1 Inner Iter. 1: Max. DC = 0.8551E-01 [K,I,J] 1 3 2
Outer Iter. 1 Inner Iter. 2: Max. DC = 0.4885E-14 [K,I,J] 1 2 3

STRESS PERIOD NO. 2

TIME STEP NO. 1
FROM TIME = 0.26280E+07 TO 0.27107E+07

ERROR READING FLOW-TRANSPORT LINK FILE.
NAME OF THE FLOW TERM REQUIRED =THKSAT
NAME OF THE FLOW TERM SAVED IN FLOW-TRANSPORT LINK FILE =÷O╟╛ N!½┐
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
Uzf_Import3.zip

Hey guys - I'm having some trouble getting a decent mass balance from a model if use RCT and the zeroth order reaction (IREACT=100). First order good, no reaction good, production ( negative rc1) good, but decay (positive rc1) bad. I created a simple 3-cell model that demonstrates the issue. Please let me know if this is some sort of user/input error.
rct_trouble.zip

After MODFLOW & More 2017, Eric Chiang with Simcore Software (PMWIN software) offered the edits contained in the attached document.
Input_Instructions_Eric_Chiang.pdf

I was getting divide by zero for itype 15 for a cell that had gone dry (volaqu = 0.0). So I added a check reset volaqu to 1.0e-5 for this case. Seems to work for me, but not sure if this has other implications...

Eric Lappala (formerly with the USGS, I believe) wrote an email addressed to Eric and Vivek on 10/24/16 that has the following recommendation:

"Why does the initial saturated thickness have to be entered separately in the dataset read by the UZT routine, which according to the published input instructions posted with MT3DUSGS have to be entered again when this information is present in the link transport file? Seems to me that the UZT routine could have code added that if the saturated thickness data is present in the link transport file that it does not have to be read again. Although moisture contents above the water table are not presently written to the link transport file by Modflow NWT, it seems that it could be and then be treated the same suggested way by the UZT package."

We should be able to accommodate this request, I would think.

In the UZT test problems involving a vertical column (in particular the problem titled 'UZT_Disp_Lamb01_TVD'), the way mass is routed/removed in the saturated zone in the bottom 2 layers of the model appears to be asymmetric even though the flow solution shows symmetry. The percent discrepancy is low enough that it isn't a problem, however, when more time is available might be worth investigating further. In looking through the source code on 9/22/16 for roughly an hour (Vivek & Eric), nothing obvious stood out, will likely need to ferret out the issue by tracking values in the A matrix.

Notes for revisiting the probelm: monitor the concentrations in cells i,j,k=(1,1,210) and (1,3,210) for example, there is not reason the concentration should be different in these cells based on what is in the flow solution. Also, the concentrations will only be different when the slug of mass moving through the system arrives at depth.

Eric,

We need to verify if the SATOLD calculation for the first time-step is correct.

Vivek

We are using xMD for solving SFT. xMD code has its own copyright info in the comments. Is that sufficient?

Using a model J White sent with 12 different instances of negative values specified for FLOW in the SFR input, only 2 are making it through to the linker file. Track down why the other 10 aren't showing up.

When users specify a negative value for the 'FLOW' variable (Dataset 4b in the Online Guide to MF-NWT, for example), the LMT package was writing this value to the linker file as a sink term, but wasn't checking to make sure it was capped by the total available flow in the network at that location.

Hey guys - I'm working nitrate transport model and I need a way to "inject" some nitrate directly into the SFR network - there are point sources of nitrate and we have estimates of the load (mass per time) for these locations along the SFR network. Is there anyway to get this done given the current version of SFT? I saw some options for specified concentration, but nothing like an itype 15. Any suggestions? I guess I could "fake" an sfr runoff and scale the concentration to get the desired mass, but I'm hoping for something cleaner...Thanks!

When READPS attempts to read the data for "UZF RECHARGE" from a formatted ftl file, it doesn't read it properly. It only read "UZF" instead of "UZF RECHARGE"
mt3d_bug.zip

From Kurt Zeiler with Brown & Caldwell on 10/30/2017:

"Hey Eric & Vivek,

I hope all is well with both of you. We stumbled across a bug in MT3D-USGS for TVD with DRYCELL. There are 2 spots where DELR is being indexed by rows instead of columns (i.e., by I instead of J) in adv1.f – both left and right face calculations of CMAS2 along x-direction (lines 2721 and 2787). Here’s the corrected version of those lines:

CMAS2=-COLD(J,I,K)*QX(J-1,I,K)/DELR(J)           !kkz - was DELR(I) previously

We haven’t fully tested yet (just got new run started 15 min ago), but certainly seems more correct. Also, it doesn’t crash for our model that has more rows than columns. 😊

Please let me know if this doesn’t make sense or the fix above isn’t correct.

Kurt

Kurt Zeiler, P.G.*
Principal Hydrogeologist
Brown and Caldwell |1527 Cole Blvd., Suite 300 | Lakewood, CO"

@vivekbedekar as a placeholder so we don't forget this, creating an issue for Daniel's email that stated:

"in a multi-component simulation, is the LKT initial concentration specified for each component? or only for the first? The I/O documentation shows multiple arrays for SFT by component, but not LKT as far as I can tell."

From: Feinstein, Daniel [mailto:[email protected]]
Sent: Sunday, May 6, 2018 2:34 PM
To: Vivek Bedekar [email protected]
Subject: ahh- the problem I think is with MT3D-USGS documentation

Hi Vivek,

Referring to my last email:

The I/O manual for UZT it seems to say that you need to make INCUZET=-1 and INCGWET=-1 even if no ET is active in UZF package. But that is not right. The "-1" values act to make UZT reuse the infiltration concentration from the last stress period - which is disastrous for my run.

When I take out the previous INCUZET and INCGWET lines:
-1
-1
then the results make sense.

Is this the correct interpretation of the manual? Maybe you have already addressed this problem in a
release?

I am sending you a new UZT package - please check the end to make sure it now
makes sense.

thanks much,
d

I hate these messages. Would be good to try and eliminate them. Easy fix is to change G10.4 to G11 or 12 .4. I don't want to do it in case it is building a table and it would affect the alignment.

.\src_temp\rct1.f(493): remark #8291: Recommended relationship between field wid
th 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+
7'.
1 /1X,'STOICHIOMETRIC RATIO (FEDEA) = ',G10.4)
---------------------------------------------------------------^

Checked to make sure the code compiles with more strict compiler settings. Code compiled with the following flags set (see attached screen shot):

1. Array bounds
2. Check uninitialize
3. Stackframes
4. Traceback

USGS does not copyright its code, so the copyright statement should be removed. Also the following needs to be updated. And what is the deal with the PERCEL bug? Does this still exist? What is sub-time stepping?

C=======================================================================
C Version history: xx-xx-2016 (1.00)
C
C  =====================================================================
C
C
C--There appears to be a bug when PERCEL (ADV input file) is set 
C--to something that causes sub-time stepping.  Mass is not      
C--conserved in this scenario and may need to be addressed before
C--code is distributed.         

SFT not properly accounting for multi-species. For example, compare the list files (*.m3d) from the two attached models, 1 and 2 species, respectively.

The following call in gnt1.f (line 39):

CALL xmdsolv(AMATSF, RHSSF, CNEWSF, hclosexmd, rrctol, IASF, 
&            JASF,NJASF,numactive, north, iter, iacl, ierr)

sends in numactive equal to the number of SFT cells. However, CNEWSF gets picked-up as x and numactive gets picked up as n. In xmdsolv.f, on line 1961, x is dimensioned as x(n), but this is a problem because CNEWSF is dimensioned as CNEWSF(1:650,1:2) when it gets sent in, and so the second species is not properly accounted for.

CrnkNic.zip

Do we want to have it write a file instead? Suggest we add something like the following to the name file

DATA              81          .\gwt\gwt.lak.out

I suspect the 64-bit exe does not have all the libraries bundled with it, and thus will not run on PC's that do not have Intel. Would be worthwhile to check this.

Hey guys. Is there any reason not to add optional use of free format lists in the ssm? I see there is support for free format arrays and free format lists would make the ssm file less shitty and easier to parameterize. I can take a shot at it but just want to make sure there is some obvious issue that I'm not thinking of???

Eric Chiang left the following comment in marked-up version of the input instruction that were attached to Issue #16:

In the RCT Package section of the input instructions, item 2c, Eric writes:
"How will SRCONC be used if ISOTHM != 4, 5, 6, or -6?"
For context, see the attachment in Issue #16 item 2c

This line needs to be changed for each release:

CHARACTER,PARAMETER :: VID*18='[Ver 1.00.00]'

The code still reports:

 MT3D-USGS - Modular 3D Multi-Species Transport Model [Ver 1.00.00] 
 and based on MT3DMS. MT3D-USGS developed in cooperation by 
 S.S. Papadopulos & Associates and the U.S. Geological Survey

With commit 0c3902f in MF-NWT, we're now using the flow out of a reach to calculate the cross-sectional area. Prior to this fix, we were using the inflow to the reach to calculate cross-sectional area. This change will likely change the solutions is some of the auto-tests that employ SFT. As this issue states, negative values for FLOW in item 4b of the SFR input (see the Online guide to MODFLOW) were never tested. However, J White was running into issues with negative values for FLOW specified in one of his models.

Consider incorporating Paul Hsieh's sorption example (emailed to Eric on Dec. 3rd) into the regression suite of tests for MT3D-USGS

I think these should be the other way around

KPERMF6=>KPER, KSTPMF6=>KSTP

added MT3D001S.ucn to UZT_NonEq .cmp directory and the nosetests fails. With this file removed, nosetests passes.

example file for testing this functionality added in commit ca030fb

While working on the Trout Lake proof-of-concept (POC) model, rather large mass balance discrepancies were noted in the MT3D-USGS listing file. Originally, an issue was reported in the Trout_Lake repo. However, after rerunning the POC model with Vivek's f4c72e9 commit, percent discrepancies are now small.

hey guys - I just ran into a situation where this number was reported for a sw concen: 0.5771964-104 (missing the e since the exponent is 3 digits and negative). I realize this is allowed in Fortran (barf!) but modern languages (haha) do not like this. Also, for single precision, this would be an underflow. Would it be ok to cast down really small numbers (< e-38) to zero prior to writing to the sft ascii output file? I can probably make this change and submit a pull request but was just checking if this was acceptable.

With commit 377efcf, an external reference to a library called IFLPORT made its way into the code. This prevents gfortran from compiling since it doesn't have access to this dependency. Commit 377efcf offered code originally written by a 3rd party for detecting whether the FTL file is written in double precision (instead of the standard single precision). This is good functionality to keep, just needs to be tweaked to avoid the use of a dependency.

The following is one of the batch files:

@ECHO OFF
REM Generate FTL file using MODFLOW-NWT ver 1.1.0 (or greater)
..\bin\MF_NWT_1.1.1_64.exe .\2ED5EAs\2ED5EA_mf.nam

REM Run MT3D-USGS only after FTL is generated
..\bin\MT3D-USGS_64.exe .\2ED5EAs\2ED5EA_mt.nam
ECHO.
ECHO Run complete. Please press enter when you want to continue.
PAUSE>NUL

Would it be better to change it so that first it changes into that subdirectory and then runs the model?

@ECHO OFF

REM Change directory
cd .\2ED5EAs

REM Generate FTL file using MODFLOW-NWT ver 1.1.0 (or greater)
..\..\bin\MF_NWT_1.1.1_64.exe 2ED5EA_mf.nam

REM Run MT3D-USGS only after FTL is generated
..\..\bin\MT3D-USGS_64.exe 2ED5EA_mt.nam

ECHO.
ECHO Run complete. Please press enter when you want to continue.
PAUSE>NUL

Vivek started getting an error with a model that works with MT3DMS 5.3. Turns out the model runs when TTSMULT=1.0 in the BTN input file; however, when equal to 1.2 in this particular model, the distributed MT3D-USGS_64.exe bombs.

Vivek and I noted that when run from Visual Studio in "Release" mode, with TTSMULT=1.2, the model will bomb when the following compile flag is set equal to "Underflow gives 0.0; Abort on other IEEE exceptions (/fpe:0)" (more to follow image):

The line of code that is bombing is shown below, DTRANS is getting really large by the time the model gets to transport step ~447 in the code below:

A model for debugging this problem was added with commit a226eee (11/8/16)

Daniel Feinstein has requested that in the mass budgets for SFT, the mass contributed to streams from rejected infiltration be separated from groundwater discharge (springs) that makes its way to streams via the IRUNBND array (note: this is not referring to groundwater discharge through the streambed, rather spring flow directed to streams). Seems like a reasonable request. A model from Daniel to test this out with is attached.

NZ_exmp_mod_split_UZT-SFT_runoff_budg_itms.zip

If the ftl file is formatted, MT3D doesn't accurately identify when the UZF package is active. I have a suggested fix.

Here are lines 127-132 of fmi1.f with the fix included. The problem is that when the ftl file is formatted, the read statement has the text followed by blanks but in the tests that follow, it is assumed that the text is preceded by blanks. The fix is to us ADJUSTR on the text.
IF(IFTLFMT.EQ.0) THEN
READ(INFTL) TEXT1
ELSEIF(IFTLFMT.EQ.1) THEN
READ(INFTL,*) TEXT1
TEXT1 = ADJUSTR(TEXT1)
ENDIF

While this corrects this particular problem, the model still won't run because of another bug.
mt3d_bug.zip

Got a question from Rumbaugh about this. Not sure why there is need to have these IREACT values, and they are not described well in the user guide.

      IF(IREACT.EQ.90.OR.IREACT.EQ.91) THEN  
        IREACTION=1                          
        IREACT=IREACT-90                     
      ENDIF                                  

It looks like all that these values do is set IREACTION to 1, which is an input variable anyway.

Any thoughts on this? Do we want to eliminate these special IREACT values for the next release? If so, I will let Jim know.