MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
tspssmmodule Module Reference

This module contains the TspSsm Module. More...

Data Types

type  tspssmtype
 Derived type for the SSM Package. More...
 

Functions/Subroutines

subroutine, public ssm_cr (ssmobj, name_model, inunit, iout, fmi, eqnsclfac, depvartype)
 @ brief Create a new SSM package More...
 
subroutine ssm_df (this)
 @ brief Define SSM Package More...
 
subroutine ssm_ar (this, dis, ibound, cnew)
 @ brief Allocate and read SSM Package More...
 
subroutine ssm_rp (this)
 @ brief Read and prepare this SSM Package More...
 
subroutine ssm_ad (this)
 @ brief Advance the SSM Package More...
 
subroutine ssm_term (this, ipackage, ientry, rrate, rhsval, hcofval, cssm, qssm)
 @ brief Calculate the SSM mass flow rate and hcof and rhs values More...
 
subroutine get_ssm_conc (this, ipackage, ientry, nbound_flow, conc, lauxmixed)
 @ brief Provide bound concentration (or temperature) and mixed flag More...
 
subroutine ssm_fc (this, matrix_sln, idxglo, rhs)
 @ brief Fill coefficients More...
 
subroutine ssm_cq (this, flowja)
 @ brief Calculate flow More...
 
subroutine ssm_bd (this, isuppress_output, model_budget)
 @ brief Calculate the global SSM budget terms More...
 
subroutine ssm_ot_flow (this, icbcfl, ibudfl, icbcun)
 @ brief Output flows More...
 
subroutine ssm_da (this)
 @ brief Deallocate More...
 
subroutine allocate_scalars (this)
 @ brief Allocate scalars More...
 
subroutine allocate_arrays (this)
 @ brief Allocate arrays More...
 
subroutine read_options (this)
 @ brief Read package options More...
 
subroutine read_data (this)
 @ brief Read package data More...
 
subroutine read_sources_aux (this)
 @ brief Read SOURCES block More...
 
subroutine read_sources_fileinput (this)
 @ brief Read FILEINPUT block More...
 
subroutine set_iauxpak (this, ip, packname)
 @ brief Set iauxpak array value for package ip More...
 
subroutine set_ssmivec (this, ip, packname)
 @ brief Set ssmivec array value for package ip More...
 
subroutine pak_setup_outputtab (this)
 @ brief Setup the output table More...
 

Variables

character(len=lenftype) ftype = 'SSM'
 
character(len=lenpackagename) text = ' SOURCE-SINK MIX'
 

Detailed Description

This module contains the code for handling sources and sinks associated with groundwater flow model stress packages.

todo: need observations for SSM terms

Function/Subroutine Documentation

◆ allocate_arrays()

