Feat (mfusg): FloPy MODFLOW-USG support (with .bas file example) #1476

* Feat(mfusg): add mfusg transport version 2.4.0 support

fixes and updates included with feature:
New Features
* MfUsgBct : Block Centered Transport (BCT)
* MfUsgPcb : Prescribed Concentration Boundary (PCB)
* MfUsgDdf: Density Driven Flow (DDF)
* MfUsgDpf: Dual Porosity Flow (DPF)
* MfUsgDpt: Dual Porosity Transport (DPT)
* MfUsgMdt: Matrix Diffusion Transport (MDT)
* MfUsgRch: Recharge (RCH)
* MfUsgLak: Lake (LAK)

Modified Packages
* MfUsgBcf : Block Centered Flow (BCF)
* MfUsgCln :  Connected Linear Network (CLN)
* MfUsgSms :  Sparse Matrix Solver (SMS)
* MfUsgWel :  Well (WEL)
* ModflowFhb : Flow and Head Boundary (FHB)

See Panday, S., 2024; USG-Transport Version 2.4.0: Transport and Other Enhancements to MODFLOW-USG, GSI Environmental, July 2024 http://www.gsi-net.com/en/software/free-software/USG-Transport.html

flopy/mfusg/__init__.py

@@ -3,10 +3,19 @@
 from .cln_dtypes import MfUsgClnDtypes
 from .mfusg import MfUsg
 from .mfusgbcf import MfUsgBcf
+from .mfusgbct import MfUsgBct
 from .mfusgcln import MfUsgCln
+from .mfusgddf import MfUsgDdf
 from .mfusgdisu import MfUsgDisU
+from .mfusgdpf import MfUsgDpf
+from .mfusgdpt import MfUsgDpt
 from .mfusggnc import MfUsgGnc
+from .mfusglak import MfUsgLak
 from .mfusglpf import MfUsgLpf
+from .mfusgmdt import MfUsgMdt
+from .mfusgoc import MfUsgOc
+from .mfusgpcb import MfUsgPcb
+from .mfusgrch import MfUsgRch
 from .mfusgsms import MfUsgSms
 from .mfusgwel import MfUsgWel
 
@@ -20,4 +29,13 @@
     "MfUsgClnDtypes",
     "MfUsgSms",
     "MfUsgGnc",
+    "MfUsgBct",
+    "MfUsgPcb",
+    "MfUsgDdf",
+    "MfUsgMdt",
+    "MfUsgDpf",
+    "MfUsgDpt",
+    "MfUsgRch",
+    "MfUsgOc",
+    "MfUsgLak",
 ]

flopy/mfusg/mfusg.py

@@ -91,6 +91,11 @@ def __init__(
             verbose=self.verbose,
             **kwargs,
         )
