From da3fe6ffb13264539523c57349003781e920c381 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 21 Dec 2022 11:17:12 -0500 Subject: [PATCH] Lots more restructuring, cleaning, refactoring, etc. Rearranged netcdf stuff to make it a bit cleaner --- src/CMakeLists.txt | 21 +- src/base/base_module.f90 | 212 +-- src/collision/collision_io.f90 | 118 +- src/collision/collision_module.f90 | 6 +- src/encounter/encounter_check.f90 | 6 +- src/encounter/encounter_io.f90 | 92 +- src/encounter/encounter_module.f90 | 37 +- src/netcdf_io/netcdf_io_implementations.f90 | 66 + src/netcdf_io/netcdf_io_module.f90 | 159 +++ src/swiftest/swiftest_io.f90 | 1278 ++++++++++++++++++- src/swiftest/swiftest_io_netcdf.f90 | 1246 ------------------ src/swiftest/swiftest_kick.f90 | 2 +- src/swiftest/swiftest_module.f90 | 200 +-- src/swiftest/swiftest_setup.f90 | 19 +- src/swiftest/swiftest_util.f90 | 14 +- src/symba/symba_kick.f90 | 2 +- src/symba/symba_util.f90 | 9 +- 17 files changed, 1756 insertions(+), 1731 deletions(-) create mode 100644 src/netcdf_io/netcdf_io_implementations.f90 create mode 100644 src/netcdf_io/netcdf_io_module.f90 delete mode 100644 src/swiftest/swiftest_io_netcdf.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 33a8f5536..b43cf0407 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,6 +24,7 @@ SET(STRICT_MATH_FILES SET(FAST_MATH_FILES ${SRC}/globals/globals_module.f90 ${SRC}/base/base_module.f90 + ${SRC}/netcdf_io/netcdf_io_module.f90 ${SRC}/misc/lambda_function_module.f90 ${SRC}/misc/io_progress_bar_module.f90 ${SRC}/encounter/encounter_module.f90 @@ -36,16 +37,6 @@ SET(FAST_MATH_FILES ${SRC}/rmvs/rmvs_module.f90 ${SRC}/helio/helio_module.f90 ${SRC}/symba/symba_module.f90 - ${SRC}/swiftest/swiftest_discard.f90 - ${SRC}/swiftest/swiftest_io.f90 - ${SRC}/swiftest/swiftest_obl.f90 - ${SRC}/swiftest/swiftest_util.f90 - ${SRC}/swiftest/swiftest_drift.f90 - ${SRC}/swiftest/swiftest_io_netcdf.f90 - ${SRC}/swiftest/swiftest_orbel.f90 - ${SRC}/swiftest/swiftest_gr.f90 - ${SRC}/swiftest/swiftest_kick.f90 - ${SRC}/swiftest/swiftest_setup.f90 ${SRC}/collision/collision_check.f90 ${SRC}/collision/collision_regime.f90 ${SRC}/collision/collision_setup.f90 @@ -67,6 +58,7 @@ SET(FAST_MATH_FILES ${SRC}/helio/helio_setup.f90 ${SRC}/helio/helio_step.f90 ${SRC}/helio/helio_util.f90 + ${SRC}/netcdf_io/netcdf_io_implementations.f90 ${SRC}/operator/operator_cross.f90 ${SRC}/operator/operator_mag.f90 ${SRC}/operator/operator_unit.f90 @@ -75,6 +67,15 @@ SET(FAST_MATH_FILES ${SRC}/rmvs/rmvs_setup.f90 ${SRC}/rmvs/rmvs_step.f90 ${SRC}/rmvs/rmvs_util.f90 + ${SRC}/swiftest/swiftest_discard.f90 + ${SRC}/swiftest/swiftest_io.f90 + ${SRC}/swiftest/swiftest_obl.f90 + ${SRC}/swiftest/swiftest_util.f90 + ${SRC}/swiftest/swiftest_drift.f90 + ${SRC}/swiftest/swiftest_orbel.f90 + ${SRC}/swiftest/swiftest_gr.f90 + ${SRC}/swiftest/swiftest_kick.f90 + ${SRC}/swiftest/swiftest_setup.f90 ${SRC}/symba/symba_discard.f90 ${SRC}/symba/symba_drift.f90 ${SRC}/symba/symba_encounter_check.f90 diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index 96e529c9d..7e19be601 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -57,6 +57,7 @@ module base real(DP) :: GMTINY = -1.0_DP !! Smallest G*mass that is fully gravitating real(DP) :: min_GMfrag = -1.0_DP !! Smallest G*mass that can be produced in a fragmentation event integer(I4B), dimension(:), allocatable :: seed !! Random seeds for fragmentation modeling + logical :: lmtiny_pl = .false. !! Include semi-interacting massive bodies logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. character(STRMAX) :: encounter_save = "NONE" !! Indicate if and how encounter data should be saved logical :: lenc_save_trajectory = .false. !! Indicates that when encounters are saved, the full trajectory through recursion steps are saved @@ -152,137 +153,6 @@ subroutine abstract_io_read_in_param(self, param_file_name) end subroutine abstract_io_read_in_param end interface - !! This derived datatype stores the NetCDF ID values for each of the variables included in the NetCDF data file. This is used as the base class defined in base - type, abstract :: base_io_netcdf_parameters - character(STRMAX) :: file_name !! Name of the output file - integer(I4B) :: out_type !! output type (will be assigned either NF90_DOUBLE or NF90_FLOAT, depending on the user parameter) - integer(I4B) :: id !! ID for the output file - integer(I4B) :: discard_body_id_varid !! ID for the id of the other body involved in the discard - integer(I4B) :: id_chunk !! Chunk size for the id dimension variables - integer(I4B) :: time_chunk !! Chunk size for the time dimension variables - logical :: lpseudo_vel_exists = .false. !! Logical flag to indicate whether or not the pseudovelocity vectors were present in an old file. - - ! Dimension ids and variable names - character(NAMELEN) :: str_dimname = "string32" !! name of the character string dimension - integer(I4B) :: str_dimid !! ID for the character string dimension - character(NAMELEN) :: time_dimname = "time" !! name of the time dimension - integer(I4B) :: time_dimid !! ID for the time dimension - integer(I4B) :: time_varid !! ID for the time variable - character(NAMELEN) :: name_dimname = "name" !! name of the particle name dimension - integer(I4B) :: name_dimid !! ID for the particle name dimension - integer(I4B) :: name_varid !! ID for the particle name variable - character(NAMELEN) :: space_dimname = "space" !! name of the space dimension - integer(I4B) :: space_dimid !! ID for the space dimension - integer(I4B) :: space_varid !! ID for the space variable - character(len=1), dimension(3) :: space_coords = ["x","y","z"] !! The space dimension coordinate labels - - ! Non-dimension ids and variable names - character(NAMELEN) :: ptype_varname = "particle_type" !! name of the particle type variable - integer(I4B) :: ptype_varid !! ID for the particle type variable - character(NAMELEN) :: id_varname = "id" !! name of the particle id variable - integer(I4B) :: id_varid !! ID for the id variable - character(NAMELEN) :: npl_varname = "npl" !! name of the number of active massive bodies variable - integer(I4B) :: npl_varid !! ID for the number of active massive bodies variable - character(NAMELEN) :: ntp_varname = "ntp" !! name of the number of active test particles variable - integer(I4B) :: ntp_varid !! ID for the number of active test particles variable - character(NAMELEN) :: nplm_varname = "nplm" !! name of the number of active fully interacting massive bodies variable (SyMBA) - integer(I4B) :: nplm_varid !! ID for the number of active fully interacting massive bodies variable (SyMBA) - character(NAMELEN) :: a_varname = "a" !! name of the semimajor axis variable - integer(I4B) :: a_varid !! ID for the semimajor axis variable - character(NAMELEN) :: e_varname = "e" !! name of the eccentricity variable - integer(I4B) :: e_varid !! ID for the eccentricity variable - character(NAMELEN) :: inc_varname = "inc" !! name of the inclination variable - integer(I4B) :: inc_varid !! ID for the inclination variable - character(NAMELEN) :: capom_varname = "capom" !! name of the long. asc. node variable - integer(I4B) :: capom_varid !! ID for the long. asc. node variable - character(NAMELEN) :: omega_varname = "omega" !! name of the arg. of periapsis variable - integer(I4B) :: omega_varid !! ID for the arg. of periapsis variable - character(NAMELEN) :: capm_varname = "capm" !! name of the mean anomaly variable - integer(I4B) :: capm_varid !! ID for the mean anomaly variable - character(NAMELEN) :: varpi_varname = "varpi" !! name of the long. of periapsis variable - integer(I4B) :: varpi_varid !! ID for the long. of periapsis variable - character(NAMELEN) :: lam_varname = "lam" !! name of the mean longitude variable - integer(I4B) :: lam_varid !! ID for the mean longitude variable - character(NAMELEN) :: f_varname = "f" !! name of the true anomaly variable - integer(I4B) :: f_varid !! ID for the true anomaly variable - character(NAMELEN) :: cape_varname = "cape" !! name of the eccentric anomaly variable - integer(I4B) :: cape_varid !! ID for the eccentric anomaly variable - character(NAMELEN) :: rh_varname = "rh" !! name of the heliocentric position vector variable - integer(I4B) :: rh_varid !! ID for the heliocentric position vector variable - character(NAMELEN) :: vh_varname = "vh" !! name of the heliocentric velocity vector variable - integer(I4B) :: vh_varid !! ID for the heliocentric velocity vector variable - character(NAMELEN) :: gr_pseudo_vh_varname = "gr_pseudo_vh" !! name of the heliocentric pseudovelocity vector variable (used in GR only) - integer(I4B) :: gr_pseudo_vh_varid !! ID for the heliocentric pseudovelocity vector variable (used in GR) - character(NAMELEN) :: gmass_varname = "Gmass" !! name of the mass variable - integer(I4B) :: Gmass_varid !! ID for the mass variable - character(NAMELEN) :: rhill_varname = "rhill" !! name of the hill radius variable - integer(I4B) :: rhill_varid !! ID for the hill radius variable - character(NAMELEN) :: radius_varname = "radius" !! name of the radius variable - integer(I4B) :: radius_varid !! ID for the radius variable - character(NAMELEN) :: Ip_varname = "Ip" !! name of the principal moment of inertial variable - integer(I4B) :: Ip_varid !! ID for the axis principal moment of inertia variable - character(NAMELEN) :: rot_varname = "rot" !! name of the rotation vector variable - integer(I4B) :: rot_varid !! ID for the rotation vector variable - character(NAMELEN) :: j2rp2_varname = "j2rp2" !! name of the j2rp2 variable - integer(I4B) :: j2rp2_varid !! ID for the j2 variable - character(NAMELEN) :: j4rp4_varname = "j4rp4" !! name of the j4pr4 variable - integer(I4B) :: j4rp4_varid !! ID for the j4 variable - character(NAMELEN) :: k2_varname = "k2" !! name of the Love number variable - integer(I4B) :: k2_varid !! ID for the Love number variable - character(NAMELEN) :: q_varname = "Q" !! name of the energy dissipation variable - integer(I4B) :: Q_varid !! ID for the energy dissipation variable - character(NAMELEN) :: ke_orb_varname = "KE_orb" !! name of the system orbital kinetic energy variable - integer(I4B) :: KE_orb_varid !! ID for the system orbital kinetic energy variable - character(NAMELEN) :: ke_spin_varname = "KE_spin" !! name of the system spin kinetic energy variable - integer(I4B) :: KE_spin_varid !! ID for the system spin kinetic energy variable - character(NAMELEN) :: pe_varname = "PE" !! name of the system potential energy variable - integer(I4B) :: PE_varid !! ID for the system potential energy variable - character(NAMELEN) :: L_orb_varname = "L_orb" !! name of the orbital angular momentum vector variable - integer(I4B) :: L_orb_varid !! ID for the system orbital angular momentum vector variable - character(NAMELEN) :: Lspin_varname = "Lspin" !! name of the spin angular momentum vector variable - integer(I4B) :: Lspin_varid !! ID for the system spin angular momentum vector variable - character(NAMELEN) :: L_escape_varname = "L_escape" !! name of the escaped angular momentum vector variable - integer(I4B) :: L_escape_varid !! ID for the escaped angular momentum vector variable - character(NAMELEN) :: Ecollisions_varname = "Ecollisions" !! name of the escaped angular momentum y variable - integer(I4B) :: Ecollisions_varid !! ID for the energy lost in collisions variable - character(NAMELEN) :: Euntracked_varname = "Euntracked" !! name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) - integer(I4B) :: Euntracked_varid !! ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) - character(NAMELEN) :: GMescape_varname = "GMescape" !! name of the G*Mass of bodies that escape the system - integer(I4B) :: GMescape_varid !! ID for the G*Mass of bodies that escape the system - character(NAMELEN) :: origin_type_varname = "origin_type" !! name of the origin type variable (Initial Conditions, Disruption, etc.) - integer(I4B) :: origin_type_varid !! ID for the origin type - character(NAMELEN) :: origin_time_varname = "origin_time" !! name of the time of origin variable - integer(I4B) :: origin_time_varid !! ID for the origin time - character(NAMELEN) :: collision_id_varname = "collision_id" !! name of the collision id variable - integer(I4B) :: collision_id_varid !! Netcdf ID for the origin collision ID - character(NAMELEN) :: origin_rh_varname = "origin_rh" !! name of the heliocentric position vector of the body at the time of origin variable - integer(I4B) :: origin_rh_varid !! ID for the origin position vector variable - character(NAMELEN) :: origin_vh_varname = "origin_vh" !! name of the heliocentric velocity vector of the body at the time of origin variable - integer(I4B) :: origin_vh_varid !! ID for the origin velocity vector component - character(NAMELEN) :: discard_time_varname = "discard_time" !! name of the time of discard variable - integer(I4B) :: discard_time_varid !! ID for the time of discard variable - character(NAMELEN) :: discard_rh_varname = "discard_rh" !! name of the heliocentric position vector of the body at the time of discard variable - integer(I4B) :: discard_rh_varid !! ID for the heliocentric position vector of the body at the time of discard variable - character(NAMELEN) :: discard_vh_varname = "discard_vh" !! name of the heliocentric velocity vector of the body at the time of discard variable - integer(I4B) :: discard_vh_varid !! ID for the heliocentric velocity vector of the body at the time of discard variable - character(NAMELEN) :: discard_body_id_varname = "discard_body_id" !! name of the id of the other body involved in the discard - contains - procedure(abstract_io_netcdf_initialize_output), deferred :: initialize - procedure(abstract_io_netcdf_open), deferred :: open - procedure :: close => base_io_netcdf_close !! Closes an open NetCDF file - procedure :: flush => base_io_netcdf_flush !! Flushes the current buffer to disk by closing and re-opening the file. - procedure :: sync => base_io_netcdf_sync !! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) - end type base_io_netcdf_parameters - - abstract interface - subroutine abstract_io_netcdf_initialize_output(self, param) - import base_io_netcdf_parameters, base_parameters - implicit none - class(base_io_netcdf_parameters), intent(inout) :: self !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine abstract_io_netcdf_initialize_output - end interface - type :: base_storage_frame class(*), allocatable :: item @@ -306,7 +176,6 @@ end subroutine abstract_io_netcdf_initialize_output integer(I4B) :: nt !! Number of unique time values in all saved snapshots real(DP), dimension(:), allocatable :: tvals !! The set of unique time values contained in the snapshots integer(I4B), dimension(:), allocatable :: tmap !! The t value -> index map - class(base_io_netcdf_parameters), allocatable :: nc !! NetCDF object attached to this storage object contains procedure :: reset => reset_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 end type base_storage @@ -338,16 +207,6 @@ end subroutine abstract_io_netcdf_initialize_output type, abstract :: base_nbody_system end type base_nbody_system - abstract interface - subroutine abstract_io_netcdf_open(self, param, readonly) - import base_io_netcdf_parameters, base_parameters - implicit none - class(base_io_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(in) :: param !! Current run configuration parameters - logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only - end subroutine abstract_io_netcdf_open - end interface - contains subroutine copy_store(self, source) @@ -377,75 +236,6 @@ subroutine final_storage_frame(self) return end subroutine final_storage_frame - subroutine netcdf_check(status, call_identifier) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Checks the status of all NetCDF operations to catch errors - use netcdf - implicit none - ! Arguments - integer, intent (in) :: status !! The status code returned by a NetCDF function - character(len=*), intent(in), optional :: call_identifier !! String that indicates which calling function caused the error for diagnostic purposes - - if(status /= nf90_noerr) then - if (present(call_identifier)) write(*,*) "NetCDF error in ",trim(call_identifier) - write(*,*) trim(nf90_strerror(status)) - call swiftest_util_exit(FAILURE) - end if - - return - end subroutine netcdf_check - - - subroutine base_io_netcdf_close(self) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Closes a NetCDF file - use netcdf - implicit none - ! Arguments - class(base_io_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - - call netcdf_check( nf90_close(self%id), "base_io_netcdf_close" ) - - return - end subroutine base_io_netcdf_close - - - subroutine base_io_netcdf_flush(self, param) - !! author: David A. Minton - !! - !! Flushes the current buffer to disk by closing and re-opening the file. - !! - implicit none - ! Arguments - class(base_io_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - - call self%close() - call self%open(param) - - return - end subroutine base_io_netcdf_flush - - - - subroutine base_io_netcdf_sync(self) - !! author: David A. Minton - !! - !! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) - !! - use netcdf - implicit none - ! Arguments - class(base_io_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - - call netcdf_check( nf90_sync(self%id), "base_io_netcdf_sync nf90_sync" ) - - return - end subroutine base_io_netcdf_sync - - subroutine base_util_final_storage(self) !! author: David A. Minton diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 864fa93f9..963d86463 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -24,7 +24,7 @@ module subroutine collision_io_dump(self, param) integer(I4B) :: i select type(nc => self%nc) - class is (collision_io_parameters) + class is (collision_netcdf_parameters) select type(param) class is (swiftest_parameters) if (self%iframe > 0) then @@ -65,7 +65,7 @@ module subroutine collision_io_initialize_output(self, param) use netcdf implicit none ! Arguments - class(collision_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(collision_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset class(base_parameters), intent(in) :: param ! Internals integer(I4B) :: nvar, varid, vartype @@ -96,93 +96,93 @@ module subroutine collision_io_initialize_output(self, param) close(unit=LUN, status="delete") end if - call netcdf_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "collision_io_initialize_output nf90_create" ) + call netcdf_io_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "collision_io_initialize_output nf90_create" ) ! Dimensions - call netcdf_check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "collision_io_initialize_output nf90_def_dim event_dimid" ) ! Dimension to store individual collision events - call netcdf_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "collision_io_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call netcdf_check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "collision_io_initialize_output nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers - call netcdf_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "collision_io_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - call netcdf_check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "collision_io_initialize_output nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" + call netcdf_io_check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "collision_io_initialize_output nf90_def_dim event_dimid" ) ! Dimension to store individual collision events + call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "collision_io_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "collision_io_initialize_output nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers + call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "collision_io_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + call netcdf_io_check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "collision_io_initialize_output nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" ! Dimension coordinates - call netcdf_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_initialize_output nf90_def_var space_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "collision_io_initialize_output nf90_def_var name_varid") - call netcdf_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "collision_io_initialize_output nf90_def_var stage_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_initialize_output nf90_def_var space_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "collision_io_initialize_output nf90_def_var name_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "collision_io_initialize_output nf90_def_var stage_varid" ) ! Variables - call netcdf_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "collision_io_initialize_output nf90_def_var id_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & + call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "collision_io_initialize_output nf90_def_var id_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & nc%event_dimid, nc%time_varid), "collision_io_initialize_output nf90_def_var time_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & + call netcdf_io_check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & [nc%str_dimid, nc%event_dimid], nc%regime_varid), "collision_io_initialize_output nf90_def_var regime_varid") - call netcdf_check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & + call netcdf_io_check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & [ nc%event_dimid], nc%Qloss_varid), "collision_io_initialize_output nf90_def_var Qloss_varid") - call netcdf_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & + call netcdf_io_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%ptype_varid), "collision_io_initialize_output nf90_def_var ptype_varid") - call netcdf_check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, & + call netcdf_io_check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, & [ nc%event_dimid], nc%loop_varid), "collision_io_initialize_output nf90_def_var loop_varid") - call netcdf_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rh_varid), "collision_io_initialize_output nf90_def_var rh_varid") - call netcdf_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%vh_varid), "collision_io_initialize_output nf90_def_var vh_varid") - call netcdf_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Gmass_varid), "collision_io_initialize_output nf90_def_var Gmass_varid") - call netcdf_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%radius_varid), "collision_io_initialize_output nf90_def_var radius_varid") - call netcdf_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Ip_varid), "collision_io_initialize_output nf90_def_var Ip_varid") - call netcdf_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rot_varid), "collision_io_initialize_output nf90_def_var rot_varid") if (param%lenergy) then - call netcdf_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "collision_io_initialize_output nf90_def_var KE_orb_varid") - call netcdf_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "collision_io_initialize_output nf90_def_var KE_spin_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "collision_io_initialize_output nf90_def_var PE_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & + call netcdf_io_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "collision_io_initialize_output nf90_def_var L_orb_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type,& + call netcdf_io_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type,& [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%Lspin_varid), "collision_io_initialize_output nf90_def_var Lspin_varid" ) end if - call netcdf_check( nf90_inquire(nc%id, nVariables=nvar), "collision_io_initialize_output nf90_inquire nVariables" ) + call netcdf_io_check( nf90_inquire(nc%id, nVariables=nvar), "collision_io_initialize_output nf90_inquire nVariables" ) do varid = 1, nvar - call netcdf_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "collision_io_initialize_output nf90_inquire_variable" ) + call netcdf_io_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "collision_io_initialize_output nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "collision_io_initialize_output nf90_def_var_fill NF90_INT" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "collision_io_initialize_output nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "collision_io_initialize_output nf90_def_var_fill NF90_FLOAT" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "collision_io_initialize_output nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "collision_io_initialize_output nf90_def_var_fill NF90_DOUBLE" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "collision_io_initialize_output nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "collision_io_initialize_output nf90_def_var_fill NF90_CHAR" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "collision_io_initialize_output nf90_def_var_fill NF90_CHAR" ) end select end do ! Take the file out of define mode - call netcdf_check( nf90_enddef(nc%id), "collision_io_initialize_output nf90_enddef" ) + call netcdf_io_check( nf90_enddef(nc%id), "collision_io_initialize_output nf90_enddef" ) ! Add in the space and stage dimension coordinates - call netcdf_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "collision_io_initialize_output nf90_put_var space" ) - call netcdf_check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(1), start=[1,1], count=[len(nc%stage_coords(1)),1]), "collision_io_initialize_output nf90_put_var stage 1" ) - call netcdf_check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(2), start=[1,2], count=[len(nc%stage_coords(2)),1]), "collision_io_initialize_output nf90_put_var stage 2" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "collision_io_initialize_output nf90_put_var space" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(1), start=[1,1], count=[len(nc%stage_coords(1)),1]), "collision_io_initialize_output nf90_put_var stage 1" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(2), start=[1,2], count=[len(nc%stage_coords(2)),1]), "collision_io_initialize_output nf90_put_var stage 2" ) end associate end select @@ -211,16 +211,16 @@ module subroutine collision_io_write_frame_snapshot(self, history, param) class(swiftest_pl), allocatable :: pl select type(nc => history%nc) - class is (collision_io_parameters) + class is (collision_netcdf_parameters) associate(system => self%collision_system, impactors => self%collision_system%impactors, fragments => self%collision_system%fragments, eslot => param%ioutput) - call netcdf_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "collision_io_write_frame_snapshot nf90_set_fill" ) + call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "collision_io_write_frame_snapshot nf90_set_fill" ) - call netcdf_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "collision_io_write_frame_snapshot nf90_put_var time_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[eslot]), "collision_io_write_frame_snapshot nf90_put_varloop_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "collision_io_write_frame_snapshot nf90_put_var time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[eslot]), "collision_io_write_frame_snapshot nf90_put_varloop_varid" ) charstring = trim(adjustl(REGIME_NAMES(impactors%regime))) - call netcdf_check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], count=[len(charstring), 1]), "collision_io_write_frame_snapshot nf90_put_var regime_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%Qloss_varid, impactors%Qloss, start=[eslot] ), "collision_io_write_frame_snapshot nf90_put_var Qloss_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], count=[len(charstring), 1]), "collision_io_write_frame_snapshot nf90_put_var regime_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Qloss_varid, impactors%Qloss, start=[eslot] ), "collision_io_write_frame_snapshot nf90_put_var Qloss_varid" ) select type(before =>self%collision_system%before) class is (swiftest_nbody_system) @@ -237,30 +237,30 @@ module subroutine collision_io_write_frame_snapshot(self, history, param) npl = pl%nbody do i = 1, npl idslot = findloc(history%idvals,pl%id(i),dim=1) - call netcdf_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), "collision_io_write_frame_snapshot nf90_put_var id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), "collision_io_write_frame_snapshot nf90_put_var id_varid" ) charstring = trim(adjustl(pl%info(i)%name)) - call netcdf_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], count=[len(charstring), 1]), "collision_io_write_frame_snapshot nf90_put_var name_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], count=[len(charstring), 1]), "collision_io_write_frame_snapshot nf90_put_var name_varid" ) charstring = trim(adjustl(pl%info(i)%particle_type)) - call netcdf_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], count=[len(charstring), 1, 1]), "collision_io_write_frame_snapshot nf90_put_var particle_type_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var rh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var vh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), "collision_io_write_frame_snapshot nf90_put_var Gmass_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), "collision_io_write_frame_snapshot nf90_put_var radius_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var Ip_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var rotx_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], count=[len(charstring), 1, 1]), "collision_io_write_frame_snapshot nf90_put_var particle_type_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var rh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var vh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), "collision_io_write_frame_snapshot nf90_put_var Gmass_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), "collision_io_write_frame_snapshot nf90_put_var radius_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var Ip_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var rotx_varid" ) end do end do end select end select if (param%lenergy) then - call netcdf_check( nf90_put_var(nc%id, nc%ke_orb_varid, system%ke_orbit(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var ke_orb_varid before" ) - call netcdf_check( nf90_put_var(nc%id, nc%ke_spin_varid, system%ke_spin(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var ke_spin_varid before" ) - call netcdf_check( nf90_put_var(nc%id, nc%pe_varid, system%pe(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var pe_varid before" ) - call netcdf_check( nf90_put_var(nc%id, nc%L_orb_varid, system%Lorbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_write_frame_snapshot nf90_put_var L_orb_varid before" ) - call netcdf_check( nf90_put_var(nc%id, nc%Lspin_varid, system%Lspin(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_write_frame_snapshot nf90_put_var Lspin_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ke_orb_varid, system%ke_orbit(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var ke_orb_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ke_spin_varid, system%ke_spin(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var ke_spin_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%pe_varid, system%pe(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var pe_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_orb_varid, system%Lorbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_write_frame_snapshot nf90_put_var L_orb_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Lspin_varid, system%Lspin(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_write_frame_snapshot nf90_put_var Lspin_varid before" ) end if - call netcdf_check( nf90_set_fill(nc%id, old_mode, old_mode) ) + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) ) end associate end select return diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index a4e6cd9bb..ccde318cc 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -137,7 +137,7 @@ module collision !! NetCDF dimension and variable names for the enounter save object - type, extends(encounter_io_parameters) :: collision_io_parameters + type, extends(encounter_netcdf_parameters) :: collision_netcdf_parameters integer(I4B) :: stage_dimid !! ID for the stage dimension integer(I4B) :: stage_varid !! ID for the stage variable character(NAMELEN) :: stage_dimname = "stage" !! name of the stage dimension (before/after) @@ -154,7 +154,7 @@ module collision integer(I4B) :: regime_varid !! ID for the collision regime variable contains procedure :: initialize => collision_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object - end type collision_io_parameters + end type collision_netcdf_parameters type, extends(encounter_snapshot) :: collision_snapshot @@ -205,7 +205,7 @@ end subroutine collision_io_dump module subroutine collision_io_initialize_output(self, param) implicit none - class(collision_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(collision_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset class(base_parameters), intent(in) :: param !! Current run configuration parameters end subroutine collision_io_initialize_output diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index cccb0bfe9..de81b75bc 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -43,7 +43,7 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, ind ! lfirst = .false. ! itimer%step_counter = INTERACTION_TIMER_CADENCE ! else - ! if (itimer%io_netcdf_check(param, nplpl)) call itimer%time_this_loop(param, nplpl) + ! if (itimer%netcdf_io_check(param, nplpl)) call itimer%time_this_loop(param, nplpl) ! end if ! else ! param%lencounter_sas_plpl = .false. @@ -117,7 +117,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, ! lfirst = .false. ! itimer%step_counter = INTERACTION_TIMER_CADENCE ! else - ! if (itimer%io_netcdf_check(param, nplplm)) call itimer%time_this_loop(param, nplplm) + ! if (itimer%netcdf_io_check(param, nplplm)) call itimer%time_this_loop(param, nplplm) ! end if ! else ! param%lencounter_sas_plpl = .false. @@ -214,7 +214,7 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, ! if (lsecond) then ! This ensures that the encounter check methods are run at least once prior to timing. Sort and sweep improves on the second pass due to the bounding box extents needing to be nearly sorted ! call itimer%time_this_loop(param, npltp) ! lsecond = .false. - ! else if (itimer%io_netcdf_check(param, npltp)) then + ! else if (itimer%netcdf_io_check(param, npltp)) then ! lsecond = .true. ! itimer%is_on = .false. ! end if diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index df82d567e..1529dd955 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -23,7 +23,7 @@ module subroutine encounter_io_dump(self, param) integer(I4B) :: i select type(nc => self%nc) - class is (encounter_io_parameters) + class is (encounter_netcdf_parameters) if (self%iframe > 0) then ! Create and save the output files for this encounter and fragmentation nc%file_number = nc%file_number + 1 @@ -62,7 +62,7 @@ module subroutine encounter_io_initialize_output(self, param) use netcdf implicit none ! Arguments - class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(encounter_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset class(base_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: nvar, varid, vartype @@ -91,54 +91,54 @@ module subroutine encounter_io_initialize_output(self, param) close(unit=LUN, status="delete") end if - call netcdf_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "encounter_io_initialize_output nf90_create" ) + call netcdf_io_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "encounter_io_initialize_output nf90_create" ) ! Dimensions - call netcdf_check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension - call netcdf_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "encounter_io_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call netcdf_check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "encounter_io_initialize_output nf90_def_dim name_dimid" ) ! dimension to store particle id numbers - call netcdf_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "encounter_io_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + call netcdf_io_check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "encounter_io_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "encounter_io_initialize_output nf90_def_dim name_dimid" ) ! dimension to store particle id numbers + call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "encounter_io_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) ! Dimension coordinates - call netcdf_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize_output nf90_def_var time_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_io_initialize_output nf90_def_var space_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "encounter_io_initialize_output nf90_def_var id_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize_output nf90_def_var time_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_io_initialize_output nf90_def_var space_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "encounter_io_initialize_output nf90_def_var id_varid" ) ! Variables - call netcdf_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "encounter_io_initialize_output nf90_def_var id_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%ptype_varid), "encounter_io_initialize_output nf90_def_var ptype_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_initialize_output nf90_def_var rh_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%vh_varid), "encounter_io_initialize_output nf90_def_var vh_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize_output nf90_def_var Gmass_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, [nc%time_dimid], nc%loop_varid), "encounter_io_initialize_output nf90_def_var loop_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "encounter_io_initialize_output nf90_def_var id_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%ptype_varid), "encounter_io_initialize_output nf90_def_var ptype_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_initialize_output nf90_def_var rh_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%vh_varid), "encounter_io_initialize_output nf90_def_var vh_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize_output nf90_def_var Gmass_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, [nc%time_dimid], nc%loop_varid), "encounter_io_initialize_output nf90_def_var loop_varid" ) if (param%lclose) then - call netcdf_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize_output nf90_def_var radius_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize_output nf90_def_var radius_varid" ) end if if (param%lrotation) then - call netcdf_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%Ip_varid), "encounter_io_initialize_output nf90_def_var Ip_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rot_varid), "encounter_io_initialize_output nf90_def_var rot_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%Ip_varid), "encounter_io_initialize_output nf90_def_var Ip_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rot_varid), "encounter_io_initialize_output nf90_def_var rot_varid" ) end if - call netcdf_check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize_output nf90_inquire nVariables" ) + call netcdf_io_check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize_output nf90_inquire nVariables" ) do varid = 1, nvar - call netcdf_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize_output nf90_inquire_variable" ) + call netcdf_io_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize_output nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "encounter_io_initialize_output nf90_def_var_fill NF90_INT" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "encounter_io_initialize_output nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "encounter_io_initialize_output nf90_def_var_fill NF90_FLOAT" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "encounter_io_initialize_output nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "encounter_io_initialize_output nf90_def_var_fill NF90_DOUBLE" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "encounter_io_initialize_output nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "encounter_io_initialize_output nf90_def_var_fill NF90_CHAR" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "encounter_io_initialize_output nf90_def_var_fill NF90_CHAR" ) end select end do ! Take the file out of define mode - call netcdf_check( nf90_enddef(nc%id), "encounter_io_initialize_output nf90_enddef" ) + call netcdf_io_check( nf90_enddef(nc%id), "encounter_io_initialize_output nf90_enddef" ) ! Add in the space dimension coordinates - call netcdf_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize_output nf90_put_var space" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize_output nf90_put_var space" ) end associate @@ -172,48 +172,48 @@ module subroutine encounter_io_write_frame_snapshot(self, history, param) select type(tp => self%tp) class is (swiftest_pl) select type (nc => history%nc) - class is (encounter_io_parameters) + class is (encounter_netcdf_parameters) associate(tslot => param%ioutput) - call netcdf_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame_snapshot nf90_set_fill" ) + call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame_snapshot nf90_set_fill" ) - call netcdf_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame_snapshot nf90_put_var time_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl loop_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame_snapshot nf90_put_var time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl loop_varid" ) npl = pl%nbody do i = 1, npl idslot = findloc(history%idvals,pl%id(i),dim=1) - call netcdf_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame_snapshot nf90_put_var pl id_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl rh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl vh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl Gmass_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame_snapshot nf90_put_var pl id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl rh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl vh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl Gmass_varid" ) - if (param%lclose) call netcdf_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl radius_varid" ) + if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl radius_varid" ) if (param%lrotation) then - call netcdf_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl Ip_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl rotx_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl Ip_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl rotx_varid" ) end if charstring = trim(adjustl(pl%info(i)%name)) - call netcdf_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var pl name_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var pl name_varid" ) charstring = trim(adjustl(pl%info(i)%particle_type)) - call netcdf_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var pl particle_type_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var pl particle_type_varid" ) end do ntp = tp%nbody do i = 1, ntp idslot = findloc(history%idvals,tp%id(i),dim=1) - call netcdf_check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame_snapshot nf90_put_var tp id_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var tp rh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var tp vh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame_snapshot nf90_put_var tp id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var tp rh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var tp vh_varid" ) charstring = trim(adjustl(tp%info(i)%name)) - call netcdf_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var tp name_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var tp name_varid" ) charstring = trim(adjustl(tp%info(i)%particle_type)) - call netcdf_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var tp particle_type_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var tp particle_type_varid" ) end do - call netcdf_check( nf90_set_fill(nc%id, old_mode, old_mode) ) + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) ) end associate end select end select diff --git a/src/encounter/encounter_module.f90 b/src/encounter/encounter_module.f90 index 04e7170c6..5c32a4c97 100644 --- a/src/encounter/encounter_module.f90 +++ b/src/encounter/encounter_module.f90 @@ -13,6 +13,7 @@ module encounter !! Definition of classes and methods used to determine close encounters use globals use base + use netcdf_io implicit none public @@ -57,8 +58,22 @@ module encounter final :: encounter_util_final_snapshot end type encounter_snapshot + !> NetCDF dimension and variable names for the enounter save object + type, extends(netcdf_parameters) :: encounter_netcdf_parameters + character(NAMELEN) :: loop_varname = "loopnum" !! Loop number for encounter + integer(I4B) :: loop_varid !! ID for the recursion level variable + integer(I4B) :: time_dimsize = 0 !! Number of time values in snapshot + integer(I4B) :: name_dimsize = 0 !! Number of potential id values in snapshot + integer(I4B) :: file_number = 1 !! The number to append on the output file + contains + procedure :: initialize => encounter_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object + procedure :: open => encounter_netcdf_io_open + end type encounter_netcdf_parameters + + !> A class that that is used to store simulation history data between file output type, extends(base_storage) :: encounter_storage + class(encounter_netcdf_parameters), allocatable :: nc !! NetCDF object attached to this storage object contains procedure :: dump => encounter_io_dump !! Dumps contents of encounter history to file procedure :: get_index_values => encounter_util_get_vals_storage !! Gets the unique values of the indices of a storage object (i.e. body id or time value) @@ -67,17 +82,7 @@ module encounter final :: encounter_util_final_storage end type encounter_storage - !> NetCDF dimension and variable names for the enounter save object - type, extends(base_io_netcdf_parameters) :: encounter_io_parameters - character(NAMELEN) :: loop_varname = "loopnum" !! Loop number for encounter - integer(I4B) :: loop_varid !! ID for the recursion level variable - integer(I4B) :: time_dimsize = 0 !! Number of time values in snapshot - integer(I4B) :: name_dimsize = 0 !! Number of potential id values in snapshot - integer(I4B) :: file_number = 1 !! The number to append on the output file - contains - procedure :: initialize => encounter_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object - procedure :: open => encounter_io_netcdf_open - end type encounter_io_parameters + type encounter_bounding_box_1D integer(I4B) :: n !! Number of bodies with extents @@ -90,6 +95,7 @@ module encounter final :: encounter_util_final_aabb !! Finalize the axis-aligned bounding box (1D) - deallocates all allocatables end type + type encounter_bounding_box type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb contains @@ -99,6 +105,7 @@ module encounter generic :: sweep => sweep_single, sweep_double end type + interface module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) use base, only: base_parameters @@ -219,16 +226,16 @@ end subroutine encounter_io_dump module subroutine encounter_io_initialize_output(self, param) implicit none - class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(encounter_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset class(base_parameters), intent(in) :: param end subroutine encounter_io_initialize_output - module subroutine encounter_io_netcdf_open(self, param, readonly) + module subroutine encounter_netcdf_io_open(self, param, readonly) implicit none - class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(encounter_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset class(base_parameters), intent(in) :: param !! Current run configuration parameters logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only - end subroutine encounter_io_netcdf_open + end subroutine encounter_netcdf_io_open module subroutine encounter_io_write_frame_snapshot(self, history, param) implicit none diff --git a/src/netcdf_io/netcdf_io_implementations.f90 b/src/netcdf_io/netcdf_io_implementations.f90 new file mode 100644 index 000000000..d08771203 --- /dev/null +++ b/src/netcdf_io/netcdf_io_implementations.f90 @@ -0,0 +1,66 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule (netcdf_io) s_netcdf_io_implementations + use netcdf +contains + + module subroutine netcdf_io_check(status, call_identifier) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Checks the status of all NetCDF operations to catch errors + use netcdf + implicit none + ! Arguments + integer, intent (in) :: status !! The status code returned by a NetCDF function + character(len=*), intent(in), optional :: call_identifier !! String that indicates which calling function caused the error for diagnostic purposes + + if(status /= nf90_noerr) then + if (present(call_identifier)) write(*,*) "NetCDF error in ",trim(call_identifier) + write(*,*) trim(nf90_strerror(status)) + call swiftest_util_exit(FAILURE) + end if + + return + end subroutine netcdf_io_check + + + module subroutine netcdf_io_close(self) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Closes a NetCDF file + use netcdf + implicit none + ! Arguments + class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + + call netcdf_io_check( nf90_close(self%id), "netcdf_io_close" ) + + return + end subroutine netcdf_io_close + + + module subroutine netcdf_io_sync(self) + !! author: David A. Minton + !! + !! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) + !! + use netcdf + implicit none + ! Arguments + class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + + call netcdf_io_check( nf90_sync(self%id), "netcdf_io_sync nf90_sync" ) + + return + end subroutine netcdf_io_sync + + + +end submodule s_netcdf_io_implementations diff --git a/src/netcdf_io/netcdf_io_module.f90 b/src/netcdf_io/netcdf_io_module.f90 new file mode 100644 index 000000000..bb3f098b8 --- /dev/null +++ b/src/netcdf_io/netcdf_io_module.f90 @@ -0,0 +1,159 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +module netcdf_io + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Base type definitions. This allows the collision and encounter modules to be defined before the swiftest module. + !! + use globals + use base + implicit none + public + + + !! This derived datatype stores the NetCDF ID values for each of the variables included in the NetCDF data file. This is used as the base class defined in base + type, abstract :: netcdf_parameters + character(STRMAX) :: file_name !! Name of the output file + integer(I4B) :: out_type !! output type (will be assigned either NF90_DOUBLE or NF90_FLOAT, depending on the user parameter) + integer(I4B) :: id !! ID for the output file + integer(I4B) :: name_chunk !! Chunk size for the id dimension variables + integer(I4B) :: time_chunk !! Chunk size for the time dimension variables + + ! Dimension ids and variable names + character(NAMELEN) :: str_dimname = "string32" !! name of the character string dimension + integer(I4B) :: str_dimid !! ID for the character string dimension + character(NAMELEN) :: time_dimname = "time" !! name of the time dimension + integer(I4B) :: time_dimid !! ID for the time dimension + integer(I4B) :: time_varid !! ID for the time variable + character(NAMELEN) :: name_dimname = "name" !! name of the particle name dimension + integer(I4B) :: name_dimid !! ID for the particle name dimension + integer(I4B) :: name_varid !! ID for the particle name variable + character(NAMELEN) :: space_dimname = "space" !! name of the space dimension + integer(I4B) :: space_dimid !! ID for the space dimension + integer(I4B) :: space_varid !! ID for the space variable + character(len=1), dimension(3) :: space_coords = ["x","y","z"] !! The space dimension coordinate labels + + ! Non-dimension ids and variable names + character(NAMELEN) :: id_varname = "id" !! name of the particle id variable + integer(I4B) :: id_varid !! ID for the id variable + character(NAMELEN) :: ptype_varname = "particle_type" !! name of the particle type variable + integer(I4B) :: ptype_varid !! ID for the particle type variable + character(NAMELEN) :: npl_varname = "npl" !! name of the number of active massive bodies variable + integer(I4B) :: npl_varid !! ID for the number of active massive bodies variable + character(NAMELEN) :: ntp_varname = "ntp" !! name of the number of active test particles variable + integer(I4B) :: ntp_varid !! ID for the number of active test particles variable + character(NAMELEN) :: nplm_varname = "nplm" !! name of the number of active fully interacting massive bodies variable (SyMBA) + integer(I4B) :: nplm_varid !! ID for the number of active fully interacting massive bodies variable (SyMBA) + character(NAMELEN) :: a_varname = "a" !! name of the semimajor axis variable + integer(I4B) :: a_varid !! ID for the semimajor axis variable + character(NAMELEN) :: e_varname = "e" !! name of the eccentricity variable + integer(I4B) :: e_varid !! ID for the eccentricity variable + character(NAMELEN) :: inc_varname = "inc" !! name of the inclination variable + integer(I4B) :: inc_varid !! ID for the inclination variable + character(NAMELEN) :: capom_varname = "capom" !! name of the long. asc. node variable + integer(I4B) :: capom_varid !! ID for the long. asc. node variable + character(NAMELEN) :: omega_varname = "omega" !! name of the arg. of periapsis variable + integer(I4B) :: omega_varid !! ID for the arg. of periapsis variable + character(NAMELEN) :: capm_varname = "capm" !! name of the mean anomaly variable + integer(I4B) :: capm_varid !! ID for the mean anomaly variable + character(NAMELEN) :: varpi_varname = "varpi" !! name of the long. of periapsis variable + integer(I4B) :: varpi_varid !! ID for the long. of periapsis variable + character(NAMELEN) :: lam_varname = "lam" !! name of the mean longitude variable + integer(I4B) :: lam_varid !! ID for the mean longitude variable + character(NAMELEN) :: f_varname = "f" !! name of the true anomaly variable + integer(I4B) :: f_varid !! ID for the true anomaly variable + character(NAMELEN) :: cape_varname = "cape" !! name of the eccentric anomaly variable + integer(I4B) :: cape_varid !! ID for the eccentric anomaly variable + character(NAMELEN) :: rh_varname = "rh" !! name of the heliocentric position vector variable + integer(I4B) :: rh_varid !! ID for the heliocentric position vector variable + character(NAMELEN) :: vh_varname = "vh" !! name of the heliocentric velocity vector variable + integer(I4B) :: vh_varid !! ID for the heliocentric velocity vector variable + character(NAMELEN) :: gr_pseudo_vh_varname = "gr_pseudo_vh" !! name of the heliocentric pseudovelocity vector variable (used in GR only) + integer(I4B) :: gr_pseudo_vh_varid !! ID for the heliocentric pseudovelocity vector variable (used in GR) + character(NAMELEN) :: gmass_varname = "Gmass" !! name of the mass variable + integer(I4B) :: Gmass_varid !! ID for the mass variable + character(NAMELEN) :: rhill_varname = "rhill" !! name of the hill radius variable + integer(I4B) :: rhill_varid !! ID for the hill radius variable + character(NAMELEN) :: radius_varname = "radius" !! name of the radius variable + integer(I4B) :: radius_varid !! ID for the radius variable + character(NAMELEN) :: Ip_varname = "Ip" !! name of the principal moment of inertial variable + integer(I4B) :: Ip_varid !! ID for the axis principal moment of inertia variable + character(NAMELEN) :: rot_varname = "rot" !! name of the rotation vector variable + integer(I4B) :: rot_varid !! ID for the rotation vector variable + character(NAMELEN) :: j2rp2_varname = "j2rp2" !! name of the j2rp2 variable + integer(I4B) :: j2rp2_varid !! ID for the j2 variable + character(NAMELEN) :: j4rp4_varname = "j4rp4" !! name of the j4pr4 variable + integer(I4B) :: j4rp4_varid !! ID for the j4 variable + character(NAMELEN) :: k2_varname = "k2" !! name of the Love number variable + integer(I4B) :: k2_varid !! ID for the Love number variable + character(NAMELEN) :: q_varname = "Q" !! name of the energy dissipation variable + integer(I4B) :: Q_varid !! ID for the energy dissipation variable + character(NAMELEN) :: ke_orb_varname = "KE_orb" !! name of the system orbital kinetic energy variable + integer(I4B) :: KE_orb_varid !! ID for the system orbital kinetic energy variable + character(NAMELEN) :: ke_spin_varname = "KE_spin" !! name of the system spin kinetic energy variable + integer(I4B) :: KE_spin_varid !! ID for the system spin kinetic energy variable + character(NAMELEN) :: pe_varname = "PE" !! name of the system potential energy variable + integer(I4B) :: PE_varid !! ID for the system potential energy variable + character(NAMELEN) :: L_orb_varname = "L_orb" !! name of the orbital angular momentum vector variable + integer(I4B) :: L_orb_varid !! ID for the system orbital angular momentum vector variable + character(NAMELEN) :: Lspin_varname = "Lspin" !! name of the spin angular momentum vector variable + integer(I4B) :: Lspin_varid !! ID for the system spin angular momentum vector variable + character(NAMELEN) :: L_escape_varname = "L_escape" !! name of the escaped angular momentum vector variable + integer(I4B) :: L_escape_varid !! ID for the escaped angular momentum vector variable + character(NAMELEN) :: Ecollisions_varname = "Ecollisions" !! name of the escaped angular momentum y variable + integer(I4B) :: Ecollisions_varid !! ID for the energy lost in collisions variable + character(NAMELEN) :: Euntracked_varname = "Euntracked" !! name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) + integer(I4B) :: Euntracked_varid !! ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) + character(NAMELEN) :: GMescape_varname = "GMescape" !! name of the G*Mass of bodies that escape the system + integer(I4B) :: GMescape_varid !! ID for the G*Mass of bodies that escape the system + character(NAMELEN) :: origin_type_varname = "origin_type" !! name of the origin type variable (Initial Conditions, Disruption, etc.) + integer(I4B) :: origin_type_varid !! ID for the origin type + character(NAMELEN) :: origin_time_varname = "origin_time" !! name of the time of origin variable + integer(I4B) :: origin_time_varid !! ID for the origin time + character(NAMELEN) :: collision_id_varname = "collision_id" !! name of the collision id variable + integer(I4B) :: collision_id_varid !! Netcdf ID for the origin collision ID + character(NAMELEN) :: origin_rh_varname = "origin_rh" !! name of the heliocentric position vector of the body at the time of origin variable + integer(I4B) :: origin_rh_varid !! ID for the origin position vector variable + character(NAMELEN) :: origin_vh_varname = "origin_vh" !! name of the heliocentric velocity vector of the body at the time of origin variable + integer(I4B) :: origin_vh_varid !! ID for the origin velocity vector component + character(NAMELEN) :: discard_time_varname = "discard_time" !! name of the time of discard variable + integer(I4B) :: discard_time_varid !! ID for the time of discard variable + character(NAMELEN) :: discard_rh_varname = "discard_rh" !! name of the heliocentric position vector of the body at the time of discard variable + integer(I4B) :: discard_rh_varid !! ID for the heliocentric position vector of the body at the time of discard variable + character(NAMELEN) :: discard_vh_varname = "discard_vh" !! name of the heliocentric velocity vector of the body at the time of discard variable + integer(I4B) :: discard_vh_varid !! ID for the heliocentric velocity vector of the body at the time of discard variable + character(NAMELEN) :: discard_body_id_varname = "discard_body_id" !! name of the id of the other body involved in the discard + integer(I4B) :: discard_body_id_varid !! ID for the id of the other body involved in the discard + logical :: lpseudo_vel_exists = .false. !! Logical flag to indicate whether or not the pseudovelocity vectors were present in an old file. + contains + procedure :: close => netcdf_io_close !! Closes an open NetCDF file + procedure :: sync => netcdf_io_sync !! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) + end type netcdf_parameters + + interface + module subroutine netcdf_io_check(status, call_identifier) + implicit none + integer, intent (in) :: status !! The status code returned by a NetCDF function + character(len=*), intent(in), optional :: call_identifier !! String that indicates which calling function caused the error for diagnostic purposes + end subroutine netcdf_io_check + + module subroutine netcdf_io_close(self) + implicit none + class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + end subroutine netcdf_io_close + + module subroutine netcdf_io_sync(self) + implicit none + class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + end subroutine netcdf_io_sync + end interface + + +end module netcdf_io diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index cf2fe02d6..9c5fa4eb2 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -9,6 +9,7 @@ submodule (swiftest) s_io use symba + use netcdf contains module subroutine swiftest_io_compact_output(self, param, timer) @@ -246,7 +247,7 @@ module subroutine swiftest_io_dump_system(self, param) dump_param%in_type = "NETCDF_DOUBLE" dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) associate(nc => dump_param%system_history%nc) - nc%id_chunk = self%pl%nbody + self%tp%nbody + nc%name_chunk = self%pl%nbody + self%tp%nbody nc%time_chunk = 1 dump_param%tstart = self%t @@ -264,11 +265,8 @@ module subroutine swiftest_io_dump_system(self, param) if (idx > NDUMPFILES) idx = 1 ! Dump the encounter history if necessary - select type(param) - class is (swiftest_parameters) - if (param%lenc_save_trajectory .or. param%lenc_save_closest) call self%encounter_history%dump(param) - call self%collision_history%dump(param) - end select + if (param%lenc_save_trajectory .or. param%lenc_save_closest) call self%encounter_history%dump(param) + call self%collision_history%dump(param) ! Dump the system history to file call param%system_history%dump(param) @@ -481,6 +479,1240 @@ module subroutine swiftest_io_log_start(param, file, header) end subroutine swiftest_io_log_start + module subroutine swiftest_io_netcdf_flush(self, param) + !! author: David A. Minton + !! + !! Flushes the current buffer to disk by closing and re-opening the file. + !! + implicit none + ! Arguments + class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + call self%close() + call self%open(param) + + return + end subroutine swiftest_io_netcdf_flush + + + module function swiftest_io_netcdf_get_old_t_final_system(self, param) result(old_t_final) + !! author: David A. Minton + !! + !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output. + !! + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self + class(swiftest_parameters), intent(inout) :: param + ! Result + real(DP) :: old_t_final + ! Internals + integer(I4B) :: itmax, idmax + real(DP), dimension(:), allocatable :: vals + real(DP), dimension(1) :: rtemp + real(DP), dimension(NDIM) :: rot0, Ip0, Lnow + real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig + + associate (nc => param%system_history%nc, cb => self%cb) + call nc%open(param) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=itmax), "netcdf_io_get_old_t_final_system time_dimid" ) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "netcdf_io_get_old_t_final_system name_dimid" ) + allocate(vals(idmax)) + call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_old_t_final_system time_varid" ) + + !old_t_final = rtemp(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 netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_old_t_final_system KE_orb_varid" ) + KE_orb_orig = rtemp(1) + + call netcdf_io_check( nf90_get_var(nc%id, nc%KE_spin_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_old_t_final_system KE_spin_varid" ) + KE_spin_orig = rtemp(1) + + call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_old_t_final_system PE_varid" ) + PE_orig = rtemp(1) + + call netcdf_io_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[1]), "netcdf_io_get_old_t_final_system Ecollisions_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[1]), "netcdf_io_get_old_t_final_system Euntracked_varid" ) + + self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + self%Ecollisions + self%Euntracked + + call netcdf_io_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_old_t_final_system L_orb_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_old_t_final_system Lspin_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_old_t_final_system L_escape_varid" ) + + self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) + + call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "netcdf_io_get_old_t_final_system Gmass_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[1]), "netcdf_io_get_old_t_final_system GMescape_varid" ) + self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + self%GMescape + + cb%GM0 = vals(1) + cb%dGM = cb%Gmass - cb%GM0 + + call netcdf_io_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1,1], count=[1,1]), "netcdf_io_get_old_t_final_system radius_varid" ) + cb%R0 = rtemp(1) + + if (param%lrotation) then + + call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_io_get_old_t_final_system rot_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_io_get_old_t_final_system Ip_varid" ) + + 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 if + + deallocate(vals) + end associate + + return + end function swiftest_io_netcdf_get_old_t_final_system + + + module subroutine swiftest_io_netcdf_initialize_output(self, param) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Initialize a NetCDF file system and defines all variables. + use, intrinsic :: ieee_arithmetic + implicit none + ! Arguments + class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: nvar, varid, vartype + real(DP) :: dfill + real(SP) :: sfill + integer(I4B), parameter :: NO_FILL = 0 + logical :: fileExists + character(len=STRMAX) :: errmsg + integer(I4B) :: ndims + + associate(nc => self) + + dfill = ieee_value(dfill, IEEE_QUIET_NAN) + sfill = ieee_value(sfill, IEEE_QUIET_NAN) + + select case (param%out_type) + case("netcdf_io_FLOAT") + nc%out_type = NF90_FLOAT + case("netcdf_io_DOUBLE") + nc%out_type = NF90_DOUBLE + end select + + ! Check if the file exists, and if it does, delete it + inquire(file=nc%file_name, exist=fileExists) + if (fileExists) then + open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg) + close(unit=LUN, status="delete") + end if + + ! Create the file + call netcdf_io_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "netcdf_io_initialize_output nf90_create" ) + + ! Dimensions + call netcdf_io_check( nf90_def_dim(nc%id, nc%time_dimname, NF90_UNLIMITED, nc%time_dimid), "netcdf_io_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "netcdf_io_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, NF90_UNLIMITED, nc%name_dimid), "netcdf_io_initialize_output nf90_def_dim name_dimid" ) ! dimension to store particle id numbers + call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "netcdf_io_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + + ! Dimension coordinates + call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "netcdf_io_initialize_output nf90_def_var time_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "netcdf_io_initialize_output nf90_def_var space_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "netcdf_io_initialize_output nf90_def_var name_varid" ) + + ! Variables + call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "netcdf_io_initialize_output nf90_def_var id_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%npl_varname, NF90_INT, nc%time_dimid, nc%npl_varid), "netcdf_io_initialize_output nf90_def_var npl_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%ntp_varname, NF90_INT, nc%time_dimid, nc%ntp_varid), "netcdf_io_initialize_output nf90_def_var ntp_varid" ) + if (param%lmtiny_pl) call netcdf_io_check( nf90_def_var(nc%id, nc%nplm_varname, NF90_INT, nc%time_dimid, nc%nplm_varid), "netcdf_io_initialize_output nf90_def_var nplm_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%ptype_varid), "netcdf_io_initialize_output nf90_def_var ptype_varid" ) + + if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then + call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rh_varid), "netcdf_io_initialize_output nf90_def_var rh_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%vh_varid), "netcdf_io_initialize_output nf90_def_var vh_varid" ) + + !! When GR is enabled, we need to save the pseudovelocity vectors in addition to the true heliocentric velocity vectors, otherwise + !! we cannnot expect bit-identical runs from restarted runs with GR enabled due to floating point errors during the conversion. + if (param%lgr) then + call netcdf_io_check( nf90_def_var(nc%id, nc%gr_pseudo_vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%gr_pseudo_vh_varid), "netcdf_io_initialize_output nf90_def_var gr_psuedo_vh_varid" ) + nc%lpseudo_vel_exists = .true. + end if + + end if + + if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then + call netcdf_io_check( nf90_def_var(nc%id, nc%a_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%a_varid), "netcdf_io_initialize_output nf90_def_var a_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%e_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%e_varid), "netcdf_io_initialize_output nf90_def_var e_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%inc_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%inc_varid), "netcdf_io_initialize_output nf90_def_var inc_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%capom_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%capom_varid), "netcdf_io_initialize_output nf90_def_var capom_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%omega_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%omega_varid), "netcdf_io_initialize_output nf90_def_var omega_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%capm_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%capm_varid), "netcdf_io_initialize_output nf90_def_var capm_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%varpi_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%varpi_varid), "netcdf_io_initialize_output nf90_def_var varpi_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%lam_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%lam_varid), "netcdf_io_initialize_output nf90_def_var lam_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%f_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%f_varid), "netcdf_io_initialize_output nf90_def_var f_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%cape_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%cape_varid), "netcdf_io_initialize_output nf90_def_var cape_varid" ) + end if + + call netcdf_io_check( nf90_def_var(nc%id, nc%gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "netcdf_io_initialize_output nf90_def_var Gmass_varid" ) + + if (param%lrhill_present) then + call netcdf_io_check( nf90_def_var(nc%id, nc%rhill_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%rhill_varid), "netcdf_io_initialize_output nf90_def_var rhill_varid" ) + end if + + if (param%lclose) then + call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%radius_varid), "netcdf_io_initialize_output nf90_def_var radius_varid" ) + + call netcdf_io_check( nf90_def_var(nc%id, nc%origin_time_varname, nc%out_type, nc%name_dimid, nc%origin_time_varid), "netcdf_io_initialize_output nf90_def_var origin_time_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%origin_type_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], & + nc%origin_type_varid), "netcdf_io_initialize_output nf90_create" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%origin_rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], nc%origin_rh_varid), "netcdf_io_initialize_output nf90_def_var origin_rh_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%origin_vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], nc%origin_vh_varid), "netcdf_io_initialize_output nf90_def_var origin_vh_varid" ) + + call netcdf_io_check( nf90_def_var(nc%id, nc%collision_id_varname, NF90_INT, nc%name_dimid, nc%collision_id_varid), "netcdf_io_initialize_output nf90_def_var collision_id_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%discard_time_varname, nc%out_type, nc%name_dimid, nc%discard_time_varid), "netcdf_io_initialize_output nf90_def_var discard_time_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%discard_rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], nc%discard_rh_varid), "netcdf_io_initialize_output nf90_def_var discard_rh_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%discard_vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], nc%discard_vh_varid), "netcdf_io_initialize_output nf90_def_var discard_vh_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%discard_body_id_varname, NF90_INT, nc%name_dimid, nc%discard_body_id_varid), "netcdf_io_initialize_output nf90_def_var discard_body_id_varid" ) + end if + + if (param%lrotation) then + call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%Ip_varid), "netcdf_io_initialize_output nf90_def_var Ip_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rot_varid), "netcdf_io_initialize_output nf90_def_var rot_varid" ) + end if + + ! if (param%ltides) then + ! call netcdf_io_check( nf90_def_var(nc%id, nc%k2_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%k2_varid), "netcdf_io_initialize_output nf90_def_var k2_varid" ) + ! call netcdf_io_check( nf90_def_var(nc%id, nc%q_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Q_varid), "netcdf_io_initialize_output nf90_def_var Q_varid" ) + ! end if + + if (param%lenergy) then + call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type, nc%time_dimid, nc%KE_orb_varid), "netcdf_io_initialize_output nf90_def_var KE_orb_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type, nc%time_dimid, nc%KE_spin_varid), "netcdf_io_initialize_output nf90_def_var KE_spin_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, nc%time_dimid, nc%PE_varid), "netcdf_io_initialize_output nf90_def_var PE_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_orb_varid), "netcdf_io_initialize_output nf90_def_var L_orb_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%Lspin_varid), "netcdf_io_initialize_output nf90_def_var Lspin_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_escape_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_escape_varid), "netcdf_io_initialize_output nf90_def_var L_escape_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%Ecollisions_varname, nc%out_type, nc%time_dimid, nc%Ecollisions_varid), "netcdf_io_initialize_output nf90_def_var Ecollisions_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%Euntracked_varname, nc%out_type, nc%time_dimid, nc%Euntracked_varid), "netcdf_io_initialize_output nf90_def_var Euntracked_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%GMescape_varname, nc%out_type, nc%time_dimid, nc%GMescape_varid), "netcdf_io_initialize_output nf90_def_var GMescape_varid" ) + end if + + call netcdf_io_check( nf90_def_var(nc%id, nc%j2rp2_varname, nc%out_type, nc%time_dimid, nc%j2rp2_varid), "netcdf_io_initialize_output nf90_def_var j2rp2_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%j4rp4_varname, nc%out_type, nc%time_dimid, nc%j4rp4_varid), "netcdf_io_initialize_output nf90_def_var j4rp4_varid" ) + + + ! Set fill mode to NaN for all variables + call netcdf_io_check( nf90_inquire(nc%id, nVariables=nvar), "netcdf_io_initialize_output nf90_inquire nVariables" ) + do varid = 1, nvar + call netcdf_io_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "netcdf_io_initialize_output nf90_inquire_variable" ) + select case(vartype) + case(NF90_INT) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "netcdf_io_initialize_output nf90_def_var_fill NF90_INT" ) + case(NF90_FLOAT) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "netcdf_io_initialize_output nf90_def_var_fill NF90_FLOAT" ) + case(NF90_DOUBLE) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "netcdf_io_initialize_output nf90_def_var_fill NF90_DOUBLE" ) + case(NF90_CHAR) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "netcdf_io_initialize_output nf90_def_var_fill NF90_CHAR" ) + end select + end do + + ! Set special fill mode for discard time so that we can make use of it for non-discarded bodies. + select case (vartype) + case(NF90_FLOAT) + call netcdf_io_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_SP)), "netcdf_io_initialize_output nf90_def_var_fill discard_time NF90_FLOAT" ) + case(NF90_DOUBLE) + call netcdf_io_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_DP)), "netcdf_io_initialize_output nf90_def_var_fill discard_time NF90_DOUBLE" ) + end select + + ! Take the file out of define mode + call netcdf_io_check( nf90_enddef(nc%id), "netcdf_io_initialize_output nf90_enddef" ) + + ! Add in the space dimension coordinates + call netcdf_io_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "netcdf_io_initialize_output nf90_put_var space" ) + + end associate + return + + 667 continue + write(*,*) "Error creating NetCDF output file. " // trim(adjustl(errmsg)) + call swiftest_util_exit(FAILURE) + end subroutine swiftest_io_netcdf_initialize_output + + + module subroutine swiftest_io_netcdf_open(self, param, readonly) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Opens a NetCDF file and does the variable inquiries to activate variable ids + implicit none + ! Arguments + class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only + ! Internals + integer(I4B) :: mode, status + character(len=STRMAX) :: errmsg + + mode = NF90_WRITE + if (present(readonly)) then + if (readonly) mode = NF90_NOWRITE + end if + + associate(nc => self) + + write(errmsg,*) "netcdf_io_open nf90_open ",trim(adjustl(nc%file_name)) + call netcdf_io_check( nf90_open(nc%file_name, mode, nc%id), errmsg) + + ! Dimensions + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), "netcdf_io_open nf90_inq_dimid time_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "netcdf_io_open nf90_inq_dimid space_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "netcdf_io_open nf90_inq_dimid name_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "netcdf_io_open nf90_inq_dimid str_dimid" ) + + ! Dimension coordinates + call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "netcdf_io_open nf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "netcdf_io_open nf90_inq_varid space_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "netcdf_io_open nf90_inq_varid name_varid" ) + + ! Required Variables + call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "netcdf_io_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%gmass_varname, nc%Gmass_varid), "netcdf_io_open nf90_inq_varid Gmass_varid" ) + + if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "netcdf_io_open nf90_inq_varid rh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "netcdf_io_open nf90_inq_varid vh_varid" ) + + if (param%lgr) then + !! check if pseudovelocity vectors exist in this file. If they are, set the correct flag so we know whe should not do the conversion. + status = nf90_inq_varid(nc%id, nc%gr_pseudo_vh_varname, nc%gr_pseudo_vh_varid) + nc%lpseudo_vel_exists = (status == nf90_noerr) + if (param%lrestart .and. .not.nc%lpseudo_vel_exists) then + write(*,*) "Warning! Pseudovelocity not found in input file for GR enabled run. If this is a restarted run, bit-identical trajectories are not guarunteed!" + end if + + end if + end if + + if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%a_varname, nc%a_varid), "netcdf_io_open nf90_inq_varid a_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%e_varname, nc%e_varid), "netcdf_io_open nf90_inq_varid e_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%inc_varname, nc%inc_varid), "netcdf_io_open nf90_inq_varid inc_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%capom_varname, nc%capom_varid), "netcdf_io_open nf90_inq_varid capom_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%omega_varname, nc%omega_varid), "netcdf_io_open nf90_inq_varid omega_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%capm_varname, nc%capm_varid), "netcdf_io_open nf90_inq_varid capm_varid" ) + end if + + if (param%lclose) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "netcdf_io_open nf90_inq_varid radius_varid" ) + end if + + if (param%lrotation) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "netcdf_io_open nf90_inq_varid Ip_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "netcdf_io_open nf90_inq_varid rot_varid" ) + end if + + ! if (param%ltides) then + ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%k2_varname, nc%k2_varid), "netcdf_io_open nf90_inq_varid k2_varid" ) + ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%q_varname, nc%Q_varid), "netcdf_io_open nf90_inq_varid Q_varid" ) + ! end if + + ! Optional Variables + if (param%lrhill_present) then + status = nf90_inq_varid(nc%id, nc%rhill_varname, nc%rhill_varid) + if (status /= nf90_noerr) write(*,*) "Warning! RHILL variable not set in input file. Calculating." + end if + + ! Optional variables The User Doesn't Need to Know About + status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) + status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid) + status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid) + status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid) + status = nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid) + status = nf90_inq_varid(nc%id, nc%varpi_varname, nc%varpi_varid) + status = nf90_inq_varid(nc%id, nc%lam_varname, nc%lam_varid) + status = nf90_inq_varid(nc%id, nc%f_varname, nc%f_varid) + status = nf90_inq_varid(nc%id, nc%cape_varname, nc%cape_varid) + + if (param%lmtiny_pl) status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid) + + if (param%lclose) then + status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid) + status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid) + status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid) + status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid) + status = nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid) + status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid) + status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid) + status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid) + status = nf90_inq_varid(nc%id, nc%discard_body_id_varname, nc%discard_body_id_varid) + end if + + if (param%lenergy) then + status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) + status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) + status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) + status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) + status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) + status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) + status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) + status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) + status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) + end if + + end associate + + return + end subroutine swiftest_io_netcdf_open + + + module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ierr) + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Read a frame (header plus records for each massive body and active test particle) from an output binary file + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Return + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + ! Internals + integer(I4B) :: i, tslot, idmax, npl_check, ntp_check, nplm_check, t_max, str_max, status + real(DP), dimension(:), allocatable :: rtemp + real(DP), dimension(:,:), allocatable :: vectemp + integer(I4B), dimension(:), allocatable :: itemp + logical, dimension(:), allocatable :: validmask, tpmask, plmask + + tslot = param%ioutput + + call nc%open(param, readonly=.true.) + call self%read_hdr(nc, param) + associate(cb => self%cb, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) + + call pl%setup(npl, param) + call tp%setup(ntp, param) + + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "netcdf_io_read_frame_system nf90_inquire_dimension name_dimid" ) + allocate(rtemp(idmax)) + allocate(vectemp(NDIM,idmax)) + allocate(itemp(idmax)) + allocate(validmask(idmax)) + allocate(tpmask(idmax)) + allocate(plmask(idmax)) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=t_max), "netcdf_io_read_frame_system nf90_inquire_dimension time_dimid" ) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%str_dimid, len=str_max), "netcdf_io_read_frame_system nf90_inquire_dimension str_dimid" ) + + ! First filter out only the id slots that contain valid bodies + if (param%in_form == "XV") then + call netcdf_io_check( nf90_get_var(nc%id, nc%rh_varid, vectemp(:,:), start=[1, 1, tslot]), "netcdf_io_read_frame_system filter pass nf90_getvar rh_varid" ) + validmask(:) = vectemp(1,:) == vectemp(1,:) + else + call netcdf_io_check( nf90_get_var(nc%id, nc%a_varid, rtemp(:), start=[1, tslot]), "netcdf_io_read_frame_system filter pass nf90_getvar a_varid" ) + validmask(:) = rtemp(:) == rtemp(:) + end if + + ! Next, filter only bodies that don't have mass (test particles) + call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, rtemp(:), start=[1, tslot]), "netcdf_io_read_frame_system nf90_getvar tp finder Gmass_varid" ) + plmask(:) = rtemp(:) == rtemp(:) .and. validmask(:) + tpmask(:) = .not. plmask(:) .and. validmask(:) + plmask(1) = .false. ! This is the central body + + ! Check to make sure the number of bodies is correct + npl_check = count(plmask(:)) + ntp_check = count(tpmask(:)) + + if (npl_check /= npl) then + write(*,*) "Error reading in NetCDF file: The recorded value of npl does not match the number of active massive bodies" + call swiftest_util_exit(failure) + end if + + if (ntp_check /= ntp) then + write(*,*) "Error reading in NetCDF file: The recorded value of ntp does not match the number of active test particles" + call swiftest_util_exit(failure) + end if + + if (param%lmtiny_pl) then + nplm_check = count(pack(rtemp,plmask) > param%GMTINY ) + 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 swiftest_util_exit(failure) + end if + end if + + ! Now read in each variable and split the outputs by body type + if ((param%in_form == "XV") .or. (param%in_form == "XVEL")) then + call netcdf_io_check( nf90_get_var(nc%id, nc%rh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_io_read_frame_system nf90_getvar rh_varid" ) + do i = 1, NDIM + if (npl > 0) pl%rh(i,:) = pack(vectemp(i,:), plmask(:)) + if (ntp > 0) tp%rh(i,:) = pack(vectemp(i,:), tpmask(:)) + end do + + if (param%lgr .and. nc%lpseudo_vel_exists) then + call netcdf_io_check( nf90_get_var(nc%id, nc%gr_pseudo_vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_io_read_frame_system nf90_getvar gr_pseudo_vh_varid" ) + do i = 1, NDIM + if (npl > 0) pl%vh(i,:) = pack(vectemp(i,:), plmask(:)) + if (ntp > 0) tp%vh(i,:) = pack(vectemp(i,:), tpmask(:)) + end do + else + call netcdf_io_check( nf90_get_var(nc%id, nc%vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_io_read_frame_system nf90_getvar vh_varid" ) + do i = 1, NDIM + if (npl > 0) pl%vh(i,:) = pack(vectemp(i,:), plmask(:)) + if (ntp > 0) tp%vh(i,:) = pack(vectemp(i,:), tpmask(:)) + end do + end if + end if + + if ((param%in_form == "EL") .or. (param%in_form == "XVEL")) then + call netcdf_io_check( nf90_get_var(nc%id, nc%a_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar a_varid" ) + if (.not.allocated(pl%a)) allocate(pl%a(npl)) + if (.not.allocated(tp%a)) allocate(tp%a(ntp)) + if (npl > 0) pl%a(:) = pack(rtemp, plmask) + if (ntp > 0) tp%a(:) = pack(rtemp, tpmask) + + call netcdf_io_check( nf90_get_var(nc%id, nc%e_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar e_varid" ) + if (.not.allocated(pl%e)) allocate(pl%e(npl)) + if (.not.allocated(tp%e)) allocate(tp%e(ntp)) + if (npl > 0) pl%e(:) = pack(rtemp, plmask) + if (ntp > 0) tp%e(:) = pack(rtemp, tpmask) + + call netcdf_io_check( nf90_get_var(nc%id, nc%inc_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar inc_varid" ) + rtemp = rtemp * DEG2RAD + if (.not.allocated(pl%inc)) allocate(pl%inc(npl)) + if (.not.allocated(tp%inc)) allocate(tp%inc(ntp)) + if (npl > 0) pl%inc(:) = pack(rtemp, plmask) + if (ntp > 0) tp%inc(:) = pack(rtemp, tpmask) + + call netcdf_io_check( nf90_get_var(nc%id, nc%capom_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar capom_varid" ) + rtemp = rtemp * DEG2RAD + if (.not.allocated(pl%capom)) allocate(pl%capom(npl)) + if (.not.allocated(tp%capom)) allocate(tp%capom(ntp)) + if (npl > 0) pl%capom(:) = pack(rtemp, plmask) + if (ntp > 0) tp%capom(:) = pack(rtemp, tpmask) + + call netcdf_io_check( nf90_get_var(nc%id, nc%omega_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar omega_varid" ) + rtemp = rtemp * DEG2RAD + if (.not.allocated(pl%omega)) allocate(pl%omega(npl)) + if (.not.allocated(tp%omega)) allocate(tp%omega(ntp)) + if (npl > 0) pl%omega(:) = pack(rtemp, plmask) + if (ntp > 0) tp%omega(:) = pack(rtemp, tpmask) + + call netcdf_io_check( nf90_get_var(nc%id, nc%capm_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar capm_varid" ) + rtemp = rtemp * DEG2RAD + if (.not.allocated(pl%capm)) allocate(pl%capm(npl)) + if (.not.allocated(tp%capm)) allocate(tp%capm(ntp)) + if (npl > 0) pl%capm(:) = pack(rtemp, plmask) + if (ntp > 0) tp%capm(:) = pack(rtemp, tpmask) + + end if + + call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar Gmass_varid" ) + cb%Gmass = rtemp(1) + cb%mass = cb%Gmass / param%GU + + ! Set initial central body mass for Helio bookkeeping + cb%GM0 = cb%Gmass + + + if (npl > 0) then + pl%Gmass(:) = pack(rtemp, plmask) + pl%mass(:) = pl%Gmass(:) / param%GU + + if (param%lrhill_present) then + call netcdf_io_check( nf90_get_var(nc%id, nc%rhill_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar rhill_varid" ) + pl%rhill(:) = pack(rtemp, plmask) + end if + end if + + if (param%lclose) then + call netcdf_io_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1, tslot], count=[idmax,1]), "netcdf_io_read_frame_system nf90_getvar radius_varid" ) + cb%radius = rtemp(1) + + ! Set initial central body radius for SyMBA bookkeeping + cb%R0 = cb%radius + if (npl > 0) pl%radius(:) = pack(rtemp, plmask) + else + cb%radius = param%rmin + if (npl > 0) pl%radius(:) = 0.0_DP + end if + + if (param%lrotation) then + call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_io_read_frame_system nf90_getvar Ip_varid" ) + cb%Ip(:) = vectemp(:,1) + do i = 1, NDIM + if (npl > 0) pl%Ip(i,:) = pack(vectemp(i,:), plmask(:)) + end do + + call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "netcdf_io_read_frame_system nf90_getvar rot_varid" ) + cb%rot(:) = vectemp(:,1) + do i = 1, NDIM + if (npl > 0) pl%rot(i,:) = pack(vectemp(i,:), plmask(:)) + end do + + ! Set initial central body angular momentum for bookkeeping + cb%L0(:) = cb%Ip(3) * cb%GM0 * cb%R0**2 * cb%rot(:) + end if + + ! if (param%ltides) then + ! call netcdf_io_check( nf90_get_var(nc%id, nc%k2_varid, rtemp, start=[1, tslot]), "netcdf_io_read_frame_system nf90_getvar k2_varid" ) + ! cb%k2 = rtemp(1) + ! if (npl > 0) pl%k2(:) = pack(rtemp, plmask) + + ! call netcdf_io_check( nf90_get_var(nc%id, nc%Q_varid, rtemp, start=[1, tslot]), "netcdf_io_read_frame_system nf90_getvar Q_varid" ) + ! cb%Q = rtemp(1) + ! if (npl > 0) pl%Q(:) = pack(rtemp, plmask) + ! end if + + status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%j2rp2_varid, cb%j2rp2, start=[tslot]), "netcdf_io_read_frame_system nf90_getvar j2rp2_varid" ) + else + cb%j2rp2 = 0.0_DP + end if + + status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%j4rp4_varid, cb%j4rp4, start=[tslot]), "netcdf_io_read_frame_system nf90_getvar j4rp4_varid" ) + else + cb%j4rp4 = 0.0_DP + end if + + call self%read_particle_info(nc, param, plmask, tpmask) + + if (param%in_form == "EL") then + call pl%el2xv(cb) + call tp%el2xv(cb) + end if + ! if this is a GR-enabled run, check to see if we got the pseudovelocities in. Otherwise, we'll need to generate them. + if (param%lgr .and. .not.(nc%lpseudo_vel_exists)) then + call pl%set_mu(cb) + call tp%set_mu(cb) + call pl%v2pv(param) + call tp%v2pv(param) + end if + + end associate + + call nc%close() + + ierr = 0 + return + + 667 continue + write(*,*) "Error reading system frame in netcdf_io_read_frame_system" + + end function swiftest_io_netcdf_read_frame_system + + + module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) + !! author: David A. Minton + !! + !! Reads header information (variables that change with time, but not particle id). + !! This subroutine swiftest_significantly improves the output over the original binary file, allowing us to track energy, momentum, and other quantities that + !! previously were handled as separate output files. + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: tslot, status, idmax + real(DP), dimension(:), allocatable :: gmtemp + logical, dimension(:), allocatable :: plmask, tpmask, plmmask + + + tslot = param%ioutput + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "netcdf_io_read_hdr_system nf90_inquire_dimension name_dimid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar time_varid" ) + + allocate(gmtemp(idmax)) + allocate(tpmask(idmax)) + allocate(plmask(idmax)) + allocate(plmmask(idmax)) + + call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, gmtemp, start=[1,1], count=[idmax,1]), "netcdf_io_read_hdr_system nf90_getvar Gmass_varid" ) + + plmask(:) = gmtemp(:) == gmtemp(:) + tpmask(:) = .not. plmask(:) + plmask(1) = .false. ! This is the central body + plmmask(:) = plmask(:) + + if (param%lmtiny_pl) then + where(plmask(:)) + plmmask(:) = gmtemp(:) > param%GMTINY + endwhere + end if + + status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar npl_varid" ) + else + self%pl%nbody = count(plmask(:)) + end if + + status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar ntp_varid" ) + else + self%tp%nbody = count(tpmask(:)) + end if + + if (param%lmtiny_pl) then + status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar nplm_varid" ) + else + self%pl%nplm = count(plmmask(:)) + end if + end if + + if (param%lenergy) then + status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar KE_orb_varid" ) + status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar KE_spin_varid" ) + status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar PE_varid" ) + status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_orb_varid" ) + status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar Lspin_varid" ) + status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1, tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_escape_varid" ) + status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar Ecollisions_varid" ) + status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar Euntracked_varid" ) + status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar GMescape_varid" ) + end if + + return + end subroutine swiftest_io_netcdf_read_hdr_system + + + module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, plmask, tpmask) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Reads particle information metadata from file + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + logical, dimension(:), intent(in) :: plmask !! Logical array indicating which index values belong to massive bodies + logical, dimension(:), intent(in) :: tpmask !! Logical array indicating which index values belong to test particles + + ! Internals + integer(I4B) :: i, idmax, status + real(DP), dimension(:), allocatable :: rtemp + real(DP), dimension(:,:), allocatable :: vectemp + integer(I4B), dimension(:), allocatable :: itemp + character(len=NAMELEN), dimension(:), allocatable :: ctemp + integer(I4B), dimension(:), allocatable :: plind, tpind + + ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables + idmax = size(plmask) + allocate(rtemp(idmax)) + allocate(vectemp(NDIM,idmax)) + allocate(itemp(idmax)) + allocate(ctemp(idmax)) + + associate(cb => self%cb, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) + + if (npl > 0) then + pl%status(:) = ACTIVE + pl%lmask(:) = .true. + do i = 1, npl + call pl%info(i)%set_value(status="ACTIVE") + end do + allocate(plind(npl)) + plind(:) = pack([(i, i = 1, idmax)], plmask(:)) + end if + if (ntp > 0) then + tp%status(:) = ACTIVE + tp%lmask(:) = .true. + do i = 1, ntp + call tp%info(i)%set_value(status="ACTIVE") + end do + allocate(tpind(ntp)) + tpind(:) = pack([(i, i = 1, idmax)], tpmask(:)) + end if + + call netcdf_io_check( nf90_get_var(nc%id, nc%id_varid, itemp), "netcdf_io_read_particle_info_system nf90_getvar id_varid" ) + cb%id = itemp(1) + pl%id(:) = pack(itemp, plmask) + tp%id(:) = pack(itemp, tpmask) + cb%id = 0 + pl%id(:) = pack([(i,i=0,idmax-1)],plmask) + tp%id(:) = pack([(i,i=0,idmax-1)],tpmask) + + call netcdf_io_check( nf90_get_var(nc%id, nc%name_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_io_read_particle_info_system nf90_getvar name_varid" ) + call cb%info%set_value(name=ctemp(1)) + do i = 1, npl + call pl%info(i)%set_value(name=ctemp(plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(name=ctemp(tpind(i))) + end do + + status = nf90_get_var(nc%id, nc%ptype_varid, ctemp, count=[NAMELEN, idmax]) + if (status /= nf90_noerr) then ! Set default particle types + call cb%info%set_value(particle_type=CB_TYPE_NAME) + + ! Handle semi-interacting bodies in SyMBA + do i = 1, npl + if (param%lmtiny_pl .and. (pl%Gmass(i) < param%GMTINY)) then + call pl%info(i)%set_value(particle_type=PL_TINY_TYPE_NAME) + else + call pl%info(i)%set_value(particle_type=PL_TYPE_NAME) + end if + end do + do i = 1, ntp + call tp%info(i)%set_value(particle_type=TP_TYPE_NAME) + end do + else ! Use particle types defined in input file + call cb%info%set_value(particle_type=ctemp(1)) + do i = 1, npl + call pl%info(i)%set_value(particle_type=ctemp(plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(particle_type=ctemp(tpind(i))) + end do + end if + + call cb%info%set_value(status="ACTIVE") + + if (param%lclose) then + + status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%origin_type_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_io_read_particle_info_system nf90_getvar origin_type_varid" ) + else + ctemp = "Initial Conditions" + end if + + call cb%info%set_value(origin_type=ctemp(1)) + do i = 1, npl + call pl%info(i)%set_value(origin_type=ctemp(plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(origin_type=ctemp(tpind(i))) + end do + + status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%origin_time_varid, rtemp), "netcdf_io_read_particle_info_system nf90_getvar origin_time_varid" ) + else + rtemp = param%t0 + end if + + call cb%info%set_value(origin_time=rtemp(1)) + do i = 1, npl + call pl%info(i)%set_value(origin_time=rtemp(plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(origin_time=rtemp(tpind(i))) + end do + + status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%origin_rh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar origin_rh_varid" ) + else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then + call netcdf_io_check( nf90_get_var(nc%id, nc%rh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar rh_varid" ) + else + vectemp(:,:) = 0._DP + end if + + do i = 1, npl + call pl%info(i)%set_value(origin_rh=vectemp(:,plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(origin_rh=vectemp(:,tpind(i))) + end do + + status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%origin_vh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar origin_vh_varid" ) + else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then + call netcdf_io_check( nf90_get_var(nc%id, nc%vh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar vh_varid" ) + else + vectemp(:,:) = 0._DP + end if + + do i = 1, npl + call pl%info(i)%set_value(origin_vh=vectemp(:,plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(origin_vh=vectemp(:,tpind(i))) + end do + + status = nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%collision_id_varid, itemp), "netcdf_io_read_particle_info_system nf90_getvar collision_id_varid" ) + else + itemp = 0 + end if + + do i = 1, npl + call pl%info(i)%set_value(collision_id=itemp(plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(collision_id=itemp(tpind(i))) + end do + + status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%discard_time_varid, rtemp), "netcdf_io_read_particle_info_system nf90_getvar discard_time_varid" ) + else + select case (param%out_type) + case("netcdf_io_FLOAT") + rtemp(:) = huge(0.0_SP) + case("netcdf_io_DOUBLE") + rtemp(:) = huge(0.0_DP) + end select + end if + + call cb%info%set_value(discard_time=rtemp(1)) + do i = 1, npl + call pl%info(i)%set_value(discard_time=rtemp(plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(discard_time=rtemp(tpind(i))) + end do + + status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%discard_rh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar discard_rh_varid" ) + else + vectemp(:,:) = 0.0_DP + end if + + do i = 1, npl + call pl%info(i)%set_value(discard_rh=vectemp(:,plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(discard_rh=vectemp(:,tpind(i))) + end do + + status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid) + if (status == nf90_noerr) then + call netcdf_io_check( nf90_get_var(nc%id, nc%discard_vh_varid, vectemp(:,:)), "netcdf_io_read_particle_info_system nf90_getvar discard_vh_varid" ) + else + vectemp(:,:) = 0.0_DP + end if + + do i = 1, npl + call pl%info(i)%set_value(discard_vh=vectemp(:,plind(i))) + end do + do i = 1, ntp + call tp%info(i)%set_value(discard_vh=vectemp(:,tpind(i))) + end do + end if + + end associate + + return + end subroutine swiftest_io_netcdf_read_particle_info_system + + + module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Write a frame of output of either test particle or massive body data to the binary output file + !! Note: If outputting to orbital elements, but sure that the conversion is done prior to calling this method + implicit none + ! Arguments + class(swiftest_body), intent(in) :: self !! Swiftest base object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, j, tslot, idslot, old_mode + integer(I4B), dimension(:), allocatable :: ind + real(DP), dimension(NDIM) :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs + real(DP) :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf + + tslot = param%ioutput + + call self%write_info(nc, param) + + call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_io_write_frame_body nf90_set_fill" ) + select type(self) + class is (swiftest_body) + select type (param) + class is (swiftest_parameters) + associate(n => self%nbody) + if (n == 0) return + + call swiftest_util_sort(self%id(1:n), ind) + + do i = 1, n + j = ind(i) + idslot = self%id(j) + 1 + + !! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity + if (param%lgr) call swiftest_gr_pseudovel2vel(param, self%mu(j), self%rh(:, j), self%vh(:, j), vh(:)) + + if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then + call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, self%rh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var rh_varid" ) + if (param%lgr) then !! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity + call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, vh(:), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var vh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%gr_pseudo_vh_varid, self%vh(:, j), start=[1,idslot, tslot],count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var gr_pseudo_vhx_varid" ) + + else + call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, self%vh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var vh_varid" ) + end if + end if + + if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then + if (param%lgr) then !! For GR-enabled runs, use the true value of velocity computed above + call swiftest_orbel_xv2el(self%mu(j), self%rh(1,j), self%rh(2,j), self%rh(3,j), & + vh(1), vh(2), vh(3), & + a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf) + else !! For non-GR runs just convert from the velocity we have + call swiftest_orbel_xv2el(self%mu(j), self%rh(1,j), self%rh(2,j), self%rh(3,j), & + self%vh(1,j), self%vh(2,j), self%vh(3,j), & + a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf) + end if + call netcdf_io_check( nf90_put_var(nc%id, nc%a_varid, a, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body a_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%e_varid, e, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body e_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%inc_varid, inc * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body inc_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%capom_varid, capom * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body capom_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%omega_varid, omega * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body omega_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%capm_varid, capm * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body capm_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%varpi_varid, varpi * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body varpi_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%lam_varid, lam * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body lam_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%f_varid, f * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body f_varid" ) + if (e < 1.0_DP) then + call netcdf_io_check( nf90_put_var(nc%id, nc%cape_varid, cape * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body cape_varid" ) + else if (e > 1.0_DP) then + call netcdf_io_check( nf90_put_var(nc%id, nc%cape_varid, capf * RAD2DEG, start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body (capf) cape_varid" ) + end if + end if + + select type(self) + class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body + call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body Gmass_varid" ) + if (param%lrhill_present) then + call netcdf_io_check( nf90_put_var(nc%id, nc%rhill_varid, self%rhill(j), start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body rhill_varid" ) + end if + if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius(j), start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body radius_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var body Ip_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_body nf90_put_var body rotx_varid" ) + end if + ! if (param%ltides) then + ! call netcdf_io_check( nf90_put_var(nc%id, nc%k2_varid, self%k2(j), start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body k2_varid" ) + ! call netcdf_io_check( nf90_put_var(nc%id, nc%Q_varid, self%Q(j), start=[idslot, tslot]), "netcdf_io_write_frame_body nf90_put_var body Q_varid" ) + ! end if + + end select + end do + end associate + end select + end select + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_io_write_frame_body nf90_set_fill old_mode" ) + + return + end subroutine swiftest_io_netcdf_write_frame_body + + + module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Write a frame of output of the central body + implicit none + ! Arguments + class(swiftest_cb), intent(in) :: self !! Swiftest base object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, j, tslot, idslot, old_mode + + tslot = param%ioutput + + call self%write_info(nc, param) + + call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_io_write_frame_cb nf90_set_fill" ) + + idslot = self%id + 1 + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "netcdf_io_write_frame_cb nf90_put_var cb id_varid" ) + + call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb Gmass_varid" ) + if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius, start=[idslot, tslot]), "netcdf_io_write_frame_cb nf90_put_var cb radius_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%j2rp2_varid, self%j2rp2, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j2rp2_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), "netcdf_io_write_frame_cb nf90_put_var cb j4rp4_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cb nf90_put_var cb Ip_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "netcdf_io_write_frame_cby nf90_put_var cb rot_varid" ) + end if + + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_io_write_frame_cb nf90_set_fill old_mode" ) + + return + end subroutine swiftest_io_netcdf_write_frame_cb + + + module subroutine swiftest_io_netcdf_write_frame_system(self, nc, param) + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Write a frame (header plus records for each massive body and active test particle) to a output binary file + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + call self%write_hdr(nc, param) + call self%cb%write_frame(nc, param) + call self%pl%write_frame(nc, param) + call self%tp%write_frame(nc, param) + + return + end subroutine swiftest_io_netcdf_write_frame_system + + + module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) + !! author: David A. Minton + !! + !! Writes header information (variables that change with time, but not particle id). + !! This subroutine swiftest_significantly improves the output over the original binary file, allowing us to track energy, momentum, and other quantities that + !! previously were handled as separate output files. + implicit none + ! Arguments + class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: tslot + + tslot = param%ioutput + + call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var npl_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var ntp_varid" ) + if (param%lmtiny_pl) call netcdf_io_check( nf90_put_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var nplm_varid" ) + + if (param%lenergy) then + call netcdf_io_check( nf90_put_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_orb_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_spin_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var PE_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_orb_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var Lspin_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_escape_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var Ecollisions_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var Euntracked_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var GMescape_varid" ) + end if + + return + end subroutine swiftest_io_netcdf_write_hdr_system + + + module subroutine swiftest_io_netcdf_write_info_body(self, nc, param) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Write all current particle to file + implicit none + ! Arguments + class(swiftest_body), intent(in) :: self !! Swiftest particle object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, j, idslot, old_mode + integer(I4B), dimension(:), allocatable :: ind + character(len=:), allocatable :: charstring + + ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables + call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_io_write_info_body nf90_set_fill nf90_nofill" ) + + select type(self) + class is (swiftest_body) + associate(n => self%nbody) + if (n == 0) return + call swiftest_util_sort(self%id(1:n), ind) + + do i = 1, n + j = ind(i) + idslot = self%id(j) + 1 + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id(j), start=[idslot]), "netcdf_io_write_info_body nf90_put_var id_varid" ) + + charstring = trim(adjustl(self%info(j)%name)) + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "netcdf_io_write_info_body nf90_put_var name_varid" ) + + charstring = trim(adjustl(self%info(j)%particle_type)) + call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "netcdf_io_write_info_body nf90_put_var particle_type_varid" ) + + if (param%lclose) then + charstring = trim(adjustl(self%info(j)%origin_type)) + call netcdf_io_check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "netcdf_io_write_info_body nf90_put_var origin_type_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%origin_time_varid, self%info(j)%origin_time, start=[idslot]), "netcdf_io_write_info_body nf90_put_var origin_time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%origin_rh_varid, self%info(j)%origin_rh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var origin_rh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%origin_vh_varid, self%info(j)%origin_vh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var origin_vh_varid" ) + + call netcdf_io_check( nf90_put_var(nc%id, nc%collision_id_varid, self%info(j)%collision_id, start=[idslot]), "netcdf_io_write_info_body nf90_put_var collision_id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%discard_time_varid, self%info(j)%discard_time, start=[idslot]), "netcdf_io_write_info_body nf90_put_var discard_time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info(j)%discard_rh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var discard_rh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info(j)%discard_vh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var discard_vh_varid" ) + end if + + end do + end associate + end select + + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) ) + return + end subroutine swiftest_io_netcdf_write_info_body + + + module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Write the central body info to file + implicit none + class(swiftest_cb), intent(in) :: self !! Swiftest particle object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: idslot, old_mode + character(len=:), allocatable :: charstring + + ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables + call netcdf_io_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "netcdf_io_write_info_body nf90_set_fill nf90_nofill" ) + + idslot = self%id + 1 + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "netcdf_io_write_info_body nf90_put_var cb id_varid" ) + + charstring = trim(adjustl(self%info%name)) + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "netcdf_io_write_info_body nf90_put_var cb name_varid" ) + + charstring = trim(adjustl(self%info%particle_type)) + call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "netcdf_io_write_info_body nf90_put_var cb ptype_varid" ) + + if (param%lclose) then + charstring = trim(adjustl(self%info%origin_type)) + call netcdf_io_check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "netcdf_io_write_info_body nf90_put_var cb origin_type_varid" ) + + call netcdf_io_check( nf90_put_var(nc%id, nc%origin_time_varid, self%info%origin_time, start=[idslot]), "netcdf_io_write_info_body nf90_put_var cb origin_time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%origin_rh_varid, self%info%origin_rh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var cb origin_rh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%origin_vh_varid, self%info%origin_vh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var cb origin_vh_varid" ) + + call netcdf_io_check( nf90_put_var(nc%id, nc%collision_id_varid, self%info%collision_id, start=[idslot]), "netcdf_io_write_info_body nf90_put_var cb collision_id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%discard_time_varid, self%info%discard_time, start=[idslot]), "netcdf_io_write_info_body nf90_put_var cb discard_time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info%discard_rh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var cb discard_rh_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info%discard_vh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var cb discard_vh_varid" ) + end if + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) ) + + return + end subroutine swiftest_io_netcdf_write_info_cb + + module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! @@ -813,21 +2045,6 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) - if (self%GMTINY < 0.0_DP) then - write(iomsg,*) "GMTINY invalid or not set: ", self%GMTINY - iostat = -1 - return - end if - - if (param%lfragmentation) then - if (seed_set) then - call random_seed(put = param%seed) - else - call random_seed(get = param%seed) - end if - if (param%min_GMfrag < 0.0_DP) param%min_GMfrag = param%GMTINY - end if - ! All reporting of collision information in SyMBA (including mergers) is now recorded in the Fraggle logfile call io_log_start(param, FRAGGLE_LOG_OUT, "Fraggle logfile") @@ -852,6 +2069,23 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i return end if end if + + param%lmtiny_pl = (integrator == INT_SYMBA) + + if (param%lmtiny_pl .and. self%GMTINY < 0.0_DP) then + write(iomsg,*) "GMTINY invalid or not set: ", self%GMTINY + iostat = -1 + return + end if + + if (param%lfragmentation) then + if (seed_set) then + call random_seed(put = param%seed) + else + call random_seed(get = param%seed) + end if + if (param%min_GMfrag < 0.0_DP) param%min_GMfrag = param%GMTINY + end if ! Determine if the GR flag is set correctly for this integrator select case(integrator) @@ -1618,7 +2852,7 @@ module subroutine swiftest_io_write_frame_system(self, param) logical :: fileExists associate (nc => param%system_history%nc, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) - nc%id_chunk = npl + ntp + nc%name_chunk = npl + ntp nc%time_chunk = max(param%dump_cadence / param%istep_out, 1) nc%file_name = param%outfile if (lfirst) then diff --git a/src/swiftest/swiftest_io_netcdf.f90 b/src/swiftest/swiftest_io_netcdf.f90 deleted file mode 100644 index 04f4b6b56..000000000 --- a/src/swiftest/swiftest_io_netcdf.f90 +++ /dev/null @@ -1,1246 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule (swiftest) s_io_netcdf - use netcdf -contains - - - - module function swiftest_io_netcdf_get_old_t_final_system(self, param) result(old_t_final) - !! author: David A. Minton - !! - !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output. - !! - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self - class(base_parameters), intent(inout) :: param - ! Result - real(DP) :: old_t_final - ! Internals - integer(I4B) :: itmax, idmax - real(DP), dimension(:), allocatable :: vals - real(DP), dimension(1) :: rtemp - real(DP), dimension(NDIM) :: rot0, Ip0, Lnow - real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig - - select type(param) - class is (swiftest_parameters) - associate (nc => param%system_history%nc, cb => self%cb) - call nc%open(param) - call netcdf_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=itmax), "swiftest_io_netcdf_get_old_t_final_system time_dimid" ) - call netcdf_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "swiftest_io_netcdf_get_old_t_final_system name_dimid" ) - allocate(vals(idmax)) - call netcdf_check( nf90_get_var(nc%id, nc%time_varid, rtemp, start=[1], count=[1]), "swiftest_io_netcdf_get_old_t_final_system time_varid" ) - - !old_t_final = rtemp(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 netcdf_check( nf90_get_var(nc%id, nc%KE_orb_varid, rtemp, start=[1], count=[1]), "swiftest_io_netcdf_get_old_t_final_system KE_orb_varid" ) - KE_orb_orig = rtemp(1) - - call netcdf_check( nf90_get_var(nc%id, nc%KE_spin_varid, rtemp, start=[1], count=[1]), "swiftest_io_netcdf_get_old_t_final_system KE_spin_varid" ) - KE_spin_orig = rtemp(1) - - call netcdf_check( nf90_get_var(nc%id, nc%PE_varid, rtemp, start=[1], count=[1]), "swiftest_io_netcdf_get_old_t_final_system PE_varid" ) - PE_orig = rtemp(1) - - call netcdf_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[1]), "swiftest_io_netcdf_get_old_t_final_system Ecollisions_varid" ) - call netcdf_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[1]), "swiftest_io_netcdf_get_old_t_final_system Euntracked_varid" ) - - self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + self%Ecollisions + self%Euntracked - - call netcdf_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "swiftest_io_netcdf_get_old_t_final_system L_orb_varid" ) - call netcdf_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "swiftest_io_netcdf_get_old_t_final_system Lspin_varid" ) - call netcdf_check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "swiftest_io_netcdf_get_old_t_final_system L_escape_varid" ) - - self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) - - call netcdf_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "swiftest_io_netcdf_get_old_t_final_system Gmass_varid" ) - call netcdf_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[1]), "swiftest_io_netcdf_get_old_t_final_system GMescape_varid" ) - self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + self%GMescape - - cb%GM0 = vals(1) - cb%dGM = cb%Gmass - cb%GM0 - - call netcdf_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1,1], count=[1,1]), "swiftest_io_netcdf_get_old_t_final_system radius_varid" ) - cb%R0 = rtemp(1) - - if (param%lrotation) then - - call netcdf_check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,1], count=[NDIM,1,1]), "swiftest_io_netcdf_get_old_t_final_system rot_varid" ) - call netcdf_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,1], count=[NDIM,1,1]), "swiftest_io_netcdf_get_old_t_final_system Ip_varid" ) - - 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 if - - deallocate(vals) - end associate - end select - - return - end function swiftest_io_netcdf_get_old_t_final_system - - - module subroutine swiftest_io_netcdf_initialize_output(self, param) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Initialize a NetCDF file system and defines all variables. - use, intrinsic :: ieee_arithmetic - implicit none - ! Arguments - class(swiftest_io_netcdf_parameters), intent(inout) :: self !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: nvar, varid, vartype - real(DP) :: dfill - real(SP) :: sfill - integer(I4B), parameter :: NO_FILL = 0 - logical :: fileExists - character(len=STRMAX) :: errmsg - integer(I4B) :: ndims - - associate(nc => self) - - dfill = ieee_value(dfill, IEEE_QUIET_NAN) - sfill = ieee_value(sfill, IEEE_QUIET_NAN) - - select case (param%out_type) - case("swiftest_io_netcdf_FLOAT") - nc%out_type = NF90_FLOAT - case("swiftest_io_netcdf_DOUBLE") - nc%out_type = NF90_DOUBLE - end select - - ! Check if the file exists, and if it does, delete it - inquire(file=nc%file_name, exist=fileExists) - if (fileExists) then - open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg) - close(unit=LUN, status="delete") - end if - - ! Create the file - call netcdf_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "swiftest_io_netcdf_initialize_output nf90_create" ) - - ! Dimensions - call netcdf_check( nf90_def_dim(nc%id, nc%time_dimname, NF90_UNLIMITED, nc%time_dimid), "swiftest_io_netcdf_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension - call netcdf_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "swiftest_io_netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call netcdf_check( nf90_def_dim(nc%id, nc%name_dimname, NF90_UNLIMITED, nc%name_dimid), "swiftest_io_netcdf_initialize_output nf90_def_dim name_dimid" ) ! dimension to store particle id numbers - call netcdf_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "swiftest_io_netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - - ! Dimension coordinates - call netcdf_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "swiftest_io_netcdf_initialize_output nf90_def_var time_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "swiftest_io_netcdf_initialize_output nf90_def_var space_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "swiftest_io_netcdf_initialize_output nf90_def_var name_varid" ) - - ! Variables - call netcdf_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "swiftest_io_netcdf_initialize_output nf90_def_var id_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%npl_varname, NF90_INT, nc%time_dimid, nc%npl_varid), "swiftest_io_netcdf_initialize_output nf90_def_var npl_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%ntp_varname, NF90_INT, nc%time_dimid, nc%ntp_varid), "swiftest_io_netcdf_initialize_output nf90_def_var ntp_varid" ) - if (param%integrator == INT_SYMBA) call netcdf_check( nf90_def_var(nc%id, nc%nplm_varname, NF90_INT, nc%time_dimid, nc%nplm_varid), "swiftest_io_netcdf_initialize_output nf90_def_var nplm_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%ptype_varid), "swiftest_io_netcdf_initialize_output nf90_def_var ptype_varid" ) - - if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call netcdf_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rh_varid), "swiftest_io_netcdf_initialize_output nf90_def_var rh_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%vh_varid), "swiftest_io_netcdf_initialize_output nf90_def_var vh_varid" ) - - !! When GR is enabled, we need to save the pseudovelocity vectors in addition to the true heliocentric velocity vectors, otherwise - !! we cannnot expect bit-identical runs from restarted runs with GR enabled due to floating point errors during the conversion. - if (param%lgr) then - call netcdf_check( nf90_def_var(nc%id, nc%gr_pseudo_vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%gr_pseudo_vh_varid), "swiftest_io_netcdf_initialize_output nf90_def_var gr_psuedo_vh_varid" ) - nc%lpseudo_vel_exists = .true. - end if - - end if - - if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then - call netcdf_check( nf90_def_var(nc%id, nc%a_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%a_varid), "swiftest_io_netcdf_initialize_output nf90_def_var a_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%e_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%e_varid), "swiftest_io_netcdf_initialize_output nf90_def_var e_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%inc_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%inc_varid), "swiftest_io_netcdf_initialize_output nf90_def_var inc_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%capom_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%capom_varid), "swiftest_io_netcdf_initialize_output nf90_def_var capom_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%omega_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%omega_varid), "swiftest_io_netcdf_initialize_output nf90_def_var omega_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%capm_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%capm_varid), "swiftest_io_netcdf_initialize_output nf90_def_var capm_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%varpi_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%varpi_varid), "swiftest_io_netcdf_initialize_output nf90_def_var varpi_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%lam_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%lam_varid), "swiftest_io_netcdf_initialize_output nf90_def_var lam_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%f_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%f_varid), "swiftest_io_netcdf_initialize_output nf90_def_var f_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%cape_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%cape_varid), "swiftest_io_netcdf_initialize_output nf90_def_var cape_varid" ) - end if - - call netcdf_check( nf90_def_var(nc%id, nc%gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "swiftest_io_netcdf_initialize_output nf90_def_var Gmass_varid" ) - - if (param%lrhill_present) then - call netcdf_check( nf90_def_var(nc%id, nc%rhill_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%rhill_varid), "swiftest_io_netcdf_initialize_output nf90_def_var rhill_varid" ) - end if - - if (param%lclose) then - call netcdf_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%radius_varid), "swiftest_io_netcdf_initialize_output nf90_def_var radius_varid" ) - - call netcdf_check( nf90_def_var(nc%id, nc%origin_time_varname, nc%out_type, nc%name_dimid, nc%origin_time_varid), "swiftest_io_netcdf_initialize_output nf90_def_var origin_time_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%origin_type_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], & - nc%origin_type_varid), "swiftest_io_netcdf_initialize_output nf90_create" ) - call netcdf_check( nf90_def_var(nc%id, nc%origin_rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], nc%origin_rh_varid), "swiftest_io_netcdf_initialize_output nf90_def_var origin_rh_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%origin_vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], nc%origin_vh_varid), "swiftest_io_netcdf_initialize_output nf90_def_var origin_vh_varid" ) - - call netcdf_check( nf90_def_var(nc%id, nc%collision_id_varname, NF90_INT, nc%name_dimid, nc%collision_id_varid), "swiftest_io_netcdf_initialize_output nf90_def_var collision_id_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%discard_time_varname, nc%out_type, nc%name_dimid, nc%discard_time_varid), "swiftest_io_netcdf_initialize_output nf90_def_var discard_time_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%discard_rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], nc%discard_rh_varid), "swiftest_io_netcdf_initialize_output nf90_def_var discard_rh_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%discard_vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], nc%discard_vh_varid), "swiftest_io_netcdf_initialize_output nf90_def_var discard_vh_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%discard_body_id_varname, NF90_INT, nc%name_dimid, nc%discard_body_id_varid), "swiftest_io_netcdf_initialize_output nf90_def_var discard_body_id_varid" ) - end if - - if (param%lrotation) then - call netcdf_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%Ip_varid), "swiftest_io_netcdf_initialize_output nf90_def_var Ip_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rot_varid), "swiftest_io_netcdf_initialize_output nf90_def_var rot_varid" ) - end if - - ! if (param%ltides) then - ! call netcdf_check( nf90_def_var(nc%id, nc%k2_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%k2_varid), "swiftest_io_netcdf_initialize_output nf90_def_var k2_varid" ) - ! call netcdf_check( nf90_def_var(nc%id, nc%q_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Q_varid), "swiftest_io_netcdf_initialize_output nf90_def_var Q_varid" ) - ! end if - - if (param%lenergy) then - call netcdf_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type, nc%time_dimid, nc%KE_orb_varid), "swiftest_io_netcdf_initialize_output nf90_def_var KE_orb_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type, nc%time_dimid, nc%KE_spin_varid), "swiftest_io_netcdf_initialize_output nf90_def_var KE_spin_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, nc%time_dimid, nc%PE_varid), "swiftest_io_netcdf_initialize_output nf90_def_var PE_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_orb_varid), "swiftest_io_netcdf_initialize_output nf90_def_var L_orb_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%Lspin_varid), "swiftest_io_netcdf_initialize_output nf90_def_var Lspin_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%L_escape_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_escape_varid), "swiftest_io_netcdf_initialize_output nf90_def_var L_escape_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%Ecollisions_varname, nc%out_type, nc%time_dimid, nc%Ecollisions_varid), "swiftest_io_netcdf_initialize_output nf90_def_var Ecollisions_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%Euntracked_varname, nc%out_type, nc%time_dimid, nc%Euntracked_varid), "swiftest_io_netcdf_initialize_output nf90_def_var Euntracked_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%GMescape_varname, nc%out_type, nc%time_dimid, nc%GMescape_varid), "swiftest_io_netcdf_initialize_output nf90_def_var GMescape_varid" ) - end if - - call netcdf_check( nf90_def_var(nc%id, nc%j2rp2_varname, nc%out_type, nc%time_dimid, nc%j2rp2_varid), "swiftest_io_netcdf_initialize_output nf90_def_var j2rp2_varid" ) - call netcdf_check( nf90_def_var(nc%id, nc%j4rp4_varname, nc%out_type, nc%time_dimid, nc%j4rp4_varid), "swiftest_io_netcdf_initialize_output nf90_def_var j4rp4_varid" ) - - - ! Set fill mode to NaN for all variables - call netcdf_check( nf90_inquire(nc%id, nVariables=nvar), "swiftest_io_netcdf_initialize_output nf90_inquire nVariables" ) - do varid = 1, nvar - call netcdf_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "swiftest_io_netcdf_initialize_output nf90_inquire_variable" ) - select case(vartype) - case(NF90_INT) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "swiftest_io_netcdf_initialize_output nf90_def_var_fill NF90_INT" ) - case(NF90_FLOAT) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "swiftest_io_netcdf_initialize_output nf90_def_var_fill NF90_FLOAT" ) - case(NF90_DOUBLE) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "swiftest_io_netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE" ) - case(NF90_CHAR) - call netcdf_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "swiftest_io_netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) - end select - end do - - ! Set special fill mode for discard time so that we can make use of it for non-discarded bodies. - select case (vartype) - case(NF90_FLOAT) - call netcdf_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_SP)), "swiftest_io_netcdf_initialize_output nf90_def_var_fill discard_time NF90_FLOAT" ) - case(NF90_DOUBLE) - call netcdf_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_DP)), "swiftest_io_netcdf_initialize_output nf90_def_var_fill discard_time NF90_DOUBLE" ) - end select - - ! Take the file out of define mode - call netcdf_check( nf90_enddef(nc%id), "swiftest_io_netcdf_initialize_output nf90_enddef" ) - - ! Add in the space dimension coordinates - call netcdf_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "swiftest_io_netcdf_initialize_output nf90_put_var space" ) - - end associate - return - - 667 continue - write(*,*) "Error creating NetCDF output file. " // trim(adjustl(errmsg)) - call swiftest_util_exit(FAILURE) - end subroutine swiftest_io_netcdf_initialize_output - - - module subroutine swiftest_io_netcdf_open(self, param, readonly) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Opens a NetCDF file and does the variable inquiries to activate variable ids - implicit none - ! Arguments - class(swiftest_io_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(in) :: param !! Current run configuration parameters - logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only - ! Internals - integer(I4B) :: mode, status - character(len=STRMAX) :: errmsg - - mode = NF90_WRITE - if (present(readonly)) then - if (readonly) mode = NF90_NOWRITE - end if - - associate(nc => self) - - write(errmsg,*) "swiftest_io_netcdf_open nf90_open ",trim(adjustl(nc%file_name)) - call netcdf_check( nf90_open(nc%file_name, mode, nc%id), errmsg) - - ! Dimensions - call netcdf_check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), "swiftest_io_netcdf_open nf90_inq_dimid time_dimid" ) - call netcdf_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "swiftest_io_netcdf_open nf90_inq_dimid space_dimid" ) - call netcdf_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "swiftest_io_netcdf_open nf90_inq_dimid name_dimid" ) - call netcdf_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "swiftest_io_netcdf_open nf90_inq_dimid str_dimid" ) - - ! Dimension coordinates - call netcdf_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "swiftest_io_netcdf_open nf90_inq_varid time_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "swiftest_io_netcdf_open nf90_inq_varid space_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "swiftest_io_netcdf_open nf90_inq_varid name_varid" ) - - ! Required Variables - call netcdf_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "swiftest_io_netcdf_open nf90_inq_varid name_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%gmass_varname, nc%Gmass_varid), "swiftest_io_netcdf_open nf90_inq_varid Gmass_varid" ) - - if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call netcdf_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "swiftest_io_netcdf_open nf90_inq_varid rh_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "swiftest_io_netcdf_open nf90_inq_varid vh_varid" ) - - if (param%lgr) then - !! check if pseudovelocity vectors exist in this file. If they are, set the correct flag so we know whe should not do the conversion. - status = nf90_inq_varid(nc%id, nc%gr_pseudo_vh_varname, nc%gr_pseudo_vh_varid) - nc%lpseudo_vel_exists = (status == nf90_noerr) - if (param%lrestart .and. .not.nc%lpseudo_vel_exists) then - write(*,*) "Warning! Pseudovelocity not found in input file for GR enabled run. If this is a restarted run, bit-identical trajectories are not guarunteed!" - end if - - end if - end if - - if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then - call netcdf_check( nf90_inq_varid(nc%id, nc%a_varname, nc%a_varid), "swiftest_io_netcdf_open nf90_inq_varid a_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%e_varname, nc%e_varid), "swiftest_io_netcdf_open nf90_inq_varid e_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%inc_varname, nc%inc_varid), "swiftest_io_netcdf_open nf90_inq_varid inc_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%capom_varname, nc%capom_varid), "swiftest_io_netcdf_open nf90_inq_varid capom_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%omega_varname, nc%omega_varid), "swiftest_io_netcdf_open nf90_inq_varid omega_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%capm_varname, nc%capm_varid), "swiftest_io_netcdf_open nf90_inq_varid capm_varid" ) - end if - - if (param%lclose) then - call netcdf_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "swiftest_io_netcdf_open nf90_inq_varid radius_varid" ) - end if - - if (param%lrotation) then - call netcdf_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "swiftest_io_netcdf_open nf90_inq_varid Ip_varid" ) - call netcdf_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "swiftest_io_netcdf_open nf90_inq_varid rot_varid" ) - end if - - ! if (param%ltides) then - ! call netcdf_check( nf90_inq_varid(nc%id, nc%k2_varname, nc%k2_varid), "swiftest_io_netcdf_open nf90_inq_varid k2_varid" ) - ! call netcdf_check( nf90_inq_varid(nc%id, nc%q_varname, nc%Q_varid), "swiftest_io_netcdf_open nf90_inq_varid Q_varid" ) - ! end if - - ! Optional Variables - if (param%lrhill_present) then - status = nf90_inq_varid(nc%id, nc%rhill_varname, nc%rhill_varid) - if (status /= nf90_noerr) write(*,*) "Warning! RHILL variable not set in input file. Calculating." - end if - - ! Optional variables The User Doesn't Need to Know About - status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) - status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid) - status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid) - status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid) - status = nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid) - status = nf90_inq_varid(nc%id, nc%varpi_varname, nc%varpi_varid) - status = nf90_inq_varid(nc%id, nc%lam_varname, nc%lam_varid) - status = nf90_inq_varid(nc%id, nc%f_varname, nc%f_varid) - status = nf90_inq_varid(nc%id, nc%cape_varname, nc%cape_varid) - - if (param%integrator == INT_SYMBA) then - status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid) - end if - - if (param%lclose) then - status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid) - status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid) - status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid) - status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid) - status = nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid) - status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid) - status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid) - status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid) - status = nf90_inq_varid(nc%id, nc%discard_body_id_varname, nc%discard_body_id_varid) - end if - - if (param%lenergy) then - status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) - status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) - status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) - status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) - status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) - status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) - status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) - status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) - status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) - end if - - end associate - - return - end subroutine swiftest_io_netcdf_open - - - module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ierr) - !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott - !! - !! Read a frame (header plus records for each massive body and active test particle) from an output binary file - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - ! Return - integer(I4B) :: ierr !! Error code: returns 0 if the read is successful - ! Internals - integer(I4B) :: i, tslot, idmax, npl_check, ntp_check, nplm_check, t_max, str_max, status - real(DP), dimension(:), allocatable :: rtemp - real(DP), dimension(:,:), allocatable :: vectemp - integer(I4B), dimension(:), allocatable :: itemp - logical, dimension(:), allocatable :: validmask, tpmask, plmask - - tslot = param%ioutput - - call nc%open(param, readonly=.true.) - call self%read_hdr(nc, param) - select type(param) - class is (swiftest_parameters) - associate(cb => self%cb, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) - - call pl%setup(npl, param) - call tp%setup(ntp, param) - - call netcdf_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "swiftest_io_netcdf_read_frame_system nf90_inquire_dimension name_dimid" ) - allocate(rtemp(idmax)) - allocate(vectemp(NDIM,idmax)) - allocate(itemp(idmax)) - allocate(validmask(idmax)) - allocate(tpmask(idmax)) - allocate(plmask(idmax)) - call netcdf_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=t_max), "swiftest_io_netcdf_read_frame_system nf90_inquire_dimension time_dimid" ) - call netcdf_check( nf90_inquire_dimension(nc%id, nc%str_dimid, len=str_max), "swiftest_io_netcdf_read_frame_system nf90_inquire_dimension str_dimid" ) - - ! First filter out only the id slots that contain valid bodies - if (param%in_form == "XV") then - call netcdf_check( nf90_get_var(nc%id, nc%rh_varid, vectemp(:,:), start=[1, 1, tslot]), "swiftest_io_netcdf_read_frame_system filter pass nf90_getvar rh_varid" ) - validmask(:) = vectemp(1,:) == vectemp(1,:) - else - call netcdf_check( nf90_get_var(nc%id, nc%a_varid, rtemp(:), start=[1, tslot]), "swiftest_io_netcdf_read_frame_system filter pass nf90_getvar a_varid" ) - validmask(:) = rtemp(:) == rtemp(:) - end if - - ! Next, filter only bodies that don't have mass (test particles) - call netcdf_check( nf90_get_var(nc%id, nc%Gmass_varid, rtemp(:), start=[1, tslot]), "swiftest_io_netcdf_read_frame_system nf90_getvar tp finder Gmass_varid" ) - plmask(:) = rtemp(:) == rtemp(:) .and. validmask(:) - tpmask(:) = .not. plmask(:) .and. validmask(:) - plmask(1) = .false. ! This is the central body - - ! Check to make sure the number of bodies is correct - npl_check = count(plmask(:)) - ntp_check = count(tpmask(:)) - - if (npl_check /= npl) then - write(*,*) "Error reading in NetCDF file: The recorded value of npl does not match the number of active massive bodies" - call swiftest_util_exit(failure) - end if - - if (ntp_check /= ntp) then - write(*,*) "Error reading in NetCDF file: The recorded value of ntp does not match the number of active test particles" - call swiftest_util_exit(failure) - end if - - if (param%integrator == INT_SYMBA) then - nplm_check = count(pack(rtemp,plmask) > param%GMTINY ) - 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 swiftest_util_exit(failure) - end if - end if - - ! Now read in each variable and split the outputs by body type - if ((param%in_form == "XV") .or. (param%in_form == "XVEL")) then - call netcdf_check( nf90_get_var(nc%id, nc%rh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar rh_varid" ) - do i = 1, NDIM - if (npl > 0) pl%rh(i,:) = pack(vectemp(i,:), plmask(:)) - if (ntp > 0) tp%rh(i,:) = pack(vectemp(i,:), tpmask(:)) - end do - - if (param%lgr .and. nc%lpseudo_vel_exists) then - call netcdf_check( nf90_get_var(nc%id, nc%gr_pseudo_vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar gr_pseudo_vh_varid" ) - do i = 1, NDIM - if (npl > 0) pl%vh(i,:) = pack(vectemp(i,:), plmask(:)) - if (ntp > 0) tp%vh(i,:) = pack(vectemp(i,:), tpmask(:)) - end do - else - call netcdf_check( nf90_get_var(nc%id, nc%vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar vh_varid" ) - do i = 1, NDIM - if (npl > 0) pl%vh(i,:) = pack(vectemp(i,:), plmask(:)) - if (ntp > 0) tp%vh(i,:) = pack(vectemp(i,:), tpmask(:)) - end do - end if - end if - - if ((param%in_form == "EL") .or. (param%in_form == "XVEL")) then - call netcdf_check( nf90_get_var(nc%id, nc%a_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar a_varid" ) - if (.not.allocated(pl%a)) allocate(pl%a(npl)) - if (.not.allocated(tp%a)) allocate(tp%a(ntp)) - if (npl > 0) pl%a(:) = pack(rtemp, plmask) - if (ntp > 0) tp%a(:) = pack(rtemp, tpmask) - - call netcdf_check( nf90_get_var(nc%id, nc%e_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar e_varid" ) - if (.not.allocated(pl%e)) allocate(pl%e(npl)) - if (.not.allocated(tp%e)) allocate(tp%e(ntp)) - if (npl > 0) pl%e(:) = pack(rtemp, plmask) - if (ntp > 0) tp%e(:) = pack(rtemp, tpmask) - - call netcdf_check( nf90_get_var(nc%id, nc%inc_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar inc_varid" ) - rtemp = rtemp * DEG2RAD - if (.not.allocated(pl%inc)) allocate(pl%inc(npl)) - if (.not.allocated(tp%inc)) allocate(tp%inc(ntp)) - if (npl > 0) pl%inc(:) = pack(rtemp, plmask) - if (ntp > 0) tp%inc(:) = pack(rtemp, tpmask) - - call netcdf_check( nf90_get_var(nc%id, nc%capom_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar capom_varid" ) - rtemp = rtemp * DEG2RAD - if (.not.allocated(pl%capom)) allocate(pl%capom(npl)) - if (.not.allocated(tp%capom)) allocate(tp%capom(ntp)) - if (npl > 0) pl%capom(:) = pack(rtemp, plmask) - if (ntp > 0) tp%capom(:) = pack(rtemp, tpmask) - - call netcdf_check( nf90_get_var(nc%id, nc%omega_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar omega_varid" ) - rtemp = rtemp * DEG2RAD - if (.not.allocated(pl%omega)) allocate(pl%omega(npl)) - if (.not.allocated(tp%omega)) allocate(tp%omega(ntp)) - if (npl > 0) pl%omega(:) = pack(rtemp, plmask) - if (ntp > 0) tp%omega(:) = pack(rtemp, tpmask) - - call netcdf_check( nf90_get_var(nc%id, nc%capm_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar capm_varid" ) - rtemp = rtemp * DEG2RAD - if (.not.allocated(pl%capm)) allocate(pl%capm(npl)) - if (.not.allocated(tp%capm)) allocate(tp%capm(ntp)) - if (npl > 0) pl%capm(:) = pack(rtemp, plmask) - if (ntp > 0) tp%capm(:) = pack(rtemp, tpmask) - - end if - - call netcdf_check( nf90_get_var(nc%id, nc%Gmass_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar Gmass_varid" ) - cb%Gmass = rtemp(1) - cb%mass = cb%Gmass / param%GU - - ! Set initial central body mass for Helio bookkeeping - cb%GM0 = cb%Gmass - - - if (npl > 0) then - pl%Gmass(:) = pack(rtemp, plmask) - pl%mass(:) = pl%Gmass(:) / param%GU - - if (param%lrhill_present) then - call netcdf_check( nf90_get_var(nc%id, nc%rhill_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar rhill_varid" ) - pl%rhill(:) = pack(rtemp, plmask) - end if - end if - - if (param%lclose) then - call netcdf_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1, tslot], count=[idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar radius_varid" ) - cb%radius = rtemp(1) - - ! Set initial central body radius for SyMBA bookkeeping - cb%R0 = cb%radius - if (npl > 0) pl%radius(:) = pack(rtemp, plmask) - else - cb%radius = param%rmin - if (npl > 0) pl%radius(:) = 0.0_DP - end if - - if (param%lrotation) then - call netcdf_check( nf90_get_var(nc%id, nc%Ip_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar Ip_varid" ) - cb%Ip(:) = vectemp(:,1) - do i = 1, NDIM - if (npl > 0) pl%Ip(i,:) = pack(vectemp(i,:), plmask(:)) - end do - - call netcdf_check( nf90_get_var(nc%id, nc%rot_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), "swiftest_io_netcdf_read_frame_system nf90_getvar rot_varid" ) - cb%rot(:) = vectemp(:,1) - do i = 1, NDIM - if (npl > 0) pl%rot(i,:) = pack(vectemp(i,:), plmask(:)) - end do - - ! Set initial central body angular momentum for bookkeeping - cb%L0(:) = cb%Ip(3) * cb%GM0 * cb%R0**2 * cb%rot(:) - end if - - ! if (param%ltides) then - ! call netcdf_check( nf90_get_var(nc%id, nc%k2_varid, rtemp, start=[1, tslot]), "swiftest_io_netcdf_read_frame_system nf90_getvar k2_varid" ) - ! cb%k2 = rtemp(1) - ! if (npl > 0) pl%k2(:) = pack(rtemp, plmask) - - ! call netcdf_check( nf90_get_var(nc%id, nc%Q_varid, rtemp, start=[1, tslot]), "swiftest_io_netcdf_read_frame_system nf90_getvar Q_varid" ) - ! cb%Q = rtemp(1) - ! if (npl > 0) pl%Q(:) = pack(rtemp, plmask) - ! end if - - status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%j2rp2_varid, cb%j2rp2, start=[tslot]), "swiftest_io_netcdf_read_frame_system nf90_getvar j2rp2_varid" ) - else - cb%j2rp2 = 0.0_DP - end if - - status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%j4rp4_varid, cb%j4rp4, start=[tslot]), "swiftest_io_netcdf_read_frame_system nf90_getvar j4rp4_varid" ) - else - cb%j4rp4 = 0.0_DP - end if - - call self%read_particle_info(nc, param, plmask, tpmask) - - if (param%in_form == "EL") then - call pl%el2xv(cb) - call tp%el2xv(cb) - end if - ! if this is a GR-enabled run, check to see if we got the pseudovelocities in. Otherwise, we'll need to generate them. - if (param%lgr .and. .not.(nc%lpseudo_vel_exists)) then - call pl%set_mu(cb) - call tp%set_mu(cb) - call pl%v2pv(param) - call tp%v2pv(param) - end if - - end associate - end select - - call nc%close() - - ierr = 0 - return - - 667 continue - write(*,*) "Error reading system frame in io_netcdf_read_frame_system" - - end function swiftest_io_netcdf_read_frame_system - - - module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) - !! author: David A. Minton - !! - !! Reads header information (variables that change with time, but not particle id). - !! This subroutine swiftest_significantly improves the output over the original binary file, allowing us to track energy, momentum, and other quantities that - !! previously were handled as separate output files. - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: tslot, status, idmax - real(DP), dimension(:), allocatable :: gmtemp - logical, dimension(:), allocatable :: plmask, tpmask, plmmask - - - tslot = param%ioutput - call netcdf_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "swiftest_io_netcdf_read_hdr_system nf90_inquire_dimension name_dimid" ) - call netcdf_check( nf90_get_var(nc%id, nc%time_varid, self%t, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar time_varid" ) - - allocate(gmtemp(idmax)) - allocate(tpmask(idmax)) - allocate(plmask(idmax)) - allocate(plmmask(idmax)) - - call netcdf_check( nf90_get_var(nc%id, nc%Gmass_varid, gmtemp, start=[1,1], count=[idmax,1]), "swiftest_io_netcdf_read_hdr_system nf90_getvar Gmass_varid" ) - - plmask(:) = gmtemp(:) == gmtemp(:) - tpmask(:) = .not. plmask(:) - plmask(1) = .false. ! This is the central body - plmmask(:) = plmask(:) - where(plmask(:)) - plmmask(:) = gmtemp(:) > param%GMTINY - endwhere - - status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar npl_varid" ) - else - self%pl%nbody = count(plmask(:)) - end if - - status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar ntp_varid" ) - else - self%tp%nbody = count(tpmask(:)) - end if - - if (param%integrator == INT_SYMBA) then - status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar nplm_varid" ) - else - self%pl%nplm = count(plmmask(:)) - end if - end if - - if (param%lenergy) then - status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar KE_orb_varid" ) - status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar KE_spin_varid" ) - status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar PE_varid" ) - status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "swiftest_io_netcdf_read_hdr_system nf90_getvar L_orb_varid" ) - status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "swiftest_io_netcdf_read_hdr_system nf90_getvar Lspin_varid" ) - status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1, tslot], count=[NDIM,1]), "swiftest_io_netcdf_read_hdr_system nf90_getvar L_escape_varid" ) - status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar Ecollisions_varid" ) - status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar Euntracked_varid" ) - status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) - if (status == nf90_noerr) call netcdf_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "swiftest_io_netcdf_read_hdr_system nf90_getvar GMescape_varid" ) - end if - - return - end subroutine swiftest_io_netcdf_read_hdr_system - - - module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, plmask, tpmask) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Reads particle information metadata from file - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - logical, dimension(:), intent(in) :: plmask !! Logical array indicating which index values belong to massive bodies - logical, dimension(:), intent(in) :: tpmask !! Logical array indicating which index values belong to test particles - - ! Internals - integer(I4B) :: i, idmax, status - real(DP), dimension(:), allocatable :: rtemp - real(DP), dimension(:,:), allocatable :: vectemp - integer(I4B), dimension(:), allocatable :: itemp - character(len=NAMELEN), dimension(:), allocatable :: ctemp - integer(I4B), dimension(:), allocatable :: plind, tpind - - ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables - idmax = size(plmask) - allocate(rtemp(idmax)) - allocate(vectemp(NDIM,idmax)) - allocate(itemp(idmax)) - allocate(ctemp(idmax)) - - associate(cb => self%cb, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) - - if (npl > 0) then - pl%status(:) = ACTIVE - pl%lmask(:) = .true. - do i = 1, npl - call pl%info(i)%set_value(status="ACTIVE") - end do - allocate(plind(npl)) - plind(:) = pack([(i, i = 1, idmax)], plmask(:)) - end if - if (ntp > 0) then - tp%status(:) = ACTIVE - tp%lmask(:) = .true. - do i = 1, ntp - call tp%info(i)%set_value(status="ACTIVE") - end do - allocate(tpind(ntp)) - tpind(:) = pack([(i, i = 1, idmax)], tpmask(:)) - end if - - call netcdf_check( nf90_get_var(nc%id, nc%id_varid, itemp), "swiftest_io_netcdf_read_particle_info_system nf90_getvar id_varid" ) - cb%id = itemp(1) - pl%id(:) = pack(itemp, plmask) - tp%id(:) = pack(itemp, tpmask) - cb%id = 0 - pl%id(:) = pack([(i,i=0,idmax-1)],plmask) - tp%id(:) = pack([(i,i=0,idmax-1)],tpmask) - - call netcdf_check( nf90_get_var(nc%id, nc%name_varid, ctemp, count=[NAMELEN, idmax]), "swiftest_io_netcdf_read_particle_info_system nf90_getvar name_varid" ) - call cb%info%set_value(name=ctemp(1)) - do i = 1, npl - call pl%info(i)%set_value(name=ctemp(plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(name=ctemp(tpind(i))) - end do - - status = nf90_get_var(nc%id, nc%ptype_varid, ctemp, count=[NAMELEN, idmax]) - if (status /= nf90_noerr) then ! Set default particle types - call cb%info%set_value(particle_type=CB_TYPE_NAME) - - ! Handle semi-interacting bodies in SyMBA - if (param%integrator == INT_SYMBA) then - select type (param) - class is (swiftest_parameters) - do i = 1, npl - if (pl%Gmass(i) < param%GMTINY) then - call pl%info(i)%set_value(particle_type=PL_TINY_TYPE_NAME) - else - call pl%info(i)%set_value(particle_type=PL_TYPE_NAME) - end if - end do - end select - else ! Non-SyMBA massive bodies - do i = 1, npl - call pl%info(i)%set_value(particle_type=PL_TYPE_NAME) - end do - end if - do i = 1, ntp - call tp%info(i)%set_value(particle_type=TP_TYPE_NAME) - end do - else ! Use particle types defined in input file - call cb%info%set_value(particle_type=ctemp(1)) - do i = 1, npl - call pl%info(i)%set_value(particle_type=ctemp(plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(particle_type=ctemp(tpind(i))) - end do - end if - - call cb%info%set_value(status="ACTIVE") - - if (param%lclose) then - - status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%origin_type_varid, ctemp, count=[NAMELEN, idmax]), "swiftest_io_netcdf_read_particle_info_system nf90_getvar origin_type_varid" ) - else - ctemp = "Initial Conditions" - end if - - call cb%info%set_value(origin_type=ctemp(1)) - do i = 1, npl - call pl%info(i)%set_value(origin_type=ctemp(plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(origin_type=ctemp(tpind(i))) - end do - - status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%origin_time_varid, rtemp), "swiftest_io_netcdf_read_particle_info_system nf90_getvar origin_time_varid" ) - else - rtemp = param%t0 - end if - - call cb%info%set_value(origin_time=rtemp(1)) - do i = 1, npl - call pl%info(i)%set_value(origin_time=rtemp(plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(origin_time=rtemp(tpind(i))) - end do - - status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%origin_rh_varid, vectemp(:,:)), "swiftest_io_netcdf_read_particle_info_system nf90_getvar origin_rh_varid" ) - else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call netcdf_check( nf90_get_var(nc%id, nc%rh_varid, vectemp(:,:)), "swiftest_io_netcdf_read_particle_info_system nf90_getvar rh_varid" ) - else - vectemp(:,:) = 0._DP - end if - - do i = 1, npl - call pl%info(i)%set_value(origin_rh=vectemp(:,plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(origin_rh=vectemp(:,tpind(i))) - end do - - status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%origin_vh_varid, vectemp(:,:)), "swiftest_io_netcdf_read_particle_info_system nf90_getvar origin_vh_varid" ) - else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call netcdf_check( nf90_get_var(nc%id, nc%vh_varid, vectemp(:,:)), "swiftest_io_netcdf_read_particle_info_system nf90_getvar vh_varid" ) - else - vectemp(:,:) = 0._DP - end if - - do i = 1, npl - call pl%info(i)%set_value(origin_vh=vectemp(:,plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(origin_vh=vectemp(:,tpind(i))) - end do - - status = nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%collision_id_varid, itemp), "swiftest_io_netcdf_read_particle_info_system nf90_getvar collision_id_varid" ) - else - itemp = 0 - end if - - do i = 1, npl - call pl%info(i)%set_value(collision_id=itemp(plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(collision_id=itemp(tpind(i))) - end do - - status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%discard_time_varid, rtemp), "swiftest_io_netcdf_read_particle_info_system nf90_getvar discard_time_varid" ) - else - select case (param%out_type) - case("swiftest_io_netcdf_FLOAT") - rtemp(:) = huge(0.0_SP) - case("swiftest_io_netcdf_DOUBLE") - rtemp(:) = huge(0.0_DP) - end select - end if - - call cb%info%set_value(discard_time=rtemp(1)) - do i = 1, npl - call pl%info(i)%set_value(discard_time=rtemp(plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(discard_time=rtemp(tpind(i))) - end do - - status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%discard_rh_varid, vectemp(:,:)), "swiftest_io_netcdf_read_particle_info_system nf90_getvar discard_rh_varid" ) - else - vectemp(:,:) = 0.0_DP - end if - - do i = 1, npl - call pl%info(i)%set_value(discard_rh=vectemp(:,plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(discard_rh=vectemp(:,tpind(i))) - end do - - status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid) - if (status == nf90_noerr) then - call netcdf_check( nf90_get_var(nc%id, nc%discard_vh_varid, vectemp(:,:)), "swiftest_io_netcdf_read_particle_info_system nf90_getvar discard_vh_varid" ) - else - vectemp(:,:) = 0.0_DP - end if - - do i = 1, npl - call pl%info(i)%set_value(discard_vh=vectemp(:,plind(i))) - end do - do i = 1, ntp - call tp%info(i)%set_value(discard_vh=vectemp(:,tpind(i))) - end do - end if - - end associate - - return - end subroutine swiftest_io_netcdf_read_particle_info_system - - - module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Write a frame of output of either test particle or massive body data to the binary output file - !! Note: If outputting to orbital elements, but sure that the conversion is done prior to calling this method - implicit none - ! Arguments - class(swiftest_body), intent(in) :: self !! Swiftest base object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i, j, tslot, idslot, old_mode - integer(I4B), dimension(:), allocatable :: ind - real(DP), dimension(NDIM) :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs - real(DP) :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf - - tslot = param%ioutput - - call self%write_info(nc, param) - - call netcdf_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "swiftest_io_netcdf_write_frame_body nf90_set_fill" ) - select type(self) - class is (swiftest_body) - select type (param) - class is (swiftest_parameters) - associate(n => self%nbody) - if (n == 0) return - - call swiftest_util_sort(self%id(1:n), ind) - - do i = 1, n - j = ind(i) - idslot = self%id(j) + 1 - - !! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity - if (param%lgr) call swiftest_gr_pseudovel2vel(param, self%mu(j), self%rh(:, j), self%vh(:, j), vh(:)) - - if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call netcdf_check( nf90_put_var(nc%id, nc%rh_varid, self%rh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_body nf90_put_var rh_varid" ) - if (param%lgr) then !! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity - call netcdf_check( nf90_put_var(nc%id, nc%vh_varid, vh(:), start=[1,idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_body nf90_put_var vh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%gr_pseudo_vh_varid, self%vh(:, j), start=[1,idslot, tslot],count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_body nf90_put_var gr_pseudo_vhx_varid" ) - - else - call netcdf_check( nf90_put_var(nc%id, nc%vh_varid, self%vh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_body nf90_put_var vh_varid" ) - end if - end if - - if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then - if (param%lgr) then !! For GR-enabled runs, use the true value of velocity computed above - call swiftest_orbel_xv2el(self%mu(j), self%rh(1,j), self%rh(2,j), self%rh(3,j), & - vh(1), vh(2), vh(3), & - a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf) - else !! For non-GR runs just convert from the velocity we have - call swiftest_orbel_xv2el(self%mu(j), self%rh(1,j), self%rh(2,j), self%rh(3,j), & - self%vh(1,j), self%vh(2,j), self%vh(3,j), & - a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf) - end if - call netcdf_check( nf90_put_var(nc%id, nc%a_varid, a, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body a_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%e_varid, e, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body e_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%inc_varid, inc * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body inc_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%capom_varid, capom * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body capom_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%omega_varid, omega * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body omega_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%capm_varid, capm * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body capm_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%varpi_varid, varpi * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body varpi_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%lam_varid, lam * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body lam_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%f_varid, f * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body f_varid" ) - if (e < 1.0_DP) then - call netcdf_check( nf90_put_var(nc%id, nc%cape_varid, cape * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body cape_varid" ) - else if (e > 1.0_DP) then - call netcdf_check( nf90_put_var(nc%id, nc%cape_varid, capf * RAD2DEG, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body (capf) cape_varid" ) - end if - end if - - select type(self) - class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body - call netcdf_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body Gmass_varid" ) - if (param%lrhill_present) then - call netcdf_check( nf90_put_var(nc%id, nc%rhill_varid, self%rhill(j), start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body rhill_varid" ) - end if - if (param%lclose) call netcdf_check( nf90_put_var(nc%id, nc%radius_varid, self%radius(j), start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body radius_varid" ) - if (param%lrotation) then - call netcdf_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_body nf90_put_var body Ip_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_body nf90_put_var body rotx_varid" ) - end if - ! if (param%ltides) then - ! call netcdf_check( nf90_put_var(nc%id, nc%k2_varid, self%k2(j), start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body k2_varid" ) - ! call netcdf_check( nf90_put_var(nc%id, nc%Q_varid, self%Q(j), start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_body nf90_put_var body Q_varid" ) - ! end if - - end select - end do - end associate - end select - end select - call netcdf_check( nf90_set_fill(nc%id, old_mode, old_mode), "swiftest_io_netcdf_write_frame_body nf90_set_fill old_mode" ) - - return - end subroutine swiftest_io_netcdf_write_frame_body - - - module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Write a frame of output of the central body - implicit none - ! Arguments - class(swiftest_cb), intent(in) :: self !! Swiftest base object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i, j, tslot, idslot, old_mode - - tslot = param%ioutput - - call self%write_info(nc, param) - - call netcdf_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "swiftest_io_netcdf_write_frame_cb nf90_set_fill" ) - - idslot = self%id + 1 - call netcdf_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb id_varid" ) - - call netcdf_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb Gmass_varid" ) - if (param%lclose) call netcdf_check( nf90_put_var(nc%id, nc%radius_varid, self%radius, start=[idslot, tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb radius_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%j2rp2_varid, self%j2rp2, start=[tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb j2rp2_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb j4rp4_varid" ) - if (param%lrotation) then - call netcdf_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_cb nf90_put_var cb Ip_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:), start=[1, idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_cby nf90_put_var cb rot_varid" ) - end if - - call netcdf_check( nf90_set_fill(nc%id, old_mode, old_mode), "swiftest_io_netcdf_write_frame_cb nf90_set_fill old_mode" ) - - return - end subroutine swiftest_io_netcdf_write_frame_cb - - - module subroutine swiftest_io_netcdf_write_frame_system(self, nc, param) - !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott - !! - !! Write a frame (header plus records for each massive body and active test particle) to a output binary file - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - - call self%write_hdr(nc, param) - call self%cb%write_frame(nc, param) - call self%pl%write_frame(nc, param) - call self%tp%write_frame(nc, param) - - return - end subroutine swiftest_io_netcdf_write_frame_system - - - module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) - !! author: David A. Minton - !! - !! Writes header information (variables that change with time, but not particle id). - !! This subroutine swiftest_significantly improves the output over the original binary file, allowing us to track energy, momentum, and other quantities that - !! previously were handled as separate output files. - implicit none - ! Arguments - class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: tslot - - tslot = param%ioutput - - call netcdf_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var time_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var npl_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var ntp_varid" ) - if (param%integrator == INT_SYMBA) call netcdf_check( nf90_put_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var nplm_varid" ) - - if (param%lenergy) then - call netcdf_check( nf90_put_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var KE_orb_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var KE_spin_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var PE_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "swiftest_io_netcdf_write_hdr_system nf90_put_var L_orb_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "swiftest_io_netcdf_write_hdr_system nf90_put_var Lspin_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,tslot], count=[NDIM,1]), "swiftest_io_netcdf_write_hdr_system nf90_put_var L_escape_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var Ecollisions_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var Euntracked_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "swiftest_io_netcdf_write_hdr_system nf90_put_var GMescape_varid" ) - end if - - return - end subroutine swiftest_io_netcdf_write_hdr_system - - - module subroutine swiftest_io_netcdf_write_info_body(self, nc, param) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Write all current particle to file - implicit none - ! Arguments - class(swiftest_body), intent(in) :: self !! Swiftest particle object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i, j, idslot, old_mode - integer(I4B), dimension(:), allocatable :: ind - character(len=:), allocatable :: charstring - - ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables - call netcdf_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "swiftest_io_netcdf_write_info_body nf90_set_fill nf90_nofill" ) - - select type(self) - class is (swiftest_body) - associate(n => self%nbody) - if (n == 0) return - call swiftest_util_sort(self%id(1:n), ind) - - do i = 1, n - j = ind(i) - idslot = self%id(j) + 1 - call netcdf_check( nf90_put_var(nc%id, nc%id_varid, self%id(j), start=[idslot]), "swiftest_io_netcdf_write_info_body nf90_put_var id_varid" ) - - charstring = trim(adjustl(self%info(j)%name)) - call netcdf_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "swiftest_io_netcdf_write_info_body nf90_put_var name_varid" ) - - charstring = trim(adjustl(self%info(j)%particle_type)) - call netcdf_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "swiftest_io_netcdf_write_info_body nf90_put_var particle_type_varid" ) - - if (param%lclose) then - charstring = trim(adjustl(self%info(j)%origin_type)) - call netcdf_check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "swiftest_io_netcdf_write_info_body nf90_put_var origin_type_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%origin_time_varid, self%info(j)%origin_time, start=[idslot]), "swiftest_io_netcdf_write_info_body nf90_put_var origin_time_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%origin_rh_varid, self%info(j)%origin_rh(:), start=[1,idslot], count=[NDIM,1]), "swiftest_io_netcdf_write_info_body nf90_put_var origin_rh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%origin_vh_varid, self%info(j)%origin_vh(:), start=[1,idslot], count=[NDIM,1]), "swiftest_io_netcdf_write_info_body nf90_put_var origin_vh_varid" ) - - call netcdf_check( nf90_put_var(nc%id, nc%collision_id_varid, self%info(j)%collision_id, start=[idslot]), "swiftest_io_netcdf_write_info_body nf90_put_var collision_id_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%discard_time_varid, self%info(j)%discard_time, start=[idslot]), "swiftest_io_netcdf_write_info_body nf90_put_var discard_time_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info(j)%discard_rh(:), start=[1,idslot], count=[NDIM,1]), "swiftest_io_netcdf_write_info_body nf90_put_var discard_rh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info(j)%discard_vh(:), start=[1,idslot], count=[NDIM,1]), "swiftest_io_netcdf_write_info_body nf90_put_var discard_vh_varid" ) - end if - - end do - end associate - end select - - call netcdf_check( nf90_set_fill(nc%id, old_mode, old_mode) ) - return - end subroutine swiftest_io_netcdf_write_info_body - - - module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param) - !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton - !! - !! Write the central body info to file - implicit none - class(swiftest_cb), intent(in) :: self !! Swiftest particle object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: idslot, old_mode - character(len=:), allocatable :: charstring - - ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables - call netcdf_check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "swiftest_io_netcdf_write_info_body nf90_set_fill nf90_nofill" ) - - idslot = self%id + 1 - call netcdf_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), "swiftest_io_netcdf_write_info_body nf90_put_var cb id_varid" ) - - charstring = trim(adjustl(self%info%name)) - call netcdf_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "swiftest_io_netcdf_write_info_body nf90_put_var cb name_varid" ) - - charstring = trim(adjustl(self%info%particle_type)) - call netcdf_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "swiftest_io_netcdf_write_info_body nf90_put_var cb ptype_varid" ) - - if (param%lclose) then - charstring = trim(adjustl(self%info%origin_type)) - call netcdf_check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "swiftest_io_netcdf_write_info_body nf90_put_var cb origin_type_varid" ) - - call netcdf_check( nf90_put_var(nc%id, nc%origin_time_varid, self%info%origin_time, start=[idslot]), "swiftest_io_netcdf_write_info_body nf90_put_var cb origin_time_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%origin_rh_varid, self%info%origin_rh(:), start=[1, idslot], count=[NDIM,1]), "swiftest_io_netcdf_write_info_body nf90_put_var cb origin_rh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%origin_vh_varid, self%info%origin_vh(:), start=[1, idslot], count=[NDIM,1]), "swiftest_io_netcdf_write_info_body nf90_put_var cb origin_vh_varid" ) - - call netcdf_check( nf90_put_var(nc%id, nc%collision_id_varid, self%info%collision_id, start=[idslot]), "swiftest_io_netcdf_write_info_body nf90_put_var cb collision_id_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%discard_time_varid, self%info%discard_time, start=[idslot]), "swiftest_io_netcdf_write_info_body nf90_put_var cb discard_time_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info%discard_rh(:), start=[1, idslot], count=[NDIM,1]), "swiftest_io_netcdf_write_info_body nf90_put_var cb discard_rh_varid" ) - call netcdf_check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info%discard_vh(:), start=[1, idslot], count=[NDIM,1]), "swiftest_io_netcdf_write_info_body nf90_put_var cb discard_vh_varid" ) - end if - call netcdf_check( nf90_set_fill(nc%id, old_mode, old_mode) ) - - return - end subroutine swiftest_io_netcdf_write_info_cb - -end submodule s_io_netcdf diff --git a/src/swiftest/swiftest_kick.f90 b/src/swiftest/swiftest_kick.f90 index 584750098..9e9576321 100644 --- a/src/swiftest/swiftest_kick.f90 +++ b/src/swiftest/swiftest_kick.f90 @@ -32,7 +32,7 @@ module subroutine swiftest_kick_getacch_int_pl(self, param) ! call itimer%time_this_loop(param, self%nplpl, self) ! lfirst = .false. ! else - ! if (itimer%io_netcdf_check(param, self%nplpl)) call itimer%time_this_loop(param, self%nplpl, self) + ! if (itimer%netcdf_io_check(param, self%nplpl)) call itimer%time_this_loop(param, self%nplpl, self) ! end if ! else ! param%lflatten_interactions = .false. diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index a7bae5336..a46f66ce2 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -42,19 +42,22 @@ module swiftest use fraggle use walltime use io_progress_bar + use netcdf_io !use advisor_annotate !$ use omp_lib implicit none public - type, extends(base_io_netcdf_parameters) :: swiftest_io_netcdf_parameters + type, extends(netcdf_parameters) :: swiftest_netcdf_parameters contains procedure :: initialize => swiftest_io_netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object - procedure :: open => swiftest_io_netcdf_open !! Opens a NetCDF file and does the variable inquiries to activate variable ids - end type swiftest_io_netcdf_parameters + procedure :: open => swiftest_io_netcdf_open !! Opens a NetCDF file and does the variable inquiries to activate variable ids + procedure :: flush => swiftest_io_netcdf_flush !! Flushes a NetCDF file by closing it then opening it again + end type swiftest_netcdf_parameters type, extends(base_storage) :: swiftest_storage + class(swiftest_netcdf_parameters), allocatable :: nc !! NetCDF object attached to this storage object contains procedure :: dump => swiftest_io_dump_storage !! Dumps storage object contents to file procedure :: get_index_values => swiftest_util_get_vals_storage !! Gets the unique values of the indices of a storage object (i.e. body id or time value) @@ -614,6 +617,104 @@ module subroutine swiftest_io_log_start(param, file, header) character(len=*), intent(in) :: header !! Header to print at top of log file end subroutine swiftest_io_log_start + module subroutine swiftest_io_netcdf_flush(self, param) + implicit none + class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_flush + + module function swiftest_io_netcdf_get_old_t_final_system(self, param) result(old_t_final) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP) :: old_t_final !! Final time from last run + end function swiftest_io_netcdf_get_old_t_final_system + + module subroutine swiftest_io_netcdf_initialize_output(self, param) + implicit none + class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_initialize_output + + module subroutine swiftest_io_netcdf_open(self, param, readonly) + implicit none + class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only + end subroutine swiftest_io_netcdf_open + + module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ierr) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + end function swiftest_io_netcdf_read_frame_system + + module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_read_hdr_system + + module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, plmask, tpmask) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + logical, dimension(:), intent(in) :: plmask !! Logical array indicating which index values belong to massive bodies + logical, dimension(:), intent(in) :: tpmask !! Logical array indicating which index values belong to test particles + end subroutine swiftest_io_netcdf_read_particle_info_system + + module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param) + implicit none + class(swiftest_body), intent(in) :: self !! Swiftest base object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_write_frame_body + + module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) + implicit none + class(swiftest_cb), intent(in) :: self !! Swiftest base object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_write_frame_cb + + module subroutine swiftest_io_netcdf_write_frame_system(self, nc, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_write_frame_system + + module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) + implicit none + class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_write_hdr_system + + module subroutine swiftest_io_netcdf_write_info_body(self, nc, param) + implicit none + class(swiftest_body), intent(in) :: self !! Swiftest particle object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_write_info_body + + module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param) + implicit none + class(swiftest_cb), intent(in) :: self !! Swiftest particle object + class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_netcdf_write_info_cb + + module subroutine swiftest_io_write_discard(self, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_io_write_discard + module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none class(swiftest_parameters), intent(inout) :: self !! Collection of parameters @@ -821,98 +922,6 @@ pure module subroutine swiftest_kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle end subroutine swiftest_kick_getacch_int_one_tp - module function swiftest_io_netcdf_get_old_t_final_system(self, param) result(old_t_final) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP) :: old_t_final !! Final time from last run - end function swiftest_io_netcdf_get_old_t_final_system - - module subroutine swiftest_io_netcdf_initialize_output(self, param) - implicit none - class(swiftest_io_netcdf_parameters), intent(inout) :: self !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine swiftest_io_netcdf_initialize_output - - module subroutine swiftest_io_netcdf_open(self, param, readonly) - implicit none - class(swiftest_io_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(in) :: param !! Current run configuration parameters - logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only - end subroutine swiftest_io_netcdf_open - - module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ierr) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B) :: ierr !! Error code: returns 0 if the read is successful - end function swiftest_io_netcdf_read_frame_system - - module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_io_netcdf_read_hdr_system - - module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, plmask, tpmask) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - logical, dimension(:), intent(in) :: plmask !! Logical array indicating which index values belong to massive bodies - logical, dimension(:), intent(in) :: tpmask !! Logical array indicating which index values belong to test particles - end subroutine swiftest_io_netcdf_read_particle_info_system - - module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param) - implicit none - class(swiftest_body), intent(in) :: self !! Swiftest base object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_io_netcdf_write_frame_body - - module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) - implicit none - class(swiftest_cb), intent(in) :: self !! Swiftest base object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_io_netcdf_write_frame_cb - - module subroutine swiftest_io_netcdf_write_frame_system(self, nc, param) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_io_netcdf_write_frame_system - - module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) - implicit none - class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_io_netcdf_write_hdr_system - - module subroutine swiftest_io_netcdf_write_info_body(self, nc, param) - implicit none - class(swiftest_body), intent(in) :: self !! Swiftest particle object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_io_netcdf_write_info_body - - module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param) - implicit none - class(swiftest_cb), intent(in) :: self !! Swiftest particle object - class(base_io_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_io_netcdf_write_info_cb - - module subroutine swiftest_io_write_discard(self, param) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_io_write_discard - module subroutine swiftest_obl_acc_body(self, system) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object @@ -1892,4 +1901,5 @@ subroutine make_impactors_pl(self, idx) return end subroutine make_impactors_pl + end module swiftest diff --git a/src/swiftest/swiftest_setup.f90 b/src/swiftest/swiftest_setup.f90 index 5429681c6..f3c604ff8 100644 --- a/src/swiftest/swiftest_setup.f90 +++ b/src/swiftest/swiftest_setup.f90 @@ -30,7 +30,7 @@ module subroutine swiftest_setup_construct_system(system, param) select type(param) class is (swiftest_parameters) allocate(swiftest_storage(param%dump_cadence) :: param%system_history) - allocate(swiftest_io_netcdf_parameters :: param%system_history%nc) + allocate(swiftest_netcdf_parameters :: param%system_history%nc) call param%system_history%reset() select case(param%integrator) @@ -89,25 +89,23 @@ module subroutine swiftest_setup_construct_system(system, param) end if if (param%lenc_save_trajectory .or. param%lenc_save_closest) then - allocate(encounter_io_parameters :: encounter_history%nc) + allocate(encounter_netcdf_parameters :: encounter_history%nc) call encounter_history%reset() select type(nc => encounter_history%nc) - class is (encounter_io_parameters) + class is (encounter_netcdf_parameters) nc%file_number = param%iloop / param%dump_cadence end select allocate(system%encounter_history, source=encounter_history) end if - allocate(collision_io_parameters :: collision_history%nc) + allocate(collision_netcdf_parameters :: collision_history%nc) call collision_history%reset() select type(nc => collision_history%nc) - class is (collision_io_parameters) + class is (collision_netcdf_parameters) nc%file_number = param%iloop / param%dump_cadence end select allocate(system%collision_history, source=collision_history) - - end select case (INT_RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' @@ -304,13 +302,11 @@ module subroutine swiftest_setup_pl(self, n, param) self%nplpl = 0 if (param%lclose) then - allocate(self%lmtiny(n)) allocate(self%nplenc(n)) allocate(self%ntpenc(n)) allocate(self%radius(n)) allocate(self%density(n)) - self%lmtiny(:) = .false. self%nplenc(:) = 0 self%ntpenc(:) = 0 self%radius(:) = 0.0_DP @@ -318,6 +314,11 @@ module subroutine swiftest_setup_pl(self, n, param) end if + if (param%lmtiny_pl) then + allocate(self%lmtiny(n)) + self%lmtiny(:) = .false. + end if + if (param%lrotation) then allocate(self%rot(NDIM, n)) allocate(self%Ip(NDIM, n)) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 7e303b214..045215d87 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -2310,12 +2310,14 @@ module subroutine swiftest_util_rearray_pl(self, system, param) pl%lcollision(1:npl) = .false. pl%lmask(1:npl) = .true. - pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY - where(pl%lmtiny(1:npl)) - pl%info(1:npl)%particle_type = PL_TINY_TYPE_NAME - elsewhere - pl%info(1:npl)%particle_type = PL_TYPE_NAME - end where + if (param%lmtiny_pl) then + pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY + where(pl%lmtiny(1:npl)) + pl%info(1:npl)%particle_type = PL_TINY_TYPE_NAME + elsewhere + pl%info(1:npl)%particle_type = PL_TYPE_NAME + end where + end if call pl%write_info(param%system_history%nc, param) deallocate(ldump_mask) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 80d3ccbc2..4e9cdd1d4 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -34,7 +34,7 @@ module subroutine symba_kick_getacch_int_pl(self, param) ! call itimer%time_this_loop(param, self%nplplm, self) ! lfirst = .false. ! else - ! if (itimer%io_netcdf_check(param, self%nplplm)) call itimer%time_this_loop(param, self%nplplm, self) + ! if (itimer%netcdf_io_check(param, self%nplplm)) call itimer%time_this_loop(param, self%nplplm, self) ! end if ! else ! param%lflatten_interactions = .false. diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index bf749a35e..6b91857a5 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -276,11 +276,12 @@ module subroutine symba_util_flatten_eucl_plpl(self, param) associate(pl => self, nplplm => self%nplplm) npl = int(self%nbody, kind=I8B) - select type(param) - class is (swiftest_parameters) + if (param%lmtiny_pl) then pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY - end select - nplm = count(.not. pl%lmtiny(1:npl)) + nplm = count(.not. pl%lmtiny(1:npl)) + else + nplm = npl + end if pl%nplm = int(nplm, kind=I4B) nplplm = nplm * npl - nplm * (nplm + 1_I8B) / 2_I8B ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies