From 2babaa8be4f58da3de0b79dcc48d744fc89c848d Mon Sep 17 00:00:00 2001
From: Eric Morway <emorway@usgs.gov>
Date: Thu, 2 Nov 2023 11:26:56 -0700
Subject: [PATCH] Add new function that accounts for auxiliary multiplier > 0
 in response to
 https://github.com/MODFLOW-USGS/modflow6/pull/1413#discussion_r1380567761

---
 src/Model/GroundWaterTransport/gwt1cnc1.f90 | 49 +++++++++++++++------
 1 file changed, 36 insertions(+), 13 deletions(-)

diff --git a/src/Model/GroundWaterTransport/gwt1cnc1.f90 b/src/Model/GroundWaterTransport/gwt1cnc1.f90
index c295fdf2723..526b32e67f3 100644
--- a/src/Model/GroundWaterTransport/gwt1cnc1.f90
+++ b/src/Model/GroundWaterTransport/gwt1cnc1.f90
@@ -39,6 +39,7 @@ module GwtCncModule
     procedure :: allocate_arrays => cnc_allocate_arrays
     procedure :: define_listlabel
     procedure :: bound_value => cnc_bound_value
+    procedure :: conc_mult
     ! -- methods for observations
     procedure, public :: bnd_obs_supported => cnc_obs_supported
     procedure, public :: bnd_df_obs => cnc_df_obs
@@ -92,7 +93,7 @@ subroutine cnc_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
     ! -- Store the appropriate label based on the dependent variable
     cncobj%depvartype = depvartype
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_create
 
@@ -128,7 +129,7 @@ subroutine cnc_allocate_arrays(this, nodelist, auxvar)
                      'TSPVAR', this%input_mempath)
     !
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_allocate_arrays
 
@@ -179,7 +180,7 @@ subroutine cnc_rp(this)
       call this%write_list()
     end if
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_rp
 
@@ -216,7 +217,7 @@ subroutine cnc_ad(this)
     !    "current" value.
     call this%obs%obs_ad()
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_ad
 
@@ -251,14 +252,14 @@ subroutine cnc_ck(this)
       call store_error_filename(this%input_fname)
     end if
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_ck
 
   !> @brief Override bnd_fc and do nothing
   !!
   !! For constant concentration/temperature boundary type, the call to bnd_fc
-  !! needs to be overwritten to prevent logic found therein from being applied
+  !! needs to be overwritten to prevent logic found in bnd from being executed
   !<
   subroutine cnc_fc(this, rhs, ia, idxglo, matrix_sln)
     ! -- dummy
@@ -269,7 +270,7 @@ subroutine cnc_fc(this, rhs, ia, idxglo, matrix_sln)
     class(MatrixBaseType), pointer :: matrix_sln
     ! -- local
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_fc
 
@@ -338,7 +339,7 @@ subroutine cnc_cq(this, x, flowja, iadv)
       !
     end if
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_cq
 
@@ -385,7 +386,7 @@ subroutine cnc_da(this)
     call mem_deallocate(this%ratecncout)
     call mem_deallocate(this%tspvar, 'TSPVAR', this%memoryPath)
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_da
 
@@ -416,7 +417,7 @@ subroutine define_listlabel(this)
       write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
     end if
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine define_listlabel
 
@@ -452,7 +453,7 @@ subroutine cnc_df_obs(this)
     call this%obs%StoreObsType(this%filtyp, .true., indx)
     this%obs%obsData(indx)%ProcessIdPtr => DefaultObsIdProcessor
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_df_obs
 
@@ -483,9 +484,31 @@ subroutine cnc_rp_ts(this)
       end if
     end do
     !
-    ! -- return
+    ! -- Return
     return
   end subroutine cnc_rp_ts
+  
+  !> @brief Apply auxiliary multiplier to specified concentration if 
+  !< appropriate
+  function conc_mult(this, row) result(conc)
+    ! -- modules
+    use ConstantsModule, only: DZERO
+    ! -- dummy variables
+    class(GwtCncType), intent(inout) :: this !< BndExtType object
+    integer(I4B), intent(in) :: row
+    ! -- result
+    real(DP) :: conc
+    !
+    if (this%iauxmultcol > 0) then
+      conc = this%tspvar(row) * this%auxvar(this%iauxmultcol, row)
+    else
+      conc = this%tspvar(row)
+    end if
+    !
+    ! -- Return
+    return
+  end function conc_mult
+
 
   !> @ brief Return a bound value
   !!
@@ -503,7 +526,7 @@ function cnc_bound_value(this, col, row) result(bndval)
     !
     select case (col)
     case (1)
-      bndval = this%tspvar(row)
+      bndval = this%conc_mult(row)
     case default
       write (errmsg, '(3a)') 'Programming error. ', &
                & adjustl(trim(this%filtyp)), ' bound value requested column '&