diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index e012ac389..ae7364ba4 100644 --- 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