From beba891bbdecbc6cd09499e4093d1021f0ba3237 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Fri, 17 Nov 2023 11:01:45 -0500 Subject: [PATCH] chore(discretization): cleanup docstrings in DIS/DISV/DISU --- src/Model/GroundWaterFlow/gwf3dis8.f90 | 470 +++------------ src/Model/GroundWaterFlow/gwf3disu8.f90 | 492 +++------------ src/Model/GroundWaterFlow/gwf3disv8.f90 | 568 ++++-------------- .../ModelUtilities/DiscretizationBase.f90 | 87 +-- 4 files changed, 334 insertions(+), 1283 deletions(-) diff --git a/src/Model/GroundWaterFlow/gwf3dis8.f90 b/src/Model/GroundWaterFlow/gwf3dis8.f90 index e661aa2ab7c..2f4772464d3 100644 --- a/src/Model/GroundWaterFlow/gwf3dis8.f90 +++ b/src/Model/GroundWaterFlow/gwf3dis8.f90 @@ -2,15 +2,18 @@ module GwfDisModule use ArrayReadersModule, only: ReadArray use KindModule, only: DP, I4B, LGP - use ConstantsModule, only: LINELENGTH, DHALF, DZERO, LENMEMPATH, LENVARNAME + use ConstantsModule, only: LINELENGTH, DHALF, DONE, DZERO, & + LENMEMPATH, LENVARNAME use BaseDisModule, only: DisBaseType - use InputOutputModule, only: get_node, URWORD, ulasav, ulaprufw, ubdsv1, & - ubdsv06 + use InputOutputModule, only: get_node, get_ijk, URWORD, ulasav, ulaprufw, & + ubdsv1, ubdsv06, urword, getunit, openfile use SimModule, only: count_errors, store_error, store_error_unit, & store_error_filename - use SimVariablesModule, only: errmsg - use MemoryManagerModule, only: mem_allocate + use SimVariablesModule, only: errmsg, idm_context + use MemoryManagerModule, only: mem_allocate, mem_deallocate + use MemoryManagerExtModule, only: mem_set_value, memorylist_remove use TdisModule, only: kstp, kper, pertim, totim, delt + use GwfDisInputModule, only: GwfDisParamFoundType implicit none private @@ -67,16 +70,8 @@ module GwfDisModule contains + !> @brief Create a new structured discretization object subroutine dis_cr(dis, name_model, input_mempath, inunit, iout) -! ****************************************************************************** -! dis_cr -- Create a new discretization 3d object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use KindModule, only: LGP - use MemoryManagerExtModule, only: mem_set_value ! -- dummy class(DisBaseType), pointer :: dis character(len=*), intent(in) :: name_model @@ -85,10 +80,11 @@ subroutine dis_cr(dis, name_model, input_mempath, inunit, iout) integer(I4B), intent(in) :: iout ! -- locals type(GwfDisType), pointer :: disnew + ! -- formats character(len=*), parameter :: fmtheader = & "(1X, /1X, 'DIS -- STRUCTURED GRID DISCRETIZATION PACKAGE,', & &' VERSION 2 : 3/27/2014 - INPUT READ FROM MEMPATH: ', A, /)" -! ------------------------------------------------------------------------------ + allocate (disnew) dis => disnew call disnew%allocate_scalars(name_model, input_mempath) @@ -103,23 +99,14 @@ subroutine dis_cr(dis, name_model, input_mempath, inunit, iout) write (iout, fmtheader) dis%input_mempath end if end if - ! - ! -- Return - return + end subroutine dis_cr + !> @brief Define the discretization subroutine dis3d_df(this) -! ****************************************************************************** -! dis3d_df -- Define -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisType) :: this - ! -- locals -! ------------------------------------------------------------------------------ + ! -- Transfer the data from the memory manager into this package object if (this%inunit /= 0) then ! @@ -135,26 +122,13 @@ subroutine dis3d_df(this) ! ! -- Final grid initialization call this%grid_finalize() - ! - ! -- Return - return + end subroutine dis3d_df + !> @brief Deallocate variables subroutine dis3d_da(this) -! ****************************************************************************** -! dis3d_da -- Deallocate discretization data -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_deallocate - use MemoryManagerExtModule, only: memorylist_remove - use SimVariablesModule, only: idm_context ! -- dummy class(GwfDisType) :: this - ! -- locals -! ------------------------------------------------------------------------------ ! ! -- Deallocate idm memory call memorylist_remove(this%name_model, 'DIS', idm_context) @@ -177,17 +151,11 @@ subroutine dis3d_da(this) call mem_deallocate(this%top2d) call mem_deallocate(this%bot3d) call mem_deallocate(this%idomain) - ! - ! -- Return - return + end subroutine dis3d_da !> @brief Copy options from IDM into package - !< subroutine source_options(this) - ! -- modules - use MemoryManagerExtModule, only: mem_set_value - use GwfDisInputModule, only: GwfDisParamFoundType ! -- dummy class(GwfDisType) :: this ! -- locals @@ -207,15 +175,11 @@ subroutine source_options(this) if (this%iout > 0) then call this%log_options(found) end if - ! - ! -- Return - return + end subroutine source_options !> @brief Write user options to list file - !< subroutine log_options(this, found) - use GwfDisInputModule, only: GwfDisParamFoundType class(GwfDisType) :: this type(GwfDisParamFoundType), intent(in) :: found @@ -248,10 +212,7 @@ subroutine log_options(this, found) end subroutine log_options !> @brief Copy dimensions from IDM into package - !< subroutine source_dimensions(this) - use MemoryManagerExtModule, only: mem_set_value - use GwfDisInputModule, only: GwfDisParamFoundType ! -- dummy class(GwfDisType) :: this ! -- locals @@ -307,15 +268,11 @@ subroutine source_dimensions(this) end do end do end do - ! - ! -- Return - return + end subroutine source_dimensions !> @brief Write dimensions to list file - !< subroutine log_dimensions(this, found) - use GwfDisInputModule, only: GwfDisParamFoundType class(GwfDisType) :: this type(GwfDisParamFoundType), intent(in) :: found @@ -337,22 +294,10 @@ subroutine log_dimensions(this, found) end subroutine log_dimensions + !> @brief Copy grid data from IDM into package subroutine source_griddata(this) -! ****************************************************************************** -! source_griddata -- update simulation mempath griddata -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerExtModule, only: mem_set_value - use GwfDisInputModule, only: GwfDisParamFoundType - ! -- dummy class(GwfDisType) :: this - ! -- locals type(GwfDisParamFoundType) :: found - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%delr, 'DELR', this%input_mempath, found%delr) @@ -365,15 +310,11 @@ subroutine source_griddata(this) if (this%iout > 0) then call this%log_griddata(found) end if - ! - ! -- Return - return + end subroutine source_griddata !> @brief Write dimensions to list file - !< subroutine log_griddata(this, found) - use GwfDisInputModule, only: GwfDisParamFoundType class(GwfDisType) :: this type(GwfDisParamFoundType), intent(in) :: found @@ -403,15 +344,9 @@ subroutine log_griddata(this, found) end subroutine log_griddata + !> @brief Finalize grid (check properties, allocate arrays, compute connections) subroutine grid_finalize(this) -! ****************************************************************************** -! grid_finalize -- Finalize grid -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules - use ConstantsModule, only: LINELENGTH, DZERO use MemoryManagerModule, only: mem_allocate ! -- dummy class(GwfDisType) :: this @@ -430,7 +365,6 @@ subroutine grid_finalize(this) "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',& &/1x, 'Number of user nodes: ',I0,& &/1X, 'Number of nodes in solution: ', I0, //)" -! ------------------------------------------------------------------------------ ! ! -- count active cells this%nodes = 0 @@ -570,22 +504,13 @@ subroutine grid_finalize(this) this%nodeuser) this%nja = this%con%nja this%njas = this%con%njas - ! - ! -- Return - return + end subroutine grid_finalize + !> @brief Write a binary grid file subroutine write_grb(this, icelltype) -! ****************************************************************************** -! write_grb -- Write the binary grid file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules - use InputOutputModule, only: getunit, openfile use OpenSpecModule, only: access, form - use ConstantsModule, only: DZERO ! -- dummy class(GwfDisType) :: this integer(I4B), dimension(:), intent(in) :: icelltype @@ -598,7 +523,6 @@ subroutine write_grb(this, icelltype) character(len=*), parameter :: fmtgrdsave = & "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', & &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)" -! ------------------------------------------------------------------------------ ! ! -- Initialize ntxt = 16 @@ -695,28 +619,17 @@ subroutine write_grb(this, icelltype) ! ! -- Close the file close (iunit) - ! - ! -- return - return + end subroutine write_grb + !> @brief Convert a user nodenumber to a string (nodenumber) or (k,i,j) 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: -! ------------------------------------------------------------------------------ - use InputOutputModule, only: get_ijk - implicit none class(GwfDisType) :: this integer(I4B), intent(in) :: nodeu character(len=*), intent(inout) :: str ! -- local integer(I4B) :: i, j, k character(len=10) :: kstr, istr, jstr -! ------------------------------------------------------------------------------ ! call get_ijk(nodeu, this%nrow, this%ncol, this%nlay, i, j, k) write (kstr, '(i10)') k @@ -725,28 +638,17 @@ subroutine nodeu_to_string(this, nodeu, str) str = '('//trim(adjustl(kstr))//','// & trim(adjustl(istr))//','// & trim(adjustl(jstr))//')' - ! - ! -- return - return + end subroutine nodeu_to_string + !> @brief Convert a user nodenumber to an array (nodenumber) or (k,i,j) 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: -! ------------------------------------------------------------------------------ - use InputOutputModule, only: get_ijk - implicit none class(GwfDisType) :: this integer(I4B), intent(in) :: nodeu integer(I4B), dimension(:), intent(inout) :: arr ! -- local integer(I4B) :: isize integer(I4B) :: i, j, k -! ------------------------------------------------------------------------------ ! ! -- check the size of arr isize = size(arr) @@ -764,29 +666,17 @@ subroutine nodeu_to_array(this, nodeu, arr) arr(1) = k arr(2) = i arr(3) = j - ! - ! -- return - return + end subroutine nodeu_to_array + !> @brief Get reduced node number from user node number 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. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use ConstantsModule, only: LINELENGTH ! -- return integer(I4B) :: nodenumber ! -- dummy class(GwfDisType), intent(in) :: this integer(I4B), intent(in) :: nodeu integer(I4B), intent(in) :: icheck - ! -- local -! ------------------------------------------------------------------------------ ! ! -- check the node number if requested if (icheck /= 0) then @@ -806,23 +696,11 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) nodenumber = nodeu if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu) end if - ! - ! -- return - return + end function get_nodenumber_idx1 - function get_nodenumber_idx3(this, k, i, j, icheck) & - result(nodenumber) -! ****************************************************************************** -! get_nodenumber_idx3 -- Return a nodenumber from the user specified layer, row, -! and column with an option to perform a check. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - use ConstantsModule, only: LINELENGTH - use InputOutputModule, only: get_node - implicit none + !> @brief Get reduced node number from layer, row and column indices + function get_nodenumber_idx3(this, k, i, j, icheck) result(nodenumber) ! -- return integer(I4B) :: nodenumber ! -- dummy @@ -835,7 +713,6 @@ function get_nodenumber_idx3(this, k, i, j, icheck) & character(len=*), parameter :: fmterr = & "('Error in structured-grid cell indices: layer = ',i0,', & &row = ',i0,', column = ',i0)" -! ------------------------------------------------------------------------------ ! nodeu = get_node(k, i, j, this%nlay, this%nrow, this%ncol) if (nodeu < 1) then @@ -862,24 +739,15 @@ function get_nodenumber_idx3(this, k, i, j, icheck) & call store_error(errmsg) end if end if - ! - ! -- return - return + end function get_nodenumber_idx3 + !> @brief Allocate and initialize scalar variables subroutine allocate_scalars(this, name_model, input_mempath) -! ****************************************************************************** -! allocate_scalars -- Allocate and initialize scalars -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisType) :: this character(len=*), intent(in) :: name_model character(len=*), intent(in) :: input_mempath -! ------------------------------------------------------------------------------ ! ! -- Allocate parent scalars call this%DisBaseType%allocate_scalars(name_model, input_mempath) @@ -894,23 +762,13 @@ subroutine allocate_scalars(this, name_model, input_mempath) this%nrow = 0 this%ncol = 0 this%ndim = 3 - ! - ! -- Return - return + end subroutine allocate_scalars + !> @brief Allocate and initialize arrays subroutine allocate_arrays(this) -! ****************************************************************************** -! allocate_arrays -- Allocate arrays -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_allocate ! -- dummy class(GwfDisType) :: this -! ------------------------------------------------------------------------------ ! ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) call this%DisBaseType%allocate_arrays() @@ -929,23 +787,17 @@ subroutine allocate_arrays(this) this%mshape(1) = this%nlay this%mshape(2) = this%nrow this%mshape(3) = this%ncol - ! - ! -- Return - return + end subroutine allocate_arrays + !> @brief Convert a string to a user nodenumber. + !! + !! Parse layer, row and column and return user nodenumber. + !! If flag_string is present and true, the first token may be + !! non-numeric (e.g. boundary name). In this case, return -2. + !< 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 discretization is DIS; read layer, row, and column. -! 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(GwfDisType) :: this integer(I4B), intent(inout) :: lloc @@ -961,7 +813,6 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & integer(I4B) :: k, i, j, nlay, nrow, ncol integer(I4B) :: lloclocal, ndum, istat, n real(DP) :: r -! ------------------------------------------------------------------------------ ! if (present(flag_string)) then if (flag_string) then @@ -1025,27 +876,19 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & call store_error(errmsg) call store_error_unit(in) end if - ! - ! -- return - return end function nodeu_from_string + !> @brief Convert a cellid string to a user nodenumber. + !! + !! If flag_string is present and true, the first token may be + !! non-numeric (e.g. boundary name). In this case, return -2. + !! + !! If allow_zero is present and true, and all indices are zero, the + !! result can be zero. If allow_zero is false, a zero in any index is an error. + !< 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: -! ------------------------------------------------------------------------------ ! -- return integer(I4B) :: nodeu ! -- dummy @@ -1060,7 +903,6 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & integer(I4B) :: k, i, j, nlay, nrow, ncol integer(I4B) :: istat real(DP) :: r -! ------------------------------------------------------------------------------ ! if (present(flag_string)) then if (flag_string) then @@ -1125,60 +967,31 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & call store_error(errmsg) call store_error_unit(inunit) end if - ! - ! -- return - return + end function nodeu_from_cellid + !> @brief Indicates whether the grid discretization supports layers. logical function supports_layers(this) - implicit none - ! -- dummy class(GwfDisType) :: this - ! supports_layers = .true. - return end function supports_layers + !> @brief Return number of cells per layer (nrow * ncol). function get_ncpl(this) -! ****************************************************************************** -! get_ncpl -- Return number of cells per layer. This is nrow * ncol -! for a DIS3D grid. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - ! -- return integer(I4B) :: get_ncpl - ! -- dummy class(GwfDisType) :: this -! ------------------------------------------------------------------------------ - ! get_ncpl = this%nrow * this%ncol - ! - ! -- Return - return end function get_ncpl + !> @brief Get normal vector components between the cell and a given neighbor. + !! The normal points outward from the shared face between noden and nodem. 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. Connection normal is a normal vector pointing -! outward from the shared face between noden and nodem. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use ConstantsModule, only: DZERO, DONE - use InputOutputModule, only: get_ijk ! -- dummy class(GwfDisType) :: this - integer(I4B), intent(in) :: noden - integer(I4B), intent(in) :: nodem - integer(I4B), intent(in) :: ihc + integer(I4B), intent(in) :: noden !< cell (reduced nn) + integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) + integer(I4B), intent(in) :: ihc !< horizontal connection flag real(DP), intent(inout) :: xcomp real(DP), intent(inout) :: ycomp real(DP), intent(inout) :: zcomp @@ -1186,7 +999,6 @@ subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, & ! -- local integer(I4B) :: nodeu1, i1, j1, k1 integer(I4B) :: nodeu2, i2, j2, k2 -! ------------------------------------------------------------------------------ ! ! -- Set vector components based on ihc if (ihc == 0) then @@ -1220,34 +1032,24 @@ subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, & end if ! end if - ! - ! -- return - return + end subroutine connection_normal + !> @brief Get unit vector components between the cell and a given neighbor. + !! Saturation must be provided to compute cell center vertical coordinates. + !! Also return the straight-line connection length. 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. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules - use ConstantsModule, only: DZERO, DONE, DHALF use DisvGeom, only: line_unit_vector - use InputOutputModule, only: get_ijk ! -- dummy class(GwfDisType) :: this - integer(I4B), intent(in) :: noden - integer(I4B), intent(in) :: nodem + integer(I4B), intent(in) :: noden !< cell (reduced nn) + integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) logical, intent(in) :: nozee real(DP), intent(in) :: satn real(DP), intent(in) :: satm - integer(I4B), intent(in) :: ihc + integer(I4B), intent(in) :: ihc !< horizontal connection flag real(DP), intent(inout) :: xcomp real(DP), intent(inout) :: ycomp real(DP), intent(inout) :: zcomp @@ -1258,7 +1060,6 @@ subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & real(DP) :: ds integer(I4B) :: i1, i2, j1, j2, k1, k2 integer(I4B) :: nodeu1, nodeu2, ipos -! ------------------------------------------------------------------------------ ! ! -- Set vector components based on ihc if (ihc == 0) then @@ -1304,24 +1105,19 @@ subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & end if call line_unit_vector(x1, y1, z1, x2, y2, z2, xcomp, ycomp, zcomp, conlen) end if - ! - ! -- return - return + end subroutine - !> @brief return discretization type + !> @brief Get the discretization type subroutine get_dis_type(this, dis_type) class(GwfDisType), intent(in) :: this character(len=*), intent(out) :: dis_type - dis_type = "DIS" end subroutine get_dis_type !> @brief Get a 2D array of polygon vertices, listed in !! clockwise order beginning with the lower left corner subroutine get_polyverts(this, ic, polyverts, closed) - ! -- modules - use InputOutputModule, only: get_ijk ! -- dummy class(GwfDisType), intent(inout) :: this integer(I4B), intent(in) :: ic !< cell number (reduced) @@ -1366,17 +1162,9 @@ subroutine get_polyverts(this, ic, polyverts, closed) end subroutine + !> @brief Read an integer array subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & iarray, aname) -! ****************************************************************************** -! read_int_array -- Read a GWF integer array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use InputOutputModule, only: urword - use ConstantsModule, only: LINELENGTH ! -- dummy class(GwfDisType), intent(inout) :: this character(len=*), intent(inout) :: line @@ -1395,7 +1183,6 @@ subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & integer(I4B) :: ncol integer(I4B) :: nval integer(I4B), dimension(:), pointer, contiguous :: itemp -! ------------------------------------------------------------------------------ ! ! -- Point the temporary pointer array, which is passed to the reading ! subroutine. The temporary array will point to ibuff if it is a @@ -1430,22 +1217,12 @@ subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & if (this%nodes < this%nodesuser) then call this%fill_grid_array(itemp, iarray) end if - ! - ! -- return - return + end subroutine read_int_array + !> @brief Read a double precision array subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & darray, aname) -! ****************************************************************************** -! read_dbl_array -- Read a GWF double precision array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use InputOutputModule, only: urword - use ConstantsModule, only: LINELENGTH ! -- dummy class(GwfDisType), intent(inout) :: this character(len=*), intent(inout) :: line @@ -1464,7 +1241,6 @@ subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & integer(I4B) :: ncol integer(I4B) :: nval real(DP), dimension(:), pointer, contiguous :: dtemp -! ------------------------------------------------------------------------------ ! ! -- Point the temporary pointer array, which is passed to the reading ! subroutine. The temporary array will point to dbuff if it is a @@ -1499,23 +1275,16 @@ subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & if (this%nodes < this%nodesuser) then call this%fill_grid_array(dtemp, darray) end if - ! - ! -- return - return + end subroutine read_dbl_array + !> @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) -! ****************************************************************************** -! 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 - use InputOutputModule, only: get_node ! -- dummy class(GwfDisType) :: this integer(I4B), intent(in) :: maxbnd @@ -1528,7 +1297,6 @@ subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, & integer(I4B), intent(in) :: iout ! -- local integer(I4B) :: ir, ic, ncol, nrow, nlay, nval, ipos, nodeu -! ------------------------------------------------------------------------------ ! ! -- set variables nlay = this%mshape(1) @@ -1553,45 +1321,28 @@ subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, & ! end do end do - ! - ! -- return + end subroutine read_layer_array + !> @brief Record a double precision array. + !! + !! The array is written to a formatted or unformatted external file + !! depending on the arguments. + !< 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 -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisType), 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 + real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record + integer(I4B), intent(in) :: iout !< ascii output unit number + integer(I4B), intent(in) :: iprint !< whether to print the array + integer(I4B), intent(in) :: idataun !< binary output unit number, if negative don't write by layers, write entire array + character(len=*), intent(in) :: aname !< text descriptor + character(len=*), intent(in) :: cdatafmp !< write format + integer(I4B), intent(in) :: nvaluesp !< values per line + integer(I4B), intent(in) :: nwidthp !< number width + character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E) + real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain ! -- local integer(I4B) :: k, ifirst integer(I4B) :: nlay @@ -1605,7 +1356,6 @@ subroutine record_array(this, darray, iout, iprint, idataun, aname, & character(len=*), parameter :: fmthsv = & "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, & &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)" -! ------------------------------------------------------------------------------ ! ! -- set variables nlay = this%mshape(1) @@ -1662,20 +1412,13 @@ subroutine record_array(this, darray, iout, iprint, idataun, aname, & call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, & iout, delt, pertim, totim) end if - ! - ! -- return - return + end subroutine record_array + !> @brief Record list header for imeth=6 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(GwfDisType) :: this character(len=16), intent(in) :: text @@ -1690,7 +1433,6 @@ subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & integer(I4B), intent(in) :: iout ! -- local integer(I4B) :: nlay, nrow, ncol -! ------------------------------------------------------------------------------ ! nlay = this%mshape(1) nrow = this%mshape(2) @@ -1700,23 +1442,11 @@ subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, & ibdchn, naux, auxtxt, ncol, nrow, nlay, & nlist, iout, delt, pertim, totim) - ! - ! -- return - return + end subroutine record_srcdst_list_header + !> @brief Convert an integer array (layer numbers) to nodelist 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 InputOutputModule, only: get_node - use ConstantsModule, only: LINELENGTH ! -- dummy class(GwfDisType) :: this integer(I4B), intent(in) :: maxbnd @@ -1726,7 +1456,6 @@ subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname) character(len=*), intent(in) :: aname ! -- local integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr -! ------------------------------------------------------------------------------ ! ! -- set variables nlay = this%mshape(1) @@ -1788,8 +1517,7 @@ subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname) nbound = maxbnd ! end if - ! - ! -- return + end subroutine nlarray_to_nodelist end module GwfDisModule diff --git a/src/Model/GroundWaterFlow/gwf3disu8.f90 b/src/Model/GroundWaterFlow/gwf3disu8.f90 index d33bf55eb2a..52eee31e738 100644 --- a/src/Model/GroundWaterFlow/gwf3disu8.f90 +++ b/src/Model/GroundWaterFlow/gwf3disu8.f90 @@ -3,15 +3,20 @@ module GwfDisuModule use ArrayReadersModule, only: ReadArray use KindModule, only: DP, I4B, LGP use ConstantsModule, only: LINELENGTH, LENMEMPATH, LENVARNAME, & - DZERO, DONE + DZERO, DONE, DHALF use ConnectionsModule, only: iac_to_ia - use InputOutputModule, only: URWORD, ulasav, ulaprufw, ubdsv1, ubdsv06 + use InputOutputModule, only: URWORD, ulasav, ulaprufw, ubdsv1, ubdsv06, & + getunit, openfile use SimModule, only: count_errors, store_error, store_error_unit, & store_error_filename - use SimVariablesModule, only: errmsg + use SimVariablesModule, only: errmsg, idm_context use BaseDisModule, only: DisBaseType - use MemoryManagerModule, only: mem_allocate + use MemoryManagerModule, only: mem_allocate, mem_deallocate, & + mem_reallocate, mem_setptr + use MemoryManagerExtModule, only: mem_set_value, memorylist_remove use TdisModule, only: kstp, kper, pertim, totim, delt + use GwfDisuInputModule, only: GwfDisuParamFoundType + use DisvGeom, only: line_unit_vector implicit none @@ -82,16 +87,8 @@ module GwfDisuModule contains + !> @brief Create a new unstructured discretization object subroutine disu_cr(dis, name_model, input_mempath, inunit, iout) -! ****************************************************************************** -! disu_cr -- Create discretization object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use KindModule, only: LGP - use MemoryManagerExtModule, only: mem_set_value ! -- dummy class(DisBaseType), pointer :: dis character(len=*), intent(in) :: name_model @@ -103,7 +100,6 @@ subroutine disu_cr(dis, name_model, input_mempath, inunit, iout) character(len=*), parameter :: fmtheader = & "(1X, /1X, 'DISU -- UNSTRUCTURED GRID DISCRETIZATION PACKAGE,', & &' VERSION 2 : 3/27/2014 - INPUT READ FROM MEMPATH: ', A, //)" -! ------------------------------------------------------------------------------ ! ! -- Create a new discretization object allocate (disnew) @@ -125,22 +121,13 @@ subroutine disu_cr(dis, name_model, input_mempath, inunit, iout) ! -- load disu call disnew%disu_load() end if - ! - ! -- Return - return + end subroutine disu_cr + !> @brief Transfer IDM data into this discretization object subroutine disu_load(this) -! ****************************************************************************** -! disu_load -- transfer data into this discretization object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisuType) :: this -! ------------------------------------------------------------------------------ ! ! -- source input data call this%source_options() @@ -160,39 +147,17 @@ subroutine disu_load(this) ! -- Make some final disu checks on the non-reduced user-provided ! input call this%disu_ck() - ! - ! -- Return - return + end subroutine disu_load + !> @brief Define the discretization subroutine disu_df(this) -! ****************************************************************************** -! disu_df -- Read discretization information from DISU input file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- dummy class(GwfDisuType) :: this -! ------------------------------------------------------------------------------ - ! - ! -- Finalize the grid by creating the connection object and reducing the - ! grid using IDOMAIN, if necessary call this%grid_finalize() - ! - ! -- Return - return end subroutine disu_df + !> @brief Finalize the grid subroutine grid_finalize(this) -! ****************************************************************************** -! grid_finalize -- Finalize grid -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_allocate, mem_reallocate ! -- dummy class(GwfDisuType) :: this ! -- locals @@ -208,7 +173,6 @@ subroutine grid_finalize(this) "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',& &/1x, 'Number of user nodes: ',I0,& &/1X, 'Number of nodes in solution: ', I0, //)" -! ------------------------------------------------------------------------------ ! ! -- count active cells this%nodes = 0 @@ -298,19 +262,11 @@ subroutine grid_finalize(this) this%iangledegx) this%nja = this%con%nja this%njas = this%con%njas - ! - ! -- Return - return + end subroutine grid_finalize + !> @brief Check discretization info subroutine disu_ck(this) -! ****************************************************************************** -! disu_ck -- Check the discretization information -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisuType) :: this ! -- local @@ -333,7 +289,6 @@ subroutine disu_ck(this) "('Top elevation (', 1pg15.6, ') for cell ', i0, ' is above bottom & &elevation (', 1pg15.6, ') for cell ', i0, '. Based on node numbering & &rules cell ', i0, ' must be below cell ', i0, '.')" -! ------------------------------------------------------------------------------ ! ! -- Check connectivity do n = 1, this%nodesuser @@ -426,25 +381,13 @@ subroutine disu_ck(this) call store_error_filename(this%input_fname) end if end if - ! - ! -- Return - return + end subroutine disu_ck + !> @brief Deallocate variables subroutine disu_da(this) -! ****************************************************************************** -! disu_da -- Deallocate discretization object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_deallocate - use MemoryManagerExtModule, only: memorylist_remove - use SimVariablesModule, only: idm_context ! -- dummy class(GwfDisuType) :: this -! ------------------------------------------------------------------------------ ! ! -- Deallocate idm memory call memorylist_remove(this%name_model, 'DISU', idm_context) @@ -483,50 +426,29 @@ subroutine disu_da(this) ! ! -- DisBaseType deallocate call this%DisBaseType%dis_da() - ! - ! -- Return - return + end subroutine disu_da + !> @brief Convert a user nodenumber to a string (nodenumber) 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(GwfDisuType) :: this integer(I4B), intent(in) :: nodeu character(len=*), intent(inout) :: str ! -- local character(len=10) :: nstr -! ------------------------------------------------------------------------------ ! write (nstr, '(i0)') nodeu str = '('//trim(adjustl(nstr))//')' - ! - ! -- return - return end subroutine nodeu_to_string + !> @brief Convert a user nodenumber to an array (nodenumber) 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: -! ------------------------------------------------------------------------------ - use InputOutputModule, only: get_ijk - implicit none class(GwfDisuType) :: this integer(I4B), intent(in) :: nodeu integer(I4B), dimension(:), intent(inout) :: arr ! -- local integer(I4B) :: isize -! ------------------------------------------------------------------------------ ! ! -- check the size of arr isize = size(arr) @@ -539,15 +461,10 @@ subroutine nodeu_to_array(this, nodeu, arr) ! ! -- fill array arr(1) = nodeu - ! - ! -- return - return end subroutine nodeu_to_array !> @brief Write user options to list file - !< subroutine log_options(this, found) - use GwfDisuInputModule, only: GwfDisuParamFoundType class(GwfDisuType) :: this type(GwfDisuParamFoundType), intent(in) :: found @@ -585,24 +502,13 @@ subroutine log_options(this, found) end subroutine log_options !> @brief Copy options from IDM into package - !< subroutine source_options(this) -! ****************************************************************************** -! source_options -- source options from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerExtModule, only: mem_set_value - use GwfDisuInputModule, only: GwfDisuParamFoundType ! -- dummy class(GwfDisuType) :: this ! -- locals character(len=LENVARNAME), dimension(3) :: lenunits = & &[character(len=LENVARNAME) :: 'FEET', 'METERS', 'CENTIMETERS'] type(GwfDisuParamFoundType) :: found -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%lenuni, 'LENGTH_UNITS', this%input_mempath, & @@ -618,15 +524,11 @@ subroutine source_options(this) if (this%iout > 0) then call this%log_options(found) end if - ! - ! -- Return - return + end subroutine source_options !> @brief Write dimensions to list file - !< subroutine log_dimensions(this, found) - use GwfDisuInputModule, only: GwfDisuParamFoundType class(GwfDisuType) :: this type(GwfDisuParamFoundType), intent(in) :: found @@ -649,22 +551,12 @@ subroutine log_dimensions(this, found) end subroutine log_dimensions !> @brief Copy dimensions from IDM into package - !< subroutine source_dimensions(this) -! ****************************************************************************** -! source_dimensions -- source dimensions from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - use MemoryManagerExtModule, only: mem_set_value - use GwfDisuInputModule, only: GwfDisuParamFoundType ! -- dummy class(GwfDisuType) :: this ! -- locals integer(I4B) :: n type(GwfDisuParamFoundType) :: found -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%nodesuser, 'NODES', this%input_mempath, found%nodes) @@ -715,15 +607,11 @@ subroutine source_dimensions(this) do n = 1, this%nodesuser this%idomain(n) = 1 end do - ! - ! -- Return - return + end subroutine source_dimensions !> @brief Write griddata found to list file - !< subroutine log_griddata(this, found) - use GwfDisuInputModule, only: GwfDisuParamFoundType class(GwfDisuType) :: this type(GwfDisuParamFoundType), intent(in) :: found @@ -749,22 +637,12 @@ subroutine log_griddata(this, found) end subroutine log_griddata + !> @brief Copy grid data from IDM into package subroutine source_griddata(this) -! ****************************************************************************** -! source_griddata -- source griddata from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerExtModule, only: mem_set_value - use GwfDisuInputModule, only: GwfDisuParamFoundType ! -- dummy class(GwfDisuType) :: this ! -- locals type(GwfDisuParamFoundType) :: found - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%top1d, 'TOP', this%input_mempath, found%top) @@ -776,15 +654,11 @@ subroutine source_griddata(this) if (this%iout > 0) then call this%log_griddata(found) end if - ! - ! -- Return - return + end subroutine source_griddata !> @brief Write griddata found to list file - !< subroutine log_connectivity(this, found, iac) - use GwfDisuInputModule, only: GwfDisuParamFoundType class(GwfDisuType) :: this type(GwfDisuParamFoundType), intent(in) :: found integer(I4B), dimension(:), contiguous, pointer, intent(in) :: iac @@ -819,24 +693,14 @@ subroutine log_connectivity(this, found, iac) end subroutine log_connectivity + !> @brief Copy grid connectivity info from IDM into package subroutine source_connectivity(this) -! ****************************************************************************** -! source_connectivity -- source connection data from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_setptr - use MemoryManagerExtModule, only: mem_set_value - use GwfDisuInputModule, only: GwfDisuParamFoundType ! -- dummy class(GwfDisuType) :: this ! -- locals type(GwfDisuParamFoundType) :: found integer(I4B), dimension(:), contiguous, pointer :: iac => null() ! -- formats -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%jainp, 'JA', this%input_mempath, found%ja) @@ -859,20 +723,11 @@ subroutine source_connectivity(this) if (this%iout > 0) then call this%log_connectivity(found, iac) end if - ! - ! -- Return - return + end subroutine source_connectivity + !> @brief Copy grid vertex data from IDM into package subroutine source_vertices(this) -! ****************************************************************************** -! source_vertices -- source vertex data from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_setptr ! -- dummy class(GwfDisuType) :: this ! -- local @@ -880,7 +735,6 @@ subroutine source_vertices(this) real(DP), dimension(:), contiguous, pointer :: vert_x => null() real(DP), dimension(:), contiguous, pointer :: vert_y => null() ! -- formats -! ------------------------------------------------------------------------------ ! ! -- set pointers to memory manager input arrays call mem_setptr(vert_x, 'XV', this%input_mempath) @@ -900,11 +754,9 @@ subroutine source_vertices(this) if (this%iout > 0) then write (this%iout, '(1x,a)') 'Discretization Vertex data loaded' end if - ! - ! -- Return - return end subroutine source_vertices + !> @brief Build data structures to hold cell vertex info subroutine define_cellverts(this, icell2d, ncvert, icvert) ! -- modules use SparseModule, only: sparsematrix @@ -917,7 +769,6 @@ subroutine define_cellverts(this, icell2d, ncvert, icvert) type(sparsematrix) :: vert_spm integer(I4B) :: i, j, ierr integer(I4B) :: icv_idx, startvert, maxnnz = 5 -! ------------------------------------------------------------------------------ ! ! -- initialize sparse matrix call vert_spm%init(this%nodesuser, this%nvert, maxnnz) @@ -942,20 +793,11 @@ subroutine define_cellverts(this, icell2d, ncvert, icvert) call mem_allocate(this%javert, vert_spm%nnz, 'JAVERT', this%memoryPath) call vert_spm%filliaja(this%iavert, this%javert, ierr) call vert_spm%destroy() - ! - ! -- Return - return + end subroutine define_cellverts + !> @brief Copy cell2d data from IDM into package subroutine source_cell2d(this) -! ****************************************************************************** -! source_cell2d -- source cell2d data from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_setptr ! -- dummy class(GwfDisuType) :: this ! -- locals @@ -965,8 +807,6 @@ subroutine source_cell2d(this) real(DP), dimension(:), contiguous, pointer :: cell_x => null() real(DP), dimension(:), contiguous, pointer :: cell_y => null() integer(I4B) :: i - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- set pointers to input path ncvert and icvert call mem_setptr(icell2d, 'ICELL2D', this%input_mempath) @@ -999,20 +839,12 @@ subroutine source_cell2d(this) if (this%iout > 0) then write (this%iout, '(1x,a)') 'Discretization Cell2d data loaded' end if - ! - ! -- Return - return + end subroutine source_cell2d + !> @brief Write a binary grid file subroutine write_grb(this, icelltype) -! ****************************************************************************** -! write_grb -- Write the binary grid file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules - use InputOutputModule, only: getunit, openfile use OpenSpecModule, only: access, form ! -- dummy class(GwfDisuType) :: this @@ -1023,10 +855,10 @@ subroutine write_grb(this, icelltype) character(len=50) :: txthdr character(len=lentxt) :: txt character(len=LINELENGTH) :: fname + ! -- formats character(len=*), parameter :: fmtgrdsave = & "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', & &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)" -! ------------------------------------------------------------------------------ ! ! -- Initialize ntxt = 10 @@ -1127,26 +959,15 @@ subroutine write_grb(this, icelltype) ! ! -- Close the file close (iunit) - ! - ! -- return - return + end subroutine write_grb + !> @brief Get reduced node number from user node number 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: -! ------------------------------------------------------------------------------ class(GwfDisuType), intent(in) :: this integer(I4B), intent(in) :: nodeu integer(I4B), intent(in) :: icheck integer(I4B) :: nodenumber -! ------------------------------------------------------------------------------ ! if (icheck /= 0) then if (nodeu < 1 .or. nodeu > this%nodes) then @@ -1164,34 +985,24 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) else nodenumber = this%nodereduced(nodeu) end if - ! - ! -- return - return + end function get_nodenumber_idx1 + !> @brief Get normal vector components between the cell and a given neighbor. + !! The normal points outward from the shared face between noden and nodem. 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 ! -- dummy class(GwfDisuType) :: this - integer(I4B), intent(in) :: noden - integer(I4B), intent(in) :: nodem - integer(I4B), intent(in) :: ihc + integer(I4B), intent(in) :: noden !< cell (reduced nn) + integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) + integer(I4B), intent(in) :: ihc !< horizontal connection flag real(DP), intent(inout) :: xcomp real(DP), intent(inout) :: ycomp real(DP), intent(inout) :: zcomp integer(I4B), intent(in) :: ipos ! -- local real(DP) :: angle, dmult -! ------------------------------------------------------------------------------ ! ! -- Set vector components based on ihc if (ihc == 0) then @@ -1218,26 +1029,14 @@ subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, & ycomp = sin(angle) * dmult zcomp = DZERO end if - ! - ! -- return - return + end subroutine connection_normal + !> @brief Get unit vector components between the cell and a given neighbor. + !! Saturation must be provided to compute cell center vertical coordinates. + !! Also return the straight-line connection length. 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 ConstantsModule, only: DHALF - use DisvGeom, only: line_unit_vector ! -- dummy class(GwfDisuType) :: this integer(I4B), intent(in) :: noden @@ -1252,7 +1051,6 @@ subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & real(DP), intent(inout) :: conlen ! -- local real(DP) :: xn, xm, yn, ym, zn, zm -! ------------------------------------------------------------------------------ ! ! -- Terminate with error if requesting unit vector components for problems ! without cell data @@ -1291,35 +1089,22 @@ subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & ! -- Use coords to find vector components and connection length call line_unit_vector(xn, yn, zn, xm, ym, zm, xcomp, ycomp, zcomp, & conlen) - ! - ! -- return - return + end subroutine connection_vector - ! return discretization type + !> @brief Get the discretization type subroutine get_dis_type(this, dis_type) class(GwfDisuType), intent(in) :: this character(len=*), intent(out) :: dis_type - dis_type = "DISU" - end subroutine get_dis_type + !> @brief Allocate and initialize scalar variables subroutine allocate_scalars(this, name_model, input_mempath) -! ****************************************************************************** -! allocate_scalars -- Allocate and initialize scalar variables in this class -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_allocate ! -- dummy class(GwfDisuType) :: this character(len=*), intent(in) :: name_model character(len=*), intent(in) :: input_mempath - ! -- local -! ------------------------------------------------------------------------------ ! ! -- Allocate parent scalars call this%DisBaseType%allocate_scalars(name_model, input_mempath) @@ -1337,24 +1122,13 @@ subroutine allocate_scalars(this, name_model, input_mempath) this%voffsettol = DZERO this%iangledegx = 0 this%readFromFile = .false. - ! - ! -- Return - return + end subroutine allocate_scalars + !> @brief Allocate and initialize arrays subroutine allocate_arrays(this) -! ****************************************************************************** -! allocate_arrays -- Read discretization information from file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_allocate ! -- dummy class(GwfDisuType) :: this - ! -- local -! ------------------------------------------------------------------------------ ! ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) call this%DisBaseType%allocate_arrays() @@ -1371,11 +1145,10 @@ subroutine allocate_arrays(this) ! ! -- Initialize this%mshape(1) = this%nodesuser - ! - ! -- Return - return + end subroutine allocate_arrays + !> @brief Allocate arrays in memory manager subroutine allocate_arrays_mem(this) use MemoryManagerModule, only: mem_allocate class(GwfDisuType) :: this @@ -1386,18 +1159,14 @@ subroutine allocate_arrays_mem(this) end subroutine allocate_arrays_mem + !> @brief Convert a string to a user nodenumber. + !! + !! Parse and return user nodenumber. + !! If flag_string is present and true, the first token may be + !! non-numeric (e.g. boundary name). In this case, return -2. + !< 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(GwfDisuType) :: this integer(I4B), intent(inout) :: lloc @@ -1412,7 +1181,6 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & ! -- local integer(I4B) :: lloclocal, ndum, istat, n real(DP) :: r -! ------------------------------------------------------------------------------ ! if (present(flag_string)) then if (flag_string) then @@ -1446,27 +1214,19 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & call store_error(errmsg) call store_error_unit(in) end if - ! - ! -- return - return end function nodeu_from_string + !> @brief Convert a cellid string to a user nodenumber. + !! + !! If flag_string is present and true, the first token may be + !! non-numeric (e.g. boundary name). In this case, return -2. + !! + !! If allow_zero is present and true, and all indices are zero, the + !! result can be zero. If allow_zero is false, a zero in any index is an error. + !< 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: -! ------------------------------------------------------------------------------ ! -- return integer(I4B) :: nodeu ! -- dummy @@ -1480,7 +1240,6 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & integer(I4B) :: lloclocal, istart, istop, ndum, n integer(I4B) :: istat real(DP) :: r -! ------------------------------------------------------------------------------ ! if (present(flag_string)) then if (flag_string) then @@ -1515,50 +1274,25 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & call store_error(errmsg) call store_error_unit(inunit) end if - ! - ! -- return - return + end function nodeu_from_cellid + !> @brief Indicates whether the grid discretization supports layers. logical function supports_layers(this) - implicit none - ! -- dummy class(GwfDisuType) :: this - ! supports_layers = .false. - return end function supports_layers + !> @brief Get number of cells per layer (total nodes since DISU isn't layered) 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(GwfDisuType) :: this -! ------------------------------------------------------------------------------ - ! get_ncpl = this%nodesuser - ! - ! -- Return - return end function get_ncpl + !> @brief Read an integer array subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & iarray, aname) -! ****************************************************************************** -! read_int_array -- Read a GWF integer array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisuType), intent(inout) :: this character(len=*), intent(inout) :: line @@ -1572,7 +1306,6 @@ subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & ! -- local integer(I4B) :: nval integer(I4B), dimension(:), pointer, contiguous :: itemp -! ------------------------------------------------------------------------------ ! ! -- Point the temporary pointer array, which is passed to the reading ! subroutine. The temporary array will point to ibuff if it is a @@ -1594,20 +1327,12 @@ subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & if (this%nodes < this%nodesuser) then call this%fill_grid_array(itemp, iarray) end if - ! - ! -- return - return + end subroutine read_int_array + !> @brief Read a double precision array subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & darray, aname) -! ****************************************************************************** -! read_dbl_array -- Read a GWF double precision array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisuType), intent(inout) :: this character(len=*), intent(inout) :: line @@ -1621,7 +1346,6 @@ subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & ! -- local integer(I4B) :: nval real(DP), dimension(:), pointer, contiguous :: dtemp -! ------------------------------------------------------------------------------ ! ! -- Point the temporary pointer array, which is passed to the reading ! subroutine. The temporary array will point to dbuff if it is a @@ -1642,46 +1366,28 @@ subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & if (this%nodes < this%nodesuser) then call this%fill_grid_array(dtemp, darray) end if - ! - ! -- return - return + end subroutine read_dbl_array + !> @brief Record a double precision array. + !! + !! The array is written to a formatted or unformatted external file + !! depending on the arguments. + !< 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 -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisuType), 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 + real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record + integer(I4B), intent(in) :: iout !< ascii output unit number + integer(I4B), intent(in) :: iprint !< whether to print the array + integer(I4B), intent(in) :: idataun !< binary output unit number + character(len=*), intent(in) :: aname !< text descriptor + character(len=*), intent(in) :: cdatafmp ! write format + integer(I4B), intent(in) :: nvaluesp !< values per line + integer(I4B), intent(in) :: nwidthp !< number width + character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E) + real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain ! -- local integer(I4B) :: k, ifirst integer(I4B) :: nlay @@ -1695,7 +1401,6 @@ subroutine record_array(this, darray, iout, iprint, idataun, aname, & character(len=*), parameter :: fmthsv = & "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, & &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)" -! ------------------------------------------------------------------------------ ! ! -- set variables nlay = 1 @@ -1752,20 +1457,13 @@ subroutine record_array(this, darray, iout, iprint, idataun, aname, & call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, & iout, delt, pertim, totim) end if - ! - ! -- return - return + end subroutine record_array + !> @brief Record list header for imeth=6 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(GwfDisuType) :: this character(len=16), intent(in) :: text @@ -1780,7 +1478,6 @@ subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & integer(I4B), intent(in) :: iout ! -- local integer(I4B) :: nlay, nrow, ncol -! ------------------------------------------------------------------------------ ! nlay = 1 nrow = 1 @@ -1790,13 +1487,10 @@ subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, & ibdchn, naux, auxtxt, ncol, nrow, nlay, & nlist, iout, delt, pertim, totim) - ! - ! -- return - return + end subroutine record_srcdst_list_header !> @brief Cast base to DISU - !< function CastAsDisuType(dis) result(disu) class(*), pointer :: dis !< base pointer to DISU object class(GwfDisuType), pointer :: disu !< the resulting DISU pointer diff --git a/src/Model/GroundWaterFlow/gwf3disv8.f90 b/src/Model/GroundWaterFlow/gwf3disv8.f90 index ac9fd2bcf9a..b0f8ac8feb9 100644 --- a/src/Model/GroundWaterFlow/gwf3disv8.f90 +++ b/src/Model/GroundWaterFlow/gwf3disv8.f90 @@ -2,17 +2,19 @@ module GwfDisvModule use ArrayReadersModule, only: ReadArray use KindModule, only: DP, I4B, LGP - use ConstantsModule, only: LINELENGTH, LENMEMPATH, LENVARNAME, DZERO, DONE, & - DHALF + use ConstantsModule, only: LINELENGTH, LENMEMPATH, LENVARNAME, & + DZERO, DONE, DHALF use BaseDisModule, only: DisBaseType - use InputOutputModule, only: get_node, URWORD, ulasav, ulaprufw, ubdsv1, & - ubdsv06 + use InputOutputModule, only: get_node, get_ijk, get_jk, URWORD, ulasav, & + ulaprufw, ubdsv1, ubdsv06, getunit, openfile use SimModule, only: count_errors, store_error, store_error_unit, & store_error_filename - use SimVariablesModule, only: errmsg - use DisvGeom, only: DisvGeomType - use MemoryManagerModule, only: mem_allocate + use SimVariablesModule, only: errmsg, idm_context + use DisvGeom, only: DisvGeomType, line_unit_vector + use MemoryManagerModule, only: mem_allocate, mem_deallocate, mem_setptr + use MemoryManagerExtModule, only: mem_set_value, memorylist_remove use TdisModule, only: kstp, kper, pertim, totim, delt + use GwfDisvInputModule, only: GwfDisvParamFoundType implicit none private @@ -67,32 +69,29 @@ module GwfDisvModule procedure :: allocate_arrays procedure :: get_cell2d_area ! + ! -- Read a node-sized model array (reduced or not) procedure :: read_int_array procedure :: read_dbl_array - ! + end type GwfDisvType contains + !> @brief Create a new discretization by vertices object subroutine disv_cr(dis, name_model, input_mempath, inunit, iout) -! ****************************************************************************** -! disv_cr -- Create a new discretization by vertices object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - use KindModule, only: LGP - use MemoryManagerExtModule, only: mem_set_value + ! -- dummy class(DisBaseType), pointer :: dis character(len=*), intent(in) :: name_model character(len=*), intent(in) :: input_mempath integer(I4B), intent(in) :: inunit integer(I4B), intent(in) :: iout + ! -- local type(GwfDisvType), pointer :: disnew + ! -- formats character(len=*), parameter :: fmtheader = & "(1X, /1X, 'DISV -- VERTEX GRID DISCRETIZATION PACKAGE,', & &' VERSION 1 : 12/23/2015 - INPUT READ FROM MEMPATH: ', A, //)" -! ------------------------------------------------------------------------------ + allocate (disnew) dis => disnew call disnew%allocate_scalars(name_model, input_mempath) @@ -110,70 +109,31 @@ subroutine disv_cr(dis, name_model, input_mempath, inunit, iout) ! -- load disv call disnew%disv_load() end if - ! - ! -- Return - return + end subroutine disv_cr + !> @brief Transfer IDM data into this discretization object subroutine disv_load(this) -! ****************************************************************************** -! disv_load -- transfer data into this discretization object -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - ! -- dummy class(GwfDisvType) :: this - ! -- locals -! ------------------------------------------------------------------------------ - ! + ! -- source input data call this%source_options() call this%source_dimensions() call this%source_griddata() call this%source_vertices() call this%source_cell2d() - ! - ! -- Return - return end subroutine disv_load + !> @brief Define the discretization subroutine disv_df(this) -! ****************************************************************************** -! disv_df -- Define -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - ! -- dummy class(GwfDisvType) :: this - ! -- locals -! ------------------------------------------------------------------------------ - ! - ! -- Final grid initialization call this%grid_finalize() - ! - ! -- Return - return end subroutine disv_df + !> @brief Deallocate variables subroutine disv_da(this) -! ****************************************************************************** -! disv_da -- Deallocate discretization data -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_deallocate - use MemoryManagerExtModule, only: memorylist_remove - use SimVariablesModule, only: idm_context ! -- dummy class(GwfDisvType) :: this - ! -- locals -! ------------------------------------------------------------------------------ ! ! -- Deallocate idm memory call memorylist_remove(this%name_model, 'DISV', idm_context) @@ -198,30 +158,17 @@ subroutine disv_da(this) call mem_deallocate(this%top1d) call mem_deallocate(this%bot2d) call mem_deallocate(this%idomain) - ! - ! -- Return - return + end subroutine disv_da !> @brief Copy options from IDM into package - !< subroutine source_options(this) -! ****************************************************************************** -! source_options -- source options from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerExtModule, only: mem_set_value - use GwfDisvInputModule, only: GwfDisvParamFoundType ! -- dummy class(GwfDisvType) :: this ! -- locals character(len=LENVARNAME), dimension(3) :: lenunits = & &[character(len=LENVARNAME) :: 'FEET', 'METERS', 'CENTIMETERS'] type(GwfDisvParamFoundType) :: found -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%lenuni, 'LENGTH_UNITS', this%input_mempath, & @@ -235,15 +182,11 @@ subroutine source_options(this) if (this%iout > 0) then call this%log_options(found) end if - ! - ! -- Return - return + end subroutine source_options !> @brief Write user options to list file - !< subroutine log_options(this, found) - use GwfDisvInputModule, only: GwfDisvParamFoundType class(GwfDisvType) :: this type(GwfDisvParamFoundType), intent(in) :: found @@ -276,22 +219,12 @@ subroutine log_options(this, found) end subroutine log_options !> @brief Copy dimensions from IDM into package - !< subroutine source_dimensions(this) -! ****************************************************************************** -! source_dimensions -- source dimensions from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - use MemoryManagerExtModule, only: mem_set_value - use GwfDisvInputModule, only: GwfDisvParamFoundType ! -- dummy class(GwfDisvType) :: this ! -- locals integer(I4B) :: j, k type(GwfDisvParamFoundType) :: found -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%nlay, 'NLAY', this%input_mempath, found%nlay) @@ -340,15 +273,11 @@ subroutine source_dimensions(this) this%idomain(j, k) = 1 end do end do - ! - ! -- Return - return + end subroutine source_dimensions !> @brief Write dimensions to list file - !< subroutine log_dimensions(this, found) - use GwfDisvInputModule, only: GwfDisvParamFoundType class(GwfDisvType) :: this type(GwfDisvParamFoundType), intent(in) :: found @@ -370,22 +299,13 @@ subroutine log_dimensions(this, found) end subroutine log_dimensions + !> @brief Copy grid data from IDM into package subroutine source_griddata(this) -! ****************************************************************************** -! source_griddata -- source griddata from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerExtModule, only: mem_set_value - use GwfDisvInputModule, only: GwfDisvParamFoundType ! -- dummy class(GwfDisvType) :: this ! -- locals type(GwfDisvParamFoundType) :: found ! -- formats -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%top1d, 'TOP', this%input_mempath, found%top) @@ -396,15 +316,11 @@ subroutine source_griddata(this) if (this%iout > 0) then call this%log_griddata(found) end if - ! - ! -- Return - return + end subroutine source_griddata !> @brief Write griddata found to list file - !< subroutine log_griddata(this, found) - use GwfDisvInputModule, only: GwfDisvParamFoundType class(GwfDisvType) :: this type(GwfDisvParamFoundType), intent(in) :: found @@ -426,14 +342,8 @@ subroutine log_griddata(this, found) end subroutine log_griddata + !> @brief Finalize grid (check properties, allocate arrays, compute connections) subroutine grid_finalize(this) -! ****************************************************************************** -! grid_finalize -- Finalize grid -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisvType) :: this ! -- locals @@ -448,8 +358,6 @@ subroutine grid_finalize(this) "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',& &/1x, 'Number of user nodes: ',I0,& &/1X, 'Number of nodes in solution: ', I0, //)" - ! -- data -! ------------------------------------------------------------------------------ ! ! -- count active cells this%nodes = 0 @@ -557,28 +465,17 @@ subroutine grid_finalize(this) ! ! -- Build connections call this%connect() - ! - ! -- Return - return + end subroutine grid_finalize + !> @brief Load grid vertices from IDM into package subroutine source_vertices(this) -! ****************************************************************************** -! source_vertices -- source vertex data from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_setptr ! -- dummy class(GwfDisvType) :: this ! -- local integer(I4B) :: i real(DP), dimension(:), contiguous, pointer :: vert_x => null() real(DP), dimension(:), contiguous, pointer :: vert_y => null() - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- set pointers to memory manager input arrays call mem_setptr(vert_x, 'XV', this%input_mempath) @@ -598,11 +495,10 @@ subroutine source_vertices(this) if (this%iout > 0) then write (this%iout, '(1x,a)') 'Discretization Vertex data loaded' end if - ! - ! -- Return - return + end subroutine source_vertices + !> @brief Build data structures to hold cell vertex info subroutine define_cellverts(this, icell2d, ncvert, icvert) ! -- modules use SparseModule, only: sparsematrix @@ -615,7 +511,6 @@ subroutine define_cellverts(this, icell2d, ncvert, icvert) type(sparsematrix) :: vert_spm integer(I4B) :: i, j, ierr integer(I4B) :: icv_idx, startvert, maxnnz = 5 -! ------------------------------------------------------------------------------ ! ! -- initialize sparse matrix call vert_spm%init(this%ncpl, this%nvert, maxnnz) @@ -640,20 +535,11 @@ subroutine define_cellverts(this, icell2d, ncvert, icvert) call mem_allocate(this%javert, vert_spm%nnz, 'JAVERT', this%memoryPath) call vert_spm%filliaja(this%iavert, this%javert, ierr) call vert_spm%destroy() - ! - ! -- Return - return + end subroutine define_cellverts + !> @brief Copy cell2d data from IDM into package subroutine source_cell2d(this) -! ****************************************************************************** -! source_cell2d -- source cell2d data from memory manager input path -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_setptr ! -- dummy class(GwfDisvType) :: this ! -- locals @@ -663,8 +549,6 @@ subroutine source_cell2d(this) real(DP), dimension(:), contiguous, pointer :: cell_x => null() real(DP), dimension(:), contiguous, pointer :: cell_y => null() integer(I4B) :: i - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- set pointers to input path ncvert and icvert call mem_setptr(icell2d, 'ICELL2D', this%input_mempath) @@ -698,19 +582,11 @@ subroutine source_cell2d(this) if (this%iout > 0) then write (this%iout, '(1x,a)') 'Discretization Cell2d data loaded' end if - ! - ! -- Return - return + end subroutine source_cell2d + !> @brief Build grid connections subroutine connect(this) -! ****************************************************************************** -! connect -- Build grid connections -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisvType) :: this ! -- local @@ -719,7 +595,6 @@ subroutine connect(this) integer(I4B) :: narea_eq_zero integer(I4B) :: narea_lt_zero real(DP) :: area -! ------------------------------------------------------------------------------ ! ! -- Initialize narea_eq_zero = 0 @@ -778,21 +653,12 @@ subroutine connect(this) this%nodereduced, this%nodeuser) this%nja = this%con%nja this%njas = this%con%njas - ! - ! - ! -- return - return + end subroutine connect + !> @brief Write a binary grid file subroutine write_grb(this, icelltype) -! ****************************************************************************** -! write_grb -- Write the binary grid file -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules - use InputOutputModule, only: getunit, openfile use OpenSpecModule, only: access, form ! -- dummy class(GwfDisvType) :: this @@ -803,10 +669,10 @@ subroutine write_grb(this, icelltype) character(len=50) :: txthdr character(len=lentxt) :: txt character(len=LINELENGTH) :: fname + ! -- formats character(len=*), parameter :: fmtgrdsave = & "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', & &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)" -! ------------------------------------------------------------------------------ ! ! -- Initialize ntxt = 20 @@ -920,21 +786,11 @@ subroutine write_grb(this, icelltype) ! ! -- Close the file close (iunit) - ! - ! -- return - return + end subroutine write_grb + !> @brief Convert a user nodenumber to a string (nodenumber) or (k,j) subroutine nodeu_to_string(this, nodeu, str) -! ****************************************************************************** -! nodeu_to_string -- Convert user node number to a string in the form of -! (nodenumber) or (k,j) -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use InputOutputModule, only: get_ijk ! -- dummy class(GwfDisvType) :: this integer(I4B), intent(in) :: nodeu @@ -942,35 +798,23 @@ subroutine nodeu_to_string(this, nodeu, str) ! -- local integer(I4B) :: i, j, k character(len=10) :: kstr, jstr -! ------------------------------------------------------------------------------ ! call get_ijk(nodeu, 1, this%ncpl, this%nlay, i, j, k) write (kstr, '(i10)') k write (jstr, '(i10)') j str = '('//trim(adjustl(kstr))//','// & trim(adjustl(jstr))//')' - ! - ! -- return - return + end subroutine nodeu_to_string + !> @brief Convert a user nodenumber to an array (nodenumber) or (k,j) 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: -! ------------------------------------------------------------------------------ - use InputOutputModule, only: get_ijk - implicit none class(GwfDisvType) :: this integer(I4B), intent(in) :: nodeu integer(I4B), dimension(:), intent(inout) :: arr ! -- local integer(I4B) :: isize integer(I4B) :: i, j, k -! ------------------------------------------------------------------------------ ! ! -- check the size of arr isize = size(arr) @@ -987,20 +831,11 @@ subroutine nodeu_to_array(this, nodeu, arr) ! -- fill array arr(1) = k arr(2) = j - ! - ! -- return - return + end subroutine nodeu_to_array + !> @brief Get reduced node number from user node number 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. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- return integer(I4B) :: nodenumber ! -- dummy @@ -1008,7 +843,6 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) integer(I4B), intent(in) :: nodeu integer(I4B), intent(in) :: icheck ! -- local -! ------------------------------------------------------------------------------ ! ! -- check the node number if requested if (icheck /= 0) then @@ -1028,22 +862,11 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) nodenumber = nodeu if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu) end if - ! - ! -- return - return + end function get_nodenumber_idx1 - function get_nodenumber_idx2(this, k, j, icheck) & - result(nodenumber) -! ****************************************************************************** -! get_nodenumber_idx2 -- Return a nodenumber from the user specified layer and -! column with an option to perform a check. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - use InputOutputModule, only: get_node - implicit none + !> @brief Get reduced node number from layer and within-layer node indices + function get_nodenumber_idx2(this, k, j, icheck) result(nodenumber) ! -- return integer(I4B) :: nodenumber ! -- dummy @@ -1055,7 +878,6 @@ function get_nodenumber_idx2(this, k, j, icheck) & ! formats character(len=*), parameter :: fmterr = & &"('Error in disv grid cell indices: layer = ',i0,', node = ',i0)" -! ------------------------------------------------------------------------------ ! nodeu = get_node(k, 1, j, this%nlay, 1, this%ncpl) if (nodeu < 1) then @@ -1093,36 +915,24 @@ function get_nodenumber_idx2(this, k, j, icheck) & end if ! end if - ! - ! -- return - return + end function get_nodenumber_idx2 + !> @brief Get normal vector components between the cell and a given neighbor. + !! The normal points outward from the shared face between noden and nodem. 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 ! -- dummy class(GwfDisvType) :: this - integer(I4B), intent(in) :: noden - integer(I4B), intent(in) :: nodem - integer(I4B), intent(in) :: ihc + integer(I4B), intent(in) :: noden !< cell (reduced nn) + integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) + integer(I4B), intent(in) :: ihc !< horizontal connection flag real(DP), intent(inout) :: xcomp real(DP), intent(inout) :: ycomp real(DP), intent(inout) :: zcomp integer(I4B), intent(in) :: ipos ! -- local - !integer(I4B) :: ipos - !integer(I4B) :: ncell3d, mcell3d real(DP) :: angle, dmult -! ------------------------------------------------------------------------------ ! ! -- Set vector components based on ihc if (ihc == 0) then @@ -1148,34 +958,22 @@ subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, & ycomp = sin(angle) * dmult zcomp = DZERO end if - ! - ! -- return - return + end subroutine connection_normal + !> @brief Get unit vector components between the cell and a given neighbor. + !! Saturation must be provided to compute cell center vertical coordinates. + !! Also return the straight-line connection length. 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 InputOutputModule, only: get_jk - use DisvGeom, only: line_unit_vector ! -- dummy class(GwfDisvType) :: this - integer(I4B), intent(in) :: noden - integer(I4B), intent(in) :: nodem + integer(I4B), intent(in) :: noden !< cell (reduced nn) + integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) logical, intent(in) :: nozee real(DP), intent(in) :: satn real(DP), intent(in) :: satm - integer(I4B), intent(in) :: ihc + integer(I4B), intent(in) :: ihc !< horizontal connection flag real(DP), intent(inout) :: xcomp real(DP), intent(inout) :: ycomp real(DP), intent(inout) :: zcomp @@ -1183,7 +981,6 @@ subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & ! -- local integer(I4B) :: nodeu, ncell2d, mcell2d, k real(DP) :: xn, xm, yn, ym, zn, zm -! ------------------------------------------------------------------------------ ! ! -- Set vector components based on ihc if (ihc == 0) then @@ -1221,34 +1018,22 @@ subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & call line_unit_vector(xn, yn, zn, xm, ym, zm, xcomp, ycomp, zcomp, & conlen) end if - ! - ! -- return - return + end subroutine connection_vector - ! return discretization type + !> @brief Get the discretization type subroutine get_dis_type(this, dis_type) class(GwfDisvType), intent(in) :: this character(len=*), intent(out) :: dis_type - dis_type = "DISV" - end subroutine get_dis_type + !> @brief Allocate and initialize scalars subroutine allocate_scalars(this, name_model, input_mempath) -! ****************************************************************************** -! allocate_scalars -- Allocate and initialize scalars -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_allocate ! -- dummy class(GwfDisvType) :: this character(len=*), intent(in) :: name_model character(len=*), intent(in) :: input_mempath -! ------------------------------------------------------------------------------ ! ! -- Allocate parent scalars call this%DisBaseType%allocate_scalars(name_model, input_mempath) @@ -1263,23 +1048,13 @@ subroutine allocate_scalars(this, name_model, input_mempath) this%ncpl = 0 this%nvert = 0 this%ndim = 2 - ! - ! -- Return - return + end subroutine allocate_scalars + !> @brief Allocate and initialize arrays subroutine allocate_arrays(this) -! ****************************************************************************** -! allocate_arrays -- Allocate arrays -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - use MemoryManagerModule, only: mem_allocate ! -- dummy class(GwfDisvType) :: this -! ------------------------------------------------------------------------------ ! ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) call this%DisBaseType%allocate_arrays() @@ -1296,22 +1071,17 @@ subroutine allocate_arrays(this) ! -- Initialize this%mshape(1) = this%nlay this%mshape(2) = this%ncpl - ! - ! -- Return - return + end subroutine allocate_arrays + !> @brief Get the signed area of the cell. + !! + !! A negative result means points are in counter-clockwise orientation. + !! Area is computed from the formula: + !! a = 1/2 *[(x1*y2 + x2*y3 + x3*y4 + ... + xn*y1) - + !! (x2*y1 + x3*y2 + x4*y3 + ... + x1*yn)] + !< function get_cell2d_area(this, icell2d) result(area) -! ****************************************************************************** -! get_cell2d_area -- Calculate and return the signed area of the cell. A -! negative area means the points are in counter clockwise orientation. -! a = 1/2 *[(x1*y2 + x2*y3 + x3*y4 + ... + xn*y1) - -! (x2*y1 + x3*y2 + x4*y3 + ... + x1*yn)] -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- module ! -- dummy class(GwfDisvType) :: this integer(I4B), intent(in) :: icell2d @@ -1326,7 +1096,6 @@ function get_cell2d_area(this, icell2d) result(area) real(DP) :: y real(DP) :: x1 real(DP) :: y1 -! ------------------------------------------------------------------------------ ! area = DZERO nvert = this%iavert(icell2d + 1) - this%iavert(icell2d) @@ -1358,23 +1127,17 @@ function get_cell2d_area(this, icell2d) result(area) end do ! area = -DONE * area * DHALF - ! - ! -- return - return + end function get_cell2d_area + !> @brief Convert a string to a user nodenumber. + !! + !! Parse layer and within-layer cell number and return user nodenumber. + !! If flag_string is present and true, the first token may be + !! non-numeric (e.g. boundary name). In this case, return -2. + !< 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 discretization is DISV; read layer and cell number. -! 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(GwfDisvType) :: this integer(I4B), intent(inout) :: lloc @@ -1390,7 +1153,6 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & integer(I4B) :: j, k, nlay, nrow, ncpl integer(I4B) :: lloclocal, ndum, istat, n real(DP) :: r -! ------------------------------------------------------------------------------ ! if (present(flag_string)) then if (flag_string) then @@ -1448,26 +1210,19 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & call store_error(errmsg) call store_error_unit(in) end if - ! - ! -- return - return + end function nodeu_from_string + !> @brief Convert a cellid string to a user nodenumber. + !! + !! If flag_string is present and true, the first token may be + !! non-numeric (e.g. boundary name). In this case, return -2. + !! + !! If allow_zero is present and true, and all indices are zero, the + !! result can be zero. If allow_zero is false, a zero in any index is an error. + !< 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: -! ------------------------------------------------------------------------------ ! -- return integer(I4B) :: nodeu ! -- dummy @@ -1482,7 +1237,6 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & integer(I4B) :: lloclocal, ndum, istat, n integer(I4B) :: istart, istop real(DP) :: r -! ------------------------------------------------------------------------------ ! if (present(flag_string)) then if (flag_string) then @@ -1541,39 +1295,20 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & call store_error(errmsg) call store_error_unit(inunit) end if - ! - ! -- return - return + end function nodeu_from_cellid + !> @brief Indicates whether the grid discretization supports layers. logical function supports_layers(this) - implicit none - ! -- dummy class(GwfDisvType) :: this - ! supports_layers = .true. - return end function supports_layers + !> @brief Get number of cells per layer (ncpl). function get_ncpl(this) -! ****************************************************************************** -! get_ncpl -- Return number of cells per layer. This is ncpl -! for a DISV grid. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules - ! -- return integer(I4B) :: get_ncpl - ! -- dummy class(GwfDisvType) :: this -! ------------------------------------------------------------------------------ - ! get_ncpl = this%ncpl - ! - ! -- Return - return end function get_ncpl !> @brief Get a 2D array of polygon vertices, listed in @@ -1622,15 +1357,9 @@ subroutine get_polyverts(this, ic, polyverts, closed) end subroutine + !> @brief Read an integer array subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & iarray, aname) -! ****************************************************************************** -! read_int_array -- Read a GWF integer array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisvType), intent(inout) :: this character(len=*), intent(inout) :: line @@ -1649,7 +1378,6 @@ subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & integer(I4B) :: ncol integer(I4B) :: nval integer(I4B), dimension(:), pointer, contiguous :: itemp -! ------------------------------------------------------------------------------ ! ! -- Point the temporary pointer array, which is passed to the reading ! subroutine. The temporary array will point to ibuff if it is a @@ -1684,20 +1412,12 @@ subroutine read_int_array(this, line, lloc, istart, istop, iout, in, & if (this%nodes < this%nodesuser) then call this%fill_grid_array(itemp, iarray) end if - ! - ! -- return - return + end subroutine read_int_array + !> @brief Read a double precision array subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & darray, aname) -! ****************************************************************************** -! read_dbl_array -- Read a GWF double precision array -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisvType), intent(inout) :: this character(len=*), intent(inout) :: line @@ -1716,7 +1436,6 @@ subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & integer(I4B) :: ncol integer(I4B) :: nval real(DP), dimension(:), pointer, contiguous :: dtemp -! ------------------------------------------------------------------------------ ! ! -- Point the temporary pointer array, which is passed to the reading ! subroutine. The temporary array will point to dbuff if it is a @@ -1756,18 +1475,13 @@ subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, & return end subroutine read_dbl_array + !> @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) -! ****************************************************************************** -! 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 - use InputOutputModule, only: get_node ! -- dummy class(GwfDisvType) :: this integer(I4B), intent(in) :: ncolbnd @@ -1780,7 +1494,6 @@ subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, & integer(I4B), intent(in) :: iout ! -- local integer(I4B) :: ir, ic, ncol, nrow, nlay, nval, ipos, nodeu -! ------------------------------------------------------------------------------ ! ! -- set variables nlay = this%mshape(1) @@ -1803,45 +1516,28 @@ subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, & ipos = ipos + 1 end do end do - ! - ! -- return + end subroutine read_layer_array + !> @brief Record a double precision array. + !! + !! The array is written to a formatted or unformatted external file + !! depending on the arguments. + !< 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 -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfDisvType), 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 + real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record + integer(I4B), intent(in) :: iout !< ascii output unit number + integer(I4B), intent(in) :: iprint !< whether to print the array + integer(I4B), intent(in) :: idataun !< binary output unit number, if negative don't write by layers, write entire array + character(len=*), intent(in) :: aname !< text descriptor + character(len=*), intent(in) :: cdatafmp !< write format + integer(I4B), intent(in) :: nvaluesp !< values per line + integer(I4B), intent(in) :: nwidthp !< number width + character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E) + real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain ! -- local integer(I4B) :: k, ifirst integer(I4B) :: nlay @@ -1855,7 +1551,6 @@ subroutine record_array(this, darray, iout, iprint, idataun, aname, & character(len=*), parameter :: fmthsv = & "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, & &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)" -! ------------------------------------------------------------------------------ ! ! -- set variables nlay = this%mshape(1) @@ -1912,20 +1607,13 @@ subroutine record_array(this, darray, iout, iprint, idataun, aname, & call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, & iout, delt, pertim, totim) end if - ! - ! -- return - return + end subroutine record_array + !> @brief Record list header for imeth=6 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(GwfDisvType) :: this character(len=16), intent(in) :: text @@ -1940,7 +1628,6 @@ subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & integer(I4B), intent(in) :: iout ! -- local integer(I4B) :: nlay, nrow, ncol -! ------------------------------------------------------------------------------ ! nlay = this%mshape(1) nrow = 1 @@ -1950,22 +1637,11 @@ subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, & ibdchn, naux, auxtxt, ncol, nrow, nlay, & nlist, iout, delt, pertim, totim) - ! - ! -- return - return + end subroutine record_srcdst_list_header + !> @brief Convert an integer array (layer numbers) to nodelist 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 InputOutputModule, only: get_node ! -- dummy class(GwfDisvType) :: this integer(I4B), intent(in) :: maxbnd @@ -1975,7 +1651,6 @@ subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname) character(len=*), intent(in) :: aname ! -- local integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr -! ------------------------------------------------------------------------------ ! ! -- set variables nlay = this%mshape(1) @@ -2021,8 +1696,7 @@ subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname) nodelist(ipos) = 0 end do end if - ! - ! -- return + end subroutine nlarray_to_nodelist end module GwfDisvModule diff --git a/src/Model/ModelUtilities/DiscretizationBase.f90 b/src/Model/ModelUtilities/DiscretizationBase.f90 index 673de9e07b1..927b49c6926 100644 --- a/src/Model/ModelUtilities/DiscretizationBase.f90 +++ b/src/Model/ModelUtilities/DiscretizationBase.f90 @@ -5,12 +5,13 @@ module BaseDisModule DZERO, LENMEMPATH, DPIO180 use SmoothingModule, only: sQuadraticSaturation use ConnectionsModule, only: ConnectionsType - use InputOutputModule, only: URWORD, ubdsv1 + use InputOutputModule, only: URWORD, ubdsv1, ubdsvd use SimVariablesModule, only: errmsg use SimModule, only: count_errors, store_error, & store_error_unit use BlockParserModule, only: BlockParserType use MemoryManagerModule, only: mem_allocate + use MemoryManagerExtModule, only: mem_set_value use MemoryHelperModule, only: create_mem_path use TdisModule, only: kstp, kper, pertim, totim, delt use TimeSeriesManagerModule, only: TimeSeriesManagerType @@ -116,9 +117,7 @@ module BaseDisModule !> @brief Define the discretization subroutine dis_df(this) - ! -- dummy class(DisBaseType) :: this - call store_error('Programmer error: dis_df must be overridden', & terminate=.true.) end subroutine dis_df @@ -189,10 +188,8 @@ end subroutine dis_ar !> @brief Write a binary grid file subroutine write_grb(this, icelltype) - ! -- dummy class(DisBaseType) :: this integer(I4B), dimension(:), intent(in) :: icelltype - ! call store_error('Programmer error: write_grb must be overridden', & terminate=.true.) end subroutine write_grb @@ -240,34 +237,30 @@ end subroutine dis_da !> @brief Convert a user nodenumber to a string (nodenumber), (k,j), or (k,i,j) subroutine nodeu_to_string(this, nodeu, str) - ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: nodeu character(len=*), intent(inout) :: str - ! + call store_error('Programmer error: nodeu_to_string must be overridden', & terminate=.true.) end subroutine nodeu_to_string !> @brief Convert a user nodenumber to an array (nodenumber), (k,j), or (k,i,j) subroutine nodeu_to_array(this, nodeu, arr) - ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: nodeu integer(I4B), dimension(:), intent(inout) :: arr - ! + call store_error('Programmer error: nodeu_to_array must be overridden', & terminate=.true.) end subroutine nodeu_to_array !> @brief Convert a reduced nodenumber to a user node number function get_nodeuser(this, noder) result(nodenumber) - ! -- return - integer(I4B) :: nodenumber - ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: noder - ! + integer(I4B) :: nodenumber + if (this%nodes < this%nodesuser) then nodenumber = this%nodeuser(noder) else @@ -276,14 +269,9 @@ function get_nodeuser(this, noder) result(nodenumber) end function get_nodeuser function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) - ! -- modules - use ConstantsModule, only: LINELENGTH - use SimModule, only: store_error - ! -- dummy class(DisBaseType), intent(in) :: this integer(I4B), intent(in) :: nodeu integer(I4B), intent(in) :: icheck - ! -- local integer(I4B) :: nodenumber nodenumber = 0 @@ -292,8 +280,6 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) end function get_nodenumber_idx1 function get_nodenumber_idx2(this, k, j, icheck) result(nodenumber) - use SimModule, only: store_error - ! -- dummy class(DisBaseType), intent(in) :: this integer(I4B), intent(in) :: k, j integer(I4B), intent(in) :: icheck @@ -305,8 +291,6 @@ function get_nodenumber_idx2(this, k, j, icheck) result(nodenumber) end function get_nodenumber_idx2 function get_nodenumber_idx3(this, k, i, j, icheck) result(nodenumber) - use SimModule, only: store_error - ! -- dummy class(DisBaseType), intent(in) :: this integer(I4B), intent(in) :: k, i, j integer(I4B), intent(in) :: icheck @@ -318,11 +302,9 @@ function get_nodenumber_idx3(this, k, i, j, icheck) result(nodenumber) end function get_nodenumber_idx3 !> @brief Get normal vector components between the cell and a given neighbor. + !! The normal points outward from the shared face between noden and nodem. subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, & ipos) - ! -- modules - use SimModule, only: store_error - ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: noden !< cell (reduced nn) integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) @@ -341,9 +323,6 @@ end subroutine connection_normal !! Also return the straight-line connection length. subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, & xcomp, ycomp, zcomp, conlen) - ! -- modules - use SimModule, only: store_error - ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: noden !< cell (reduced nn) integer(I4B), intent(in) :: nodem !< neighbor (reduced nn) @@ -392,7 +371,6 @@ 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('Programmer error: get_dis_type must be overridden', & terminate=.true.) @@ -400,9 +378,6 @@ end subroutine get_dis_type !> @brief Allocate and initialize scalar variables subroutine allocate_scalars(this, name_model, input_mempath) - ! -- modules - use MemoryManagerModule, only: mem_allocate - use MemoryManagerExtModule, only: mem_set_value ! -- dummy class(DisBaseType) :: this character(len=*), intent(in) :: name_model @@ -455,9 +430,6 @@ end subroutine allocate_scalars !> @brief Allocate and initialize arrays subroutine allocate_arrays(this) - ! -- modules - use MemoryManagerModule, only: mem_allocate - ! -- dummy class(DisBaseType) :: this integer :: isize ! @@ -486,11 +458,10 @@ end subroutine allocate_arrays !> @brief Convert a string to a user nodenumber. !! - !! If 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. - ! + !! If DIS or DISV, read indices. If DISU, read user node number directly. + !! If flag_string is present and true, the first token may be + !! non-numeric (e.g. boundary name). In this case, return -2. + !< function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & flag_string, allow_zero) result(nodeu) ! -- dummy @@ -512,10 +483,10 @@ end function nodeu_from_string !> @brief Convert a cellid 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 + !! If flag_string is present and true, the first token may be + !! non-numeric (e.g. boundary name). In this case, return -2. + !! + !! If allow_zero is present and true, and all indices are zero, the !! result can be zero. If allow_zero is false, a zero in any index is an error. !< function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & @@ -637,9 +608,7 @@ end function noder_from_cellid !> @brief Indicates whether the grid discretization supports layers. logical function supports_layers(this) - ! -- dummy class(DisBaseType) :: this - supports_layers = .false. call store_error('Programmer error: supports_layers must be overridden', & terminate=.true.) @@ -648,11 +617,8 @@ end function supports_layers !> @brief Return number of cells per layer. !! This is nodes for a DISU grid, as there are no layers. function get_ncpl(this) - ! -- return integer(I4B) :: get_ncpl - ! -- dummy class(DisBaseType) :: this - get_ncpl = 0 call store_error('Programmer error: get_ncpl must be overridden', & terminate=.true.) @@ -1004,7 +970,6 @@ end subroutine noder_to_array subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & dstmodel, dstpackage, naux, auxtxt, & ibdchn, nlist, iout) - ! -- dummy class(DisBaseType) :: this character(len=16), intent(in) :: text character(len=16), intent(in) :: textmodel @@ -1024,8 +989,6 @@ end subroutine record_srcdst_list_header !> @brief Record list header subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, & naux, aux, olconv, olconv2) - ! -- modules - use InputOutputModule, only: ubdsvd ! -- dummy class(DisBaseType) :: this integer(I4B), intent(in) :: ibdchn @@ -1066,16 +1029,11 @@ subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, & call ubdsvd(ibdchn, nodeu, nodeu2, q, naux, aux) end subroutine record_srcdst_list_entry - !> @brief Convert an integer array into nodelist + !> @brief Convert an integer array to nodelist. !! - !! For structured model, integer array is layer number. - !! For unstructured model, integer array is node number. + !! For DIS/DISV, the array is layer number, for DISU it's node number. !< subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname) - ! -- 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 @@ -1125,15 +1083,12 @@ subroutine highest_active(this, n, ibound) end do end subroutine highest_active - !> @brief Return the cell area for this node + !> @brief Return the cell area for the given node function get_area(this, node) result(area) - ! -- return - real(DP) :: area - ! -- dummy class(DisBaseType) :: this - integer(I4B), intent(in) :: node - ! - ! -- Return the cell area + integer(I4B), intent(in) :: node !< reduced node number + real(DP) :: area + area = this%area(node) end function get_area