+
+        self.mcomp = 0  # number of chemical components
+        self.iheat = 0  # flag for heat transport
+        self.idpf = 0  # flag for dual porosity flow
+
         # Create a dictionary to map package with package object.
         # This is used for loading models.
         self.mfnam_packages = {
@@ -104,15 +109,12 @@ def __init__(
             "fhb": flopy.modflow.ModflowFhb,
             "drn": flopy.modflow.ModflowDrn,
             "drt": flopy.modflow.ModflowDrt,
-            "rch": flopy.modflow.ModflowRch,
             "evt": flopy.modflow.ModflowEvt,
             "ghb": flopy.modflow.ModflowGhb,
             "riv": flopy.modflow.ModflowRiv,
             "str": flopy.modflow.ModflowStr,
             "sfr": flopy.modflow.ModflowSfr2,
-            "lak": flopy.modflow.ModflowLak,
             "gage": flopy.modflow.ModflowGage,
-            "oc": flopy.modflow.ModflowOc,
             "sub": flopy.modflow.ModflowSub,
             "swt": flopy.modflow.ModflowSwt,
             "disu": flopy.mfusg.MfUsgDisU,
@@ -122,6 +124,15 @@ def __init__(
             "lpf": flopy.mfusg.MfUsgLpf,
             "cln": flopy.mfusg.MfUsgCln,
             "gnc": flopy.mfusg.MfUsgGnc,
+            "bct": flopy.mfusg.MfUsgBct,
+            "pcb": flopy.mfusg.MfUsgPcb,
+            "ddf": flopy.mfusg.MfUsgDdf,
+            "mdt": flopy.mfusg.MfUsgMdt,
+            "dpf": flopy.mfusg.MfUsgDpf,
+            "dpt": flopy.mfusg.MfUsgDpt,
+            "rch": flopy.mfusg.MfUsgRch,
+            "oc": flopy.mfusg.MfUsgOc,
+            "lak": flopy.mfusg.MfUsgLak,
         }
 
     def __repr__(self):
@@ -529,7 +540,7 @@ def fmt_string(array):
         if vtype in ("i", "b"):
             fmts.append("%10d")
         elif vtype == "f":
-            fmts.append("%14.6g")
+            fmts.append("%10.2g")
         elif vtype == "o":
             fmts.append("%10s")
         elif vtype == "s":

flopy/mfusg/mfusgbcf.py

@@ -138,6 +138,12 @@ def __init__(
         ihdwet=0,
         ikvflag=0,
         ikcflag=0,
+        tabrich=False,
+        nuzones=0,
+        nutabrows=0,
+        bubblept=False,
+        fullydry=False,
+        altsto=False,
         tran=1.0,
         hy=1.0,
         vcont=1.0,
@@ -194,6 +200,13 @@ def __init__(
 
         self.ikvflag = ikvflag
         self.ikcflag = ikcflag
+        self.tabrich = tabrich
+        self.nuzones = nuzones
+        self.nutabrows = nutabrows
+        self.bubblept = bubblept
+        self.fullydry = fullydry
+        self.altsto = altsto
+
         self.kv = kv
         self.anglex = anglex
         self.ksat = ksat
@@ -255,9 +268,19 @@ def write_file(self, f=None):
         f_obj.write(
             f" {self.ipakcb:9d} {self.hdry:9.3G} {self.iwdflg:9d}"
             f" {self.wetfct:9.3G} {self.iwetit:9d} {self.ihdwet:9d}"
-            f" {self.ikvflag:9d} {self.ikcflag:9d}\n"
+            f" {self.ikvflag:9d} {self.ikcflag:9d}"
         )
 
+        if self.tabrich:
+            f_obj.write(f" TABRICH {self.nuzones:9d} {self.nutabrows:9d}")
+        if self.bubblept:
+            f_obj.write(" BUBBLEPT")
+        if self.fullydry:
+            f_obj.write(" FULLYDRY")
+        if self.altsto:
+            f_obj.write(" ALTSTO")
+        f_obj.write("\n")
+
         # LAYCON array
         for layer in range(nlay):
             if self.intercellt[layer] > 0:
@@ -386,6 +409,37 @@ def load(cls, f, model, ext_unit_dict=None):
         ikvflag = type_from_iterable(text_list, index=6, _type=int, default_val=0)
         ikcflag = type_from_iterable(text_list, index=7, _type=int, default_val=0)
 
+        # options
+        if "TABRICH" in text_list:
+            idx = text_list.index("TABRICH")
+            tabrich = True
+            nuzones = float(text_list[idx + 1])
+            nutabrows = float(text_list[idx + 2])
+        else:
+            tabrich = False
+            nuzones = None
+            nutabrows = None
+
+        if "BUBBLEPT" in text_list:
+            bubblept = True
+        else:
+            bubblept = False
+
+        if "FULLYDRY" in text_list:
+            fullydry = True
+        else:
+            fullydry = False
+
+        if "ALTSTO" in text_list:
+            altsto = True
+        else:
+            altsto = False
+
+        # item 1c and d  -- Not implemented
+        if tabrich:
+            iuzontab = []
+            retcrvs = []
+
         # LAYCON array
         laycon, intercellt = cls._load_laycon(f_obj, model)
 
@@ -446,6 +500,12 @@ def load(cls, f, model, ext_unit_dict=None):
             ihdwet=ihdwet,
             ikvflag=ikvflag,
             ikcflag=ikcflag,
+            tabrich=tabrich,
+            nuzones=nuzones,
+            nutabrows=nutabrows,
+            bubblept=bubblept,
+            fullydry=fullydry,
+            altsto=altsto,
             tran=tran,
             hy=hy,
             vcont=vcont,