From 1a8d18f5121fdcbe2b62741e9ac76ca6d6080764 Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Mon, 15 Nov 2021 15:46:56 -0500
Subject: [PATCH] Fixedd calculation of initial angular momentum and change of
 angular momentum of central body for a restarted run

---
 src/io/io.f90         |  5 ++++-
 src/netcdf/netcdf.f90 | 29 +++++++++++++++++++++--------
 2 files changed, 25 insertions(+), 9 deletions(-)

diff --git a/src/io/io.f90 b/src/io/io.f90
index 0266af0bc..b84569cca 100644
--- a/src/io/io.f90
+++ b/src/io/io.f90
@@ -1360,7 +1360,10 @@ subroutine io_read_in_cb(self, param)
             cb%GM0 = cb%Gmass
             cb%dGM = 0.0_DP
             cb%R0 = cb%radius
-            cb%L0(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:)
+            if (param%lrotation) then
+               cb%L0(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:)
+               cb%dL(:) = 0.0_DP
+            end if
          end select
       end if
       return
diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90
index fd2406720..d28fcc788 100644
--- a/src/netcdf/netcdf.f90
+++ b/src/netcdf/netcdf.f90
@@ -66,7 +66,7 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final)
       integer(I4B)                              :: itmax, idmax
       real(DP), dimension(:), allocatable       :: vals
       real(DP), dimension(1)                    :: val
-      real(DP), dimension(NDIM)                 :: rot0
+      real(DP), dimension(NDIM)                 :: rot0, Ip0, Lnow
       real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, Ltmp
 
       call param%nciu%open(param)
@@ -117,14 +117,27 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final)
             call check( nf90_get_var(param%nciu%ncid, param%nciu%radius_varid, val, start=[1,1], count=[1,1]) )
             cb%R0 = val(1) 
 
-            call check( nf90_get_var(param%nciu%ncid, param%nciu%rotx_varid, val, start=[1,1], count=[1,1]) )
-            rot0(1) = val(1)
-            call check( nf90_get_var(param%nciu%ncid, param%nciu%roty_varid, val, start=[1,1], count=[1,1]) )
-            rot0(2) = val(1)
-            call check( nf90_get_var(param%nciu%ncid, param%nciu%rotz_varid, val, start=[1,1], count=[1,1]) )
-            rot0(3) = val(1)
+            if (param%lrotation) then
 
-            cb%L0(:) = cb%Ip(3) * cb%GM0 * cb%R0**2 * rot0(:)
+               call check( nf90_get_var(param%nciu%ncid, param%nciu%rotx_varid, val, start=[1,1], count=[1,1]) )
+               rot0(1) = val(1)
+               call check( nf90_get_var(param%nciu%ncid, param%nciu%roty_varid, val, start=[1,1], count=[1,1]) )
+               rot0(2) = val(1)
+               call check( nf90_get_var(param%nciu%ncid, param%nciu%rotz_varid, val, start=[1,1], count=[1,1]) )
+               rot0(3) = val(1)
+
+               call check( nf90_get_var(param%nciu%ncid, param%nciu%Ip1_varid, val, start=[1,1], count=[1,1]) )
+               Ip0(1) = val(1)
+               call check( nf90_get_var(param%nciu%ncid, param%nciu%Ip2_varid, val, start=[1,1], count=[1,1]) )
+               Ip0(2) = val(1)
+               call check( nf90_get_var(param%nciu%ncid, param%nciu%Ip3_varid, val, start=[1,1], count=[1,1]) )
+               Ip0(3) = val(1)
+
+               cb%L0(:) = Ip0(3) * cb%GM0 * cb%R0**2 * rot0(:)
+
+               Lnow(:) = cb%Ip(3) * cb%Gmass * cb%radius**2 * cb%rot(:)
+               cb%dL(:) = Lnow(:) - cb%L0(:)
+            end if
          end select
 
       end if