subroutine tspssmmodule::allocate_arrays ( class(tspssmtype this)

Allocate array variables for this derived type

Parameters
thisTspSsmType object

Definition at line 712 of file tsp-ssm.f90.

713  ! -- modules
715  ! -- dummy
716  class(TspSsmType) :: this !< TspSsmType object
717  ! -- local
718  integer(I4B) :: nflowpack
719  integer(I4B) :: i
720  !
721  ! -- Allocate
722  nflowpack = this%fmi%nflowpack
723  call mem_allocate(this%iauxpak, nflowpack, 'IAUXPAK', this%memoryPath)
724  call mem_allocate(this%isrctype, nflowpack, 'ISRCTYPE', this%memoryPath)
725  !
726  ! -- Initialize
727  do i = 1, nflowpack
728  this%iauxpak(i) = 0
729  this%isrctype(i) = 0
730  end do
731  !
732  ! -- Allocate the ssmivec array
733  allocate (this%ssmivec(nflowpack))

◆ allocate_scalars()

subroutine tspssmmodule::allocate_scalars ( class(tspssmtype this)

Allocate scalar variables for this derived type

Parameters
thisTspSsmType object

Definition at line 692 of file tsp-ssm.f90.

693  ! -- modules
695  ! -- dummy
696  class(TspSsmType) :: this !< TspSsmType object
697  !
698  ! -- allocate scalars in NumericalPackageType
699  call this%NumericalPackageType%allocate_scalars()
700  !
701  ! -- Allocate
702  call mem_allocate(this%nbound, 'NBOUND', this%memoryPath)
703  !
704  ! -- Initialize
705  this%nbound = 0

◆ get_ssm_conc()

subroutine tspssmmodule::get_ssm_conc ( class(tspssmtype this,
integer(i4b), intent(in)  ipackage,
integer(i4b), intent(in)  ientry,
integer(i4b), intent(in)  nbound_flow,
real(dp), intent(out)  conc,
logical(lgp), intent(out)  lauxmixed 
)

SSM concentrations and temperatures can be provided in auxiliary variables or through separate SPC files. If not provided, the default concentration (or temperature) is zero. This single routine provides the SSM bound concentration (or temperature) based on these different approaches. The mixed flag indicates whether or not the boundary as a mixed type.

Parameters
thisTspSsmType
[in]ipackagepackage number
[in]ientrybound number
[in]nbound_flowsize of flow package bound list
[out]concuser-specified concentration/temperature for this bound
[out]lauxmixeduser-specified flag for marking this as a mixed boundary

Definition at line 356 of file tsp-ssm.f90.

358  ! -- dummy
359  class(TspSsmType) :: this !< TspSsmType
360  integer(I4B), intent(in) :: ipackage !< package number
361  integer(I4B), intent(in) :: ientry !< bound number
362  integer(I4B), intent(in) :: nbound_flow !< size of flow package bound list
363  real(DP), intent(out) :: conc !< user-specified concentration/temperature for this bound
364  logical(LGP), intent(out) :: lauxmixed !< user-specified flag for marking this as a mixed boundary
365  ! -- local
366  integer(I4B) :: isrctype
367  integer(I4B) :: iauxpos
368  !
369  conc = dzero
370  lauxmixed = .false.
371  isrctype = this%isrctype(ipackage)
372  !
373  select case (isrctype)
374  case (1, 2)
375  iauxpos = this%iauxpak(ipackage)
376  conc = this%fmi%gwfpackages(ipackage)%auxvar(iauxpos, ientry)
377  if (isrctype == 2) lauxmixed = .true.
378  case (3, 4)
379  conc = this%ssmivec(ipackage)%get_value(ientry, nbound_flow)
380  if (isrctype == 4) lauxmixed = .true.
381  end select

◆ pak_setup_outputtab()

subroutine tspssmmodule::pak_setup_outputtab ( class(tspssmtype), intent(inout)  this)

Setup the output table by creating the column headers.

Definition at line 1085 of file tsp-ssm.f90.

1086  ! -- dummy
1087  class(TspSsmType), intent(inout) :: this
1088  ! -- local
1089  character(len=LINELENGTH) :: title
1090  character(len=LINELENGTH) :: text
1091  integer(I4B) :: ntabcol
1092  !
1093  ! -- allocate and initialize the output table
1094  if (this%iprflow /= 0) then
1095  !
1096  ! -- dimension table
1097  ntabcol = 6
1098  !if (this%inamedbound > 0) then
1099  ! ntabcol = ntabcol + 1
1100  !end if
1101  !
1102  ! -- initialize the output table object
1103  title = 'SSM PACKAGE ('//trim(this%packName)// &
1104  ') FLOW RATES'
1105  call table_cr(this%outputtab, this%packName, title)
1106  call this%outputtab%table_df(1, ntabcol, this%iout, transient=.true.)
1107  text = 'NUMBER'
1108  call this%outputtab%initialize_column(text, 10, alignment=tabcenter)
1109  text = 'CELLID'
1110  call this%outputtab%initialize_column(text, 20, alignment=tableft)
1111  text = 'BOUND Q'
1112  call this%outputtab%initialize_column(text, 15, alignment=tabcenter)
1113  text = 'SSM CONC'
1114  call this%outputtab%initialize_column(text, 15, alignment=tabcenter)
1115  text = 'RATE'
1116  call this%outputtab%initialize_column(text, 15, alignment=tabcenter)
1117  text = 'PACKAGE NAME'
1118  call this%outputtab%initialize_column(text, 16, alignment=tabcenter)
1119  !if (this%inamedbound > 0) then
1120  ! text = 'NAME'
1121  ! call this%outputtab%initialize_column(text, 20, alignment=TABLEFT)
1122  !end if
1123  end if
Here is the call graph for this function:

◆ read_data()

subroutine tspssmmodule::read_data ( class(tspssmtype this)

Read and set the SSM Package data

Parameters
thisTspSsmType object

Definition at line 787 of file tsp-ssm.f90.

788  ! -- dummy
789  class(TspSsmType) :: this !< TspSsmType object
790  !
791  ! -- read and process required SOURCES block
792  call this%read_sources_aux()
793  !
794  ! -- read and process optional FILEINPUT block
795  call this%read_sources_fileinput()

◆ read_options()

subroutine tspssmmodule::read_options ( class(tspssmtype this)

Read and set the SSM Package options

Parameters
thisTspSsmType object

Definition at line 740 of file tsp-ssm.f90.

741  ! -- dummy
742  class(TspSsmType) :: this !< TspSsmType object
743  ! -- local
744  character(len=LINELENGTH) :: keyword
745  integer(I4B) :: ierr
746  logical :: isfound, endOfBlock
747  ! -- formats
748  character(len=*), parameter :: fmtiprflow = &
749  "(4x,'SSM FLOW INFORMATION WILL BE PRINTED TO LISTING FILE &
750  &WHENEVER ICBCFL IS NOT ZERO.')"
751  character(len=*), parameter :: fmtisvflow = &
752  "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
753  &WHENEVER ICBCFL IS NOT ZERO.')"
754  !
755  ! -- get options block
756  call this%parser%GetBlock('OPTIONS', isfound, ierr, blockrequired=.false., &
757  supportopenclose=.true.)
758  !
759  ! -- parse options block if detected
760  if (isfound) then
761  write (this%iout, '(1x,a)') 'PROCESSING SSM OPTIONS'
762  do
763  call this%parser%GetNextLine(endofblock)
764  if (endofblock) exit
765  call this%parser%GetStringCaps(keyword)
766  select case (keyword)
767  case ('PRINT_FLOWS')
768  this%iprflow = 1
769  write (this%iout, fmtiprflow)
770  case ('SAVE_FLOWS')
771  this%ipakcb = -1
772  write (this%iout, fmtisvflow)
773  case default
774  write (errmsg, '(a,a)') 'Unknown SSM option: ', trim(keyword)
775  call store_error(errmsg)
776  call this%parser%StoreErrorUnit()
777  end select
778  end do
779  write (this%iout, '(1x,a)') 'END OF SSM OPTIONS'
780  end if
Here is the call graph for this function:

◆ read_sources_aux()

subroutine tspssmmodule::read_sources_aux ( class(tspssmtype this)

Read SOURCES block and look for auxiliary columns in corresponding flow data.

Parameters
thisTspSsmType object

Definition at line 803 of file tsp-ssm.f90.

804  ! -- dummy
805  class(TspSsmType) :: this !< TspSsmType object
806  ! -- local
807  character(len=LINELENGTH) :: keyword
808  character(len=20) :: srctype
809  integer(I4B) :: ierr
810  integer(I4B) :: ip
811  integer(I4B) :: nflowpack
812  integer(I4B) :: isrctype
813  logical :: isfound, endOfBlock
814  logical :: pakfound
815  logical :: lauxmixed
816  !
817  ! -- initialize
818  isfound = .false.
819  lauxmixed = .false.
820  nflowpack = this%fmi%nflowpack
821  !
822  ! -- get sources block
823  call this%parser%GetBlock('SOURCES', isfound, ierr, &
824  supportopenclose=.true., &
825  blockrequired=.true.)
826  if (isfound) then
827  write (this%iout, '(1x,a)') 'PROCESSING SOURCES'
828  do
829  call this%parser%GetNextLine(endofblock)
830  if (endofblock) exit
831  !
832  ! -- read package name and make sure it can be found
833  call this%parser%GetStringCaps(keyword)
834  pakfound = .false.
835  do ip = 1, nflowpack
836  if (trim(adjustl(this%fmi%gwfpackages(ip)%name)) == keyword) then
837  pakfound = .true.
838  exit
839  end if
840  end do
841  if (.not. pakfound) then
842  write (errmsg, '(a,a)') 'Flow package cannot be found: ', &
843  trim(keyword)
844  call store_error(errmsg)
845  call this%parser%StoreErrorUnit()
846  end if
847  !
848  ! -- Ensure package was not specified more than once in SOURCES block
849  if (this%isrctype(ip) /= 0) then
850  write (errmsg, '(a, a)') &
851  'A package cannot be specified more than once in the SSM SOURCES &
852  &block. The following package was specified more than once: ', &
853  trim(keyword)
854  call store_error(errmsg)
855  call this%parser%StoreErrorUnit()
856  end if
857  !
858  ! -- read the source type
859  call this%parser%GetStringCaps(srctype)
860  select case (srctype)
861  case ('AUX')
862  write (this%iout, '(1x,a)') 'AUX SOURCE DETECTED.'
863  isrctype = 1
864  case ('AUXMIXED')
865  write (this%iout, '(1x,a)') 'AUXMIXED SOURCE DETECTED.'
866  lauxmixed = .true.
867  isrctype = 2
868  case default
869  write (errmsg, '(a, a)') &
870  'SRCTYPE must be AUX or AUXMIXED. Found: ', trim(srctype)
871  call store_error(errmsg)
872  call this%parser%StoreErrorUnit()
873  end select
874  !
875  ! -- Store the source type (1 or 2)
876  this%isrctype(ip) = isrctype
877  !
878  ! -- Find and store the auxiliary column
879  call this%set_iauxpak(ip, trim(keyword))
880 
881  end do
882  write (this%iout, '(1x,a)') 'END PROCESSING SOURCES'
883  else
884  write (errmsg, '(a)') 'Required SOURCES block not found.'
885  call store_error(errmsg)
886  call this%parser%StoreErrorUnit()
887  end if
888  !
889  ! -- terminate if errors
890  if (count_errors() > 0) then
891  call this%parser%StoreErrorUnit()
892  end if
Here is the call graph for this function:

◆ read_sources_fileinput()

subroutine tspssmmodule::read_sources_fileinput ( class(tspssmtype this)

Read optional FILEINPUT block and initialize an SPC input file reader for each entry.

Parameters
thisTspSsmType object

Definition at line 900 of file tsp-ssm.f90.

901  ! -- dummy
902  class(TspSsmType) :: this !< TspSsmType object
903  ! -- local
904  character(len=LINELENGTH) :: keyword
905  character(len=LINELENGTH) :: keyword2
906  character(len=20) :: srctype
907  integer(I4B) :: ierr
908  integer(I4B) :: ip
909  integer(I4B) :: nflowpack
910  integer(I4B) :: isrctype
911  logical :: isfound, endOfBlock
912  logical :: pakfound
913  logical :: lauxmixed
914  !
915  ! -- initialize
916  isfound = .false.
917  lauxmixed = .false.
918  nflowpack = this%fmi%nflowpack
919  !
920  ! -- get sources_file block
921  call this%parser%GetBlock('FILEINPUT', isfound, ierr, &
922  supportopenclose=.true., &
923  blockrequired=.false.)
924  if (isfound) then
925  write (this%iout, '(1x,a)') 'PROCESSING FILEINPUT'
926  do
927  call this%parser%GetNextLine(endofblock)
928  if (endofblock) exit
929  !
930  ! -- read package name and make sure it can be found
931  call this%parser%GetStringCaps(keyword)
932  pakfound = .false.
933  do ip = 1, nflowpack
934  if (trim(adjustl(this%fmi%gwfpackages(ip)%name)) == keyword) then
935  pakfound = .true.
936  exit
937  end if
938  end do
939  if (.not. pakfound) then
940  write (errmsg, '(a,a)') 'Flow package cannot be found: ', &
941  trim(keyword)
942  call store_error(errmsg)
943  call this%parser%StoreErrorUnit()
944  end if
945  !
946  ! -- Ensure package was not specified more than once in SOURCES block
947  if (this%isrctype(ip) /= 0) then
948  write (errmsg, '(a, a)') &
949  'A package cannot be specified more than once in the SSM SOURCES &
950  &and SOURCES_FILES blocks. The following package was specified &
951  &more than once: ', &
952  trim(keyword)
953  call store_error(errmsg)
954  call this%parser%StoreErrorUnit()
955  end if
956  !
957  ! -- read the source type
958  call this%parser%GetStringCaps(srctype)
959  select case (srctype)
960  case ('SPC6')
961  write (this%iout, '(1x,a)') 'SPC6 SOURCE DETECTED.'
962  isrctype = 3
963  !
964  ! verify filein is next
965  call this%parser%GetStringCaps(keyword2)
966  if (trim(adjustl(keyword2)) /= 'FILEIN') then
967  errmsg = 'SPC6 keyword must be followed by "FILEIN" '// &
968  'then by filename and optionally by <MIXED>.'
969  call store_error(errmsg)
970  call this%parser%StoreErrorUnit()
971  end if
972  !
973  ! -- Use set_ssmivec to read file name and set up
974  ! ssmi file object
975  call this%set_ssmivec(ip, trim(keyword))
976  !
977  ! -- check for optional MIXED keyword and set isrctype to 4 if found
978  call this%parser%GetStringCaps(keyword2)
979  if (trim(keyword2) == 'MIXED') then
980  isrctype = 4
981  write (this%iout, '(1x,a,a)') 'ASSIGNED MIXED SSM TYPE TO PACKAGE ', &
982  trim(keyword)
983  end if
984  case default
985  write (errmsg, '(a,a)') &
986  'SRCTYPE must be SPC6. Found: ', trim(srctype)
987  call store_error(errmsg)
988  call this%parser%StoreErrorUnit()
989  end select
990  !
991  ! -- Store the source type (3 or 4)
992  this%isrctype(ip) = isrctype
993  !
994  end do
995  write (this%iout, '(1x,a)') 'END PROCESSING FILEINPUT'
996  else
997  write (this%iout, '(1x,a)') &
998  'OPTIONAL FILEINPUT BLOCK NOT FOUND. CONTINUING.'
999  end if
1000  !
1001  ! -- terminate if errors
1002  if (count_errors() > 0) then
1003  call this%parser%StoreErrorUnit()
1004  end if
Here is the call graph for this function:

◆ set_iauxpak()

subroutine tspssmmodule::set_iauxpak ( class(tspssmtype), intent(inout)  this,
integer(i4b), intent(in)  ip,
character(len=*), intent(in)  packname 
)

The next call to parser will return the auxiliary name for package ip in the SSM SOURCES block. The routine searches through the auxiliary names in package ip and sets iauxpak to the column number corresponding to the correct auxiliary column.

Parameters
[in,out]thisTspSsmType
[in]ippackage number
[in]packnamename of package

Definition at line 1015 of file tsp-ssm.f90.

1016  ! -- dummy
1017  class(TspSsmType), intent(inout) :: this !< TspSsmType
1018  integer(I4B), intent(in) :: ip !< package number
1019  character(len=*), intent(in) :: packname !< name of package
1020  ! -- local
1021  character(len=LENAUXNAME) :: auxname
1022  logical :: auxfound
1023  integer(I4B) :: iaux
1024  !
1025  ! -- read name of auxiliary column
1026  call this%parser%GetStringCaps(auxname)
1027  auxfound = .false.
1028  do iaux = 1, this%fmi%gwfpackages(ip)%naux
1029  if (trim(this%fmi%gwfpackages(ip)%auxname(iaux)) == &
1030  trim(auxname)) then
1031  auxfound = .true.
1032  exit
1033  end if
1034  end do
1035  if (.not. auxfound) then
1036  write (errmsg, '(a, a)') &
1037  'Auxiliary name cannot be found: ', trim(auxname)
1038  call store_error(errmsg)
1039  call this%parser%StoreErrorUnit()
1040  end if
1041  !
1042  ! -- set iauxpak and write message
1043  this%iauxpak(ip) = iaux
1044  write (this%iout, '(4x, a, i0, a, a)') 'USING AUX COLUMN ', &
1045  iaux, ' IN PACKAGE ', trim(packname)
Here is the call graph for this function:

◆ set_ssmivec()

subroutine tspssmmodule::set_ssmivec ( class(tspssmtype), intent(inout)  this,
integer(i4b), intent(in)  ip,
character(len=*), intent(in)  packname 
)

The next call to parser will return the input file name for package ip in the SSM SOURCES block. The routine then initializes the SPC input file.

Parameters
[in,out]thisTspSsmType
[in]ippackage number
[in]packnamename of package

Definition at line 1054 of file tsp-ssm.f90.

1055  ! -- module
1056  use inputoutputmodule, only: openfile, getunit
1057  ! -- dummy
1058  class(TspSsmType), intent(inout) :: this !< TspSsmType
1059  integer(I4B), intent(in) :: ip !< package number
1060  character(len=*), intent(in) :: packname !< name of package
1061  ! -- local
1062  character(len=LINELENGTH) :: filename
1063  type(TspSpcType), pointer :: ssmiptr
1064  integer(I4B) :: inunit
1065  !
1066  ! -- read file name
1067  call this%parser%GetString(filename)
1068  inunit = getunit()
1069  call openfile(inunit, this%iout, filename, 'SPC', filstat_opt='OLD')
1070  !
1071  ! -- Create the SPC file object
1072  ssmiptr => this%ssmivec(ip)
1073  call ssmiptr%initialize(this%dis, ip, inunit, this%iout, this%name_model, &
1074  trim(packname), this%depvartype)
1075  !
1076  write (this%iout, '(4x, a, a, a, a, a)') 'USING SPC INPUT FILE ', &
1077  trim(filename), ' TO SET ', trim(this%depvartype), &
1078  'S FOR PACKAGE ', trim(packname)
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
Here is the call graph for this function:

◆ ssm_ad()

subroutine tspssmmodule::ssm_ad ( class(tspssmtype this)

This routine is called from gwt_ad(). It is called at the beginning of each time step. The total number of flow boundaries is counted and stored in thisnbound. Also, if any SPC input files are used to provide source and sink concentrations (or temperatures) and time series are referenced in those files, then ssm concenrations must be interpolated for the time step.

Parameters
thisTspSsmType object

Definition at line 214 of file tsp-ssm.f90.

215  ! -- dummy
216  class(TspSsmType) :: this !< TspSsmType object
217  ! -- local
218  integer(I4B) :: ip
219  type(TspSpcType), pointer :: ssmiptr
220  integer(I4B) :: i
221  integer(I4B) :: node
222  !
223  ! -- Calculate total number of existing flow boundaries. It is possible
224  ! that a node may equal zero. In this case, the bound should be
225  ! skipped and not written to ssm output.
226  this%nbound = 0
227  do ip = 1, this%fmi%nflowpack
228  if (this%fmi%iatp(ip) /= 0) cycle
229  do i = 1, this%fmi%gwfpackages(ip)%nbound
230  node = this%fmi%gwfpackages(ip)%nodelist(i)
231  if (node > 0) then
232  this%nbound = this%nbound + 1
233  end if
234  end do
235  end do
236  !
237  ! -- Call the ad method on any ssm input files so that values for
238  ! time-series are interpolated
239  do ip = 1, this%fmi%nflowpack
240  if (this%fmi%iatp(ip) /= 0) cycle
241  if (this%isrctype(ip) == 3 .or. this%isrctype(ip) == 4) then
242  ssmiptr => this%ssmivec(ip)
243  call ssmiptr%spc_ad(this%fmi%gwfpackages(ip)%nbound, &
244  this%fmi%gwfpackages(ip)%budtxt)
245  end if
246  end do

◆ ssm_ar()

subroutine tspssmmodule::ssm_ar ( class(tspssmtype this,
class(disbasetype), intent(in), pointer  dis,
integer(i4b), dimension(:), pointer, contiguous  ibound,
real(dp), dimension(:), pointer, contiguous  cnew 
)

This routine is called from gwt_ar(). It allocates arrays, reads options and data, and sets up the output table.

Parameters
thisTspSsmType object
[in]disdiscretization package
iboundGWT model ibound
cnewGWT model dependent variable

Definition at line 134 of file tsp-ssm.f90.

135  ! -- modules
137  ! -- dummy
138  class(TspSsmType) :: this !< TspSsmType object
139  class(DisBaseType), pointer, intent(in) :: dis !< discretization package
140  integer(I4B), dimension(:), pointer, contiguous :: ibound !< GWT model ibound
141  real(DP), dimension(:), pointer, contiguous :: cnew !< GWT model dependent variable
142  ! -- formats
143  character(len=*), parameter :: fmtssm = &
144  "(1x,/1x,'SSM -- SOURCE-SINK MIXING PACKAGE, VERSION 1, 8/25/2017', &
145  &' INPUT READ FROM UNIT ', i0, //)"
146  !
147  ! -- print a message identifying the storage package.
148  write (this%iout, fmtssm) this%inunit
149  !
150  ! -- store pointers to arguments that were passed in
151  this%dis => dis
152  this%ibound => ibound
153  this%cnew => cnew
154  !
155  ! -- Check to make sure that there are flow packages
156  if (this%fmi%nflowpack == 0) then
157  write (errmsg, '(a)') 'SSM package does not detect any boundary flows &
158  &that require SSM terms. Activate GWF-GWT &
159  &exchange or activate FMI package and provide a &
160  &budget file that contains boundary flows. If no &
161  &boundary flows are present in corresponding GWF &
162  &model then this SSM package should be removed.'
163  call store_error(errmsg)
164  call this%parser%StoreErrorUnit()
165  end if
166  !
167  ! -- Allocate arrays
168  call this%allocate_arrays()
169  !
170  ! -- Read ssm options
171  call this%read_options()
172  !
173  ! -- read the data block
174  call this%read_data()
175  !
176  ! -- setup the output table
177  call this%pak_setup_outputtab()
Here is the call graph for this function:

◆ ssm_bd()

subroutine tspssmmodule::ssm_bd ( class(tspssmtype this,
integer(i4b), intent(in)  isuppress_output,
type(budgettype), intent(inout)  model_budget 
)

Calculate the global SSM budget terms using separate in and out entries for each flow package.

Parameters
thisTspSsmType object
[in]isuppress_outputflag to suppress output
[in,out]model_budgetbudget object for the GWT model

Definition at line 466 of file tsp-ssm.f90.

467  ! -- modules
468  use tdismodule, only: delt
469  use budgetmodule, only: budgettype
470  ! -- dummy
471  class(TspSsmType) :: this !< TspSsmType object
472  integer(I4B), intent(in) :: isuppress_output !< flag to suppress output
473  type(BudgetType), intent(inout) :: model_budget !< budget object for the GWT model
474  ! -- local
475  character(len=LENBUDROWLABEL) :: rowlabel
476  integer(I4B) :: ip
477  integer(I4B) :: i
478  integer(I4B) :: n
479  real(DP) :: rate
480  real(DP) :: rin
481  real(DP) :: rout
482  !
483  ! -- do for each flow package, unless it is being handled by an advanced
484  ! transport package
485  do ip = 1, this%fmi%nflowpack
486  !
487  ! -- cycle if package is being managed as an advanced package
488  if (this%fmi%iatp(ip) /= 0) cycle
489  !
490  ! -- Initialize the rate accumulators
491  rin = dzero
492  rout = dzero
493  !
494  ! -- do for each boundary
495  do i = 1, this%fmi%gwfpackages(ip)%nbound
496  n = this%fmi%gwfpackages(ip)%nodelist(i)
497  if (n <= 0) cycle
498  call this%ssm_term(ip, i, rrate=rate)
499  if (rate < dzero) then
500  rout = rout - rate
501  else
502  rin = rin + rate
503  end if
504  !
505  end do
506  !
507  rowlabel = 'SSM_'//adjustl(trim(this%fmi%flowpacknamearray(ip)))
508  call model_budget%addentry(rin, rout, delt, &
509  this%fmi%gwfpackages(ip)%budtxt, &
510  isuppress_output, rowlabel=rowlabel)
511  end do
This module contains the BudgetModule.
Definition: Budget.f90:20
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29
Derived type for the Budget object.
Definition: Budget.f90:39

◆ ssm_cq()

subroutine tspssmmodule::ssm_cq ( class(tspssmtype this,
real(dp), dimension(:), intent(inout), contiguous  flowja 
)

Calculate the resulting mass flow between the boundary and the connected GWT/GWE model cell. Update the diagonal position of the flowja array so that it ultimately contains the solute balance residual.

Parameters
thisTspSsmType object
[in,out]flowjaflow across each face in the model grid

Definition at line 431 of file tsp-ssm.f90.

432  ! -- dummy
433  class(TspSsmType) :: this !< TspSsmType object
434  real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow across each face in the model grid
435  ! -- local
436  integer(I4B) :: ip
437  integer(I4B) :: i
438  integer(I4B) :: n
439  integer(I4B) :: idiag
440  real(DP) :: rate
441  !
442  ! -- do for each flow package
443  do ip = 1, this%fmi%nflowpack
444  !
445  ! -- cycle if package is being managed as an advanced package
446  if (this%fmi%iatp(ip) /= 0) cycle
447  !
448  ! -- do for each boundary
449  do i = 1, this%fmi%gwfpackages(ip)%nbound
450  n = this%fmi%gwfpackages(ip)%nodelist(i)
451  if (n <= 0) cycle
452  call this%ssm_term(ip, i, rrate=rate)
453  idiag = this%dis%con%ia(n)
454  flowja(idiag) = flowja(idiag) + rate
455  !
456  end do
457  !
458  end do

◆ ssm_cr()

subroutine, public tspssmmodule::ssm_cr ( type(tspssmtype), pointer  ssmobj,
character(len=*), intent(in)  name_model,
integer(i4b), intent(in)  inunit,
integer(i4b), intent(in)  iout,
type(tspfmitype), intent(in), target  fmi,
real(dp), intent(in), pointer  eqnsclfac,
character(len=lenvarname), intent(in)  depvartype 
)

Create a new SSM package by defining names, allocating scalars and initializing the parser.

Parameters
ssmobjTspSsmType object
[in]name_modelname of the model
[in]inunitfortran unit for input
[in]ioutfortran unit for output
[in]fmiTransport FMI package
[in]eqnsclfacgoverning equation scale factor
[in]depvartypedependent variable type ('concentration' or 'temperature')

Definition at line 82 of file tsp-ssm.f90.

84  ! -- dummy
85  type(TspSsmType), pointer :: ssmobj !< TspSsmType object
86  character(len=*), intent(in) :: name_model !< name of the model
87  integer(I4B), intent(in) :: inunit !< fortran unit for input
88  integer(I4B), intent(in) :: iout !< fortran unit for output
89  type(TspFmiType), intent(in), target :: fmi !< Transport FMI package
90  real(DP), intent(in), pointer :: eqnsclfac !< governing equation scale factor
91  character(len=LENVARNAME), intent(in) :: depvartype !< dependent variable type ('concentration' or 'temperature')
92  !
93  ! -- Create the object
94  allocate (ssmobj)
95  !
96  ! -- create name and memory path
97  call ssmobj%set_names(1, name_model, 'SSM', 'SSM')
98  !
99  ! -- Allocate scalars
100  call ssmobj%allocate_scalars()
101  !
102  ! -- Set variables
103  ssmobj%inunit = inunit
104  ssmobj%iout = iout
105  ssmobj%fmi => fmi
106  ssmobj%eqnsclfac => eqnsclfac
107  !
108  ! -- Initialize block parser
109  call ssmobj%parser%Initialize(ssmobj%inunit, ssmobj%iout)
110  !
111  ! -- Store pointer to labels associated with the current model so that the
112  ! package has access to the corresponding dependent variable type
113  ssmobj%depvartype = depvartype
Here is the caller graph for this function:

◆ ssm_da()

subroutine tspssmmodule::ssm_da ( class(tspssmtype this)

Deallocate the memory associated with this derived type

Parameters
thisTspSsmType object

Definition at line 646 of file tsp-ssm.f90.

647  ! -- modules
649  ! -- dummy
650  class(TspSsmType) :: this !< TspSsmType object
651  ! -- local
652  integer(I4B) :: ip
653  type(TspSpcType), pointer :: ssmiptr
654  !
655  ! -- Deallocate the ssmi objects if package was active
656  if (this%inunit > 0) then
657  do ip = 1, size(this%ssmivec)
658  if (this%isrctype(ip) == 3 .or. this%isrctype(ip) == 4) then
659  ssmiptr => this%ssmivec(ip)
660  call ssmiptr%spc_da()
661  end if
662  end do
663  deallocate (this%ssmivec)
664  end if
665  !
666  ! -- Deallocate arrays if package was active
667  if (this%inunit > 0) then
668  call mem_deallocate(this%iauxpak)
669  call mem_deallocate(this%isrctype)
670  this%ibound => null()
671  this%fmi => null()
672  end if
673  !
674  ! -- output table object
675  if (associated(this%outputtab)) then
676  call this%outputtab%table_da()
677  deallocate (this%outputtab)
678  nullify (this%outputtab)
679  end if
680  !
681  ! -- Scalars
682  call mem_deallocate(this%nbound)
683  !
684  ! -- deallocate parent
685  call this%NumericalPackageType%da()

◆ ssm_df()

subroutine tspssmmodule::ssm_df ( class(tspssmtype this)

This routine is called from gwt_df(), but does not do anything because df is typically used to set up dimensions. For the ssm package, the total number of ssm entries is defined by the flow model.

Parameters
thisTspSsmType object

Definition at line 122 of file tsp-ssm.f90.

123  ! -- modules
125  ! -- dummy
126  class(TspSsmType) :: this !< TspSsmType object

◆ ssm_fc()

subroutine tspssmmodule::ssm_fc ( class(tspssmtype this,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(:), intent(in)  idxglo,
real(dp), dimension(:), intent(inout)  rhs 
)

This routine adds the effects of the SSM to the matrix equations by updating the a matrix and right-hand side vector.

Definition at line 389 of file tsp-ssm.f90.

390  ! -- dummy
391  class(TspSsmType) :: this
392  class(MatrixBaseType), pointer :: matrix_sln
393  integer(I4B), intent(in), dimension(:) :: idxglo
394  real(DP), intent(inout), dimension(:) :: rhs
395  ! -- local
396  integer(I4B) :: ip
397  integer(I4B) :: i
398  integer(I4B) :: n
399  integer(I4B) :: idiag
400  integer(I4B) :: nflowpack
401  integer(I4B) :: nbound
402  real(DP) :: hcofval
403  real(DP) :: rhsval
404  !
405  ! -- do for each flow package
406  nflowpack = this%fmi%nflowpack
407  do ip = 1, nflowpack
408  if (this%fmi%iatp(ip) /= 0) cycle
409  !
410  ! -- do for each entry in package (ip)
411  nbound = this%fmi%gwfpackages(ip)%nbound
412  do i = 1, nbound
413  n = this%fmi%gwfpackages(ip)%nodelist(i)
414  if (n <= 0) cycle
415  call this%ssm_term(ip, i, rhsval=rhsval, hcofval=hcofval)
416  idiag = idxglo(this%dis%con%ia(n))
417  call matrix_sln%add_value_pos(idiag, hcofval)
418  rhs(n) = rhs(n) + rhsval
419  !
420  end do
421  !
422  end do

◆ ssm_ot_flow()

subroutine tspssmmodule::ssm_ot_flow ( class(tspssmtype this,
integer(i4b), intent(in)  icbcfl,
integer(i4b), intent(in)  ibudfl,
integer(i4b), intent(in)  icbcun 
)

Based on user-specified controls, print SSM mass flow rates to the GWT listing file and/or write the SSM mass flow rates to the GWT binary budget file.

Parameters
thisTspSsmType object
[in]icbcflflag for writing binary budget terms
[in]ibudflflag for printing budget terms to list file
[in]icbcunfortran unit number for binary budget file

Definition at line 520 of file tsp-ssm.f90.

521  ! -- modules
522  use tdismodule, only: kstp, kper
524  ! -- dummy
525  class(TspSsmType) :: this !< TspSsmType object
526  integer(I4B), intent(in) :: icbcfl !< flag for writing binary budget terms
527  integer(I4B), intent(in) :: ibudfl !< flag for printing budget terms to list file
528  integer(I4B), intent(in) :: icbcun !< fortran unit number for binary budget file
529  ! -- local
530  character(len=LINELENGTH) :: title
531  integer(I4B) :: node, nodeu
532  character(len=20) :: nodestr
533  integer(I4B) :: maxrows
534  integer(I4B) :: ip
535  integer(I4B) :: i, n2, ibinun
536  real(DP) :: rrate
537  real(DP) :: qssm
538  real(DP) :: cssm
539  integer(I4B) :: naux
540  real(DP), dimension(0) :: auxrow
541  character(len=LENAUXNAME), dimension(0) :: auxname
542  ! -- for observations
543  character(len=LENBOUNDNAME) :: bname
544  ! -- formats
545  character(len=*), parameter :: fmttkk = &
546  &"(1X,/1X,A,' PERIOD ',I0,' STEP ',I0)"
547  !
548  ! -- initialize
549  naux = 0
550  maxrows = 0
551  if (ibudfl /= 0 .and. this%iprflow /= 0) then
552  call this%outputtab%set_kstpkper(kstp, kper)
553  do ip = 1, this%fmi%nflowpack
554  if (this%fmi%iatp(ip) /= 0) cycle
555  !
556  ! -- do for each boundary
557  do i = 1, this%fmi%gwfpackages(ip)%nbound
558  node = this%fmi%gwfpackages(ip)%nodelist(i)
559  if (node > 0) then
560  maxrows = maxrows + 1
561  end if
562  end do
563  end do
564  if (maxrows > 0) then
565  call this%outputtab%set_maxbound(maxrows)
566  end if
567  title = 'SSM PACKAGE ('//trim(this%packName)// &
568  ') FLOW RATES'
569  call this%outputtab%set_title(title)
570  end if
571  !
572  ! -- Set unit number for binary output
573  if (this%ipakcb < 0) then
574  ibinun = icbcun
575  else if (this%ipakcb == 0) then
576  ibinun = 0
577  else
578  ibinun = this%ipakcb
579  end if
580  if (icbcfl == 0) ibinun = 0
581  !
582  ! -- If cell-by-cell flows will be saved as a list, write header.
583  if (ibinun /= 0) then
584  call this%dis%record_srcdst_list_header(text, this%name_model, &
585  this%name_model, this%name_model, &
586  this%packName, naux, auxname, &
587  ibinun, this%nbound, this%iout)
588  end if
589  !
590  ! -- If no boundaries, skip flow calculations.
591  if (this%nbound > 0) then
592  !
593  ! -- Loop through each boundary calculating flow.
594  do ip = 1, this%fmi%nflowpack
595  if (this%fmi%iatp(ip) /= 0) cycle
596  !
597  ! -- do for each boundary
598  do i = 1, this%fmi%gwfpackages(ip)%nbound
599  !
600  ! -- Calculate rate for this entry
601  node = this%fmi%gwfpackages(ip)%nodelist(i)
602  if (node <= 0) cycle
603  call this%ssm_term(ip, i, rrate=rrate, qssm=qssm, cssm=cssm)
604  !
605  ! -- Print the individual rates if the budget is being printed
606  ! and PRINT_FLOWS was specified (this%iprflow<0)
607  if (ibudfl /= 0) then
608  if (this%iprflow /= 0) then
609  !
610  ! -- set nodestr and write outputtab table
611  nodeu = this%dis%get_nodeuser(node)
612  call this%dis%nodeu_to_string(nodeu, nodestr)
613  bname = this%fmi%gwfpackages(ip)%name
614  call this%outputtab%add_term(i)
615  call this%outputtab%add_term(trim(adjustl(nodestr)))
616  call this%outputtab%add_term(qssm)
617  call this%outputtab%add_term(cssm)
618  call this%outputtab%add_term(rrate)
619  call this%outputtab%add_term(bname)
620  end if
621  end if
622  !
623  ! -- If saving cell-by-cell flows in list, write flow
624  if (ibinun /= 0) then
625  n2 = i
626  call this%dis%record_mf6_list_entry(ibinun, node, n2, rrate, &
627  naux, auxrow, &
628  olconv2=.false.)
629  end if
630  !
631  end do
632  !
633  end do
634  end if
635  if (ibudfl /= 0) then
636  if (this%iprflow /= 0) then
637  write (this%iout, '(1x)')
638  end if
639  end if
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter lenpackagename
maximum length of the package name
Definition: Constants.f90:23
integer(i4b), parameter lenauxname
maximum length of a aux variable
Definition: Constants.f90:35
integer(i4b), parameter lenboundname
maximum length of a bound name
Definition: Constants.f90:36
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23

◆ ssm_rp()

subroutine tspssmmodule::ssm_rp ( class(tspssmtype this)

This routine is called from gwt_rp(). It is called at the beginning of each stress period. If any SPC input files are used to provide source and sink concentrations (or temperatures), then period blocks for the current stress period are read.

Parameters
thisTspSsmType object

Definition at line 187 of file tsp-ssm.f90.

188  ! -- dummy
189  class(TspSsmType) :: this !< TspSsmType object
190  ! -- local
191  integer(I4B) :: ip
192  type(TspSpcType), pointer :: ssmiptr
193  ! -- formats
194  !
195  ! -- Call the rp method on any ssm input files
196  do ip = 1, this%fmi%nflowpack
197  if (this%fmi%iatp(ip) /= 0) cycle
198  if (this%isrctype(ip) == 3 .or. this%isrctype(ip) == 4) then
199  ssmiptr => this%ssmivec(ip)
200  call ssmiptr%spc_rp()
201  end if
202  end do

◆ ssm_term()

subroutine tspssmmodule::ssm_term ( class(tspssmtype this,
integer(i4b), intent(in)  ipackage,
integer(i4b), intent(in)  ientry,
real(dp), intent(out), optional  rrate,
real(dp), intent(out), optional  rhsval,
real(dp), intent(out), optional  hcofval,
real(dp), intent(out), optional  cssm,
real(dp), intent(out), optional  qssm 
)

This is the primary SSM routine that calculates the matrix coefficient and right-hand-side value for any package and package entry. It returns several different optional variables that are used throughout this package to update matrix terms, budget calculations, and output tables.

Parameters
thisTspSsmType
[in]ipackagepackage number
[in]ientrybound number
[out]rratecalculated mass flow rate
[out]rhsvalcalculated rhs value
[out]hcofvalcalculated hcof value
[out]cssmcalculated source concentration/temperature depending on flow direction
[out]qssmwater flow rate into model cell from boundary package

Definition at line 256 of file tsp-ssm.f90.

258  ! -- dummy
259  class(TspSsmType) :: this !< TspSsmType
260  integer(I4B), intent(in) :: ipackage !< package number
261  integer(I4B), intent(in) :: ientry !< bound number
262  real(DP), intent(out), optional :: rrate !< calculated mass flow rate
263  real(DP), intent(out), optional :: rhsval !< calculated rhs value
264  real(DP), intent(out), optional :: hcofval !< calculated hcof value
265  real(DP), intent(out), optional :: cssm !< calculated source concentration/temperature depending on flow direction
266  real(DP), intent(out), optional :: qssm !< water flow rate into model cell from boundary package
267  ! -- local
268  logical(LGP) :: lauxmixed
269  integer(I4B) :: n
270  integer(I4B) :: nbound_flow
271  real(DP) :: qbnd
272  real(DP) :: ctmp
273  real(DP) :: omega
274  real(DP) :: hcoftmp
275  real(DP) :: rhstmp
276  !
277  ! -- retrieve node number, qbnd and iauxpos
278  hcoftmp = dzero
279  rhstmp = dzero
280  ctmp = dzero
281  qbnd = dzero
282  nbound_flow = this%fmi%gwfpackages(ipackage)%nbound
283  n = this%fmi%gwfpackages(ipackage)%nodelist(ientry)
284  !
285  ! -- If cell is active (ibound > 0) then calculate values
286  if (this%ibound(n) > 0) then
287  !
288  ! -- retrieve qbnd and iauxpos
289  qbnd = this%fmi%gwfpackages(ipackage)%get_flow(ientry)
290  call this%get_ssm_conc(ipackage, ientry, nbound_flow, ctmp, lauxmixed)
291  !
292  ! -- assign values for hcoftmp, rhstmp, and ctmp for subsequent assignment
293  ! of hcof, rhs, and rate
294  if (.not. lauxmixed) then
295  !
296  ! -- If qbnd is positive, then concentration (or temperature) represents
297  ! the inflow concentration (or temperature). If qbnd is negative,
298  ! then the outflow concentration (or temperature) is set equal to the
299  ! simulated cell's concentration (or temperature).
300  if (qbnd >= dzero) then
301  omega = dzero ! rhs
302  else
303  ctmp = this%cnew(n)
304  omega = done ! lhs
305  if (ctmp < dzero) then
306  omega = dzero ! concentration/temperature is negative, so set mass flux to zero
307  end if
308  end if
309  else
310  !
311  ! -- lauxmixed value indicates that this is a mixed sink type where
312  ! the concentration (or temperature) value represents the injected
313  ! concentration (or temperature) if qbnd is positive. If qbnd is
314  ! negative, then the withdrawn water is equal to the minimum of the
315  ! aux concentration (or temperature) and the cell concentration
316  ! (or temperature).
317  if (qbnd >= dzero) then
318  omega = dzero ! rhs (ctmp is aux value)
319  else
320  if (ctmp < this%cnew(n)) then
321  omega = dzero ! rhs (ctmp is aux value)
322  else
323  omega = done ! lhs (ctmp is cell concentration)
324  ctmp = this%cnew(n)
325  end if
326  end if
327  end if
328  !
329  ! -- Add terms based on qbnd sign
330  if (qbnd <= dzero) then
331  hcoftmp = qbnd * omega * this%eqnsclfac
332  else
333  rhstmp = -qbnd * ctmp * (done - omega) * this%eqnsclfac
334  end if
335  !
336  ! -- end of active ibound
337  end if
338  !
339  ! -- set requested values
340  if (present(hcofval)) hcofval = hcoftmp
341  if (present(rhsval)) rhsval = rhstmp
342  if (present(rrate)) rrate = hcoftmp * ctmp - rhstmp
343  if (present(cssm)) cssm = ctmp
344  if (present(qssm)) qssm = qbnd

Variable Documentation

◆ ftype

character(len=lenftype) tspssmmodule::ftype = 'SSM'

Definition at line 28 of file tsp-ssm.f90.

28  character(len=LENFTYPE) :: ftype = 'SSM'

◆ text

character(len=lenpackagename) tspssmmodule::text = ' SOURCE-SINK MIX'

Definition at line 29 of file tsp-ssm.f90.

29  character(len=LENPACKAGENAME) :: text = ' SOURCE-SINK MIX'