diff --git a/src/Model/GroundWaterFlow/gwf3npf8.f90 b/src/Model/GroundWaterFlow/gwf3npf8.f90 index 7872edcd030..0e564046049 100644 --- a/src/Model/GroundWaterFlow/gwf3npf8.f90 +++ b/src/Model/GroundWaterFlow/gwf3npf8.f90 @@ -103,8 +103,9 @@ module GwfNpfModule integer(I4B), pointer :: kchangeper => null() ! last stress period in which any node K (or K22, or K33) values were changed (0 if unchanged from start of simulation) integer(I4B), pointer :: kchangestp => null() ! last time step in which any node K (or K22, or K33) values were changed (0 if unchanged from start of simulation) integer(I4B), dimension(:), pointer, contiguous :: nodekchange => null() ! grid array of flags indicating for each node whether its K (or K22, or K33) value changed (1) at (kchangeper, kchangestp) or not (0) - ! + contains + procedure :: npf_df procedure :: npf_ac procedure :: npf_mc @@ -144,18 +145,15 @@ module GwfNpfModule procedure, public :: increase_edge_count procedure, public :: set_edge_properties procedure, public :: calcSatThickness + end type contains + !> @brief Create a new NPF object. Pass a inunit value of 0 if npf data will + !! initialized from memory + !< subroutine npf_cr(npfobj, name_model, input_mempath, inunit, iout) -! ****************************************************************************** -! npf_cr -- Create a new NPF object. Pass a inunit value of 0 if npf data will -! initialized from memory -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use KindModule, only: LGP use MemoryManagerExtModule, only: mem_set_value @@ -165,12 +163,10 @@ subroutine npf_cr(npfobj, name_model, input_mempath, inunit, iout) character(len=*), intent(in) :: input_mempath integer(I4B), intent(in) :: inunit integer(I4B), intent(in) :: iout - ! -- locals ! -- formats character(len=*), parameter :: fmtheader = & "(1x, /1x, 'NPF -- NODE PROPERTY FLOW PACKAGE, VERSION 1, 3/30/2015', & &' INPUT READ FROM MEMPATH: ', A, /)" -! ------------------------------------------------------------------------------ ! ! -- Create the object allocate (npfobj) @@ -196,7 +192,7 @@ subroutine npf_cr(npfobj, name_model, input_mempath, inunit, iout) return end subroutine npf_cr - !> @brief define the NPF package instance + !> @brief Define the NPF package instance !! !! This is a hybrid routine: it either reads the options for this package !! from the input file, or the optional argument @param npf_options @@ -204,12 +200,6 @@ end subroutine npf_cr !! xt3d_df is called, when enabled. !< subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options) -! ****************************************************************************** -! npf_df -- Define -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SimModule, only: store_error use Xt3dModule, only: xt3d_cr @@ -220,9 +210,6 @@ subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options) integer(I4B), intent(in) :: ingnc !< ghostnodes enabled? (>0 means yes) integer(I4B), intent(in) :: invsc !< viscosity enabled? (>0 means yes) type(GwfNpfOptionsType), optional, intent(in) :: npf_options !< the optional options, for when not constructing from file - ! -- local - ! -- data -! ------------------------------------------------------------------------------ ! ! -- Set a pointer to dis this%dis => dis @@ -247,7 +234,7 @@ subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options) ! -- allocate arrays call this%allocate_arrays(this%dis%nodes, this%dis%njas) end if - + ! call this%check_options() ! ! -- Save pointer to xt3d object @@ -266,21 +253,15 @@ subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options) return end subroutine npf_df + !> @brief Add connections for extended neighbors to the sparse matrix + !< subroutine npf_ac(this, moffset, sparse) -! ****************************************************************************** -! npf_ac -- Add connections for extended neighbors to the sparse matrix -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SparseModule, only: sparsematrix ! -- dummy class(GwfNpftype) :: this integer(I4B), intent(in) :: moffset type(sparsematrix), intent(inout) :: sparse - ! -- local -! ------------------------------------------------------------------------------ ! ! -- Add extended neighbors (neighbors of neighbors) if (this%ixt3d /= 0) call this%xt3d%xt3d_ac(moffset, sparse) @@ -289,20 +270,13 @@ subroutine npf_ac(this, moffset, sparse) return end subroutine npf_ac + !> @brief Map connections and construct iax, jax, and idxglox + !< subroutine npf_mc(this, moffset, matrix_sln) -! ****************************************************************************** -! npf_mc -- Map connections and construct iax, jax, and idxglox -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfNpftype) :: this integer(I4B), intent(in) :: moffset class(MatrixBaseType), pointer :: matrix_sln - ! -- local -! ------------------------------------------------------------------------------ ! if (this%ixt3d /= 0) call this%xt3d%xt3d_mc(moffset, matrix_sln) ! @@ -310,18 +284,12 @@ subroutine npf_mc(this, moffset, matrix_sln) return end subroutine npf_mc - !> @brief allocate and read this NPF instance + !> @brief Allocate and read this NPF instance !! !! Allocate remaining package arrays, preprocess the input data and !! call *_ar on xt3d, when active. !< subroutine npf_ar(this, ic, vsc, ibound, hnew) -! ****************************************************************************** -! npf_ar -- Allocate and Read -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use MemoryManagerModule, only: mem_reallocate ! -- dummy @@ -332,9 +300,6 @@ subroutine npf_ar(this, ic, vsc, ibound, hnew) real(DP), dimension(:), pointer, contiguous, intent(inout) :: hnew !< pointer to model head array ! -- local integer(I4B) :: n - ! -- formats - ! -- data -! ------------------------------------------------------------------------------ ! ! -- Store pointers to arguments that were passed in this%ic => ic @@ -356,7 +321,6 @@ subroutine npf_ar(this, ic, vsc, ibound, hnew) if (this%invsc /= 0) then this%vsc => vsc end if - ! ! -- allocate arrays to store original user input in case TVK/VSC modify them if (this%invsc > 0) then @@ -398,10 +362,10 @@ end subroutine npf_ar !> @brief Read and prepare method for package !! !! Read and prepare NPF stress period data. - !! !< subroutine npf_rp(this) implicit none + ! -- dummy class(GwfNpfType) :: this ! ! -- TVK @@ -409,27 +373,27 @@ subroutine npf_rp(this) call this%tvk%rp() end if ! + ! -- Return return end subroutine npf_rp + !> @brief Advance + !! + !! Sets hold (head old) to bot whenever a wettable cell is dry + !< subroutine npf_ad(this, nodes, hold, hnew, irestore) -! ****************************************************************************** -! npf_ad -- Advance -! Subroutine (1) Sets hold to bot whenever a wettable cell is dry -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + ! -- modules use TdisModule, only: kper, kstp ! implicit none + ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: nodes real(DP), dimension(nodes), intent(inout) :: hold real(DP), dimension(nodes), intent(inout) :: hnew integer(I4B), intent(in) :: irestore + ! -- local integer(I4B) :: n -! ------------------------------------------------------------------------------ ! ! -- loop through all cells and set hold=bot if wettable cell is dry if (this%irewet > 0) then @@ -482,13 +446,9 @@ subroutine npf_ad(this, nodes, hold, hnew, irestore) return end subroutine npf_ad + !> @brief Routines associated fill coefficients + !< subroutine npf_cf(this, kiter, nodes, hnew) -! ****************************************************************************** -! npf_cf -- Formulate -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwfNpfType) :: this integer(I4B) :: kiter @@ -497,7 +457,6 @@ subroutine npf_cf(this, kiter, nodes, hnew) ! -- local integer(I4B) :: n real(DP) :: satn -! ------------------------------------------------------------------------------ ! ! -- Perform wetting and drying if (this%inewton /= 1) then @@ -520,13 +479,9 @@ subroutine npf_cf(this, kiter, nodes, hnew) return end subroutine npf_cf + !> @brief Formulate coefficients + !< subroutine npf_fc(this, kiter, matrix_sln, idxglo, rhs, hnew) -! ****************************************************************************** -! npf_fc -- Formulate -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: DONE ! -- dummy @@ -541,7 +496,6 @@ subroutine npf_fc(this, kiter, matrix_sln, idxglo, rhs, hnew) integer(I4B) :: isymcon, idiagm real(DP) :: hyn, hym real(DP) :: cond -! ------------------------------------------------------------------------------ ! ! -- Calculate conductance and put into amat ! @@ -634,13 +588,9 @@ subroutine npf_fc(this, kiter, matrix_sln, idxglo, rhs, hnew) return end subroutine npf_fc + !> @brief Fill newton terms + !< subroutine npf_fn(this, kiter, matrix_sln, idxglo, rhs, hnew) -! ****************************************************************************** -! npf_fn -- Fill newton terms -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwfNpfType) :: this integer(I4B) :: kiter @@ -665,10 +615,8 @@ subroutine npf_fn(this, kiter, matrix_sln, idxglo, rhs, hnew) real(DP) :: botup real(DP) :: topdn real(DP) :: botdn -! ------------------------------------------------------------------------------ ! ! -- add newton terms to solution matrix - ! nodes = this%dis%nodes nja = this%dis%con%nja if (this%ixt3d /= 0) then @@ -779,16 +727,12 @@ subroutine npf_fn(this, kiter, matrix_sln, idxglo, rhs, hnew) return end subroutine npf_fn + !> @brief Under-relaxation + !! + !! Under-relaxation of Groundwater Flow Model Heads for current outer + !! iteration using the cell bottoms at the bottom of the model + !< subroutine npf_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax) -! ****************************************************************************** -! bnd_nur -- under-relaxation -! Subroutine: (1) Under-relaxation of Groundwater Flow Model Heads for current -! outer iteration using the cell bottoms at the bottom of the -! model -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: neqmod @@ -803,8 +747,6 @@ subroutine npf_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax) real(DP) :: botm real(DP) :: xx real(DP) :: dxx -! ------------------------------------------------------------------------------ - ! ! -- Newton-Raphson under-relaxation do n = 1, this%dis%nodes @@ -827,17 +769,13 @@ subroutine npf_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax) end if end do ! - ! -- return + ! -- Return return end subroutine npf_nur + !> @brief Calculate flowja + !< subroutine npf_cq(this, hnew, flowja) -! ****************************************************************************** -! npf_cq -- Calculate flowja -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwfNpfType) :: this real(DP), intent(inout), dimension(:) :: hnew @@ -845,7 +783,6 @@ subroutine npf_cq(this, hnew, flowja) ! -- local integer(I4B) :: n, ipos, m real(DP) :: qnm -! ------------------------------------------------------------------------------ ! ! -- Calculate the flow across each cell face and store in flowja ! @@ -869,19 +806,14 @@ subroutine npf_cq(this, hnew, flowja) return end subroutine npf_cq + !> @brief Fractional cell saturation + !< subroutine sgwf_npf_thksat(this, n, hn, thksat) -! ****************************************************************************** -! sgwf_npf_thksat -- Fractional cell saturation -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: n real(DP), intent(in) :: hn real(DP), intent(inout) :: thksat -! ------------------------------------------------------------------------------ ! ! -- Standard Formulation if (hn >= this%dis%top(n)) then @@ -901,13 +833,9 @@ subroutine sgwf_npf_thksat(this, n, hn, thksat) return end subroutine sgwf_npf_thksat + !> @brief Flow between two cells + !< subroutine sgwf_npf_qcalc(this, n, m, hn, hm, icon, qnm) -! ****************************************************************************** -! sgwf_npf_qcalc -- Flow between two cells -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: n @@ -921,7 +849,6 @@ subroutine sgwf_npf_qcalc(this, n, m, hn, hm, icon, qnm) real(DP) :: condnm real(DP) :: hntemp, hmtemp integer(I4B) :: ihc -! ------------------------------------------------------------------------------ ! ! -- Initialize ihc = this%dis%con%ihc(this%dis%con%jas(icon)) @@ -981,13 +908,9 @@ subroutine sgwf_npf_qcalc(this, n, m, hn, hm, icon, qnm) return end subroutine sgwf_npf_qcalc + !> @brief Record flowja and calculate specific discharge if requested + !< subroutine npf_save_model_flows(this, flowja, icbcfl, icbcun) -! ****************************************************************************** -! npf_save_model_flows -- Record flowja and calculate specific discharge if requested -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- dummy class(GwfNpfType) :: this real(DP), dimension(:), intent(in) :: flowja @@ -995,9 +918,6 @@ subroutine npf_save_model_flows(this, flowja, icbcfl, icbcun) integer(I4B), intent(in) :: icbcun ! -- local integer(I4B) :: ibinun - !data - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- Set unit number for binary output if (this%ipakcb < 0) then @@ -1028,13 +948,9 @@ subroutine npf_save_model_flows(this, flowja, icbcfl, icbcun) return end subroutine npf_save_model_flows + !> @brief Print budget + !< subroutine npf_print_model_flows(this, ibudfl, flowja) -! ****************************************************************************** -! npf_ot -- Budget -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use TdisModule, only: kper, kstp use ConstantsModule, only: LENBIGLINE @@ -1050,7 +966,6 @@ subroutine npf_print_model_flows(this, ibudfl, flowja) ! -- formats character(len=*), parameter :: fmtiprflow = & &"(/,4x,'CALCULATED INTERCELL FLOW FOR PERIOD ', i0, ' STEP ', i0)" -! ------------------------------------------------------------------------------ ! ! -- Write flowja to list file if requested if (ibudfl /= 0 .and. this%iprflow > 0) then @@ -1075,19 +990,14 @@ subroutine npf_print_model_flows(this, ibudfl, flowja) return end subroutine npf_print_model_flows + !> @brief Deallocate variables + !< subroutine npf_da(this) -! ****************************************************************************** -! npf_da -- Deallocate variables -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use MemoryManagerExtModule, only: memorylist_remove use SimVariablesModule, only: idm_context ! -- dummy class(GwfNpftype) :: this -! ------------------------------------------------------------------------------ ! ! -- Deallocate input memory call memorylist_remove(this%name_model, 'NPF', idm_context) @@ -1177,14 +1087,12 @@ end subroutine npf_da !! !! Allocate and initialize scalars for the VSC package. The base model !! allocate scalars method is also called. - !! !< subroutine allocate_scalars(this) ! -- modules use MemoryHelperModule, only: create_mem_path ! -- dummy class(GwfNpftype) :: this -! ------------------------------------------------------------------------------ ! ! -- allocate scalars in NumericalPackageType call this%NumericalPackageType%allocate_scalars() @@ -1279,14 +1187,12 @@ end subroutine allocate_scalars !> @ brief Store backup copy of hydraulic conductivity when the VSC !! package is activate !! - !! The K arrays (K11, etc.) get multiplied by the viscosity ratio - !! so that subsequent uses of K already take into account the effect - !! of viscosity. Thus the original user-specified K array values are - !! lost unless they are backed up in k11input, for example. In a new - !! stress period/time step, the values in k11input are multiplied by - !! the viscosity ratio, not k11 since it contains viscosity-adjusted - !! hydraulic conductivity values. - !! + !! The K arrays (K11, etc.) get multiplied by the viscosity ratio so that + !! subsequent uses of K already take into account the effect of viscosity. + !! Thus the original user-specified K array values are lost unless they are + !! backed up in k11input, for example. In a new stress period/time step, + !! the values in k11input are multiplied by the viscosity ratio, not k11 + !! since it contains viscosity-adjusted hydraulic conductivity values. !< subroutine store_original_k_arrays(this, ncells, njas) ! -- modules @@ -1297,7 +1203,6 @@ subroutine store_original_k_arrays(this, ncells, njas) integer(I4B), intent(in) :: njas ! -- local integer(I4B) :: n -! ------------------------------------------------------------------------------ ! ! -- Retain copy of user-specified K arrays do n = 1, ncells @@ -1310,21 +1215,15 @@ subroutine store_original_k_arrays(this, ncells, njas) return end subroutine store_original_k_arrays + !> @brief Allocate npf arrays + !< subroutine allocate_arrays(this, ncells, njas) -! ****************************************************************************** -! allocate_arrays -- Allocate npf arrays -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfNpftype) :: this integer(I4B), intent(in) :: ncells integer(I4B), intent(in) :: njas ! -- local integer(I4B) :: n -! ------------------------------------------------------------------------------ ! call mem_allocate(this%ithickstartflag, ncells, 'ITHICKSTARTFLAG', & this%memoryPath) @@ -1374,17 +1273,13 @@ subroutine allocate_arrays(this, ncells, njas) ' WETDRY', ' ANGLE1', & ' ANGLE2', ' ANGLE3'] ! - ! -- return + ! -- Return return end subroutine allocate_arrays + !> @brief Log npf options sourced from the input mempath + !< subroutine log_options(this, found) -! ****************************************************************************** -! log_options -- log npf options sourced from the input mempath -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use KindModule, only: LGP use GwfNpfInputModule, only: GwfNpfParamFoundType @@ -1392,7 +1287,6 @@ subroutine log_options(this, found) class(GwfNpftype) :: this ! -- locals type(GwfNpfParamFoundType), intent(in) :: found -! ------------------------------------------------------------------------------ ! write (this%iout, '(1x,a)') 'Setting NPF Options' if (found%iprflow) & @@ -1459,16 +1353,14 @@ subroutine log_options(this, found) write (this%iout, '(4x,a,i5)') & 'Head rewet equation (IHDWET) has been set to: ', this%ihdwet write (this%iout, '(1x,a,/)') 'End Setting NPF Options' - + ! + ! -- Return + return end subroutine log_options + !> @brief Update simulation options from input mempath + !< subroutine source_options(this) -! ****************************************************************************** -! source_options -- update simulation options from input mempath -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SimModule, only: store_error, store_error_filename use MemoryManagerModule, only: mem_setptr, get_isize @@ -1483,7 +1375,6 @@ subroutine source_options(this) &[character(len=LENVARNAME) :: 'LOGARITHMIC', 'AMT-LMK', 'AMT-HMK'] type(GwfNpfParamFoundType) :: found character(len=LINELENGTH) :: tvk6_filename -! ------------------------------------------------------------------------------ ! ! -- update defaults with idm sourced values call mem_set_value(this%iprflow, 'IPRFLOW', this%input_mempath, found%iprflow) @@ -1550,10 +1441,13 @@ subroutine source_options(this) return end subroutine source_options + !> @brief Set options in the NPF object + !< subroutine set_options(this, options) + ! -- dummy class(GwfNpftype) :: this type(GwfNpfOptionsType), intent(in) :: options - + ! this%icellavg = options%icellavg this%ithickstrt = options%ithickstrt this%iperched = options%iperched @@ -1563,23 +1457,20 @@ subroutine set_options(this, options) this%wetfct = options%wetfct this%iwetit = options%iwetit this%ihdwet = options%ihdwet - + ! + ! -- Return + return end subroutine set_options + !> @brief Check for conflicting NPF options + !< subroutine check_options(this) -! ****************************************************************************** -! check_options -- Check for conflicting NPF options -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SimModule, only: store_error, count_errors, store_error_filename use ConstantsModule, only: LINELENGTH ! -- dummy class(GwfNpftype) :: this - ! -- local -! ------------------------------------------------------------------------------ + ! ! -- check if this%iusgnrhc has been enabled for a model that is not using ! the Newton-Raphson formulation if (this%iusgnrhc > 0 .and. this%inewton == 0) then @@ -1680,59 +1571,59 @@ end subroutine check_options !> @brief Write dimensions to list file !< subroutine log_griddata(this, found) + ! -- modules use GwfNpfInputModule, only: GwfNpfParamFoundType + ! -- dummy class(GwfNpfType) :: this type(GwfNpfParamFoundType), intent(in) :: found - + ! write (this%iout, '(1x,a)') 'Setting NPF Griddata' - + ! if (found%icelltype) then write (this%iout, '(4x,a)') 'ICELLTYPE set from input file' end if - + ! if (found%k) then write (this%iout, '(4x,a)') 'K set from input file' end if - + ! if (found%k33) then write (this%iout, '(4x,a)') 'K33 set from input file' else write (this%iout, '(4x,a)') 'K33 not provided. Setting K33 = K.' end if - + ! if (found%k22) then write (this%iout, '(4x,a)') 'K22 set from input file' else write (this%iout, '(4x,a)') 'K22 not provided. Setting K22 = K.' end if - + ! if (found%wetdry) then write (this%iout, '(4x,a)') 'WETDRY set from input file' end if - + ! if (found%angle1) then write (this%iout, '(4x,a)') 'ANGLE1 set from input file' end if - + ! if (found%angle2) then write (this%iout, '(4x,a)') 'ANGLE2 set from input file' end if - + ! if (found%angle3) then write (this%iout, '(4x,a)') 'ANGLE3 set from input file' end if - + ! write (this%iout, '(1x,a,/)') 'End Setting NPF Griddata' - + ! + ! -- Return + return end subroutine log_griddata + !> @brief Update simulation griddata from input mempath + !< subroutine source_griddata(this) -! ****************************************************************************** -! source_griddata -- update simulation griddata from input mempath -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SimModule, only: count_errors, store_error use MemoryManagerModule, only: mem_reallocate @@ -1745,8 +1636,6 @@ subroutine source_griddata(this) type(GwfNpfParamFoundType) :: found logical, dimension(2) :: afound integer(I4B), dimension(:), pointer, contiguous :: map - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- set map to convert user input data into reduced data map => null() @@ -1824,13 +1713,10 @@ subroutine source_griddata(this) return end subroutine source_griddata + !> @brief Initialize and check NPF data + !< subroutine prepcheck(this) -! ****************************************************************************** -! prepcheck -- Initialize and check NPF data -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ + ! -- modules use ConstantsModule, only: LINELENGTH, DPIO180 use SimModule, only: store_error, count_errors, store_error_filename ! -- dummy @@ -1844,7 +1730,6 @@ subroutine prepcheck(this) &"(1x, 'Hydraulic property ',a,' is <= 0 for cell ',a, ' ', 1pg15.6)" character(len=*), parameter :: fmtkerr2 = & &"(1x, '... ', i0,' additional errors not shown for ',a)" -! ------------------------------------------------------------------------------ ! ! -- initialize aname => this%aname @@ -1979,7 +1864,8 @@ subroutine prepcheck(this) if (count_errors() > 0) then call store_error_filename(this%input_fname) end if - + ! + ! -- Return return end subroutine prepcheck @@ -1994,10 +1880,12 @@ end subroutine prepcheck !! 5. If NEWTON under-relaxation, determine lower most node !< subroutine preprocess_input(this) + ! -- modules use ConstantsModule, only: LINELENGTH use SimModule, only: store_error, count_errors, store_error_filename + ! -- dummy class(GwfNpfType) :: this !< the instance of the NPF package - ! local + ! -- local integer(I4B) :: n, m, ii, nn real(DP) :: hyn, hym real(DP) :: satn, topn, botn @@ -2005,7 +1893,7 @@ subroutine preprocess_input(this) real(DP) :: minbot, botm logical :: finished character(len=LINELENGTH) :: cellstr, errmsg - ! format strings + ! -- format character(len=*), parameter :: fmtcnv = & "(1X,'CELL ', A, & &' ELIMINATED BECAUSE ALL HYDRAULIC CONDUCTIVITIES TO NODE ARE 0.')" @@ -2189,14 +2077,12 @@ subroutine preprocess_input(this) ! ! -- Return return - end subroutine preprocess_input !> @brief Calculate CONDSAT array entries for the given node !! !! Calculate saturated conductances for all connections of the given node, !! or optionally for the upper portion of the matrix only. - !! !< subroutine calc_condsat(this, node, upperOnly) ! -- dummy variables @@ -2286,6 +2172,7 @@ subroutine calc_condsat(this, node, upperOnly) this%condsat(jj) = csat end do ! + ! -- Return return end subroutine calc_condsat @@ -2300,7 +2187,7 @@ function calc_initial_sat(this, n) result(satn) ! -- dummy variables class(GwfNpfType) :: this integer(I4B), intent(in) :: n - ! -- return + ! -- Return real(DP) :: satn ! satn = DONE @@ -2311,13 +2198,9 @@ function calc_initial_sat(this, n) result(satn) return end function calc_initial_sat + !> @brief Perform wetting and drying + !< subroutine sgwf_npf_wetdry(this, kiter, hnew) -! ****************************************************************************** -! sgwf_npf_wetdry -- Perform wetting and drying -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use TdisModule, only: kstp, kper use SimModule, only: store_error, store_error_filename @@ -2347,7 +2230,7 @@ subroutine sgwf_npf_wetdry(this, kiter, hnew) &"(1X,/1X,'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED')" character(len=*), parameter :: fmtni = & &"(1X,'CELLID=',a,' ITERATION=',I0,' TIME STEP=',I0,' STRESS PERIOD=',I0)" -! ------------------------------------------------------------------------------ + ! ! -- Initialize ncnvrt = 0 ihdcnv = 0 @@ -2420,18 +2303,14 @@ subroutine sgwf_npf_wetdry(this, kiter, hnew) return end subroutine sgwf_npf_wetdry + !> @brief Determine if a cell should rewet + !! + !! This method can be called from any external object that has a head that + !! can be used to rewet the GWF cell node. The ihc value is used to + !! determine if it is a vertical or horizontal connection, which can operate + !! differently depending on user settings. + !< subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet) -! ****************************************************************************** -! rewet_check -- Determine if a cell should rewet. This method can -! be called from any external object that has a head that can be used to -! rewet the GWF cell node. The ihc value is used to determine if it is a -! vertical or horizontal connection, which can operate differently depending -! on user settings. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: kiter @@ -2444,8 +2323,6 @@ subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet) ! -- local integer(I4B) :: itflg real(DP) :: wd, awd, turnon, bbot - ! -- formats -! ------------------------------------------------------------------------------ ! irewet = 0 ! @@ -2494,14 +2371,10 @@ subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet) return end subroutine rewet_check + !> @brief Print wet/dry message + !< subroutine sgwf_npf_wdmsg(this, icode, ncnvrt, nodcnvrt, acnvrt, ihdcnv, & kiter, n) -! ****************************************************************************** -! sgwf_npf_wdmsg -- Print wet/dry message -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use TdisModule, only: kstp, kper ! -- dummy @@ -2520,7 +2393,7 @@ subroutine sgwf_npf_wdmsg(this, icode, ncnvrt, nodcnvrt, acnvrt, ihdcnv, & "(1X,/1X,'CELL CONVERSIONS FOR ITER.=',I0, & &' STEP=',I0,' PERIOD=',I0,' (NODE or LRC)')" character(len=*), parameter :: fmtnode = "(1X,3X,5(A4, A20))" -! ------------------------------------------------------------------------------ + ! ! -- Keep track of cell conversions if (icode > 0) then ncnvrt = ncnvrt + 1 @@ -2546,22 +2419,17 @@ subroutine sgwf_npf_wdmsg(this, icode, ncnvrt, nodcnvrt, acnvrt, ihdcnv, & return end subroutine sgwf_npf_wdmsg + !> @brief Calculate the effective hydraulic conductivity for the n-m connection + !! + !! n is primary node node number + !! m is connected node (not used if vg is provided) + !! ihc is horizontal indicator (0 vertical, 1 horizontal, 2 vertically + !! staggered) + !! ipos_opt is position of connection in ja array + !! vg is the global unit vector that expresses the direction from which to + !! calculate an effective hydraulic conductivity. + !< function hy_eff(this, n, m, ihc, ipos, vg) result(hy) -! ****************************************************************************** -! hy_eff -- Calculate the effective hydraulic conductivity for the n-m -! connection. -! n is primary node node number -! m is connected node (not used if vg is provided) -! ihc is horizontal indicator (0 vertical, 1 horizontal, 2 vertically -! staggered) -! ipos_opt is position of connection in ja array -! vg is the global unit vector that expresses the direction from which to -! calculate an effective hydraulic conductivity. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- return real(DP) :: hy ! -- dummy @@ -2576,8 +2444,6 @@ function hy_eff(this, n, m, ihc, ipos, vg) result(hy) real(DP) :: hy11, hy22, hy33 real(DP) :: ang1, ang2, ang3 real(DP) :: vg1, vg2, vg3 - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- Initialize iipos = 0 @@ -2643,22 +2509,19 @@ function hy_eff(this, n, m, ihc, ipos, vg) result(hy) return end function hy_eff + !> @brief Horizontal conductance between two cells + !! + !! inwtup: if 1, then upstream-weight condsat, otherwise recalculate + !! + !! This function uses a weighted transmissivity in the harmonic mean + !! conductance calculations. This differs from the MODFLOW-NWT and + !! MODFLOW-USG conductance calculations for the Newton-Raphson formulation + !! which use a weighted hydraulic conductivity. + !< function hcond(ibdn, ibdm, ictn, ictm, inewton, inwtup, ihc, icellavg, iusg, & iupw, condsat, hn, hm, satn, satm, hkn, hkm, topn, topm, & botn, botm, cln, clm, fawidth, satomega, satminopt) & result(condnm) -! ****************************************************************************** -! hcond -- Horizontal conductance between two cells -! inwtup: if 1, then upstream-weight condsat, otherwise recalculate -! -! hcond function uses a weighted transmissivity in the harmonic mean -! conductance calculations. This differs from the MODFLOW-NWT and MODFLOW-USG -! conductance calculations for the Newton-Raphson formulation which use a -! weighted hydraulic conductivity. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- return real(DP) :: condnm ! -- dummy @@ -2700,7 +2563,7 @@ function hcond(ibdn, ibdm, ictn, ictm, inewton, inwtup, ihc, icellavg, iusg, & real(DP) :: top, bot real(DP) :: athk real(DP) :: afac -! ------------------------------------------------------------------------------ + ! if (present(satminopt)) then satmin = satminopt else @@ -2748,7 +2611,7 @@ function hcond(ibdn, ibdm, ictn, ictm, inewton, inwtup, ihc, icellavg, iusg, & sn = sQuadraticSaturation(topn, botn, hn, satomega, satmin) sm = sQuadraticSaturation(topm, botm, hm, satomega, satmin) end if - + ! if (hn > hm) then condnm = sn else @@ -2789,7 +2652,7 @@ function hcond(ibdn, ibdm, ictn, ictm, inewton, inwtup, ihc, icellavg, iusg, & thksatn = max(min(tpn, sill_top) - sill_bot, DZERO) thksatm = max(min(tpm, sill_top) - sill_bot, DZERO) end if - + ! athk = DONE if (iusg == 1) then if (ihc == 2) then @@ -2810,15 +2673,11 @@ function hcond(ibdn, ibdm, ictn, ictm, inewton, inwtup, ihc, icellavg, iusg, & return end function hcond + !> @brief Vertical conductance between two cells + !< function vcond(ibdn, ibdm, ictn, ictm, inewton, ivarcv, idewatcv, & condsat, hn, hm, vkn, vkm, satn, satm, topn, topm, botn, & botm, flowarea) result(condnm) -! ****************************************************************************** -! vcond -- Vertical conductance between two cells -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- return real(DP) :: condnm ! -- dummy @@ -2846,7 +2705,6 @@ function vcond(ibdn, ibdm, ictn, ictm, inewton, ivarcv, idewatcv, & real(DP) :: bovk1 real(DP) :: bovk2 real(DP) :: denom -! ------------------------------------------------------------------------------ ! ! -- If either n or m is inactive then conductance is zero if (ibdn == 0 .or. ibdm == 0) then @@ -2896,30 +2754,26 @@ function vcond(ibdn, ibdm, ictn, ictm, inewton, ivarcv, idewatcv, & return end function vcond + !> @brief Calculate the conductance between two cells + !! + !! k1 is hydraulic conductivity for cell 1 (in the direction of cell2) + !! k2 is hydraulic conductivity for cell 2 (in the direction of cell1) + !! thick1 is the saturated thickness for cell 1 + !! thick2 is the saturated thickness for cell 2 + !! cl1 is the distance from the center of cell1 to the shared face with cell2 + !! cl2 is the distance from the center of cell2 to the shared face with cell1 + !! h1 is the head for cell1 + !! h2 is the head for cell2 + !! width is the width perpendicular to flow + !! iavgmeth is the averaging method: + !! 0 is harmonic averaging + !! 1 is logarithmic averaging + !! 2 is arithmetic averaging of sat thickness and logarithmic averaging of + !! hydraulic conductivity + !! 3 is arithmetic averaging of sat thickness and harmonic averaging of + !! hydraulic conductivity + !< function condmean(k1, k2, thick1, thick2, cl1, cl2, width, iavgmeth) -! ****************************************************************************** -! condmean -- Calculate the conductance between two cells -! -! k1 is hydraulic conductivity for cell 1 (in the direction of cell2) -! k2 is hydraulic conductivity for cell 2 (in the direction of cell1) -! thick1 is the saturated thickness for cell 1 -! thick2 is the saturated thickness for cell 2 -! cl1 is the distance from the center of cell1 to the shared face with cell2 -! cl2 is the distance from the center of cell2 to the shared face with cell1 -! h1 is the head for cell1 -! h2 is the head for cell2 -! width is the width perpendicular to flow -! iavgmeth is the averaging method: -! 0 is harmonic averaging -! 1 is logarithmic averaging -! 2 is arithmetic averaging of sat thickness and logarithmic averaging of -! hydraulic conductivity -! 3 is arithmetic averaging of sat thickness and harmonic averaging of -! hydraulic conductivity -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- return real(DP) :: condmean ! -- dummy @@ -2935,7 +2789,6 @@ function condmean(k1, k2, thick1, thick2, cl1, cl2, width, iavgmeth) real(DP) :: t1 real(DP) :: t2 real(DP) :: tmean, kmean, denom -! ------------------------------------------------------------------------------ ! ! -- Initialize t1 = k1 * thick1 @@ -2986,14 +2839,11 @@ function condmean(k1, k2, thick1, thick2, cl1, cl2, width, iavgmeth) return end function condmean + !> @brief Calculate the the logarithmic mean of two double precision numbers + !! + !! Use an approximation if the ratio is near 1 + !< function logmean(d1, d2) -! ****************************************************************************** -! logmean -- Calculate the the logarithmic mean of two double precision -! numbers. Use an approximation if the ratio is near 1. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- return real(DP) :: logmean ! -- dummy @@ -3001,7 +2851,6 @@ function logmean(d1, d2) real(DP), intent(in) :: d2 ! -- local real(DP) :: drat -! ------------------------------------------------------------------------------ ! drat = d2 / d1 if (drat <= DLNLOW .or. drat >= DLNHIGH) then @@ -3014,31 +2863,28 @@ function logmean(d1, d2) return end function logmean + !> @brief Calculate the effective horizontal hydraulic conductivity from an + !! ellipse using a specified direction (unit vector vg1, vg2, vg3) + !! + !! k11 is the hydraulic conductivity of the major ellipse axis + !! k22 is the hydraulic conductivity of first minor axis + !! k33 is the hydraulic conductivity of the second minor axis + !! ang1 is the counter-clockwise rotation (radians) of the ellipse in + !! the (x, y) plane + !! ang2 is the rotation of the conductivity ellipsoid upward or + !! downward from the (x, y) plane + !! ang3 is the rotation of the conductivity ellipsoid about the major + !! axis + !! vg1, vg2, and vg3 are the components of a unit vector in model coordinates + !! in the direction of the connection between cell n and m + !!iavgmeth is the averaging method. If zero, then use harmonic averaging. + !! if one, then use arithmetic averaging. + !< function hyeff_calc(k11, k22, k33, ang1, ang2, ang3, vg1, vg2, vg3, & iavgmeth) result(hyeff) -! ****************************************************************************** -! hyeff_calc -- Calculate the effective horizontal hydraulic conductivity from -! an ellipse using a specified direction (unit vector vg1, vg2, vg3). -! k11 is the hydraulic conductivity of the major ellipse axis -! k22 is the hydraulic conductivity of first minor axis -! k33 is the hydraulic conductivity of the second minor axis -! ang1 is the counter-clockwise rotation (radians) of the ellipse in -! the (x, y) plane -! ang2 is the rotation of the conductivity ellipsoid upward or -! downward from the (x, y) plane -! ang3 is the rotation of the conductivity ellipsoid about the major -! axis -! vg1, vg2, and vg3 are the components of a unit vector in model coordinates -! in the direction of the connection between cell n and m -! iavgmeth is the averaging method. If zero, then use harmonic averaging. -! if one, then use arithmetic averaging. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: DONE - ! -- result + ! -- return real(DP) :: hyeff ! -- dummy real(DP), intent(in) :: k11 @@ -3056,7 +2902,6 @@ function hyeff_calc(k11, k22, k33, ang1, ang2, ang3, vg1, vg2, vg3, & real(DP), dimension(3, 3) :: r real(DP) :: ve1, ve2, ve3 real(DP) :: denom, dnum, d1, d2, d3 -! ------------------------------------------------------------------------------ ! ! -- Sin and cos of angles s1 = sin(ang1) @@ -3120,14 +2965,9 @@ function hyeff_calc(k11, k22, k33, ang1, ang2, ang3, vg1, vg2, vg3, & return end function hyeff_calc + !> @brief Calculate the 3 conmponents of specific discharge at the cell center + !< subroutine calc_spdis(this, flowja) -! ****************************************************************************** -! calc_spdis -- Calculate the 3 conmponents of specific discharge -! at the cell center. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use SimModule, only: store_error ! -- dummy @@ -3177,7 +3017,6 @@ subroutine calc_spdis(this, flowja) real(DP), allocatable, dimension(:) :: bix real(DP), allocatable, dimension(:) :: biy logical :: nozee = .true. -! ------------------------------------------------------------------------------ ! ! -- Ensure dis has necessary information if (this%icalcspdis /= 0 .and. this%dis%con%ianglex == 0) then @@ -3419,18 +3258,13 @@ subroutine calc_spdis(this, flowja) deallocate (bix) deallocate (biy) ! - ! -- return + ! -- Return return end subroutine calc_spdis + !> @brief Save specific discharge in binary format to ibinun + !< subroutine sav_spdis(this, ibinun) -! ****************************************************************************** -! sav_spdis -- save specific discharge in binary format to ibinun -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: ibinun @@ -3439,7 +3273,6 @@ subroutine sav_spdis(this, ibinun) character(len=16), dimension(3) :: auxtxt integer(I4B) :: n integer(I4B) :: naux -! ------------------------------------------------------------------------------ ! ! -- Write the header text = ' DATA-SPDIS' @@ -3456,18 +3289,13 @@ subroutine sav_spdis(this, ibinun) this%spdis(:, n)) end do ! - ! -- return + ! -- Return return end subroutine sav_spdis + !> @brief Save saturation in binary format to ibinun + !< subroutine sav_sat(this, ibinun) -! ****************************************************************************** -! sav_sat -- save saturation in binary format to ibinun -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: ibinun @@ -3477,7 +3305,6 @@ subroutine sav_sat(this, ibinun) real(DP), dimension(1) :: a integer(I4B) :: n integer(I4B) :: naux -! ------------------------------------------------------------------------------ ! ! -- Write the header text = ' DATA-SAT' @@ -3494,41 +3321,30 @@ subroutine sav_sat(this, ibinun) call this%dis%record_mf6_list_entry(ibinun, n, n, DZERO, naux, a) end do ! - ! -- return + ! -- Return return end subroutine sav_sat + !> @brief Reserve space for nedges cells that have an edge on them. + !! + !! This must be called before the npf%allocate_arrays routine, which is + !! called from npf%ar. + !< subroutine increase_edge_count(this, nedges) -! ****************************************************************************** -! increase_edge_count -- reserve space for nedges cells that have an edge on them. -! This must be called before the npf%allocate_arrays routine, which is called -! from npf%ar. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: nedges - ! -- local -! ------------------------------------------------------------------------------ ! this%nedges = this%nedges + nedges ! - ! -- return + ! -- Return return end subroutine increase_edge_count + !> @brief Provide the npf package with edge properties + !< subroutine set_edge_properties(this, nodedge, ihcedge, q, area, nx, ny, & distance) -! ****************************************************************************** -! edge_count -- provide the npf package with edge properties. -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- dummy class(GwfNpfType) :: this integer(I4B), intent(in) :: nodedge @@ -3540,7 +3356,6 @@ subroutine set_edge_properties(this, nodedge, ihcedge, q, area, nx, ny, & real(DP), intent(in) :: distance ! -- local integer(I4B) :: lastedge -! ------------------------------------------------------------------------------ ! this%lastedge = this%lastedge + 1 lastedge = this%lastedge @@ -3556,19 +3371,21 @@ subroutine set_edge_properties(this, nodedge, ihcedge, q, area, nx, ny, & ! edge properties assignment loop, so need to reset lastedge to 0 if (this%lastedge == this%nedges) this%lastedge = 0 ! - ! -- return + ! -- Return return end subroutine set_edge_properties !> Calculate saturated thickness between cell n and m !< function calcSatThickness(this, n, m, ihc) result(satThickness) + ! -- dummy class(GwfNpfType) :: this !< this NPF instance integer(I4B) :: n !< node n integer(I4B) :: m !< node m integer(I4B) :: ihc !< 1 = horizonal connection, 0 for vertical + ! -- return real(DP) :: satThickness !< saturated thickness - + ! satThickness = thksatnm(this%ibound(n), & this%ibound(m), & this%icelltype(n), & @@ -3586,19 +3403,16 @@ function calcSatThickness(this, n, m, ihc) result(satThickness) this%dis%bot(m), & this%satomega, & this%satmin) - + ! + ! -- Return + return end function calcSatThickness + !> @brief Calculate saturated thickness at interface between two cells + !< function thksatnm(ibdn, ibdm, ictn, ictm, inwtup, ihc, iusg, & hn, hm, satn, satm, topn, topm, botn, botm, & satomega, satminopt) result(res) -! ****************************************************************************** -! thksatnm -- calculate saturated thickness at interface between two cells -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ - ! -- modules ! -- return real(DP) :: res ! -- dummy @@ -3629,7 +3443,7 @@ function thksatnm(ibdn, ibdm, ictn, ictm, inwtup, ihc, iusg, & real(DP) :: sill_top, sill_bot real(DP) :: tpn, tpm real(DP) :: top, bot -! ------------------------------------------------------------------------------ + ! if (present(satminopt)) then satmin = satminopt else