Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Added nplm as a NetCDF variable for SyMBA runs
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Oct 11, 2022
1 parent 31a2e3a commit da76873
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 6 deletions.
4 changes: 2 additions & 2 deletions src/modules/swiftest_globals.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
30 changes: 26 additions & 4 deletions src/netcdf/netcdf.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(:)

Expand Down Expand Up @@ -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) )
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]) )
Expand Down Expand Up @@ -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]) )
Expand Down Expand Up @@ -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]) )
Expand Down

0 comments on commit da76873

Please sign in to comment.