diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 7a66bf9fe..640503986 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -145,7 +145,7 @@ module swiftest_globals character(*), parameter :: NAME_VARNAME = "name" !! NetCDF name of the particle name variable character(*), parameter :: NPL_VARNAME = "npl" !! NetCDF name of the number of active massive bodies variable character(*), parameter :: NTP_VARNAME = "ntp" !! NetCDF name of the number of active test particles variable - character(*), parameter :: NPLM_VARNAME = "ntpm" !! NetCDF name of the number of active fully interacting massive bodies variable (SyMBA) + character(*), parameter :: NPLM_VARNAME = "nplm" !! NetCDF name of the number of active fully interacting massive bodies variable (SyMBA) character(*), parameter :: A_VARNAME = "a" !! NetCDF name of the semimajor axis variable character(*), parameter :: E_VARNAME = "e" !! NetCDF name of the eccentricity variable character(*), parameter :: INC_VARNAME = "inc" !! NetCDF name of the inclination variable @@ -222,7 +222,7 @@ module swiftest_globals integer(I4B) :: ptype_varid !! NetCDF ID for the particle type variable integer(I4B) :: npl_varid !! NetCDF ID for the number of active massive bodies variable integer(I4B) :: ntp_varid !! NetCDF ID for the number of active test particles variable - integer(I4B) :: nplm_varid !! NetCDF ID for the number of massive bodies variable (SyMBA) + integer(I4B) :: nplm_varid !! NetCDF ID for the number of active fully interacting massive bodies variable (SyMBA) integer(I4B) :: a_varid !! NetCDF ID for the semimajor axis variable integer(I4B) :: e_varid !! NetCDF ID for the eccentricity variable integer(I4B) :: inc_varid !! NetCDF ID for the inclination variable diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index ae7364ba4..d787bcaa6 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -107,9 +107,9 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) call check( nf90_get_var(param%nciu%ncid, param%nciu%L_spinz_varid, val, start=[1], count=[1]) ) self%Lspin_orig(3) = val(1) - call check( nf90_get_var(param%nciu%ncid, param%nciu%L_escapex_varid, self%Lescape(1), start=[1]) ) - call check( nf90_get_var(param%nciu%ncid, param%nciu%L_escapey_varid, self%Lescape(2), start=[1]) ) - call check( nf90_get_var(param%nciu%ncid, param%nciu%L_escapez_varid, self%Lescape(3), start=[1]) ) + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_escapex_varid, self%Lescape(1), start=[1]) ) + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_escapey_varid, self%Lescape(2), start=[1]) ) + call check( nf90_get_var(param%nciu%ncid, param%nciu%L_escapez_varid, self%Lescape(3), start=[1]) ) self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) @@ -202,6 +202,7 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_def_var(self%ncid, ID_DIMNAME, NF90_INT, self%id_dimid, self%id_varid) ) call check( nf90_def_var(self%ncid, NPL_VARNAME, NF90_INT, self%time_dimid, self%npl_varid) ) call check( nf90_def_var(self%ncid, NTP_VARNAME, NF90_INT, self%time_dimid, self%ntp_varid) ) + if (param%integrator == SYMBA) call check( nf90_def_var(self%ncid, NPLM_VARNAME, NF90_INT, self%time_dimid, self%nplm_varid) ) call check( nf90_def_var(self%ncid, NAME_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%name_varid) ) call check( nf90_def_var(self%ncid, PTYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%ptype_varid) ) call check( nf90_def_var(self%ncid, STATUS_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%status_varid) ) @@ -358,6 +359,7 @@ module subroutine netcdf_open(self, param, readonly) call check( nf90_inq_varid(self%ncid, ID_DIMNAME, self%id_varid)) call check( nf90_inq_varid(self%ncid, NPL_VARNAME, self%npl_varid)) call check( nf90_inq_varid(self%ncid, NTP_VARNAME, self%ntp_varid)) + if (param%integrator == SYMBA) call check( nf90_inq_varid(self%ncid, NPLM_VARNAME, self%nplm_varid)) call check( nf90_inq_varid(self%ncid, NAME_VARNAME, self%name_varid)) call check( nf90_inq_varid(self%ncid, PTYPE_VARNAME, self%ptype_varid)) call check( nf90_inq_varid(self%ncid, STATUS_VARNAME, self%status_varid)) @@ -476,7 +478,7 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) ! Return integer(I4B) :: ierr !! Error code: returns 0 if the read is successful ! Internals - integer(I4B) :: dim, i, j, tslot, idmax, npl_check, ntp_check, t_max, str_max + integer(I4B) :: dim, i, j, tslot, idmax, npl_check, ntp_check, nplm_check, t_max, str_max real(DP), dimension(:), allocatable :: rtemp integer(I4B), dimension(:), allocatable :: itemp logical, dimension(:), allocatable :: validmask, tpmask, plmask @@ -529,6 +531,18 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) call util_exit(failure) end if + select type (pl) + class is (symba_pl) + select type (param) + class is (symba_parameters) + nplm_check = count(rtemp(:) > param%GMTINY .and. plmask(:)) + if (nplm_check /= pl%nplm) then + write(*,*) "Error reading in NetCDF file: The recorded value of nplm does not match the number of active fully interacting massive bodies" + call util_exit(failure) + end if + end select + end select + ! Now read in each variable and split the outputs by body type if ((param%in_form == XV) .or. (param%in_form == XVEL)) then call check( nf90_get_var(iu%ncid, iu%xhx_varid, rtemp, start=[1, tslot]) ) @@ -721,6 +735,10 @@ module subroutine netcdf_read_hdr_system(self, iu, param) call check( nf90_get_var(iu%ncid, iu%time_varid, param%t, start=[tslot]) ) call check( nf90_get_var(iu%ncid, iu%npl_varid, self%pl%nbody, start=[tslot]) ) call check( nf90_get_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]) ) + select type(pl => self%pl) + class is (symba_pl) + call check( nf90_get_var(iu%ncid, iu%nplm_varid, pl%nplm, start=[tslot]) ) + end select if (param%lenergy) then call check( nf90_get_var(iu%ncid, iu%KE_orb_varid, self%ke_orbit, start=[tslot]) ) @@ -1185,6 +1203,10 @@ module subroutine netcdf_write_hdr_system(self, iu, param) call check( nf90_put_var(iu%ncid, iu%time_varid, param%t, start=[tslot]) ) call check( nf90_put_var(iu%ncid, iu%npl_varid, self%pl%nbody, start=[tslot]) ) call check( nf90_put_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]) ) + select type(pl => self%pl) + class is (symba_pl) + call check( nf90_put_var(iu%ncid, iu%nplm_varid, pl%nplm, start=[tslot]) ) + end select if (param%lenergy) then call check( nf90_put_var(iu%ncid, iu%KE_orb_varid, self%ke_orbit, start=[tslot]) )