diff --git a/src/netcdf_io/netcdf_io_module.f90 b/src/netcdf_io/netcdf_io_module.f90 index 2c072186f..34196dc50 100644 --- a/src/netcdf_io/netcdf_io_module.f90 +++ b/src/netcdf_io/netcdf_io_module.f90 @@ -64,25 +64,6 @@ module netcdf_io !! ID for the space variable character(len=1), dimension(3) :: space_coords = ["x","y","z"] !! The space dimension coordinate labels - character(NAMELEN) :: sign_dimname = "sign" - !! name of the sign dimension for c_lm - integer(I4B) :: sign_dimid - !! ID for sign dimension - integer(I4B) :: sign_varid - !! ID for sign variable - character(NAMELEN) :: l_dimname = "l" - !! name of l dimension for c_lm - integer(I4B) :: l_dimid - !! ID for the l dimension for c_lm - integer(I4B) :: l_varid - !! ID for the l variable - character(NAMELEN) :: m_dimname = "m" - !! name of m dimension for c_lm - integer(I4B) :: m_dimid - !! ID for the m dimension for c_lm - integer(I4B) :: m_varid - !! ID for the m variable - ! Non-dimension ids and variable names character(NAMELEN) :: id_varname = "id" @@ -291,8 +272,35 @@ module netcdf_io !! 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. + + ! Gravitational harmonics ids and variable names logical :: lc_lm_exists = .false. !! Logical flag to indicate whether or not the c_lm array was present in an old file. + character(NAMELEN) :: sign_dimname = "sign" + !! name of the sign dimension for c_lm + integer(I4B) :: sign_dimid + !! ID for sign dimension + integer(I4B) :: sign_varid + !! ID for sign variable + integer(I4B), dimension(2) :: sign_coords = [-1,1] + !! The sign dimension coordinate labels + character(NAMELEN) :: l_dimname = "l" + !! name of l dimension for c_lm + integer(I4B) :: l_dimid + !! ID for the l dimension for c_lm + integer(I4B) :: l_varid + !! ID for the l variable + character(NAMELEN) :: m_dimname = "m" + !! name of m dimension for c_lm + integer(I4B) :: m_dimid + !! ID for the m dimension for c_lm + integer(I4B) :: m_varid + !! ID for the m variable + integer(I4B) :: m_dim_max + !! Maximum value of the m dimension + integer(I4B) :: l_dim_max + !! Maximum value of the l dimension + contains procedure :: close => netcdf_io_close !! Closes an open NetCDF file diff --git a/src/swiftest/swiftest_gr.f90 b/src/swiftest/swiftest_gr.f90 index a24ff2507..8de7e3dc8 100644 --- a/src/swiftest/swiftest_gr.f90 +++ b/src/swiftest/swiftest_gr.f90 @@ -23,9 +23,12 @@ pure module subroutine swiftest_gr_kick_getaccb_ns_body(self, nbody_system, para !! Adapted from David A. Minton's Swifter routine routine gr_kick_getaccb_ns.f90 implicit none ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest generic body object - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_body), intent(inout) :: self + !! Swiftest generic body object + class(swiftest_nbody_system), intent(inout) :: nbody_system + !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param + !! Current run configuration parameters ! Internals real(DP) :: rmag, rdotv, vmag2 integer(I4B) :: i @@ -62,12 +65,18 @@ pure module subroutine swiftest_gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) !! Adapted from David A. Minton's Swifter routine routine gr_whm_kick_getacch.f90 implicit none ! Arguments - real(DP), dimension(:), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:,:), intent(in) :: x !! Position vectors - logical, dimension(:), intent(in) :: lmask !! Logical mask indicating which bodies to compute - integer(I4B), intent(in) :: n !! Total number of bodies - real(DP), intent(in) :: inv_c2 !! Inverse speed of light squared: 1 / c**2 - real(DP), dimension(:,:), intent(out) :: agr !! Accelerations + real(DP), dimension(:), intent(in) :: mu + !! Gravitational constant + real(DP), dimension(:,:), intent(in) :: x + !! Position vectors + logical, dimension(:), intent(in) :: lmask + !! Logical mask indicating which bodies to compute + integer(I4B), intent(in) :: n + !! Total number of bodies + real(DP), intent(in) :: inv_c2 + !! Inverse speed of light squared: 1 / c**2 + real(DP), dimension(:,:), intent(out) :: agr + !! Accelerations ! Internals integer(I4B) :: i real(DP) :: beta, rjmag4 @@ -100,10 +109,14 @@ pure elemental module subroutine swiftest_gr_p4_pos_kick(inv_c2, rx, ry, rz, vx, !! Adapted from David A. Minton's Swifter routine gr_whm_p4.f90 implicit none ! Arguments - real(DP), intent(in) :: inv_c2 !! One over speed of light squared (1/c**2) - real(DP), intent(inout) :: rx, ry, rz !! Position vector - real(DP), intent(in) :: vx, vy, vz !! Velocity vector - real(DP), intent(in) :: dt !! Step size + real(DP), intent(in) :: inv_c2 + !! One over speed of light squared (1/c**2) + real(DP), intent(inout) :: rx, ry, rz + !! Position vector + real(DP), intent(in) :: vx, vy, vz + !! Velocity vector + real(DP), intent(in) :: dt + !! Step size ! Internals real(DP) :: drx, dry, drz real(DP) :: vmag2 @@ -133,11 +146,16 @@ pure module subroutine swiftest_gr_pseudovel2vel(param, mu, rh, pv, vh) !! Adapted from David A. Minton's Swifter routine gr_pseudovel2vel.f90 implicit none ! Arguments - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body - real(DP), dimension(:), intent(in) :: rh !! Heliocentric position vector - real(DP), dimension(:), intent(in) :: pv !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32) - real(DP), dimension(:), intent(out) :: vh !! Heliocentric velocity vector + class(swiftest_parameters), intent(in) :: param + !! Current run configuration parameters + real(DP), intent(in) :: mu + !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body + real(DP), dimension(:), intent(in) :: rh + !! Heliocentric position vector + real(DP), dimension(:), intent(in) :: pv + !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32) + real(DP), dimension(:), intent(out) :: vh + !! Heliocentric velocity vector ! Internals real(DP) :: vmag2, rmag, grterm @@ -191,11 +209,16 @@ pure module subroutine swiftest_gr_vel2pseudovel(param, mu, rh, vh, pv) !! Adapted from David A. Minton's Swifter routine gr_vel2pseudovel.f90 implicit none ! Arguments - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body - real(DP), dimension(:), intent(in) :: rh !! Heliocentric position vector - real(DP), dimension(:), intent(in) :: vh !! Heliocentric velocity vector - real(DP), dimension(:), intent(out) :: pv !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32) + class(swiftest_parameters), intent(in) :: param + !! Current run configuration parameters + real(DP), intent(in) :: mu + !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body + real(DP), dimension(:), intent(in) :: rh + !! Heliocentric position vector + real(DP), dimension(:), intent(in) :: vh + !! Heliocentric velocity vector + real(DP), dimension(:), intent(out) :: pv + !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32) ! Internals real(DP) :: v2, G, pv2, rterm, det real(DP), dimension(NDIM,NDIM) :: J,Jinv @@ -264,11 +287,14 @@ pure module subroutine swiftest_gr_vh2pv_body(self, param) !! Wrapper function that converts from heliocentric velocity to pseudovelocity for Swiftest bodies implicit none ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_body), intent(inout) :: self + !! Swiftest particle object + class(swiftest_parameters), intent(in) :: param + !! Current run configuration parameters ! Internals integer(I4B) :: i - real(DP), dimension(:,:), allocatable :: pv !! Temporary holder of pseudovelocity for in-place conversion + real(DP), dimension(:,:), allocatable :: pv + !! Temporary holder of pseudovelocity for in-place conversion associate(n => self%nbody) if (n == 0) return diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index c50345cd3..0dabd4f99 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -686,7 +686,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param) class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param ! Internals - integer(I4B) :: itmax, idmax, tslot + integer(I4B) :: itmax, idmax, tslot, status real(DP), dimension(:), allocatable :: vals logical, dimension(:), allocatable :: plmask, tpmask real(DP), dimension(1) :: rtemp @@ -754,6 +754,23 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param) cb%L0(:) = Ip0(3) * mass0 * cb%R0**2 * rot0(:) L(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) cb%dL(:) = L(:) - cb%L0 + + 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_get_t0_values_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_get_t0_values_system nf90_getvar j4rp4_varid" ) + else + cb%j4rp4 = 0.0_DP + end if + end if ! Retrieve the current bookkeeping variables @@ -767,15 +784,6 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param) call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[tslot]), & "netcdf_io_get_t0_values_system E_untracked_varid" ) - ! ! SH gravity variable dimensions - - ! call netcdf_io_check( nf90_get_var(nc%id, nc%sign_dimname, nc%sign_dimid), & - ! "swiftest_io_netcdf_open nf90_inq_dimid sign_dimid") - ! call netcdf_io_check( nf90_inq_dimid(nc%id, nc%l_dimname, nc%l_dimid), & - ! "swiftest_io_netcdf_open nf90_inq_dimid l_dimid") - ! call netcdf_io_check( nf90_inq_dimid(nc%id, nc%m_dimname, nc%m_dimid), & - ! "swiftest_io_netcdf_open nf90_inq_dimid m_dimid") - end if deallocate(vals) @@ -803,7 +811,7 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) integer(I4B), parameter :: NO_FILL = 0 logical :: fileExists character(len=STRMAX) :: errmsg - integer(I4B) :: ndims + integer(I4B) :: i, ndims associate(nc => self) @@ -839,13 +847,7 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) "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 - - call netcdf_io_check( nf90_def_dim(nc%id, nc%sign_dimname, NF90_UNLIMITED, nc%sign_dimid), & - "swiftest_io_netcdf_open nf90_def_dim sign_dimid") - call netcdf_io_check( nf90_def_dim(nc%id, nc%l_dimname, NF90_UNLIMITED, nc%l_dimid), & - "swiftest_io_netcdf_open nf90_def_dim l_dimid") - call netcdf_io_check( nf90_def_dim(nc%id, nc%m_dimname, NF90_UNLIMITED, nc%m_dimid), & - "swiftest_io_netcdf_open nf90_def_dim m_dimid") + ! Dimension coordinates call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), & @@ -855,13 +857,6 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) 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" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%sign_dimname, nc%out_type, nc%sign_dimid, nc%sign_varid), & - "swiftest_io_netcdf_open nf90_def_var sign_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%l_dimname, nc%out_type, nc%l_dimid, nc%l_varid), & - "swiftest_io_netcdf_open nf90_def_var l_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%m_dimname, nc%out_type, nc%m_dimid, nc%m_varid), & - "swiftest_io_netcdf_open nf90_def_var m_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" ) @@ -962,13 +957,37 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) 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" ) + 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" ) + nc%rot_varid), "netcdf_io_initialize_output nf90_def_var rot_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%rotphase_varname, nc%out_type, nc%time_dimid, nc%rotphase_varid), & "netcdf_io_initialize_output nf90_def_var rotphase_varid" ) + + + if (nc%lc_lm_exists) then + call netcdf_io_check( nf90_def_dim(nc%id, nc%sign_dimname, NF90_UNLIMITED, nc%sign_dimid), & + "swiftest_io_netcdf_open nf90_def_dim sign_dimid") + call netcdf_io_check( nf90_def_dim(nc%id, nc%l_dimname, NF90_UNLIMITED, nc%l_dimid), & + "swiftest_io_netcdf_open nf90_def_dim l_dimid") + call netcdf_io_check( nf90_def_dim(nc%id, nc%m_dimname, NF90_UNLIMITED, nc%m_dimid), & + "swiftest_io_netcdf_open nf90_def_dim m_dimid") + + call netcdf_io_check( nf90_def_var(nc%id, nc%sign_dimname, nc%out_type, nc%sign_dimid, nc%sign_varid), & + "swiftest_io_netcdf_open nf90_def_var sign_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%l_dimname, nc%out_type, nc%l_dimid, nc%l_varid), & + "swiftest_io_netcdf_open nf90_def_var l_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%m_dimname, nc%out_type, nc%m_dimid, nc%m_varid), & + "swiftest_io_netcdf_open nf90_def_var m_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%c_lm_varname, nc%out_type, [nc%m_dimid, nc%l_dimid, nc%sign_dimid], & + nc%c_lm_varid), "netcdf_io_initialize_output nf90_def_var c_lm_varid" ) + + else + 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" ) + end if + end if ! if (param%ltides) then @@ -1006,16 +1025,6 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) "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" ) - - if (nc%lc_lm_exists) then - call netcdf_io_check( nf90_def_var(nc%id, nc%c_lm_varname, nc%out_type, [nc%m_dimid, nc%l_dimid, nc%sign_dimid], & - nc%c_lm_varid), "netcdf_io_initialize_output nf90_def_var c_lm_varid" ) - end if - ! 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 @@ -1062,6 +1071,16 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) 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" ) + if (param%lrotation .and. nc%lc_lm_exists) then + ! Populate coordinate values for l, m, and sign + call netcdf_io_check( nf90_put_var(nc%id, nc%l_varid, [(i, i=0, nc%l_dim_max-1)], start=[1], count=[nc%l_dim_max]), & + "netcdf_io_netcdf_initialize_output nf90_put_var l_varid") + call netcdf_io_check( nf90_put_var(nc%id, nc%m_varid, [(i, i=0, nc%m_dim_max-1)], start=[1], count=[nc%m_dim_max]), & + "netcdf_io_netcdf_initialize_output nf90_put_var m_varid") + call netcdf_io_check( nf90_put_var(nc%id, nc%sign_varid, nc%sign_coords, start=[1], count=[2] ), & + "netcdf_io_netcdf_initialize_output nf90_put_var sign_varid") + end if + end associate return @@ -1354,7 +1373,6 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier integer(I4B) :: ierr !! Error code: returns 0 if the read is successful ! Internals integer(I4B) :: i, idmax, npl_check, ntp_check, str_max, status, npl, ntp - integer(I4B) :: l_dim_max, m_dim_max ! dimensions for c_lm array real(DP), dimension(:), allocatable :: rtemp real(DP), dimension(:,:), allocatable :: vectemp integer(I4B), dimension(:), allocatable :: itemp @@ -1527,8 +1545,8 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier ! Set initial central body angular momentum for bookkeeping cb%L0(:) = cb%Ip(3) * cb%mass * cb%R0**2 * cb%rot(:) - - ! rotphase may not be input by the user + + ! Check if rotphase is input by user. If not, set to 0 status = nf90_inq_varid(nc%id, nc%rotphase_varname, nc%rotphase_varid) if (status == NF90_NOERR) then call netcdf_io_check( nf90_get_var(nc%id, nc%rotphase_varid, cb%rotphase, start=[tslot]), & @@ -1537,6 +1555,53 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier else cb%rotphase = 0.0_DP end if + + ! Check if the J2 and J4 terms are input by user. If not, set them to 0 + 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 + + ! Check if gravitational harmonics coefficients, c_lm, are set by the user. Also ensure that + ! c_lm and j2rp2/j4rp4 are not both set. + status = nf90_inq_varid(nc%id, nc%c_lm_varname, nc%c_lm_varid) + if (status == NF90_NOERR) then + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%l_dimid, len = nc%l_dim_max),& + "netcdf_io_read_frame_system nf90_inquire_dimension l_dimid" ) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%m_dimid, len = nc%m_dim_max), & + "netcdf_io_read_frame_system nf90_inquire_dimension m_dimid") + + if (nc%l_dim_max /= nc%m_dim_max) then + write(*,*) "Error reading in NetCDF file: c_lm requires l_dim_max = m_dim_max" + call base_util_exit(FAILURE,param%display_unit) + end if + + if(.not. allocated(cb%c_lm)) allocate(cb%c_lm(nc%m_dim_max, nc%l_dim_max, 2)) + call netcdf_io_check( nf90_get_var(nc%id, nc%c_lm_varid, cb%c_lm, count = [nc%m_dim_max, nc%l_dim_max, 2]), & + "netcdf_io_read_frame_system nf90_getvar c_lm_varid") + + if (abs(cb%j2rp2) > tiny(1.0_DP) .or. (abs(cb%j4rp4) > tiny(1.0_DP))) then + write(*,*) "Error reading in NetCDF file: cannot use both c_lm and j2rp2/j4rp4" + call base_util_exit(FAILURE,param%display_unit) + end if + + nc%lc_lm_exists = .true. + else + if (allocated(cb%c_lm)) deallocate(cb%c_lm) + nc%lc_lm_exists = .false. + end if + end if ! if (param%ltides) then @@ -1551,39 +1616,6 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier ! 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 - - status = nf90_inq_varid(nc%id, nc%c_lm_varname, nc%c_lm_varid) - if (status == NF90_NOERR) then - call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%l_dimid, len = l_dim_max), "netcdf_io_read_frame_system nf90_inquire_dimension l_dimid" ) - call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%m_dimid, len = m_dim_max), "netcdf_io_read_frame_system nf90_inquire_dimension m_dimid") - - if(.not. allocated(cb%c_lm)) then - allocate(cb%c_lm(m_dim_max, l_dim_max, 2)) - end if - call netcdf_io_check( nf90_get_var(nc%id, nc%c_lm_varid, cb%c_lm, count = [m_dim_max, l_dim_max, 2]), "netcdf_io_read_frame_system nf90_getvar c_lm_varid") - - ! ordering of dimensions above seen to stackoverflow to prevent error 'NetCDF: Start + count exceeds dimension bound' - nc%lc_lm_exists = .true. - else - if (allocated(cb%c_lm)) deallocate(cb%c_lm) - nc%lc_lm_exists = .false. - end if - call self%read_particle_info(nc, param, plmask, tpmask) if (param%in_form == "EL") then @@ -2097,7 +2129,6 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: idslot, old_mode, tmp, i - integer(I4B) :: l_dim_max, m_dim_max integer(I4B), dimension(:), allocatable :: lm_coords integer(I4B) :: status @@ -2119,10 +2150,7 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) "swiftest_io_netcdf_write_frame_cb nf90_put_var cb mass_varid" ) if (param%lclose) call netcdf_io_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_io_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_io_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_io_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" ) @@ -2130,30 +2158,17 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) count=[NDIM,1,1]), & "swiftest_io_netcdf_write_frame_cb nf90_put_var cb rot_varid" ) - ! Following the template of j2rp2 call netcdf_io_check( nf90_put_var(nc%id, nc%rotphase_varid, self%rotphase * RAD2DEG, start = [tslot]), & "swiftest_io_netcdf_write_frame_cb nf90_put_var cb rotphase") - end if - - if (allocated(self%c_lm)) then - status = nf90_inq_varid(nc%id, nc%c_lm_varname, nc%c_lm_varid) - if (status == NF90_NOERR) then - m_dim_max = size(self%c_lm, 1) - l_dim_max = size(self%c_lm, 2) - - ! Populate coordinate values for l and m and export to hdf file - allocate(lm_coords(l_dim_max)) - do i = 0, l_dim_max - 1 - lm_coords(i + 1) = i - end do - - call netcdf_io_check( nf90_put_var(nc%id, nc%l_varid, lm_coords), "netcdf_io_write_frame_cb nf90_put_var l_varid") - call netcdf_io_check( nf90_put_var(nc%id, nc%m_varid, lm_coords), "netcdf_io_write_frame_cb nf90_put_var m_varid") - call netcdf_io_check( nf90_put_var(nc%id, nc%sign_varid, [1,-1]), "netcdf_io_write_frame_cb nf90_put_var sign_varid") - - ! Write dimension-coordinates to file - call netcdf_io_check( nf90_put_var(nc%id, nc%c_lm_varid, self%c_lm, count = [m_dim_max, l_dim_max, 2]), & - "netcdf_io_write_frame_cb nf90_put_var c_lm_varid") + + if (nc%lc_lm_exists) then + call netcdf_io_check( nf90_put_var(nc%id, nc%c_lm_varid, self%c_lm, count = [nc%m_dim_max, nc%l_dim_max, 2]), & + "netcdf_io_write_frame_cb nf90_put_var c_lm_varid") + else + call netcdf_io_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_io_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" ) end if end if