From a79319aa9a5034b7e00bec52f7561357879fb607 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 1 Mar 2024 11:44:15 -0500 Subject: [PATCH] Rearranged file io to better manage how and when gravitational harmonics dimensions, coordinates, and variables are defined and stored --- src/netcdf_io/netcdf_io_module.f90 | 46 +++++---- src/swiftest/swiftest_io.f90 | 156 +++++++++++++++-------------- 2 files changed, 110 insertions(+), 92 deletions(-) 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_io.f90 b/src/swiftest/swiftest_io.f90 index c50345cd3..8209f2a34 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,17 @@ 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 and m and export to hdf file + call netcdf_io_check( nf90_put_var(nc%id, nc%l_varid, [(i, i=0, nc%l_dim_max-1)]), & + "netcdf_io_write_frame_cb 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)]), & + "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") + end if + end associate return @@ -1354,7 +1374,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 @@ -1569,16 +1588,20 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier 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") + 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(.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") + 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") - ! ordering of dimensions above seen to stackoverflow to prevent error 'NetCDF: Start + count exceeds dimension bound' nc%lc_lm_exists = .true. + 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 else if (allocated(cb%c_lm)) deallocate(cb%c_lm) nc%lc_lm_exists = .false. @@ -2097,7 +2120,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 +2141,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 +2149,21 @@ 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") + + if (nc%lc_lm_exists) then + status = nf90_inq_varid(nc%id, nc%c_lm_varname, nc%c_lm_varid) + if (status == NF90_NOERR) then - ! 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") + 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 end if