Merge branch 'debug'

MintonGroup · Sep 22, 2022 · 35c946f · 35c946f
2 parents 9f4c5f3 + 887d0ef
commit 35c946f
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 1 deletion.
diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90
@@ -106,7 +106,7 @@ program swiftest_driver
                   write(*, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody
                end select
                if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.)
-               call integration_timer%report(message="Integration steps:", nsubsteps=istep_dump)
+               call integration_timer%report(message="Integration steps:", nsubsteps=istep_out)
                call integration_timer%reset()
 
                iout = istep_out

diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90
@@ -602,6 +602,13 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr)
          cb%Gmass = rtemp(1)
          cb%mass = cb%Gmass / param%GU
 
+         ! Set initial central body mass for Helio bookkeeping
+         select type(cb)
+            class is (symba_cb)
+               cb%GM0 = cb%Gmass
+         end select
+
+
          if (npl > 0) then
             pl%Gmass(:) = pack(rtemp, plmask)
             pl%mass(:) = pl%Gmass(:) / param%GU
@@ -615,6 +622,12 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr)
          if (param%lclose) then
             call check( nf90_get_var(iu%ncid, iu%radius_varid, rtemp, start=[1, tslot]) )
             cb%radius = rtemp(1)
+
+            ! Set initial central body radius for SyMBA bookkeeping
+            select type(cb)
+               class is (symba_cb)
+                  cb%R0 = cb%radius
+            end select
             if (npl > 0) pl%radius(:) = pack(rtemp, plmask)
          else
             cb%radius = param%rmin
@@ -645,6 +658,12 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr)
             call check( nf90_get_var(iu%ncid, iu%rotz_varid, rtemp, start=[1, tslot]) )
             cb%rot(3) = rtemp(1)
             if (npl > 0) pl%rot(3,:) = pack(rtemp, plmask)
+
+            ! Set initial central body angular momentum for Helio bookkeeping
+            select type(cb)
+               class is (symba_cb)
+                  cb%L0(:) = cb%Ip(3) * cb%GM0 * cb%R0**2 * cb%rot(:)         
+            end select
          end if
 
          if (param%ltides) then