diff --git a/src/io/io.f90 b/src/io/io.f90 index ebdd68531..7ec71d9ae 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -72,11 +72,20 @@ module subroutine io_conservation_report(self, param, lterminal) if (lterminal) write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror if (abs(Merror) > 100 * epsilon(Merror)) then write(*,*) "Severe error! Mass not conserved! Halting!" - call pl%xv2el(cb) - call self%write_hdr(param%nciu, param) - call cb%write_frame(param%nciu, param) - call pl%write_frame(param%nciu, param) - call param%nciu%close() + if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then + write(*,*) "Merror = ", Merror + write(*,*) "GMtot_now : ",GMtot_now + write(*,*) "GMtot_orig: ",system%GMtot_orig + write(*,*) "Difference: ",GMtot_now - system%GMtot_orig + else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + ! Save the frame of data to the bin file in the slot just after the present one for diagnostics + param%ioutput = param%ioutput + 1_I8B + call pl%xv2el(cb) + call self%write_hdr(param%nciu, param) + call cb%write_frame(param%nciu, param) + call pl%write_frame(param%nciu, param) + call param%nciu%close() + end if call util_exit(FAILURE) end if end if @@ -1351,7 +1360,6 @@ subroutine io_read_in_cb(self, param) close(iu, err = 667, iomsg = errmsg) if (ierr == 0) then - if (self%j2rp2 /= 0.0_DP) param%loblatecb = .true. if (param%rmin < 0.0) param%rmin = self%radius @@ -1361,7 +1369,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/kick/kick.f90 b/src/kick/kick.f90 index 891936f20..c3e37d927 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -189,7 +189,7 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius end do !$omp end parallel do else - !$omp parallel do default(private) schedule(static)& + !$omp parallel do default(private) schedule(static)& !$omp shared(npl, nplm, x, Gmass, radius) & !$omp lastprivate(rji2, xr, yr, zr) & !$omp reduction(+:ahi) & diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index b99188855..d28fcc788 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -66,6 +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, Ip0, Lnow real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, Ltmp call param%nciu%open(param) @@ -74,7 +75,8 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) allocate(vals(idmax)) call check( nf90_get_var(param%nciu%ncid, param%nciu%time_varid, val, start=[1], count=[1]) ) - old_t_final = val(1) + !old_t_final = val(1) + old_t_final = param%t0 ! For NetCDF it is safe to overwrite the final t value on a restart if (param%lenergy) then call check( nf90_get_var(param%nciu%ncid, param%nciu%KE_orb_varid, val, start=[1], count=[1]) ) @@ -107,6 +109,37 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) call check( nf90_get_var(param%nciu%ncid, param%nciu%Gmass_varid, vals, start=[1,1], count=[idmax,1]) ) self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + select type(cb => self%cb) + class is (symba_cb) + cb%GM0 = vals(1) + cb%dGM = cb%Gmass - cb%GM0 + + call check( nf90_get_var(param%nciu%ncid, param%nciu%radius_varid, val, start=[1,1], count=[1,1]) ) + cb%R0 = val(1) + + if (param%lrotation) then + + 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 deallocate(vals)