Fixedd calculation of initial angular momentum and change of angular …

…momentum of central body for a restarted run
MintonGroup · Nov 15, 2021 · 1a8d18f · 1a8d18f
1 parent 4260a05
commit 1a8d18f
Show file tree

Hide file tree

Showing 2 changed files with 25 additions and 9 deletions.
diff --git 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
@@ -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