diff --git a/src/Model/ParticleTracking/prt-fmi.f90 b/src/Model/ParticleTracking/prt-fmi.f90 index ef84224a0ee..0b4b3c22bb2 100644 --- a/src/Model/ParticleTracking/prt-fmi.f90 +++ b/src/Model/ParticleTracking/prt-fmi.f90 @@ -35,41 +35,41 @@ module PrtFmiModule !> @brief Create a new PrtFmi object subroutine fmi_cr(fmiobj, name_model, inunit, iout) - ! -- dummy + ! dummy type(PrtFmiType), pointer :: fmiobj character(len=*), intent(in) :: name_model integer(I4B), intent(inout) :: inunit integer(I4B), intent(in) :: iout ! - ! -- Create the object + ! Create the object allocate (fmiobj) ! - ! -- create name and memory path + ! create name and memory path call fmiobj%set_names(1, name_model, 'FMI', 'FMI') fmiobj%text = text ! - ! -- Allocate scalars + ! Allocate scalars call fmiobj%allocate_scalars() ! - ! -- Set variables + ! Set variables fmiobj%inunit = inunit fmiobj%iout = iout ! - ! -- Initialize block parser + ! Initialize block parser call fmiobj%parser%Initialize(fmiobj%inunit, fmiobj%iout) ! - ! -- Assign dependent variable label + ! Assign dependent variable label fmiobj%depvartype = 'TRACKS ' end subroutine fmi_cr !> @brief Time step advance subroutine fmi_ad(this) - ! -- modules + ! modules use ConstantsModule, only: DHDRY - ! -- dummy + ! dummy class(PrtFmiType) :: this - ! -- local + ! local integer(I4B) :: n character(len=15) :: nodestr character(len=*), parameter :: fmtdry = & @@ -77,33 +77,33 @@ subroutine fmi_ad(this) character(len=*), parameter :: fmtrewet = & &"(/1X,'DRY CELL REACTIVATED AT ', a)" ! - ! -- Set flag to indicated that flows are being updated. For the case where + ! Set flag to indicated that flows are being updated. For the case where ! flows may be reused (only when flows are read from a file) then set ! the flag to zero to indicated that flows were not updated this%iflowsupdated = 1 ! - ! -- If reading flows from a budget file, read the next set of records + ! If reading flows from a budget file, read the next set of records if (this%iubud /= 0) then call this%advance_bfr() end if ! - ! -- If reading heads from a head file, read the next set of records + ! If reading heads from a head file, read the next set of records if (this%iuhds /= 0) then call this%advance_hfr() end if ! - ! -- If mover flows are being read from file, read the next set of records + ! If mover flows are being read from file, read the next set of records if (this%iumvr /= 0) then call this%mvrbudobj%bfr_advance(this%dis, this%iout) end if ! - ! -- Accumulate flows + ! Accumulate flows call this%accumulate_flows() ! - ! -- if flow cell is dry, then set this%ibound = 0 + ! if flow cell is dry, then set this%ibound = 0 do n = 1, this%dis%nodes ! - ! -- Calculate the ibound-like array that has 0 if saturation + ! Calculate the ibound-like array that has 0 if saturation ! is zero and 1 otherwise if (this%gwfsat(n) > DZERO) then this%ibdgwfsat0(n) = 1 @@ -111,20 +111,20 @@ subroutine fmi_ad(this) this%ibdgwfsat0(n) = 0 end if ! - ! -- Check if active model cell is inactive for flow + ! Check if active model cell is inactive for flow if (this%ibound(n) > 0) then if (this%gwfhead(n) == DHDRY) then - ! -- cell should be made inactive + ! cell should be made inactive this%ibound(n) = 0 call this%dis%noder_to_string(n, nodestr) write (this%iout, fmtdry) trim(nodestr) end if end if ! - ! -- Convert dry model cell to active if flow has rewet + ! Convert dry model cell to active if flow has rewet if (this%ibound(n) == 0) then if (this%gwfhead(n) /= DHDRY) then - ! -- cell is now wet + ! cell is now wet this%ibound(n) = 1 call this%dis%noder_to_string(n, nodestr) write (this%iout, fmtrewet) trim(nodestr) @@ -136,17 +136,17 @@ end subroutine fmi_ad !> @brief Define the flow model interface subroutine prtfmi_df(this, dis, idryinactive) - ! -- modules + ! modules use SimModule, only: store_error - ! -- dummy + ! dummy class(PrtFmiType) :: this class(DisBaseType), pointer, intent(in) :: dis integer(I4B), intent(in) :: idryinactive ! - ! -- Call parent class define + ! Call parent class define call this%FlowModelInterfaceType%fmi_df(dis, idryinactive) ! - ! -- Allocate arrays + ! Allocate arrays allocate (this%StorageFlows(this%dis%nodes)) allocate (this%SourceFlows(this%dis%nodes)) allocate (this%SinkFlows(this%dis%nodes)) @@ -157,9 +157,9 @@ end subroutine prtfmi_df !> @brief Accumulate flows subroutine accumulate_flows(this) implicit none - ! -- dummy + ! dummy class(PrtFmiType) :: this - ! -- local + ! local integer(I4B) :: j, i, ip, ib integer(I4B) :: ioffset, iflowface, iauxiflowface real(DP) :: qbnd diff --git a/src/Model/ParticleTracking/prt-mip.f90 b/src/Model/ParticleTracking/prt-mip.f90 index 6edd72a2740..03decbe2b64 100644 --- a/src/Model/ParticleTracking/prt-mip.f90 +++ b/src/Model/ParticleTracking/prt-mip.f90 @@ -31,35 +31,35 @@ module PrtMipModule !> @brief Create a model input object subroutine mip_cr(mip, name_model, input_mempath, inunit, iout, dis) - ! -- dummy + ! dummy type(PrtMipType), pointer :: mip character(len=*), intent(in) :: name_model character(len=*), intent(in) :: input_mempath integer(I4B), intent(in) :: inunit integer(I4B), intent(in) :: iout class(DisBaseType), pointer, intent(in) :: dis - ! -- formats + ! formats character(len=*), parameter :: fmtheader = & - "(1x, /1x, 'MIP -- MODEL INPUT PACKAGE', & + "(1x, /1x, 'MIP MODEL INPUT PACKAGE', & &' INPUT READ FROM MEMPATH: ', A, /)" ! - ! -- Create the object + ! Create the object allocate (mip) ! - ! -- Create name and memory path + ! Create name and memory path call mip%set_names(1, name_model, 'MIP', 'MIP', input_mempath) ! - ! -- Allocate scalars + ! Allocate scalars call mip%allocate_scalars() ! - ! -- Set variables + ! Set variables mip%inunit = inunit mip%iout = iout ! - ! -- Set pointers + ! Set pointers mip%dis => dis ! - ! -- Print a message identifying the package if enabled + ! Print a message identifying the package if enabled if (inunit > 0) & write (iout, fmtheader) input_mempath @@ -69,13 +69,13 @@ end subroutine mip_cr subroutine mip_da(this) class(PrtMipType) :: this ! - ! -- Deallocate input memory + ! Deallocate input memory call memorystore_remove(this%name_model, 'MIP', idm_context) ! - ! -- Deallocate parent package + ! Deallocate parent package call this%NumericalPackageType%da() ! - ! -- Deallocate arrays + ! Deallocate arrays call mem_deallocate(this%porosity) call mem_deallocate(this%retfactor) call mem_deallocate(this%izone) @@ -91,10 +91,10 @@ end subroutine allocate_scalars subroutine allocate_arrays(this, nodes) class(PrtMipType) :: this integer(I4B), intent(in) :: nodes - ! -- local + ! local integer(I4B) :: i ! - ! -- Allocate + ! Allocate call mem_allocate(this%porosity, nodes, 'POROSITY', this%memoryPath) call mem_allocate(this%retfactor, nodes, 'RETFACTOR', this%memoryPath) call mem_allocate(this%izone, nodes, 'IZONE', this%memoryPath) @@ -109,20 +109,20 @@ end subroutine allocate_arrays !> @ brief Initialize package inputs subroutine mip_ar(this) - ! -- dummy variables + ! dummy variables class(PrtMipType), intent(inout) :: this !< PrtMipType object - ! -- local variables + ! local variables character(len=LINELENGTH) :: errmsg type(PrtMipParamFoundType) :: found integer(I4B), dimension(:), pointer, contiguous :: map => null() ! - ! -- set map to convert user input data into reduced data + ! set map to convert user input data into reduced data if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser ! - ! -- Allocate arrays + ! Allocate arrays call this%allocate_arrays(this%dis%nodes) ! - ! -- Source array inputs from IDM + ! Source array inputs from IDM call mem_set_value(this%porosity, 'POROSITY', this%input_mempath, & map, found%porosity) call mem_set_value(this%retfactor, 'RETFACTOR', this%input_mempath, & @@ -130,7 +130,7 @@ subroutine mip_ar(this) call mem_set_value(this%izone, 'IZONE', this%input_mempath, map, & found%izone) ! - ! -- Ensure POROSITY was found + ! Ensure POROSITY was found if (.not. found%porosity) then write (errmsg, '(a)') 'Error in GRIDDATA block: POROSITY not found' call store_error(errmsg) diff --git a/src/Solution/ParticleTracker/CellRectQuad.f90 b/src/Solution/ParticleTracker/CellRectQuad.f90 index 9c034d2be2d..7e7066c9cd5 100644 --- a/src/Solution/ParticleTracker/CellRectQuad.f90 +++ b/src/Solution/ParticleTracker/CellRectQuad.f90 @@ -75,9 +75,9 @@ end subroutine init_from !! vertices and face flows of a rectangular-quad cell. !< subroutine load_rect_verts_flows(this) - ! -- dummy + ! dummy class(CellRectQuadType), intent(inout) :: this - ! -- local + ! local integer(I4B) :: n, m n = 0 @@ -107,21 +107,21 @@ end subroutine load_rect_verts_flows !! rectangle vertex of a rectangular-quad cell !< function get_rect_ivert_sw(this) result(irv1) - ! -- dummy + ! dummy class(CellRectQuadType), intent(inout) :: this integer(I4B) :: irv1 - ! -- local + ! local integer(I4B) :: irv, irv2, irv4, ipv1, ipv2, ipv4 integer(I4B), dimension(4) :: irvnxt = (/2, 3, 4, 1/) real(DP) :: x1, y1, x2, y2, x4, y4 - ! -- Find the "southwest" rectangle vertex by finding the vertex formed - ! -- either by (1) a rectangle edge over which x decreases (going - ! -- clockwise) followed by an edge over which x does not increase, or by - ! -- (2) a rectangle edge over which y does not decrease (again going - ! -- clockwise) followed by a rectangle edge over which y increases. In - ! -- the end, ipv1 is the index (1, 2, 3, or 4) of the southwest - ! -- rectangle vertex. + ! Find the "southwest" rectangle vertex by finding the vertex formed + ! either by (1) a rectangle edge over which x decreases (going + ! clockwise) followed by an edge over which x does not increase, or by + ! (2) a rectangle edge over which y does not decrease (again going + ! clockwise) followed by a rectangle edge over which y increases. In + ! the end, ipv1 is the index (1, 2, 3, or 4) of the southwest + ! rectangle vertex. do irv = 1, 4 irv4 = irv irv1 = irvnxt(irv4) @@ -152,18 +152,18 @@ end function get_rect_ivert_sw !! as the origin !< subroutine get_rect_dim_rot(this) - ! -- dummy + ! dummy class(CellRectQuadType), intent(inout) :: this - ! -- local + ! local integer(I4B) :: irv2, irv4, ipv1, ipv2, ipv4 integer(I4B), dimension(4) :: irvnxt = (/2, 3, 4, 1/) real(DP) :: x1, y1, x2, y2, x4, y4, dx2, dy2, dx4, dy4 - ! -- Get rectangle vertex neighbors irv2 and irv4 + ! Get rectangle vertex neighbors irv2 and irv4 irv2 = irvnxt(this%irvOrigin) irv4 = irvnxt(irvnxt(irv2)) - ! -- Get model coordinates at irv1, irv2, and irv4 + ! Get model coordinates at irv1, irv2, and irv4 ipv1 = this%irectvert(this%irvOrigin) x1 = this%defn%polyvert(1, ipv1) y1 = this%defn%polyvert(2, ipv1) @@ -174,7 +174,7 @@ subroutine get_rect_dim_rot(this) x4 = this%defn%polyvert(1, ipv4) y4 = this%defn%polyvert(2, ipv4) - ! -- Compute rectangle dimensions + ! Compute rectangle dimensions this%xOrigin = x1 this%yOrigin = y1 this%zOrigin = this%defn%bot @@ -186,8 +186,8 @@ subroutine get_rect_dim_rot(this) this%dy = dsqrt(dx2 * dx2 + dy2 * dy2) this%dz = this%defn%top - this%zOrigin - ! -- Compute sine and cosine of rotation angle (angle between "southern" - ! -- rectangle side irv1-irv4 and the model x axis) + ! Compute sine and cosine of rotation angle (angle between "southern" + ! rectangle side irv1-irv4 and the model x axis) this%sinrot = dy4 / this%dx this%cosrot = dx4 / this%dx end subroutine get_rect_dim_rot @@ -202,7 +202,7 @@ end function get_rect_flow !> @brief Return whether a rectangle face is refined function face_is_refined(this, i) result(is_refined) - ! -- dummy + ! dummy class(CellRectQuadType), intent(inout) :: this integer(I4B) :: i !< face index logical(LGP) :: is_refined diff --git a/src/Solution/ParticleTracker/CellUtil.f90 b/src/Solution/ParticleTracker/CellUtil.f90 index 86934982fce..fdb9f72cfb8 100644 --- a/src/Solution/ParticleTracker/CellUtil.f90 +++ b/src/Solution/ParticleTracker/CellUtil.f90 @@ -15,10 +15,10 @@ subroutine cell_poly_to_rect(poly, rect) use CellRectModule, only: CellRectType, create_cell_rect use CellPolyModule, only: CellPolyType use CellDefnModule, only: CellDefnType - ! -- dummy + ! dummy type(CellPolyType), intent(in), pointer :: poly type(CellRectType), intent(inout), pointer :: rect - ! -- local + ! local type(CellDefnType), pointer :: defn integer(I4B) :: ipv, ipv1, ipv2, ipv3, ipv4 integer(I4B), dimension(4) :: ipvnxt = (/2, 3, 4, 1/) @@ -29,19 +29,19 @@ subroutine cell_poly_to_rect(poly, rect) call create_cell_rect(rect) defn => poly%defn - ! -- Translate and rotate the rectangular cell into local coordinates - ! -- with x varying from 0 to dx and y varying from 0 to dy. Choose the - ! -- "south-west" vertex to be the local origin so that the rotation - ! -- angle is zero if the cell already aligns with the model x and y - ! -- coordinates. The "southwest" vertex is found by finding the vertex - ! -- formed either by (1) an edge over which x decreases (going - ! -- clockwise) followed by an edge over which x does not increase, or - ! -- by (2) an edge over which y does not decrease (again going - ! -- clockwise) followed by an edge over which y increases. In the end, - ! -- ipv1 is the local vertex number (within the cell, taking a value - ! -- of 1, 2, 3, or 4) of the southwest vertex, and ipv2, ipv3, and - ! -- ipv4 are the local vertex numbers of the remaining three vertices - ! -- going clockwise. + ! Translate and rotate the rectangular cell into local coordinates + ! with x varying from 0 to dx and y varying from 0 to dy. Choose the + ! "south-west" vertex to be the local origin so that the rotation + ! angle is zero if the cell already aligns with the model x and y + ! coordinates. The "southwest" vertex is found by finding the vertex + ! formed either by (1) an edge over which x decreases (going + ! clockwise) followed by an edge over which x does not increase, or + ! by (2) an edge over which y does not decrease (again going + ! clockwise) followed by an edge over which y increases. In the end, + ! ipv1 is the local vertex number (within the cell, taking a value + ! of 1, 2, 3, or 4) of the southwest vertex, and ipv2, ipv3, and + ! ipv4 are the local vertex numbers of the remaining three vertices + ! going clockwise. do ipv = 1, 4 ipv4 = ipv ipv1 = ipvnxt(ipv4) @@ -67,9 +67,9 @@ subroutine cell_poly_to_rect(poly, rect) end do ipv3 = ipvnxt(ipv2) - ! -- Compute upper bounds on the local coordinates (the rectangular - ! -- dimensions of the cell) and the sine and cosine of the rotation - ! -- angle, and store local origin information + ! Compute upper bounds on the local coordinates (the rectangular + ! dimensions of the cell) and the sine and cosine of the rotation + ! angle, and store local origin information xOrigin = x1 yOrigin = y1 zOrigin = defn%bot @@ -93,7 +93,7 @@ subroutine cell_poly_to_rect(poly, rect) rect%zOrigin = zOrigin rect%ipvOrigin = ipv1 - ! -- Compute (unscaled) cell edge velocities from face flows + ! Compute (unscaled) cell edge velocities from face flows areax = dx * dz areay = dy * dz areaz = dx * dy @@ -115,24 +115,24 @@ subroutine cell_poly_to_quad(poly, quad) use CellRectQuadModule, only: CellRectQuadType, create_cell_rect_quad use CellPolyModule, only: CellPolyType use MathUtilModule, only: mod_offset - ! -- dummy + ! dummy type(CellPolyType), intent(in), pointer :: poly type(CellRectQuadType), intent(inout), pointer :: quad - ! -- local + ! local integer(I4B) :: i, irv, isc real(DP) :: qhalf, qdisttopbot, q1, q2, q4 call create_cell_rect_quad(quad) call quad%init_from(poly%defn) - ! -- Translate and rotate the rect-quad cell into local coordinates with - ! -- x varying from 0 to dx and y varying from 0 to dy. Choose the "south- - ! -- west" rectangle vertex to be the local origin so that the rotation - ! -- angle is zero if the cell already aligns with the model x and y - ! -- coordinates. + ! Translate and rotate the rect-quad cell into local coordinates with + ! x varying from 0 to dx and y varying from 0 to dy. Choose the "south- + ! west" rectangle vertex to be the local origin so that the rotation + ! angle is zero if the cell already aligns with the model x and y + ! coordinates. quad%irvOrigin = quad%get_rect_ivert_sw() call quad%get_rect_dim_rot() - ! -- Set the external and internal face flows used for subcells + ! Set the external and internal face flows used for subcells do i = 0, 3 irv = mod_offset(i + quad%irvOrigin, 4, 1) isc = mod_offset(i + 3, 4, 1) diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index 03a401173c2..9d653250dae 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -432,15 +432,6 @@ subroutine vertvelo(this) ! TODO: can we use some of the terms of the centroid calculation ! to do a cheap point in polygon check? - ! - ! allocate(poly(2, nvert)) - ! poly(1,:) = this%xvert - ! poly(2,:) = this%yvert - ! if (.not. point_in_polygon(this%xctr, this%yctr, poly)) then - ! print *, "error -- centroid not in cell ", this%cell%defn%icell, this%xctr, this%yctr - ! call pstop(1) - ! end if - ! deallocate(poly) ! Subcell areas do i = 1, this%nverts diff --git a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 index b5c16540786..c002a567b5a 100644 --- a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 @@ -28,9 +28,9 @@ module MethodSubcellPollockModule !> @brief Create a new Pollock's subcell method subroutine create_method_subcell_pollock(method) - ! -- dummy + ! dummy type(MethodSubcellPollockType), pointer :: method - ! -- local + ! local type(SubcellRectType), pointer :: subcell allocate (method) @@ -48,11 +48,11 @@ end subroutine deallocate !> @brief Apply Pollock's method to a rectangular subcell subroutine apply_msp(this, particle, tmax) - ! -- dummy + ! dummy class(MethodSubcellPollockType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle real(DP), intent(in) :: tmax - ! -- local + ! local real(DP) :: xOrigin real(DP) :: yOrigin real(DP) :: zOrigin @@ -61,7 +61,7 @@ subroutine apply_msp(this, particle, tmax) select type (subcell => this%subcell) type is (SubcellRectType) - ! -- Transform particle position into local subcell coordinates, + ! Transform particle position into local subcell coordinates, ! track particle across subcell, convert back to model coords ! (sinrot and cosrot should be 0 and 1, respectively, i.e. no ! rotation, also no z translation; only x and y translations) @@ -302,7 +302,7 @@ function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status) real(DP) :: v1v2 logical(LGP) :: noOutflow - ! -- Initialize variables. + ! Initialize variables. status = -1 dt = 1.0d+20 v2a = v2 @@ -313,8 +313,8 @@ function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status) dva = dv if (dva .lt. DZERO) dva = -dva - ! -- Check for a uniform zero velocity in this direction. - ! -- If so, set status = 2 and return (dt = 1.0d+20). + ! Check for a uniform zero velocity in this direction. + ! If so, set status = 2 and return (dt = 1.0d+20). tol = 1.0d-15 if ((v2a .lt. tol) .and. (v1a .lt. tol)) then v = DZERO @@ -323,9 +323,9 @@ function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status) return end if - ! -- Check for uniform non-zero velocity in this direction. - ! -- If so, set compute dt using the constant velocity, - ! -- set status = 1 and return. + ! Check for uniform non-zero velocity in this direction. + ! If so, set compute dt using the constant velocity, + ! set status = 1 and return. vv = v1a if (v2a .gt. vv) vv = v2a vvv = dva / vv @@ -341,13 +341,13 @@ function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status) return end if - ! -- Velocity has a linear variation. - ! -- Compute velocity corresponding to particle position. + ! Velocity has a linear variation. + ! Compute velocity corresponding to particle position. dvdx = dv / dx v = (DONE - xL) * v1 + xL * v2 - ! -- If flow is into the cell from both sides there is no outflow. - ! -- In that case, set status = 3 and return. + ! If flow is into the cell from both sides there is no outflow. + ! In that case, set status = 3 and return. noOutflow = .true. if (v1 .lt. DZERO) noOutflow = .false. if (v2 .gt. DZERO) noOutflow = .false. @@ -356,10 +356,10 @@ function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status) return end if - ! -- If there is a divide in the cell for this flow direction, check to - ! -- see if the particle is located exactly on the divide. If it is, move - ! -- it very slightly to get it off the divide. This avoids possible - ! -- numerical problems related to stagnation points. + ! If there is a divide in the cell for this flow direction, check to + ! see if the particle is located exactly on the divide. If it is, move + ! it very slightly to get it off the divide. This avoids possible + ! numerical problems related to stagnation points. if ((v1 .le. DZERO) .and. (v2 .ge. DZERO)) then if (abs(v) .le. DZERO) then v = 1.0d-20 @@ -367,9 +367,9 @@ function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status) end if end if - ! -- If there is a flow divide, this check finds out what side of the - ! -- divide the particle is on and sets the value of vr appropriately - ! -- to reflect that location. + ! If there is a flow divide, this check finds out what side of the + ! divide the particle is on and sets the value of vr appropriately + ! to reflect that location. vr1 = v1 / v vr2 = v2 / v vr = vr1 @@ -377,17 +377,17 @@ function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status) vr = vr2 end if - ! -- If the product v1*v2 > 0, the velocity is in the same direction - ! -- throughout the cell (i.e. no flow divide). If so, set the value - ! -- of vr to reflect the appropriate direction. + ! If the product v1*v2 > 0, the velocity is in the same direction + ! throughout the cell (i.e. no flow divide). If so, set the value + ! of vr to reflect the appropriate direction. v1v2 = v1 * v2 if (v1v2 .gt. DZERO) then if (v .gt. DZERO) vr = vr2 if (v .lt. DZERO) vr = vr1 end if - ! -- Check if vr is (very close to) zero. - ! -- If so, set status = 2 and return (dt = 1.0d+20). + ! Check if vr is (very close to) zero. + ! If so, set status = 2 and return (dt = 1.0d+20). if (dabs(vr) .lt. 1.0d-10) then v = DZERO dvdx = DZERO @@ -395,7 +395,7 @@ function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status) return end if - ! -- Compute travel time to exit face. Return with status = 0. + ! Compute travel time to exit face. Return with status = 0. dt = log(vr) / dvdx status = 0 @@ -422,14 +422,14 @@ pure function new_x(v, dvdx, v1, v2, dt, x, dx, velocity_profile) result(newx) real(DP) :: newx logical(LGP) :: lprofile - ! -- process optional arguments + ! process optional arguments if (present(velocity_profile)) then lprofile = velocity_profile else lprofile = .false. end if - ! -- recompute coordinate + ! recompute coordinate newx = x if (lprofile) then newx = newx + (v1 * dt / dx) @@ -437,7 +437,7 @@ pure function new_x(v, dvdx, v1, v2, dt, x, dx, velocity_profile) result(newx) newx = newx + (v * (exp(dvdx * dt) - DONE) / dvdx / dx) end if - ! -- clamp to [0, 1] + ! clamp to [0, 1] if (newx .lt. DZERO) newx = DZERO if (newx .gt. DONE) newx = DONE diff --git a/src/Solution/ParticleTracker/TernarySolveTrack.f90 b/src/Solution/ParticleTracker/TernarySolveTrack.f90 index d5fb5318ef2..8fb97847b0d 100644 --- a/src/Solution/ParticleTracker/TernarySolveTrack.f90 +++ b/src/Solution/ParticleTracker/TernarySolveTrack.f90 @@ -33,7 +33,7 @@ subroutine traverse_triangle(isolv, tol, texit, & alpexit, betexit, & itrifaceenter, itrifaceexit, & alp1, bet1, alp2, bet2, alpi, beti) - ! -- dummy + ! dummy integer(I4B), intent(in) :: isolv !< solution method real(DP), intent(in) :: tol !< solution tolerance real(DP), intent(out) :: texit !< time particle exits the cell @@ -47,7 +47,7 @@ subroutine traverse_triangle(isolv, tol, texit, & real(DP) :: bet2 real(DP) :: alpi real(DP) :: beti !< alpha and beta coefficients - ! -- local + ! local real(DP) :: texit0 real(DP) :: alpexit0 real(DP) :: betexit0 @@ -58,13 +58,13 @@ subroutine traverse_triangle(isolv, tol, texit, & real(DP) :: alpexit2 real(DP) :: betexit2 - ! -- Compute elements of matrix W + ! Compute elements of matrix W call get_w(alp1, bet1, alp2, bet2, waa, wab, wba, wbb) - ! -- Determine alpha and beta analytically as functions of time + ! Determine alpha and beta analytically as functions of time call solve_coefs(alpi, beti) - ! -- Compute exit time (travel time to exit) and exit location + ! Compute exit time (travel time to exit) and exit location call find_exit_bary(isolv, 0, itrifaceenter, & alpi, beti, tol, & texit0, alpexit0, betexit0) @@ -76,8 +76,8 @@ subroutine traverse_triangle(isolv, tol, texit, & texit2, alpexit2, betexit2) texit = min(texit0, texit1, texit2) - ! -- Note that while the numbering of triangle faces is generally zero-based - ! -- (0, 1, 2), itrifaceexit, which gets passed out, is one-based (1, 2, 3). + ! Note that while the numbering of triangle faces is generally zero-based + ! (0, 1, 2), itrifaceexit, which gets passed out, is one-based (1, 2, 3). if (texit .eq. texit0) then alpexit = alpexit0 betexit = betexit0 @@ -102,7 +102,7 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, & rxx, rxy, ryx, ryy, & sxx, sxy, syy, & alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti) - ! -- dummy + ! dummy real(DP) :: x0 real(DP) :: y0 real(DP) :: x1 @@ -130,7 +130,7 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, & real(DP) :: bet2 real(DP) :: alpi real(DP) :: beti !< alpha and beta coefficients - ! -- local + ! local real(DP) :: baselen real(DP) :: oobaselen real(DP) :: sinomega @@ -143,7 +143,7 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, & real(DP) :: yidiff real(DP) :: rot(2, 2), res(2) - ! -- Translate and rotate coordinates to "canonical" configuration + ! Translate and rotate coordinates to "canonical" configuration x1diff = x1 - x0 y1diff = y1 - y0 x2diff = x2 - x0 @@ -196,7 +196,7 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, & subroutine get_w( & alp1, bet1, alp2, bet2, & waa, wab, wba, wbb) - ! -- dummy + ! dummy real(DP) :: alp1 real(DP) :: bet1 real(DP) :: alp2 @@ -205,12 +205,12 @@ subroutine get_w( & real(DP) :: wab real(DP) :: wba real(DP) :: wbb !< w matrix - ! -- local + ! local real(DP) :: v1alpdiff real(DP) :: v2alpdiff real(DP) :: v2betdiff - ! -- Note: wab is the "alpha,beta" entry in matrix W + ! Note: wab is the "alpha,beta" entry in matrix W ! and the alpha component of the w^(beta) vector v1alpdiff = cv1(1) - cv0(1) v2alpdiff = cv2(1) - cv0(1) @@ -224,10 +224,10 @@ subroutine get_w( & !> @brief Compute analytical solution coefficients depending on case subroutine solve_coefs(alpi, beti) - ! -- dummy + ! dummy real(DP) :: alpi real(DP) :: beti - ! -- local + ! local real(DP) :: zerotol real(DP) :: wratv real(DP) :: acoef @@ -244,28 +244,28 @@ subroutine solve_coefs(alpi, beti) bcoef = wratv + wab * beti afact = acoef / waa vfact = cv0(2) / wbb - ! -- Coefs for beta do not depend on whether waa = 0 or not + ! Coefs for beta do not depend on whether waa = 0 or not cb1 = -vfact ! const term in beta cb2 = vfact + beti ! coef for e(wbb*t) term in beta - ! -- Coefs for alpha + ! Coefs for alpha if (dabs(waa) .gt. zerotol) then - ! -- Case waa <> 0, wbb <> 0 + ! Case waa <> 0, wbb <> 0 if (dabs(wbb - waa) .gt. zerotol) then - ! -- Subcase wbb <> waa + ! Subcase wbb <> waa bfact = bcoef / (wbb - waa) ca1 = -afact ! const term in alpha ca2 = alpi + afact - bfact ! coef for exp(waa*t) term in alpha ca3 = bfact ! coef for exp(wbb*t) term in alpha icase = 1 else - ! -- Subcase wbb = waa + ! Subcase wbb = waa ca1 = -afact ! const term in alpha ca2 = alpi + afact ! coef for exp(waa*t) term in alpha ca3 = bcoef ! coef for t*exp(waa*t) term in alpha icase = -1 end if else - ! -- Case waa = 0, wbb <> 0 + ! Case waa = 0, wbb <> 0 bfact = bcoef / wbb ca1 = alpi - bfact ! const term in alpha ca2 = acoef ! coef for t term in alpha @@ -273,11 +273,11 @@ subroutine solve_coefs(alpi, beti) icase = 2 end if else - ! -- Coefs for beta do not depend on whether waa = 0 or not + ! Coefs for beta do not depend on whether waa = 0 or not cb1 = beti ! const term in beta cb2 = cv0(2) ! coef for t term in beta if (dabs(waa) .gt. zerotol) then - ! -- Case waa <> 0, wbb = 0 + ! Case waa <> 0, wbb = 0 oowaa = 1d0 / waa vfact = (wab * oowaa) * cv0(2) ca1 = -oowaa * (cv0(1) + wab * beti + vfact) ! const term in alpha @@ -285,7 +285,7 @@ subroutine solve_coefs(alpi, beti) ca3 = alpi - ca1 ! coef for exp(waa*t) term in alpha icase = 3 else - ! -- Case waa = 0, wbb = 0 + ! Case waa = 0, wbb = 0 ca1 = alpi ! const term in alpha ca2 = cv0(1) + wab * beti ! coef for t term in alpha ca3 = 5d-1 * wab * cv0(2) ! coef for t^2 term in alpha @@ -297,7 +297,7 @@ subroutine solve_coefs(alpi, beti) !> @brief Step (evaluate) analytically depending on case subroutine step_analytical(t, alp, bet) - ! -- dummy + ! dummy real(DP), intent(in) :: t real(DP) :: alp real(DP) :: bet @@ -325,7 +325,7 @@ subroutine step_analytical(t, alp, bet) subroutine find_exit_bary(isolv, itriface, itrifaceenter, & alpi, beti, tol, & texit, alpexit, betexit) - ! -- dummy + ! dummy integer(I4B) :: isolv integer(I4B) :: itriface integer(I4B) :: itrifaceenter @@ -335,7 +335,7 @@ subroutine find_exit_bary(isolv, itriface, itrifaceenter, & real(DP) :: texit real(DP) :: alpexit real(DP) :: betexit - ! -- local + ! local real(DP) :: alplo real(DP) :: alphi real(DP) :: alpt @@ -361,82 +361,82 @@ subroutine find_exit_bary(isolv, itriface, itrifaceenter, & zerotol = 1d-10 ! todo AMP: consider tolerance - ! -- Use iterative scheme or numerical integration indicated by isolv. + ! Use iterative scheme or numerical integration indicated by isolv. if (itriface .eq. 0) then - ! -- Checking for exit on canonical face 0 (beta = 0) + ! Checking for exit on canonical face 0 (beta = 0) if (itrifaceenter .eq. 0) then - ! -- Entrance face, so no exit. (Normal velocity is uniform along face 0, - ! -- so it cannot be both an entrance and an exit.) + ! Entrance face, so no exit. (Normal velocity is uniform along face 0, + ! so it cannot be both an entrance and an exit.) texit = huge(1d0) else - ! -- Not the entrance face, so check for outflow + ! Not the entrance face, so check for outflow if (cv0(2) .ge. DZERO) then ! check beta velocity along beta=0 face - ! -- Inflow or no flow, so no exit + ! Inflow or no flow, so no exit texit = huge(DONE) else - ! -- Outflow, so check beta-velocity at the initial location, - ! -- recognizing that it will never change sign along the - ! -- trajectory (and will not be blocked from zero by an asymptote) + ! Outflow, so check beta-velocity at the initial location, + ! recognizing that it will never change sign along the + ! trajectory (and will not be blocked from zero by an asymptote) vbeti = cv0(2) + wbb * beti if (vbeti .ge. DZERO) then - ! -- Can't exit along beta = 0 + ! Can't exit along beta = 0 texit = huge(DONE) else - ! -- get alpt and check it + ! get alpt and check it call get_t_alpt(DZERO, t, alpt) if ((alpt .ge. DZERO) .and. (alpt .le. DONE)) then - ! -- alpt within the edge, so exit found + ! alpt within the edge, so exit found texit = t alpexit = alpt betexit = DZERO else - ! -- alpt not within the edge, so not an exit + ! alpt not within the edge, so not an exit texit = huge(DONE) end if end if end if end if - ! -- End canonical face 0 (beta = 0) + ! End canonical face 0 (beta = 0) else - ! -- Checking for exit on canonical face 1 (gamma = 0.) or 2 (alpha = 0.) + ! Checking for exit on canonical face 1 (gamma = 0.) or 2 (alpha = 0.) if (itriface .eq. 1) then - ! -- Normal velocities (gamma components) at ends of canonical face 1 + ! Normal velocities (gamma components) at ends of canonical face 1 v1n = -cv1(1) - cv1(2) v2n = -cv2(1) - cv2(2) else - ! -- Normal velocities (alpha components) at ends of canonical face 2 + ! Normal velocities (alpha components) at ends of canonical face 2 v1n = cv0(1) v2n = cv2(1) end if if ((v1n .ge. DZERO) .and. (v2n .ge. DZERO)) then - ! -- No outflow at vn1 and vn2 corners; no outflow interval, so no exit. + ! No outflow at vn1 and vn2 corners; no outflow interval, so no exit. texit = huge(DONE) else - ! -- Find outflow interval + ! Find outflow interval call get_bet_outflow_bary(v1n, v2n, betoutlo, betouthi) - ! -- Find trend of and limits on beta from beta{t} solution + ! Find trend of and limits on beta from beta{t} solution call get_bet_soln_limits(beti, betsollo, betsolhi, ibettrend) - ! -- Look for exit + ! Look for exit if (ibettrend .eq. 0) then - ! -- Beta is constant, so check if it's within the outflow interval; - ! -- if not, no exit; if so, solve for t and alpha + ! Beta is constant, so check if it's within the outflow interval; + ! if not, no exit; if so, solve for t and alpha if ((beti .gt. betouthi) .or. (beti .lt. betoutlo)) then texit = huge(1d0) else - ! -- Check alpha-velocity at the initial location, - ! -- recognizing that it will never change sign along the - ! -- trajectory (and will not be blocked from zero by an asymptote) - ! -- in this special case + ! Check alpha-velocity at the initial location, + ! recognizing that it will never change sign along the + ! trajectory (and will not be blocked from zero by an asymptote) + ! in this special case v0alpstar = cv0(1) + wab * beti valpi = v0alpstar + waa * alpi if ((itriface .eq. 1) .and. (valpi .le. DZERO)) then - ! -- Can't exit along gamma = 0. + ! Can't exit along gamma = 0. texit = huge(DONE) else if ((itriface .eq. 2) .and. (valpi .ge. DZERO)) then - ! -- Can't exit along alpha = 0. + ! Can't exit along alpha = 0. texit = huge(DONE) else - ! -- get exit + ! get exit if (itriface .eq. 1) then alpexit = DONE - beti else @@ -451,16 +451,16 @@ subroutine find_exit_bary(isolv, itriface, itrifaceenter, & end if end if end if - ! -- End constant-beta case + ! End constant-beta case else - ! -- Beta varies along trajectory; combine outflow and soln limits on beta + ! Beta varies along trajectory; combine outflow and soln limits on beta bethi = min(betouthi, betsolhi) betlo = max(betoutlo, betsollo) if (betlo .gt. bethi) then - ! -- If bounds on bet leave no feasible interval, no exit + ! If bounds on bet leave no feasible interval, no exit texit = huge(DONE) else - ! -- Check sign of function value at beta bounds + ! Check sign of function value at beta bounds call get_t_alpt(bethi, thi, alphi) call get_t_alpt(betlo, tlo, alplo) if (itriface .eq. 1) then @@ -471,17 +471,17 @@ subroutine find_exit_bary(isolv, itriface, itrifaceenter, & fbx = alphi end if if (fax * fbx .gt. DZERO) then - ! -- Root not bracketed; no exit + ! Root not bracketed; no exit texit = huge(DONE) else if (isolv .eq. 1) then - ! -- Use Brent's method with initial bounds on beta of betlo and bethi, - ! -- assuming they bound the root + ! Use Brent's method with initial bounds on beta of betlo and bethi, + ! assuming they bound the root call soln_brent(itriface, betlo, bethi, tol, texit, & alpexit, betexit) else if (isolv .eq. 2) then - ! -- Use Chandrupatla's method with initial bounds on beta of betlo and bethi, - ! -- assuming they bound the root + ! Use Chandrupatla's method with initial bounds on beta of betlo and bethi, + ! assuming they bound the root call soln_chand(itriface, betlo, bethi, tol, texit, & alpexit, betexit) else @@ -489,10 +489,10 @@ subroutine find_exit_bary(isolv, itriface, itrifaceenter, & end if end if end if - ! -- End variable-beta case + ! End variable-beta case end if end if - ! -- End canonical face 1 (gamma = 0.) or 2 (alpha = 0.) + ! End canonical face 1 (gamma = 0.) or 2 (alpha = 0.) end if if (texit .ne. huge(DONE) .and. texit .lt. DZERO) then @@ -503,39 +503,39 @@ subroutine find_exit_bary(isolv, itriface, itrifaceenter, & !> @brief Brent's method applied to canonical face 1 (gamma = 0) function fbary1(bet) result(fb) - ! -- dummy + ! dummy real(DP), intent(in) :: bet real(DP) :: fb - ! -- local + ! local real(DP) :: t real(DP) :: alpt - ! -- Evaluate gamma{t{beta}} = 1. - alpha{t{beta}} - beta + ! Evaluate gamma{t{beta}} = 1. - alpha{t{beta}} - beta call get_t_alpt(bet, t, alpt) fb = DONE - alpt - bet end function !> @brief Brent's method applied to canonical face 2 (alpha = 0) function fbary2(bet) result(fb) - ! -- dummy + ! dummy real(DP), intent(in) :: bet real(DP) :: fb - ! -- local + ! local real(DP) :: t real(DP) :: alpt - ! -- Evaluate alpha{t{beta}} + ! Evaluate alpha{t{beta}} call get_t_alpt(bet, t, alpt) fb = alpt end function !> @brief Given beta evaluate t and alpha depending on case subroutine get_t_alpt(bet, t, alp) - ! -- dummy + ! dummy real(DP), intent(in) :: bet real(DP) :: t real(DP) :: alp - ! -- local + ! local real(DP) :: term real(DP) :: zerotol real(DP) :: waat @@ -589,33 +589,33 @@ subroutine get_t_alpt(bet, t, alp) !> @brief Find outflow interval subroutine get_bet_outflow_bary(vn1, vn2, betoutlo, betouthi) - ! -- dummy + ! dummy real(DP) :: vn1 real(DP) :: vn2 real(DP) :: betoutlo real(DP) :: betouthi - ! -- local + ! local real(DP) :: vndiff vndiff = vn2 - vn1 if (vn1 .lt. DZERO) then - ! -- Outflow at vn1 corner + ! Outflow at vn1 corner betoutlo = DZERO if (vn2 .le. DZERO) then - ! -- Outflow along entire edge (except possibly no-flow right at vn2 corner) + ! Outflow along entire edge (except possibly no-flow right at vn2 corner) betouthi = DONE else - ! -- Outflow along part of edge + ! Outflow along part of edge betouthi = -vn1 / vndiff end if else - ! -- Outflow at vn2 corner + ! Outflow at vn2 corner betouthi = DONE if (vn1 .le. DZERO) then - ! -- Outflow along entire edge (except possibly no-flow right at vn1 corner) + ! Outflow along entire edge (except possibly no-flow right at vn1 corner) betoutlo = DZERO else - ! -- Outflow along part of edge + ! Outflow along part of edge betoutlo = -vn1 / vndiff end if end if @@ -624,12 +624,12 @@ subroutine get_bet_outflow_bary(vn1, vn2, betoutlo, betouthi) !> @brief Find trend of and limits on beta from beta{t} solution subroutine get_bet_soln_limits(beti, betsollo, betsolhi, ibettrend) - ! -- dummy + ! dummy real(DP), intent(in) :: beti real(DP) :: betsollo real(DP) :: betsolhi integer(I4B), intent(inout) :: ibettrend - ! -- local + ! local real(DP) :: betlim if (icase > 2) then @@ -683,7 +683,7 @@ subroutine get_bet_soln_limits(beti, betsollo, betsolhi, ibettrend) !> @brief Use Brent's method with initial bounds on beta of betlo and bethi subroutine soln_brent(itriface, betlo, bethi, tol, & texit, alpexit, betexit) - ! -- dummy + ! dummy integer(I4B), intent(in) :: itriface real(DP) :: betlo real(DP) :: bethi @@ -691,12 +691,12 @@ subroutine soln_brent(itriface, betlo, bethi, tol, & real(DP) :: texit real(DP) :: alpexit real(DP) :: betexit - ! -- local + ! local real(DP) :: blo real(DP) :: bhi procedure(f1d), pointer :: f - ! -- assuming betlo and bethi bracket the root + ! assuming betlo and bethi bracket the root blo = betlo bhi = bethi if (itriface .eq. 1) then @@ -713,7 +713,7 @@ subroutine soln_brent(itriface, betlo, bethi, tol, & !> @brief Use Chandrupatla's method with initial bounds on beta of betlo and bethi subroutine soln_chand(itriface, betlo, bethi, tol, & texit, alpexit, betexit) - ! -- dummy + ! dummy integer(I4B), intent(in) :: itriface real(DP) :: betlo real(DP) :: bethi @@ -721,12 +721,12 @@ subroutine soln_chand(itriface, betlo, bethi, tol, & real(DP) :: texit real(DP) :: alpexit real(DP) :: betexit - ! -- local + ! local real(DP) :: blo real(DP) :: bhi procedure(f1d), pointer :: f - ! -- note: assuming betlo and bethi bracket the root + ! note: assuming betlo and bethi bracket the root blo = betlo bhi = bethi if (itriface .eq. 1) then