diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 22827f611..60748c5b0 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -112,9 +112,10 @@ module subroutine encounter_io_initialize(self, param) integer(I4B) :: nvar, varid, vartype real(DP) :: dfill real(SP) :: sfill + integer(I4B), parameter :: NO_FILL = 0 logical :: fileExists character(len=STRMAX) :: errmsg - integer(I4B) :: ndims, i + integer(I4B) :: ndims associate(nc => self) dfill = ieee_value(dfill, IEEE_QUIET_NAN) @@ -167,13 +168,13 @@ module subroutine encounter_io_initialize(self, param) call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" ) end select end do diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index e4f0c1fbc..9a1be2901 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -690,6 +690,8 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k) call orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi) snapshot%t = t + tperi + if ((snapshot%t < maxval(pl_snap%info(:)%origin_time)) .or. & + (snapshot%t > minval(pl_snap%info(:)%discard_time))) cycle ! Computer the center mass of the pair rcom(:) = (plplenc_list%r1(:,k) * pl_snap%Gmass(1) + plplenc_list%r2(:,k) * pl_snap%Gmass(2)) / Gmtot diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index f47a64047..6f75c7d01 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -26,9 +26,10 @@ module subroutine fraggle_io_initialize_output(self, param) integer(I4B) :: nvar, varid, vartype real(DP) :: dfill real(SP) :: sfill + integer(I4B), parameter :: NO_FILL = 0 logical :: fileExists character(len=STRMAX) :: errmsg - integer(I4B) :: i, ndims + integer(I4B) :: ndims select type(param) class is (symba_parameters) @@ -120,13 +121,13 @@ module subroutine fraggle_io_initialize_output(self, param) call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "fraggle_io_initialize nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "fraggle_io_initialize nf90_def_var_fill NF90_INT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "fraggle_io_initialize nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "fraggle_io_initialize nf90_def_var_fill NF90_FLOAT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "fraggle_io_initialize nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "fraggle_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "fraggle_io_initialize nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "fraggle_io_initialize nf90_def_var_fill NF90_CHAR" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "fraggle_io_initialize nf90_def_var_fill NF90_CHAR" ) end select end do ! Take the file out of define mode diff --git a/src/io/io.f90 b/src/io/io.f90 index f159e6ac7..54f99e42b 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -294,11 +294,11 @@ module subroutine io_dump_storage(self, param) integer(I8B) :: iloop_start if (self%iframe == 0) return - iloop_start = param%iloop - int(param%istep_out * param%dump_cadence, kind=I8B) + 1 + iloop_start = max(param%iloop - int(param%istep_out * param%dump_cadence, kind=I8B) + 1_I8B,0_I8B) call self%make_index_map() - do i = 1, param%dump_cadence - param%ioutput = iloop_start + self%tmap(i) + do i = 1, self%iframe if (allocated(self%frame(i)%item)) then + param%ioutput = iloop_start + self%tmap(i) select type(system => self%frame(i)%item) class is (swiftest_nbody_system) call system%write_frame(param) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 1d5a708b3..1a2c2ef8e 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -232,18 +232,20 @@ module subroutine symba_collision_make_colliders_pl(self,idx) integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision end subroutine symba_collision_make_colliders_pl - module subroutine symba_resolve_collision_fragmentations(self, system, param) + module subroutine symba_resolve_collision_fragmentations(self, system, param, t) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision end subroutine symba_resolve_collision_fragmentations - module subroutine symba_resolve_collision_mergers(self, system, param) + module subroutine symba_resolve_collision_mergers(self, system, param, t) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision end subroutine symba_resolve_collision_mergers module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, irec) @@ -343,25 +345,28 @@ pure module subroutine symba_gr_p4_tp(self, system, param, dt) real(DP), intent(in) :: dt !! Step size end subroutine symba_gr_p4_tp - module function symba_collision_casedisruption(system, param) result(status) + module function symba_collision_casedisruption(system, param, t) result(status) implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casedisruption - module function symba_collision_casehitandrun(system, param) result(status) + module function symba_collision_casehitandrun(system, param, t) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casehitandrun - module function symba_collision_casemerge(system, param) result(status) + module function symba_collision_casemerge(system, param, t) result(status) implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casemerge module subroutine symba_util_set_renc(self, scale) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 588a138d6..4ce217124 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -157,6 +157,7 @@ module subroutine netcdf_initialize_output(self, param) integer(I4B) :: nvar, varid, vartype real(DP) :: dfill real(SP) :: sfill + integer(I4B), parameter :: NO_FILL = 0 logical :: fileExists character(len=STRMAX) :: errmsg integer(I4B) :: ndims @@ -281,16 +282,24 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "netcdf_initialize_output nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "netcdf_initialize_output nf90_def_var_fill NF90_INT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "netcdf_initialize_output nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "netcdf_initialize_output nf90_def_var_fill NF90_FLOAT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "netcdf_initialize_output nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) end select end do + ! Set special fill mode for discard time so that we can make use of it for non-discarded bodies. + select case (vartype) + case(NF90_FLOAT) + call check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_SP)), "netcdf_initialize_output nf90_def_var_fill discard_time NF90_FLOAT" ) + case(NF90_DOUBLE) + call check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_DP)), "netcdf_initialize_output nf90_def_var_fill discard_time NF90_DOUBLE" ) + end select + ! Take the file out of define mode call check( nf90_enddef(nc%id), "netcdf_initialize_output nf90_enddef" ) @@ -944,7 +953,7 @@ module subroutine netcdf_read_particle_info_system(self, nc, param, plmask, tpma if (status == nf90_noerr) then call check( nf90_get_var(nc%id, nc%collision_id_varid, itemp), "netcdf_read_particle_info_system nf90_getvar collision_id_varid" ) else - itemp = 0.0_DP + itemp = 0 end if do i = 1, npl @@ -958,7 +967,12 @@ module subroutine netcdf_read_particle_info_system(self, nc, param, plmask, tpma if (status == nf90_noerr) then call check( nf90_get_var(nc%id, nc%discard_time_varid, rtemp), "netcdf_read_particle_info_system nf90_getvar discard_time_varid" ) else - rtemp = 0.0_DP + select case (param%out_type) + case("NETCDF_FLOAT") + rtemp(:) = huge(0.0_SP) + case("NETCDF_DOUBLE") + rtemp(:) = huge(0.0_DP) + end select end if call cb%info%set_value(discard_time=rtemp(1)) diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index ef8558aef..57679f622 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -248,7 +248,7 @@ module subroutine setup_body(self, n, param) origin_time = -huge(1.0_DP), & origin_rh = [0.0_DP, 0.0_DP, 0.0_DP], & origin_vh = [0.0_DP, 0.0_DP, 0.0_DP], & - discard_time = -huge(1.0_DP), & + discard_time = huge(1.0_DP), & discard_rh = [0.0_DP, 0.0_DP, 0.0_DP], & discard_vh = [0.0_DP, 0.0_DP, 0.0_DP], & discard_body_id = -1 & diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 67540bfa3..0353cc7d5 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -12,15 +12,16 @@ contains - module function symba_collision_casedisruption(system, param) result(status) + module function symba_collision_casedisruption(system, param, t) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Create the fragments resulting from a non-catastrophic disruption collision !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals @@ -73,7 +74,7 @@ module function symba_collision_casedisruption(system, param) result(status) param%maxid = fragments%id(nfrag) end select - call symba_collision_mergeaddsub(system, param, status) + call symba_collision_mergeaddsub(system, param, t, status) end if end associate @@ -81,15 +82,16 @@ module function symba_collision_casedisruption(system, param) result(status) end function symba_collision_casedisruption - module function symba_collision_casehitandrun(system, param) result(status) + module function symba_collision_casehitandrun(system, param, t) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Create the fragments resulting from a non-catastrophic hit-and-run collision !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Result integer(I4B) :: status !! Status flag assigned to this outcom ! Internals @@ -144,7 +146,7 @@ module function symba_collision_casehitandrun(system, param) result(status) fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] param%maxid = fragments%id(nfrag) status = HIT_AND_RUN_DISRUPT - call symba_collision_mergeaddsub(system, param, status) + call symba_collision_mergeaddsub(system, param, t, status) end if end associate @@ -153,7 +155,7 @@ module function symba_collision_casehitandrun(system, param) result(status) end function symba_collision_casehitandrun - module function symba_collision_casemerge(system, param) result(status) + module function symba_collision_casemerge(system, param, t) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Merge massive bodies. @@ -163,8 +165,9 @@ module function symba_collision_casemerge(system, param) result(status) !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals @@ -226,7 +229,7 @@ module function symba_collision_casemerge(system, param) result(status) status = MERGED - call symba_collision_mergeaddsub(system, param, status) + call symba_collision_mergeaddsub(system, param, t, status) end select end associate @@ -364,8 +367,8 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step pl%lcollision([i, j]) = .true. pl%status([i, j]) = COLLISION - call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i)) - call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,j), discard_vh=pl%vh(:,j)) + call pl%info(i)%set_value(status="COLLISION") + call pl%info(j)%set_value(status="COLLISION") end if else self%r2(:,k) = tp%rh(:,j) + system%cb%rb(:) @@ -719,16 +722,17 @@ module subroutine symba_collision_make_colliders_pl(self, idx) end subroutine symba_collision_make_colliders_pl - subroutine symba_collision_mergeaddsub(system, param, status) + subroutine symba_collision_mergeaddsub(system, param, t, status) !! author: David A. Minton !! !! Fills the pl_discards and pl_adds with removed and added bodies !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), intent(in) :: status !! Status flag to assign to adds + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B), intent(in) :: status !! Status flag to assign to adds ! Internals integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag logical, dimension(system%pl%nbody) :: lmask @@ -775,7 +779,7 @@ subroutine symba_collision_mergeaddsub(system, param, status) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) fragments%id(i) - call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, & + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=t, name=newname, & origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do @@ -785,7 +789,7 @@ subroutine symba_collision_mergeaddsub(system, param, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, & + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -799,14 +803,14 @@ subroutine symba_collision_mergeaddsub(system, param, status) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag write(newname, FRAGFMT) fragments%id(i) - call plnew%info(i)%set_value(origin_type=origin_type, origin_time=system%t, name=newname, & + call plnew%info(i)%set_value(origin_type=origin_type, origin_time=t, name=newname, & origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=system%t, & + call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -817,7 +821,7 @@ subroutine symba_collision_mergeaddsub(system, param, status) if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_rh=pl%rh(:,i), & + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=t, discard_rh=pl%rh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select @@ -880,7 +884,7 @@ subroutine symba_collision_mergeaddsub(system, param, status) end subroutine symba_collision_mergeaddsub - module subroutine symba_resolve_collision_fragmentations(self, system, param) + module subroutine symba_resolve_collision_fragmentations(self, system, param, t) !! author: David A. Minton !! !! Process list of collisions, determine the collisional regime, and then create fragments. @@ -890,13 +894,14 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Internals ! Internals integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision integer(I4B) :: i - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, t => system%t, collision_history => param%collision_history) + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, collision_history => param%collision_history) select type(pl => system%pl) class is (symba_pl) select type (cb => system%cb) @@ -914,11 +919,11 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "before") select case (system%fragments%regime) case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - plplcollision_list%status(i) = symba_collision_casedisruption(system, param) + plplcollision_list%status(i) = symba_collision_casedisruption(system, param, t) case (COLLRESOLVE_REGIME_HIT_AND_RUN) - plplcollision_list%status(i) = symba_collision_casehitandrun(system, param) + plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, t) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - plplcollision_list%status(i) = symba_collision_casemerge(system, param) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, t) case default write(*,*) "Error in symba_collision, unrecognized collision regime" call util_exit(FAILURE) @@ -934,7 +939,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) end subroutine symba_resolve_collision_fragmentations - module subroutine symba_resolve_collision_mergers(self, system, param) + module subroutine symba_resolve_collision_mergers(self, system, param, t) !! author: David A. Minton !! !! Process list of collisions and merge colliding bodies together. @@ -944,6 +949,7 @@ module subroutine symba_resolve_collision_mergers(self, system, param) class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Internals integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision @@ -968,7 +974,7 @@ module subroutine symba_resolve_collision_mergers(self, system, param) fragments%mass_dist(3) = 0.0_DP fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot - plplcollision_list%status(i) = symba_collision_casemerge(system, param) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, t) end do end select end select @@ -1024,9 +1030,9 @@ module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, ir "***********************************************************") allocate(tmp_param, source=param) if (param%lfragmentation) then - call plplcollision_list%resolve_fragmentations(system, param) + call plplcollision_list%resolve_fragmentations(system, param, t) else - call plplcollision_list%resolve_mergers(system, param) + call plplcollision_list%resolve_mergers(system, param, t) end if ! Destroy the collision list now that the collisions are resolved diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 157c3f5af..cf817622e 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -659,15 +659,12 @@ module subroutine symba_util_rearray_pl(self, system, param) pl%lcollision(1:npl) = .false. pl%lmask(1:npl) = .true. - select type(param) - class is (symba_parameters) - pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY - where(pl%lmtiny(1:npl)) - pl%info(1:npl)%particle_type = PL_TINY_TYPE_NAME - elsewhere - pl%info(1:npl)%particle_type = PL_TYPE_NAME - end where - end select + pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY + where(pl%lmtiny(1:npl)) + pl%info(1:npl)%particle_type = PL_TINY_TYPE_NAME + elsewhere + pl%info(1:npl)%particle_type = PL_TYPE_NAME + end where call pl%write_info(param%system_history%nc, param) deallocate(ldump_mask)