diff --git a/src/Model/GroundWaterFlow/gwf3dis8.f90 b/src/Model/GroundWaterFlow/gwf3dis8.f90 index d165c275a3f..2e1b1e08b9e 100644 --- a/src/Model/GroundWaterFlow/gwf3dis8.f90 +++ b/src/Model/GroundWaterFlow/gwf3dis8.f90 @@ -3,7 +3,7 @@ module GwfDisModule use ArrayReadersModule, only: ReadArray use KindModule, only: DP, I4B use ConstantsModule, only: LINELENGTH, DHALF, DZERO, LENMEMPATH, LENVARNAME - use BaseDisModule, only: DisBaseType + use BaseDisModule, only: DisBaseType, dis_da use InputOutputModule, only: get_node, URWORD, ulasav, ulaprufw, ubdsv1, & ubdsv06 use SimModule, only: count_errors, store_error, store_error_unit, & @@ -56,8 +56,8 @@ module GwfDisModule procedure :: log_griddata procedure :: grid_finalize procedure :: write_grb - procedure :: allocate_scalars - procedure :: allocate_arrays + procedure :: allocate_scalars => allocate_scalars_dis + procedure :: allocate_arrays => allocate_arrays_dis ! ! -- Read a node-sized model array (reduced or not) procedure :: read_int_array @@ -159,7 +159,7 @@ subroutine dis3d_da(this) call memorylist_remove(this%name_model, 'DIS', idm_context) ! ! -- DisBaseType deallocate - call this%DisBaseType%dis_da() + call dis_da(this) ! ! -- Deallocate scalars call mem_deallocate(this%nlay) @@ -866,7 +866,7 @@ function get_nodenumber_idx3(this, k, i, j, icheck) & return end function get_nodenumber_idx3 - subroutine allocate_scalars(this, name_model, input_mempath) + subroutine allocate_scalars_dis(this, name_model, input_mempath) ! ****************************************************************************** ! allocate_scalars -- Allocate and initialize scalars ! ****************************************************************************** @@ -881,7 +881,7 @@ subroutine allocate_scalars(this, name_model, input_mempath) ! ------------------------------------------------------------------------------ ! ! -- Allocate parent scalars - call this%DisBaseType%allocate_scalars(name_model, input_mempath) + call this%allocate_scalars_default(name_model, input_mempath) ! ! -- Allocate call mem_allocate(this%nlay, 'NLAY', this%memoryPath) @@ -896,9 +896,9 @@ subroutine allocate_scalars(this, name_model, input_mempath) ! ! -- Return return - end subroutine allocate_scalars + end subroutine allocate_scalars_dis - subroutine allocate_arrays(this) + subroutine allocate_arrays_dis(this) ! ****************************************************************************** ! allocate_arrays -- Allocate arrays ! ****************************************************************************** @@ -912,7 +912,7 @@ subroutine allocate_arrays(this) ! ------------------------------------------------------------------------------ ! ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) - call this%DisBaseType%allocate_arrays() + call this%allocate_arrays_default() ! ! -- Allocate arrays for GwfDisType if (this%nodes < this%nodesuser) then @@ -931,7 +931,7 @@ subroutine allocate_arrays(this) ! ! -- Return return - end subroutine allocate_arrays + end subroutine allocate_arrays_dis function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & flag_string, allow_zero) result(nodeu) diff --git a/src/Model/GroundWaterFlow/gwf3disu8.f90 b/src/Model/GroundWaterFlow/gwf3disu8.f90 index d33bf55eb2a..0fbb34eb817 100644 --- a/src/Model/GroundWaterFlow/gwf3disu8.f90 +++ b/src/Model/GroundWaterFlow/gwf3disu8.f90 @@ -9,7 +9,7 @@ module GwfDisuModule use SimModule, only: count_errors, store_error, store_error_unit, & store_error_filename use SimVariablesModule, only: errmsg - use BaseDisModule, only: DisBaseType + use BaseDisModule, only: DisBaseType, dis_da use MemoryManagerModule, only: mem_allocate use TdisModule, only: kstp, kper, pertim, totim, delt @@ -59,8 +59,8 @@ module GwfDisuModule procedure, public :: record_array procedure, public :: record_srcdst_list_header ! -- private - procedure :: allocate_scalars - procedure :: allocate_arrays + procedure :: allocate_scalars => allocate_scalars_disu + procedure :: allocate_arrays => allocate_arrays_disu procedure :: allocate_arrays_mem procedure :: source_options procedure :: source_dimensions @@ -78,6 +78,9 @@ module GwfDisuModule ! -- Read a node-sized model array (reduced or not) procedure :: read_int_array procedure :: read_dbl_array + ! + procedure :: nlarray_to_nodelist + procedure :: read_layer_array end type GwfDisuType contains @@ -482,7 +485,7 @@ subroutine disu_da(this) call mem_deallocate(this%nodereduced) ! ! -- DisBaseType deallocate - call this%DisBaseType%dis_da() + call dis_da(this) ! ! -- Return return @@ -1305,7 +1308,7 @@ subroutine get_dis_type(this, dis_type) end subroutine get_dis_type - subroutine allocate_scalars(this, name_model, input_mempath) + subroutine allocate_scalars_disu(this, name_model, input_mempath) ! ****************************************************************************** ! allocate_scalars -- Allocate and initialize scalar variables in this class ! ****************************************************************************** @@ -1322,7 +1325,7 @@ subroutine allocate_scalars(this, name_model, input_mempath) ! ------------------------------------------------------------------------------ ! ! -- Allocate parent scalars - call this%DisBaseType%allocate_scalars(name_model, input_mempath) + call this%allocate_scalars_default(name_model, input_mempath) ! ! -- Allocate variables for DISU call mem_allocate(this%njausr, 'NJAUSR', this%memoryPath) @@ -1340,9 +1343,9 @@ subroutine allocate_scalars(this, name_model, input_mempath) ! ! -- Return return - end subroutine allocate_scalars + end subroutine allocate_scalars_disu - subroutine allocate_arrays(this) + subroutine allocate_arrays_disu(this) ! ****************************************************************************** ! allocate_arrays -- Read discretization information from file ! ****************************************************************************** @@ -1357,7 +1360,7 @@ subroutine allocate_arrays(this) ! ------------------------------------------------------------------------------ ! ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) - call this%DisBaseType%allocate_arrays() + call this%allocate_arrays_default() ! ! -- Allocate arrays in DISU if (this%nodes < this%nodesuser) then @@ -1374,7 +1377,7 @@ subroutine allocate_arrays(this) ! ! -- Return return - end subroutine allocate_arrays + end subroutine allocate_arrays_disu subroutine allocate_arrays_mem(this) use MemoryManagerModule, only: mem_allocate @@ -1809,4 +1812,43 @@ function CastAsDisuType(dis) result(disu) end function CastAsDisuType + ! todo: routines below are not used for disu, remove from DisBaseType? + + subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname) + ! -- modules + use SimModule, only: store_error + use ConstantsModule, only: LINELENGTH + ! -- dummy + class(GwfDisuType) :: this + integer(I4B), intent(in) :: maxbnd + integer(I4B), dimension(:), pointer, contiguous :: darray + integer(I4B), dimension(maxbnd), intent(inout) :: nodelist + integer(I4B), intent(inout) :: nbound + character(len=*), intent(in) :: aname + ! + errmsg = 'Programmer error: nlarray_to_nodelist called for DISU grid.' + call store_error(errmsg, terminate=.TRUE.) + + end subroutine nlarray_to_nodelist + + subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, & + icolbnd, aname, inunit, iout) + ! -- dummy + class(GwfDisuType) :: this + integer(I4B), intent(in) :: ncolbnd + integer(I4B), intent(in) :: maxbnd + integer(I4B), dimension(maxbnd) :: nodelist + real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray + integer(I4B), intent(in) :: icolbnd + character(len=*), intent(in) :: aname + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + ! + ! + errmsg = 'Programmer error: read_layer_array called for DISU grid.' + call store_error(errmsg, terminate=.TRUE.) + ! + ! -- return + end subroutine read_layer_array + end module GwfDisuModule diff --git a/src/Model/GroundWaterFlow/gwf3disv8.f90 b/src/Model/GroundWaterFlow/gwf3disv8.f90 index 505cd680e89..f9ae8079b38 100644 --- a/src/Model/GroundWaterFlow/gwf3disv8.f90 +++ b/src/Model/GroundWaterFlow/gwf3disv8.f90 @@ -4,7 +4,7 @@ module GwfDisvModule use KindModule, only: DP, I4B use ConstantsModule, only: LINELENGTH, LENMEMPATH, LENVARNAME, DZERO, DONE, & DHALF - use BaseDisModule, only: DisBaseType + use BaseDisModule, only: DisBaseType, dis_da use InputOutputModule, only: get_node, URWORD, ulasav, ulaprufw, ubdsv1, & ubdsv06 use SimModule, only: count_errors, store_error, store_error_unit, & @@ -62,8 +62,8 @@ module GwfDisvModule procedure :: grid_finalize procedure :: connect procedure :: write_grb - procedure :: allocate_scalars - procedure :: allocate_arrays + procedure :: allocate_scalars => allocate_scalars_disv + procedure :: allocate_arrays => allocate_arrays_disv procedure :: get_cell2d_area ! procedure :: read_int_array @@ -180,7 +180,7 @@ subroutine disv_da(this) context=idm_context) ! ! -- DisBaseType deallocate - call this%DisBaseType%dis_da() + call dis_da(this) ! ! -- Deallocate scalars call mem_deallocate(this%nlay) @@ -1234,7 +1234,7 @@ subroutine get_dis_type(this, dis_type) end subroutine get_dis_type - subroutine allocate_scalars(this, name_model, input_mempath) + subroutine allocate_scalars_disv(this, name_model, input_mempath) ! ****************************************************************************** ! allocate_scalars -- Allocate and initialize scalars ! ****************************************************************************** @@ -1250,7 +1250,7 @@ subroutine allocate_scalars(this, name_model, input_mempath) ! ------------------------------------------------------------------------------ ! ! -- Allocate parent scalars - call this%DisBaseType%allocate_scalars(name_model, input_mempath) + call this%allocate_scalars_default(name_model, input_mempath) ! ! -- Allocate call mem_allocate(this%nlay, 'NLAY', this%memoryPath) @@ -1265,9 +1265,9 @@ subroutine allocate_scalars(this, name_model, input_mempath) ! ! -- Return return - end subroutine allocate_scalars + end subroutine allocate_scalars_disv - subroutine allocate_arrays(this) + subroutine allocate_arrays_disv(this) ! ****************************************************************************** ! allocate_arrays -- Allocate arrays ! ****************************************************************************** @@ -1281,7 +1281,7 @@ subroutine allocate_arrays(this) ! ------------------------------------------------------------------------------ ! ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) - call this%DisBaseType%allocate_arrays() + call this%allocate_arrays_default() ! ! -- Allocate arrays for GwfDisvType if (this%nodes < this%nodesuser) then @@ -1298,7 +1298,7 @@ subroutine allocate_arrays(this) ! ! -- Return return - end subroutine allocate_arrays + end subroutine allocate_arrays_disv function get_cell2d_area(this, icell2d) result(area) ! ****************************************************************************** diff --git a/src/Model/ModelUtilities/DiscretizationBase.f90 b/src/Model/ModelUtilities/DiscretizationBase.f90 index adcbe831379..19c429df825 100644 --- a/src/Model/ModelUtilities/DiscretizationBase.f90 +++ b/src/Model/ModelUtilities/DiscretizationBase.f90 @@ -21,8 +21,11 @@ module BaseDisModule private public :: DisBaseType public :: dis_transform_xy + public :: dis_da + + ! Abstract base type for grid discretizations. + type, abstract :: DisBaseType - type :: DisBaseType character(len=LENMEMPATH) :: memoryPath !< path for memory allocation character(len=LENMEMPATH) :: input_mempath = '' !< input context mempath character(len=LENMODELNAME), pointer :: name_model => null() !< name of the model @@ -53,92 +56,347 @@ module BaseDisModule integer(I4B), dimension(:), pointer, contiguous :: nodereduced => null() !< (size:nodesuser)contains reduced nodenumber (size 0 if not reduced); -1 means vertical pass through, 0 is idomain = 0 integer(I4B), dimension(:), pointer, contiguous :: nodeuser => null() !< (size:nodes) given a reduced nodenumber, provide the user nodenumber (size 0 if not reduced) contains - procedure :: dis_df + + ! lifecycle procedures + procedure(dis_df), deferred :: dis_df procedure :: dis_ac procedure :: dis_mc procedure :: dis_ar procedure :: dis_da - ! -- helper functions - ! - ! -- get_nodenumber is an overloaded integer function that will always - ! return the reduced nodenumber. For all grids, get_nodenumber can - ! be passed the user nodenumber. For some other grids, it can also - ! be passed an index. For dis3d the index is k, i, j, and for - ! disv the index is k, n. + + ! allocate scalar and array variables + procedure :: allocate_scalars_default + procedure :: allocate_scalars => allocate_scalars_default + procedure :: allocate_arrays_default + procedure :: allocate_arrays => allocate_arrays_default + + ! get_nodenumber is an overloaded integer function that will always + ! return the reduced nodenumber. For all grids, get_nodenumber can + ! be passed the user nodenumber. For some other grids, it can also + ! be passed an index. For dis3d the index is k, i, j, and for + ! disv the index is k, n. generic :: get_nodenumber => get_nodenumber_idx1, & get_nodenumber_idx2, & get_nodenumber_idx3 procedure :: get_nodenumber_idx1 procedure :: get_nodenumber_idx2 procedure :: get_nodenumber_idx3 - procedure :: get_nodeuser - procedure :: nodeu_to_string - procedure :: nodeu_to_array - procedure :: nodeu_from_string - procedure :: nodeu_from_cellid + + ! conversion between user and reduced node numbers, and utilities + ! to parse user node numbers from lines of text and cellid strings + procedure :: get_nodeuser ! todo rename to noder_from_nodeu? procedure :: noder_from_string procedure :: noder_from_cellid - procedure :: connection_normal - procedure :: connection_vector - procedure :: get_dis_type - procedure :: supports_layers - procedure :: allocate_scalars - procedure :: allocate_arrays - procedure :: get_ncpl + procedure :: noder_to_string + procedure :: noder_to_array + procedure(nodeu_from_string), deferred :: nodeu_from_string + procedure(nodeu_from_cellid), deferred :: nodeu_from_cellid + procedure(nodeu_to_string), deferred :: nodeu_to_string + procedure(nodeu_to_array), deferred :: nodeu_to_array + procedure(nlarray_to_nodelist), deferred :: nlarray_to_nodelist + + ! calculate normal and unit vector components for shared faces + procedure(connection_normal), deferred :: connection_normal + procedure(connection_vector), deferred :: connection_vector + + ! basic discretization info + procedure(get_dis_type), deferred :: get_dis_type + procedure(supports_layers), deferred :: supports_layers + + ! cell-related information + procedure(get_ncpl), deferred :: get_ncpl procedure :: get_cell_volume - procedure :: write_grb - ! + procedure :: get_area + procedure :: get_area_factor + procedure :: highest_active ! todo rename to get_highest_active? + + ! write a binary grid file + procedure(write_grb), deferred :: write_grb + + ! generic procedure to read an array from a line of text + generic :: read_grid_array => read_int_array, read_dbl_array procedure :: read_int_array procedure :: read_dbl_array - generic, public :: read_grid_array => read_int_array, read_dbl_array - procedure, public :: read_layer_array + + ! read a 2D double array from file + procedure(read_layer_array), deferred, public :: read_layer_array + + ! generic procedure to fill arrays indexed by reduced node number procedure :: fill_int_array procedure :: fill_dbl_array - generic, public :: fill_grid_array => fill_int_array, fill_dbl_array - procedure, public :: read_list - ! - procedure, public :: record_array - procedure, public :: record_connection_array - procedure, public :: noder_to_string - procedure, public :: noder_to_array - procedure, public :: record_srcdst_list_header + generic :: fill_grid_array => fill_int_array, fill_dbl_array + + ! read from input files and write to output files + procedure :: read_list + procedure(record_srcdst_list_header), deferred :: record_srcdst_list_header procedure, private :: record_srcdst_list_entry - generic, public :: record_mf6_list_entry => record_srcdst_list_entry - procedure, public :: nlarray_to_nodelist - procedure, public :: highest_active - procedure, public :: get_area - procedure, public :: get_area_factor + procedure(record_array), deferred :: record_array + procedure :: record_connection_array + generic :: record_mf6_list_entry => record_srcdst_list_entry end type DisBaseType -contains + abstract interface + !> @brief Get the discretization type. + subroutine get_dis_type(this, dis_type) + import DisBaseType + class(DisBaseType), intent(in) :: this + character(len=*), intent(out) :: dis_type + end subroutine get_dis_type + end interface - subroutine dis_df(this) -! ****************************************************************************** -! dis_df -- Read discretization information from DISU input file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- dummy - class(DisBaseType) :: this -! ------------------------------------------------------------------------------ - ! - call store_error('Program error: DisBaseType method dis_df not & - &implemented.', terminate=.TRUE.) - ! - ! -- Return - return - end subroutine dis_df + abstract interface + !> @brief Define the discretization. + subroutine dis_df(this) + import DisBaseType + class(DisBaseType) :: this + end subroutine dis_df + end interface + + abstract interface + !> @brief Write binary grid file. Called from AR procedure. + subroutine write_grb(this, icelltype) + import DisBaseType + import I4B + class(DisBaseType) :: this + integer(I4B), dimension(:), intent(in) :: icelltype + end subroutine write_grb + end interface + + abstract interface + !> @brief Convert user node number to a string in the form of (nodenumber) or (k,i,j). + subroutine nodeu_to_string(this, nodeu, str) + import DisBaseType + import I4B + class(DisBaseType) :: this + integer(I4B), intent(in) :: nodeu + character(len=*), intent(inout) :: str + end subroutine nodeu_to_string + end interface + + abstract interface + !> @brief Convert user node number to array with format (nodenumber) or (k,j) or (k,i,j). + subroutine nodeu_to_array(this, nodeu, arr) + import DisBaseType + import I4B + class(DisBaseType) :: this + integer(I4B), intent(in) :: nodeu + integer(I4B), dimension(:), intent(inout) :: arr + end subroutine nodeu_to_array + end interface + + abstract interface + !> @brief Calculate normal vector components for shared cell faces. + !! + !! Reduced nodenumber cell is noden and its shared face with is cell nodem. + !! ihc is the horizontal connection flag. + subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, & + ipos) + import DisBaseType + import I4B + import DP + class(DisBaseType) :: this + integer(I4B), intent(in) :: noden + integer(I4B), intent(in) :: nodem + integer(I4B), intent(in) :: ihc + real(DP), intent(inout) :: xcomp + real(DP), intent(inout) :: ycomp + real(DP), intent(inout) :: zcomp + integer(I4B), intent(in) :: ipos + end subroutine connection_normal + end interface + + abstract interface + !> @brief Calculate unit vector components for shared cell faces. + !! + !! Reduced nodenumber cell is noden and its neighbor cell is nodem. + !! Saturations for these cells are also required so the vertical + !! position of the cell centers can be calculated. ihc is the + !! horizontal flag. Also return the straight-line connection length. + subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & + xcomp, ycomp, zcomp, conlen) + import DisBaseType + import I4B + import DP + class(DisBaseType) :: this + integer(I4B), intent(in) :: noden + integer(I4B), intent(in) :: nodem + logical, intent(in) :: nozee + real(DP), intent(in) :: satn + real(DP), intent(in) :: satm + integer(I4B), intent(in) :: ihc + real(DP), intent(inout) :: xcomp + real(DP), intent(inout) :: ycomp + real(DP), intent(inout) :: zcomp + real(DP), intent(inout) :: conlen + end subroutine connection_vector + end interface + + abstract interface + !> @brief Parse a user nodenumber from a line of text. + !! + !! The model is unstructured; just read user nodenumber. + !! If flag_string argument is present and true, the first token in string + !! is allowed to be a string (e.g. boundary name). In this case, if a string + !! is encountered, return value as -2. + function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & + flag_string, allow_zero) result(nodeu) + import DisBaseType + import I4B + class(DisBaseType) :: this + integer(I4B), intent(inout) :: lloc + integer(I4B), intent(inout) :: istart + integer(I4B), intent(inout) :: istop + integer(I4B), intent(in) :: in + integer(I4B), intent(in) :: iout + character(len=*), intent(inout) :: line + logical, optional, intent(in) :: flag_string + logical, optional, intent(in) :: allow_zero + integer(I4B) :: nodeu + end function nodeu_from_string + end interface + + abstract interface + !> @brief Parse a user nodenumber from a cell ID string. + !! + !! If flag_string argument is present and true, the first token in string + !! is allowed to be a string (e.g. boundary name). In this case, if a string + !! is encountered, return value as -2. + !! If allow_zero argument is present and true, if all indices equal zero, the + !! result can be zero. If allow_zero is false, a zero in any index causes an + !! error. + function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & + allow_zero) result(nodeu) + import DisBaseType + import I4B + class(DisBaseType) :: this + character(len=*), intent(inout) :: cellid + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + logical, optional, intent(in) :: flag_string + logical, optional, intent(in) :: allow_zero + integer(I4B) :: nodeu + end function nodeu_from_cellid + end interface + + abstract interface + !> @brief Indicates whether the discretization supports layers. + logical function supports_layers(this) + import DisBaseType + class(DisBaseType) :: this + end function supports_layers + end interface + abstract interface + !> @brief Get the number of cells per layer (nodes for DISU, since no layers). + function get_ncpl(this) + import DisBaseType + import I4B + integer(I4B) :: get_ncpl + class(DisBaseType) :: this + end function get_ncpl + end interface + + abstract interface + !> @brief Read a 2d double array into col icolbnd of darray. + !! For cells that are outside of the active domain, + !! do not copy the array value into darray. + subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, & + icolbnd, aname, inunit, iout) + import DisBaseType + import I4B + import DP + class(DisBaseType) :: this + integer(I4B), intent(in) :: ncolbnd + integer(I4B), intent(in) :: maxbnd + integer(I4B), dimension(maxbnd) :: nodelist + real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray + integer(I4B), intent(in) :: icolbnd + character(len=*), intent(in) :: aname + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + end subroutine read_layer_array + end interface + + abstract interface + !> @brief Record a double precision array. + !! + !! The array will be printed to an external file and/or written to an + !! unformatted external file depending on the argument specifications. + !! + !! SPECIFICATIONS: + !! darray is the double precision array to record + !! iout is the unit number for ascii output + !! iprint is a flag indicating whether or not to print the array + !! idataun is the unit number to which the array will be written in binary + !! form; if negative then do not write by layers, write entire array + !! aname is the text descriptor of the array + !! cdatafmp is the fortran format for writing the array + !! nvaluesp is the number of values per line for printing + !! nwidthp is the width of the number for printing + !! editdesc is the format type (I, G, F, S, E) + !! dinact is the double precision value to use for cells that are excluded + !! from the model domain + subroutine record_array(this, darray, iout, iprint, idataun, aname, & + cdatafmp, nvaluesp, nwidthp, editdesc, dinact) + import DisBaseType + import I4B + import DP + class(DisBaseType), intent(inout) :: this + real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray + integer(I4B), intent(in) :: iout + integer(I4B), intent(in) :: iprint + integer(I4B), intent(in) :: idataun + character(len=*), intent(in) :: aname + character(len=*), intent(in) :: cdatafmp + integer(I4B), intent(in) :: nvaluesp + integer(I4B), intent(in) :: nwidthp + character(len=*), intent(in) :: editdesc + real(DP), intent(in) :: dinact + end subroutine record_array + end interface + + abstract interface + !> @brief Record list header for imeth=6 + subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & + dstmodel, dstpackage, naux, auxtxt, & + ibdchn, nlist, iout) + import DisBaseType + import I4B + class(DisBaseType) :: this + character(len=16), intent(in) :: text + character(len=16), intent(in) :: textmodel + character(len=16), intent(in) :: textpackage + character(len=16), intent(in) :: dstmodel + character(len=16), intent(in) :: dstpackage + integer(I4B), intent(in) :: naux + character(len=16), dimension(:), intent(in) :: auxtxt + integer(I4B), intent(in) :: ibdchn + integer(I4B), intent(in) :: nlist + integer(I4B), intent(in) :: iout + end subroutine record_srcdst_list_header + end interface + + abstract interface + !> @brief Convert an integer array to nodelist. + !! For structured model, integer array is layer number; for unstructured + !! model, integer array is node number. + subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname) + import DisBaseType + import I4B + class(DisBaseType) :: this + integer(I4B), intent(in) :: maxbnd + integer(I4B), dimension(:), pointer, contiguous :: darray + integer(I4B), dimension(maxbnd), intent(inout) :: nodelist + integer(I4B), intent(inout) :: nbound + character(len=*), intent(in) :: aname + end subroutine nlarray_to_nodelist + end interface + +contains + + !> @brief Add connections to sparse based on cell connectivity. subroutine dis_ac(this, moffset, sparse) -! ****************************************************************************** -! dis_ac -- Add connections to sparse based on cell connectivity -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules use SparseModule, only: sparsematrix ! -- dummy class(DisBaseType) :: this @@ -146,7 +404,6 @@ subroutine dis_ac(this, moffset, sparse) type(sparsematrix), intent(inout) :: sparse ! -- local integer(I4B) :: i, j, ipos, iglo, jglo -! ------------------------------------------------------------------------------ ! do i = 1, this%nodes do ipos = this%con%ia(i), this%con%ia(i + 1) - 1 @@ -156,20 +413,10 @@ subroutine dis_ac(this, moffset, sparse) call sparse%addconnection(iglo, jglo, 1) end do end do - ! - ! -- Return - return end subroutine dis_ac + !> @brief Map cell connection positions in numerical solution coefficient matrix. subroutine dis_mc(this, moffset, idxglo, matrix_sln) -! ****************************************************************************** -! dis_mc -- Map the positions of cell connections in the numerical solution -! coefficient matrix. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: moffset @@ -177,7 +424,6 @@ subroutine dis_mc(this, moffset, idxglo, matrix_sln) class(MatrixBaseType), pointer :: matrix_sln ! -- local integer(I4B) :: i, j, ipos, iglo, jglo -! ------------------------------------------------------------------------------ ! do i = 1, this%nodes iglo = i + moffset @@ -187,26 +433,16 @@ subroutine dis_mc(this, moffset, idxglo, matrix_sln) idxglo(ipos) = matrix_sln%get_position(iglo, jglo) end do end do - ! - ! -- Return - return end subroutine dis_mc + !> @brief Allocate and read. subroutine dis_ar(this, icelltype) -! ****************************************************************************** -! dis_ar -- Called from AR procedure. Only task is to write binary grid file. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(DisBaseType) :: this integer(I4B), dimension(:), intent(in) :: icelltype ! -- local integer(I4B), dimension(:), allocatable :: ict integer(I4B) :: nu, nr -! ------------------------------------------------------------------------------ ! ! -- Expand icelltype to full grid; fill with 0 if cell is excluded allocate (ict(this%nodesuser)) @@ -220,45 +456,14 @@ subroutine dis_ar(this, icelltype) end do ! if (this%nogrb == 0) call this%write_grb(ict) - ! - ! -- Return - return end subroutine dis_ar - subroutine write_grb(this, icelltype) -! ****************************************************************************** -! write_grb -- Called from AR procedure. Only task is to write binary grid file. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - ! -- dummy - class(DisBaseType) :: this - integer(I4B), dimension(:), intent(in) :: icelltype - ! -- local -! ------------------------------------------------------------------------------ - ! - ! - call store_error('Program error: DisBaseType method write_grb not & - &implemented.', terminate=.TRUE.) - ! - ! -- Return - return - end subroutine write_grb - + !> @brief Deallocate discretization object. subroutine dis_da(this) -! ****************************************************************************** -! dis_da -- Deallocate discretization object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use MemoryManagerModule, only: mem_deallocate ! -- dummy class(DisBaseType) :: this -! ------------------------------------------------------------------------------ ! ! -- Strings deallocate (this%name_model) @@ -292,89 +497,27 @@ subroutine dis_da(this) ! -- Connections call this%con%con_da() deallocate (this%con) - ! - ! -- Return - return end subroutine dis_da - subroutine nodeu_to_string(this, nodeu, str) -! ****************************************************************************** -! nodeu_to_string -- Convert user node number to a string in the form of -! (nodenumber) or (k,i,j) -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- dummy - class(DisBaseType) :: this - integer(I4B), intent(in) :: nodeu - character(len=*), intent(inout) :: str - ! -- local -! ------------------------------------------------------------------------------ - ! - call store_error('Program error: DisBaseType method nodeu_to_string not & - &implemented.', terminate=.TRUE.) - ! - ! -- return - return - end subroutine nodeu_to_string - - subroutine nodeu_to_array(this, nodeu, arr) -! ****************************************************************************** -! nodeu_to_array -- Convert user node number to cellid and fill array with -! (nodenumber) or (k,j) or (k,i,j) -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- dummy - class(DisBaseType) :: this - integer(I4B), intent(in) :: nodeu - integer(I4B), dimension(:), intent(inout) :: arr - ! -- local -! ------------------------------------------------------------------------------ - ! - call store_error('Program error: DisBaseType method nodeu_to_array not & - &implemented.', terminate=.TRUE.) - ! - ! -- return - return - end subroutine nodeu_to_array - + !> @brief Get user nodenumber from reduced node number. function get_nodeuser(this, noder) result(nodenumber) -! ****************************************************************************** -! get_nodeuser -- Return the user nodenumber from the reduced node number -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- return integer(I4B) :: nodenumber ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: noder -! ------------------------------------------------------------------------------ ! if (this%nodes < this%nodesuser) then nodenumber = this%nodeuser(noder) else nodenumber = noder end if - ! - ! -- return - return end function get_nodeuser + !> @brief Get nodenumber from the user specified node number. + !! Overridde in child classes to map to a model node number, + !! optionally performing a check. function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) -! ****************************************************************************** -! get_nodenumber -- Return a nodenumber from the user specified node number -! with an option to perform a check. This subroutine -! can be overridden by child classes to perform mapping -! to a model node number -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: LINELENGTH use SimModule, only: store_error @@ -384,133 +527,42 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) integer(I4B), intent(in) :: icheck ! -- local integer(I4B) :: nodenumber -! ------------------------------------------------------------------------------ ! nodenumber = 0 call store_error('Program error: get_nodenumber_idx1 not implemented.', & terminate=.TRUE.) - ! - ! -- return - return end function get_nodenumber_idx1 + !> @brief Override in child classes to map to a model node number. function get_nodenumber_idx2(this, k, j, icheck) result(nodenumber) -! ****************************************************************************** -! get_nodenumber_idx2 -- This function should never be called. It must be -! overridden by a child class. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ use SimModule, only: store_error ! -- dummy class(DisBaseType), intent(in) :: this integer(I4B), intent(in) :: k, j integer(I4B), intent(in) :: icheck integer(I4B) :: nodenumber -! ------------------------------------------------------------------------------ ! nodenumber = 0 call store_error('Program error: get_nodenumber_idx2 not implemented.', & terminate=.TRUE.) - ! - ! -- Return - return end function get_nodenumber_idx2 + !> @brief This function will not be invoked for an unstructured + !! model, but it may be from a Discretization3dType model. function get_nodenumber_idx3(this, k, i, j, icheck) result(nodenumber) -! ****************************************************************************** -! get_nodenumber_idx3 -- This function will not be invoked for an unstructured -! model, but it may be from a Discretization3dType model. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ use SimModule, only: store_error ! -- dummy class(DisBaseType), intent(in) :: this integer(I4B), intent(in) :: k, i, j integer(I4B), intent(in) :: icheck integer(I4B) :: nodenumber -! ------------------------------------------------------------------------------ ! nodenumber = 0 call store_error('Program error: get_nodenumber_idx3 not implemented.', & terminate=.TRUE.) - ! - ! -- Return - return end function get_nodenumber_idx3 - subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, & - ipos) -! ****************************************************************************** -! connection_normal -- calculate the normal vector components for reduced -! nodenumber cell (noden) and its shared face with cell nodem. ihc is the -! horizontal connection flag. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use SimModule, only: store_error - ! -- dummy - class(DisBaseType) :: this - integer(I4B), intent(in) :: noden - integer(I4B), intent(in) :: nodem - integer(I4B), intent(in) :: ihc - real(DP), intent(inout) :: xcomp - real(DP), intent(inout) :: ycomp - real(DP), intent(inout) :: zcomp - integer(I4B), intent(in) :: ipos -! ------------------------------------------------------------------------------ - ! - call store_error('Program error: connection_normal not implemented.', & - terminate=.TRUE.) - ! - ! -- return - return - end subroutine connection_normal - - subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & - xcomp, ycomp, zcomp, conlen) -! ****************************************************************************** -! connection_vector -- calculate the unit vector components from reduced -! nodenumber cell (noden) to its neighbor cell (nodem). The saturation for -! for these cells are also required so that the vertical position of the cell -! cell centers can be calculated. ihc is the horizontal flag. Also return -! the straight-line connection length. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use SimModule, only: store_error - ! -- dummy - class(DisBaseType) :: this - integer(I4B), intent(in) :: noden - integer(I4B), intent(in) :: nodem - logical, intent(in) :: nozee - real(DP), intent(in) :: satn - real(DP), intent(in) :: satm - integer(I4B), intent(in) :: ihc - real(DP), intent(inout) :: xcomp - real(DP), intent(inout) :: ycomp - real(DP), intent(inout) :: zcomp - real(DP), intent(inout) :: conlen - ! -- local -! ------------------------------------------------------------------------------ - ! - call store_error('Program error: connection_vector not implemented.', & - terminate=.TRUE.) - ! - ! -- return - return - end subroutine connection_vector - - !> @brief get the x,y for a node transformed into - !! 'global coordinates' using xorigin, yorigin, angrot, - !< analogously to how flopy does this. + !> @brief Transform local to global coords using xorigin, yorigin, angrot. subroutine dis_transform_xy(x, y, xorigin, yorigin, angrot, xglo, yglo) real(DP), intent(in) :: x !< the cell-x coordinate to transform real(DP), intent(in) :: y !< the cell-y coordinate to transform @@ -535,30 +587,10 @@ subroutine dis_transform_xy(x, y, xorigin, yorigin, angrot, xglo, yglo) ! then _translate_ xglo = xglo + xorigin yglo = yglo + yorigin - end subroutine dis_transform_xy - !> @brief return discretization type - !< - subroutine get_dis_type(this, dis_type) - class(DisBaseType), intent(in) :: this - character(len=*), intent(out) :: dis_type - - ! suppress warning - dis_type = "Not implemented" - - call store_error('Program error: get_dis_type not implemented.', & - terminate=.TRUE.) - - end subroutine get_dis_type - - subroutine allocate_scalars(this, name_model, input_mempath) -! ****************************************************************************** -! allocate_scalars -- Allocate and initialize scalar variables in this class -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + !> @brief Allocate and initialize scalar variables. + subroutine allocate_scalars_default(this, name_model, input_mempath) ! -- modules use MemoryManagerModule, only: mem_allocate use MemoryManagerExtModule, only: mem_set_value @@ -567,8 +599,6 @@ subroutine allocate_scalars(this, name_model, input_mempath) character(len=*), intent(in) :: name_model character(len=*), intent(in) :: input_mempath logical(LGP) :: found - ! -- local -! ------------------------------------------------------------------------------ ! ! -- Create memory path this%memoryPath = create_mem_path(name_model, 'DIS') @@ -612,24 +642,15 @@ subroutine allocate_scalars(this, name_model, input_mempath) ! -- update input filename call mem_set_value(this%input_fname, 'INPUT_FNAME', & this%input_mempath, found) - ! - ! -- Return - return - end subroutine allocate_scalars - - subroutine allocate_arrays(this) -! ****************************************************************************** -! allocate_arrays -- Read discretization information from file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + end subroutine allocate_scalars_default + + !> @brief Allocate and initialize array variables. + subroutine allocate_arrays_default(this) ! -- modules use MemoryManagerModule, only: mem_allocate ! -- dummy class(DisBaseType) :: this integer :: isize -! ------------------------------------------------------------------------------ ! ! -- Allocate call mem_allocate(this%mshape, this%ndim, 'MSHAPE', this%memoryPath) @@ -652,91 +673,16 @@ subroutine allocate_arrays(this) ! -- Allocate the arrays call mem_allocate(this%dbuff, isize, 'DBUFF', this%name_model) call mem_allocate(this%ibuff, isize, 'IBUFF', this%name_model) - ! - ! -- Return - return - end subroutine allocate_arrays - - function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & - flag_string, allow_zero) result(nodeu) -! ****************************************************************************** -! nodeu_from_string -- Receive a string and convert the string to a user -! nodenumber. The model is unstructured; just read user nodenumber. -! If flag_string argument is present and true, the first token in string -! is allowed to be a string (e.g. boundary name). In this case, if a string -! is encountered, return value as -2. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- dummy - class(DisBaseType) :: this - integer(I4B), intent(inout) :: lloc - integer(I4B), intent(inout) :: istart - integer(I4B), intent(inout) :: istop - integer(I4B), intent(in) :: in - integer(I4B), intent(in) :: iout - character(len=*), intent(inout) :: line - logical, optional, intent(in) :: flag_string - logical, optional, intent(in) :: allow_zero - integer(I4B) :: nodeu - ! -- local -! ------------------------------------------------------------------------------ - ! - ! - nodeu = 0 - call store_error('Program error: DisBaseType method nodeu_from_string & - ¬ implemented.', terminate=.TRUE.) - ! - ! -- return - return - end function nodeu_from_string - - function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & - allow_zero) result(nodeu) -! ****************************************************************************** -! nodeu_from_cellid -- Receive cellid as a string and convert the string to a -! user nodenumber. -! If flag_string argument is present and true, the first token in string -! is allowed to be a string (e.g. boundary name). In this case, if a string -! is encountered, return value as -2. -! If allow_zero argument is present and true, if all indices equal zero, the -! result can be zero. If allow_zero is false, a zero in any index causes an -! error. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- dummy - class(DisBaseType) :: this - character(len=*), intent(inout) :: cellid - integer(I4B), intent(in) :: inunit - integer(I4B), intent(in) :: iout - logical, optional, intent(in) :: flag_string - logical, optional, intent(in) :: allow_zero - integer(I4B) :: nodeu -! ------------------------------------------------------------------------------ - ! - nodeu = 0 - call store_error('Program error: DisBaseType method nodeu_from_cellid & - ¬ implemented.', terminate=.TRUE.) - ! - ! -- return - return - end function nodeu_from_cellid + end subroutine allocate_arrays_default + !> @brief Parse a reduced nodenumber from a line of text. + !! + !! The model is unstructured; just read user nodenumber. + !! If flag_string argument is present and true, the first token in string + !! is allowed to be a string (e.g. boundary name). In this case, if a string + !! is encountered, return value as -2. function noder_from_string(this, lloc, istart, istop, in, iout, line, & flag_string) result(noder) -! ****************************************************************************** -! noder_from_string -- Receive a string and convert the string to a reduced -! nodenumber. The model is unstructured; just read user nodenumber. -! If flag_string argument is present and true, the first token in string -! is allowed to be a string (e.g. boundary name). In this case, if a string -! is encountered, return value as -2. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(DisBaseType) :: this integer(I4B), intent(inout) :: lloc @@ -751,7 +697,6 @@ function noder_from_string(this, lloc, istart, istop, in, iout, line, & integer(I4B) :: nodeu character(len=LINELENGTH) :: nodestr logical :: flag_string_local -! ------------------------------------------------------------------------------ ! if (present(flag_string)) then flag_string_local = flag_string @@ -774,27 +719,17 @@ function noder_from_string(this, lloc, istart, istop, in, iout, line, & trim(adjustl(nodestr)) call store_error(errmsg) end if - ! - ! -- return - return end function noder_from_string + !> @brief Parse a reduced nodenumber from a cell ID string. + !! + !! If flag_string argument is present and true, the first token in string + !! is allowed to be a string (e.g. boundary name). In this case, if a string + !! is encountered, return value as -2. + !! If allow_zero argument is present and true, if all indices equal zero, the + !! result can be zero. If allow_zero is false, a zero in any index is an error. function noder_from_cellid(this, cellid, inunit, iout, flag_string, & allow_zero) result(noder) -! ****************************************************************************** -! noder_from_cellid -- Receive cellid as a string and convert it to a reduced -! nodenumber. -! If flag_string argument is present and true, the first token in string -! is allowed to be a string (e.g. boundary name). In this case, if a string -! is encountered, return value as -2. -! If allow_zero argument is present and true, if all indices equal zero, the -! result can be zero. If allow_zero is false, a zero in any index causes an -! error. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- return integer(I4B) :: noder ! -- dummy class(DisBaseType) :: this @@ -808,7 +743,6 @@ function noder_from_cellid(this, cellid, inunit, iout, flag_string, & logical :: allowzerolocal character(len=LINELENGTH) :: nodestr logical :: flag_string_local -! ------------------------------------------------------------------------------ ! if (present(flag_string)) then flag_string_local = flag_string @@ -837,60 +771,10 @@ function noder_from_cellid(this, cellid, inunit, iout, flag_string, & trim(adjustl(nodestr)) call store_error(errmsg) end if - ! - ! -- return - return end function noder_from_cellid - logical function supports_layers(this) -! ****************************************************************************** -! supports_layers -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- dummy - class(DisBaseType) :: this -! ------------------------------------------------------------------------------ - ! - ! - supports_layers = .false. - call store_error('Program error: DisBaseType method supports_layers not & - &implemented.', terminate=.TRUE.) - return - end function supports_layers - - function get_ncpl(this) -! ****************************************************************************** -! get_ncpl -- Return number of cells per layer. This is nodes -! for a DISU grid, as there are no layers. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - ! -- return - integer(I4B) :: get_ncpl - ! -- dummy - class(DisBaseType) :: this -! ------------------------------------------------------------------------------ - ! - ! - get_ncpl = 0 - call store_error('Program error: DisBaseType method get_ncpl not & - &implemented.', terminate=.TRUE.) - ! - ! -- Return - return - end function get_ncpl - + !> @brief Get the volume of cell n based on x value passed. function get_cell_volume(this, n, x) -! ****************************************************************************** -! get_cell_volume -- Return volume of cell n based on x value passed. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules ! -- return real(DP) :: get_cell_volume @@ -903,7 +787,6 @@ function get_cell_volume(this, n, x) real(DP) :: bt real(DP) :: sat real(DP) :: thick -! ------------------------------------------------------------------------------ ! get_cell_volume = DZERO tp = this%top(n) @@ -911,19 +794,11 @@ function get_cell_volume(this, n, x) sat = sQuadraticSaturation(tp, bt, x) thick = (tp - bt) * sat get_cell_volume = this%area(n) * thick - ! - ! -- Return - return end function get_cell_volume + !> @brief Read an integer array from a line of text. subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & iarray, aname) -! ****************************************************************************** -! read_int_array -- Read a GWF integer array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(DisBaseType), intent(inout) :: this character(len=*), intent(inout) :: line @@ -939,19 +814,11 @@ subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & errmsg = 'Programmer error: read_int_array needs to be overridden & &in any DIS type that extends DisBaseType' call store_error(errmsg, terminate=.TRUE.) - ! - ! -- return - return end subroutine read_int_array + !> @brief Read a double precision array from a line of text. subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & darray, aname) -! ****************************************************************************** -! read_dbl_array -- Read a GWF double precision array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(DisBaseType), intent(inout) :: this character(len=*), intent(inout) :: line @@ -967,18 +834,10 @@ subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & errmsg = 'Programmer error: read_dbl_array needs to be overridden & &in any DIS type that extends DisBaseType' call store_error(errmsg, terminate=.TRUE.) - ! - ! -- return - return end subroutine read_dbl_array + !> @brief Fill an integer array indexed by reduced node number. subroutine fill_int_array(this, ibuff1, ibuff2) -! ****************************************************************************** -! fill_dbl_array -- Fill a GWF integer array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(DisBaseType), intent(inout) :: this integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibuff1 @@ -986,24 +845,16 @@ subroutine fill_int_array(this, ibuff1, ibuff2) ! -- local integer(I4B) :: nodeu integer(I4B) :: noder -! ------------------------------------------------------------------------------ + ! do nodeu = 1, this%nodesuser noder = this%get_nodenumber(nodeu, 0) if (noder <= 0) cycle ibuff2(noder) = ibuff1(nodeu) end do - ! - ! -- return - return end subroutine fill_int_array + !> @brief Fill a double precision array indexed by reduced node number. subroutine fill_dbl_array(this, buff1, buff2) -! ****************************************************************************** -! fill_dbl_array -- Fill a GWF double precision array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(DisBaseType), intent(inout) :: this real(DP), dimension(:), pointer, contiguous, intent(in) :: buff1 @@ -1011,31 +862,24 @@ subroutine fill_dbl_array(this, buff1, buff2) ! -- local integer(I4B) :: nodeu integer(I4B) :: noder -! ------------------------------------------------------------------------------ + ! do nodeu = 1, this%nodesuser noder = this%get_nodenumber(nodeu, 0) if (noder <= 0) cycle buff2(noder) = buff1(nodeu) end do - ! - ! -- return - return end subroutine fill_dbl_array + !> @brief Read a list using the list reader object. + !! + !! Convert user node numbers to reduced numbers. + !! Terminate if any nodenumbers are within an inactive domain. + !! Set up time series and multiply by iauxmultcol if it exists. + !! Write the list to iout if iprpak is set. subroutine read_list(this, line_reader, in, iout, iprpak, nlist, & inamedbound, iauxmultcol, nodelist, rlist, auxvar, & auxname, boundname, label, pkgname, tsManager, iscloc, & indxconvertflux) -! ****************************************************************************** -! read_list -- Read a list using the list reader object. -! Convert user node numbers to reduced numbers. -! Terminate if any nodenumbers are within an inactive domain. -! Set up time series and multiply by iauxmultcol if it exists. -! Write the list to iout if iprpak is set. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: LENBOUNDNAME, LINELENGTH use LongLineReaderModule, only: LongLineReaderType @@ -1059,8 +903,6 @@ subroutine read_list(this, line_reader, in, iout, iprpak, nlist, & character(len=LENAUXNAME), dimension(:), intent(inout) :: auxname character(len=LENBOUNDNAME), dimension(:), pointer, contiguous, & intent(inout) :: boundname - !character(len=:), dimension(:), pointer, contiguous, intent(inout) :: auxname - !character(len=:), dimension(:), pointer, contiguous, intent(inout) :: boundname character(len=*), intent(in) :: label character(len=*), intent(in) :: pkgName type(TimeSeriesManagerType) :: tsManager @@ -1075,7 +917,6 @@ subroutine read_list(this, line_reader, in, iout, iprpak, nlist, & type(ListReaderType) :: lstrdobj type(TimeSeriesLinkType), pointer :: tsLinkBnd => null() type(TimeSeriesLinkType), pointer :: tsLinkAux => null() -! ------------------------------------------------------------------------------ ! ! -- Read the list call lstrdobj%read_list(line_reader, in, iout, nlist, inamedbound, & @@ -1176,90 +1017,10 @@ subroutine read_list(this, line_reader, in, iout, iprpak, nlist, & call store_error_unit(in) end if end if - ! - ! -- return end subroutine read_list - subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, & - icolbnd, aname, inunit, iout) -! ****************************************************************************** -! read_layer_array -- Read a 2d double array into col icolbnd of darray. -! For cells that are outside of the active domain, -! do not copy the array value into darray. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - ! -- dummy - class(DisBaseType) :: this - integer(I4B), intent(in) :: ncolbnd - integer(I4B), intent(in) :: maxbnd - integer(I4B), dimension(maxbnd) :: nodelist - real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray - integer(I4B), intent(in) :: icolbnd - character(len=*), intent(in) :: aname - integer(I4B), intent(in) :: inunit - integer(I4B), intent(in) :: iout - ! - ! - errmsg = 'Programmer error: read_layer_array needs to be overridden & - &in any DIS type that extends DisBaseType' - call store_error(errmsg, terminate=.TRUE.) - ! - ! -- return - end subroutine read_layer_array - - subroutine record_array(this, darray, iout, iprint, idataun, aname, & - cdatafmp, nvaluesp, nwidthp, editdesc, dinact) -! ****************************************************************************** -! record_array -- Record a double precision array. The array will be -! printed to an external file and/or written to an unformatted external file -! depending on the argument specifications. -! ****************************************************************************** -! -! SPECIFICATIONS: -! darray is the double precision array to record -! iout is the unit number for ascii output -! iprint is a flag indicating whether or not to print the array -! idataun is the unit number to which the array will be written in binary -! form; if negative then do not write by layers, write entire array -! aname is the text descriptor of the array -! cdatafmp is the fortran format for writing the array -! nvaluesp is the number of values per line for printing -! nwidthp is the width of the number for printing -! editdesc is the format type (I, G, F, S, E) -! dinact is the double precision value to use for cells that are excluded -! from the model domain -! ------------------------------------------------------------------------------ - ! -- dummy - class(DisBaseType), intent(inout) :: this - real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray - integer(I4B), intent(in) :: iout - integer(I4B), intent(in) :: iprint - integer(I4B), intent(in) :: idataun - character(len=*), intent(in) :: aname - character(len=*), intent(in) :: cdatafmp - integer(I4B), intent(in) :: nvaluesp - integer(I4B), intent(in) :: nwidthp - character(len=*), intent(in) :: editdesc - real(DP), intent(in) :: dinact - ! - ! -- - errmsg = 'Programmer error: record_array needs to be overridden & - &in any DIS type that extends DisBaseType' - call store_error(errmsg, terminate=.TRUE.) - ! - end subroutine record_array - + !> @brief Record a connection-based double precision array. subroutine record_connection_array(this, flowja, ibinun, iout) -! ****************************************************************************** -! record_connection_array -- Record a connection-based double precision -! array for either a structured or an unstructured grid. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(DisBaseType) :: this real(DP), dimension(:), intent(in) :: flowja @@ -1269,103 +1030,42 @@ subroutine record_connection_array(this, flowja, ibinun, iout) character(len=16), dimension(1) :: text ! -- data data text(1)/' FLOW-JA-FACE'/ -! ------------------------------------------------------------------------------ ! ! -- write full ja array call ubdsv1(kstp, kper, text(1), ibinun, flowja, size(flowja), 1, 1, & iout, delt, pertim, totim) - ! - ! -- return - return end subroutine record_connection_array + !> @brief Convert reduced node number to string in form (nodenumber) or (k,i,j). subroutine noder_to_string(this, noder, str) -! ****************************************************************************** -! noder_to_string -- Convert reduced node number to a string in the form of -! (nodenumber) or (k,i,j) -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: noder character(len=*), intent(inout) :: str ! -- local integer(I4B) :: nodeu -! ------------------------------------------------------------------------------ ! nodeu = this%get_nodeuser(noder) call this%nodeu_to_string(nodeu, str) - ! - ! -- return - return end subroutine noder_to_string + !> @brief Convert reduced node number to cellid and fill array with (nodenumber) + !! or (k,j) or (k,i,j) subroutine noder_to_array(this, noder, arr) -! ****************************************************************************** -! noder_to_array -- Convert reduced node number to cellid and fill array with -! (nodenumber) or (k,j) or (k,i,j) -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: noder integer(I4B), dimension(:), intent(inout) :: arr ! -- local integer(I4B) :: nodeu -! ------------------------------------------------------------------------------ ! nodeu = this%get_nodeuser(noder) call this%nodeu_to_array(nodeu, arr) - ! - ! -- return - return end subroutine noder_to_array - subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & - dstmodel, dstpackage, naux, auxtxt, & - ibdchn, nlist, iout) -! ****************************************************************************** -! record_srcdst_list_header -- Record list header for imeth=6 -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- dummy - class(DisBaseType) :: this - character(len=16), intent(in) :: text - character(len=16), intent(in) :: textmodel - character(len=16), intent(in) :: textpackage - character(len=16), intent(in) :: dstmodel - character(len=16), intent(in) :: dstpackage - integer(I4B), intent(in) :: naux - character(len=16), dimension(:), intent(in) :: auxtxt - integer(I4B), intent(in) :: ibdchn - integer(I4B), intent(in) :: nlist - integer(I4B), intent(in) :: iout - ! - ! -- - errmsg = 'Programmer error: record_srcdst_list_header needs to be & - &overridden in any DIS type that extends DisBaseType' - call store_error(errmsg, terminate=.TRUE.) - ! - ! -- return - return - end subroutine record_srcdst_list_header - + !> @brief Record list header subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, & naux, aux, olconv, olconv2) -! ****************************************************************************** -! record_srcdst_list_header -- Record list header -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use InputOutputModule, only: ubdsvd ! -- dummy @@ -1383,7 +1083,6 @@ subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, & logical :: lconv2 integer(I4B) :: nodeu integer(I4B) :: nodeu2 -! ------------------------------------------------------------------------------ ! ! -- Use ubdsvb to write list header if (present(olconv)) then @@ -1407,47 +1106,10 @@ subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, & nodeu2 = noder2 end if call ubdsvd(ibdchn, nodeu, nodeu2, q, naux, aux) - ! - ! -- return - return end subroutine record_srcdst_list_entry - subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname) -! ****************************************************************************** -! nlarray_to_nodelist -- Convert an integer array into nodelist. For structured -! model, integer array is layer number; for unstructured -! model, integer array is node number. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use SimModule, only: store_error - use ConstantsModule, only: LINELENGTH - ! -- dummy - class(DisBaseType) :: this - integer(I4B), intent(in) :: maxbnd - integer(I4B), dimension(:), pointer, contiguous :: darray - integer(I4B), dimension(maxbnd), intent(inout) :: nodelist - integer(I4B), intent(inout) :: nbound - character(len=*), intent(in) :: aname - ! - ! -- - errmsg = 'Programmer error: nlarray_to_nodelist needs to be & - &overridden in any DIS type that extends DisBaseType' - call store_error(errmsg, terminate=.TRUE.) - ! - ! -- return - return - end subroutine nlarray_to_nodelist - + !> @brief Find the first highest active cell beneath cell n. subroutine highest_active(this, n, ibound) -! ****************************************************************************** -! highest_active -- Find the first highest active cell beneath cell n -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(DisBaseType) :: this integer(I4B), intent(inout) :: n @@ -1455,7 +1117,6 @@ subroutine highest_active(this, n, ibound) ! -- locals integer(I4B) :: m, ii, iis logical done, bottomcell -! ------------------------------------------------------------------------------ ! ! -- Loop through connected cells until the highest active one (including a ! constant head cell) is found. Return that cell as n. @@ -1483,44 +1144,28 @@ subroutine highest_active(this, n, ibound) end do cloop if (bottomcell) done = .true. end do - ! - ! -- return - return end subroutine highest_active + !> @brief Get the cell area for this node. function get_area(this, node) result(area) -! ****************************************************************************** -! get_area -- Return the cell area for this node -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- return real(DP) :: area ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: node -! ------------------------------------------------------------------------------ ! ! -- Return the cell area area = this%area(node) - ! - ! -- return - return end function get_area - !> @ brief Calculate the area factor for the cell connection - !! - !! Function calculates the area factor for the cell connection. The sum of - !! all area factors for all cell connections to overlying or underlying - !! cells cells will be 1. + !> @brief Calculate the area factor for the cell connection !! - !! TODO: confirm that this works for cells that are only partially covered - !! by overlying or underlying cells. + !! Function calculates the area factor for the cell connection. The sum of + !! all area factors for all cell connections to overlying or underlying + !! cells cells will be 1. !! - !< + !! TODO: confirm that this works for cells that are only partially covered + !! by overlying or underlying cells. function get_area_factor(this, node, idx_conn) result(area_factor) - ! -- return real(DP) :: area_factor !< connection cell area factor ! -- dummy class(DisBaseType) :: this @@ -1529,7 +1174,6 @@ function get_area_factor(this, node, idx_conn) result(area_factor) ! -- local real(DP) :: area_node real(DP) :: area_conn - ! ------------------------------------------------------------------------------ ! ! -- calculate the cell area fraction area_node = this%area(node) @@ -1537,9 +1181,6 @@ function get_area_factor(this, node, idx_conn) result(area_factor) ! ! -- return the cell area factor area_factor = area_conn / area_node - ! - ! -- return - return end function get_area_factor end module BaseDisModule