From 184c70ad03e3c97d7a3ebb6046b79d582aa5938e Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 31 Jan 2024 10:34:59 -0700 Subject: [PATCH] First version with PIO domain write working: - Generate decompositions. - Transform domain_write routines to PIO. - Trasform more netcdf_write routines. - fixes in pio namelist handling. --- fms2_pio_io/fms_netcdf_domain_io.F90 | 228 ++++++ fms2_pio_io/include/domain_write.inc | 833 +++++----------------- fms2_pio_io/include/netcdf_write_data.inc | 124 ++-- fms2_pio_io/netcdf_io.F90 | 80 ++- 4 files changed, 517 insertions(+), 748 deletions(-) diff --git a/fms2_pio_io/fms_netcdf_domain_io.F90 b/fms2_pio_io/fms_netcdf_domain_io.F90 index cc437808b..8d8464539 100644 --- a/fms2_pio_io/fms_netcdf_domain_io.F90 +++ b/fms2_pio_io/fms_netcdf_domain_io.F90 @@ -29,6 +29,12 @@ module fms_netcdf_domain_io_mod use fms_io_utils_mod use netcdf_io_mod use platform_mod + +use pio, only : File_desc_t, var_desc_t, iosystem_desc_t, IO_desc_t +use pio, only : PIO_initdecomp, PIO_write_darray, pio_inq_varndims +use pio, only : PIO_INT, PIO_REAL, PIO_DOUBLE, PIO_OFFSET_KIND +use pio, only : PIO_setframe + implicit none private @@ -70,6 +76,19 @@ module fms_netcdf_domain_io_mod !! for domain-decomposed read. endtype FmsNetcdfDomainFile_t +! A type for encapsulating IOdesc instances and their distinct attributes +type PIO_decomp + type(IO_desc_t) :: iodesc + integer :: basetype_nf + integer :: unlim_dim_index + integer :: ndims + integer, dimension(:), allocatable :: dshape + type(PIO_decomp), pointer :: next => null() +end type PIO_decomp + +! head of PIO_decomp linked list +type(PIO_decomp), pointer :: decomps_head => null() + public :: open_domain_file public :: close_domain_file @@ -900,6 +919,215 @@ subroutine get_mosaic_tile_grid(grid_file,mosaic_file, domain, tile_count) end subroutine get_mosaic_tile_grid +function get_decomp(fileobj, variable_name, vardesc, vdata_shape, basetype_nf, ndims, xdim_index, ydim_index, unlim_dim_level) & + result(decomp) + + type(FmsNetcdfDomainFile_t), intent(in) :: fileobj !< File object. + character(len=*), intent(in) :: variable_name !< Variable name. + type(var_desc_t), intent(inout) :: vardesc + integer, dimension(:), intent(in) :: vdata_shape + integer, intent(in) :: basetype_nf + integer, intent(in) :: ndims + integer, intent(in) :: xdim_index + integer, intent(in) :: ydim_index + integer, intent(in), optional :: unlim_dim_level + type(PIO_decomp), pointer :: decomp + ! local + logical :: decomp_found + type(PIO_decomp), pointer :: decomps_tail => null() + integer :: global_size + integer :: local_size + integer :: ioff, joff + integer :: n + integer, dimension(:), allocatable :: dof + integer :: isc, iec, jsc, jec + integer :: ni, nj, nig, njg + integer :: i,j,k,l + integer :: unlim_dim_index + + !print *, "dbg-pio: looking for decomp ", basetype_nf, ndims + + ! TODO, calling get_variable_unlimited_dimension_index everytime get_decomp gets called + ! is inefficient, find a more efficient workaround. + unlim_dim_index = get_variable_unlimited_dimension_index(fileobj, variable_name) + + decomp => decomps_head + decomp_found = .false. + do while ( associated(decomp) ) + if (basetype_nf == decomp%basetype_nf .and. & + unlim_dim_index == decomp%unlim_dim_index .and. & + ndims == decomp%ndims) then + if (all(vdata_shape==decomp%dshape)) then + !print *, " dbg-pio: found decomp" + decomp_found = .true. + exit + endif + endif + !print *, " iter..." + decomps_tail => decomp + decomp => decomp%next + enddo + + ! If the compatible decomposition is found, return it + if (decomp_found) then + ! print *, "dbg-pio: found the decomp ", basetype_nf, ndims + else ! decomp not found + ! instantiate a new decomp and return it: + !print *, "dbg-pio: couldn't find the decomposition, creating a new one" + decomp => create_decomp(fileobj, variable_name, vdata_shape, basetype_nf, ndims, xdim_index, ydim_index) + + ! add new decomp to the end of decomp linked list: + if (.not. associated(decomps_head)) then + ! inialize the empty linked list + decomps_head => decomp + else + ! add the new decomp to the end of the list + decomps_tail%next => decomp + endif + endif + + ! set frame + if (unlim_dim_index /= -1) then + if (.not. present(unlim_dim_level)) call error("Must provide unlim_dim_level when calling get_decomp for a variable with unlimited dimension.") + call PIO_setframe(fileobj%file_desc, vardesc, int(unlim_dim_level, kind=PIO_OFFSET_KIND)) + endif + +end function get_decomp + +function create_decomp(fileobj, variable_name, vdata_shape, basetype_nf, ndims, xdim_index, ydim_index) & + result(decomp) + + type(FmsNetcdfDomainFile_t), intent(in) :: fileobj !< File object. + character(len=*), intent(in) :: variable_name !< Variable name. + integer, dimension(:), intent(in) :: vdata_shape + integer, intent(in) :: basetype_nf + integer, intent(in) :: ndims + integer, intent(in) :: xdim_index + integer, intent(in) :: ydim_index + type(PIO_decomp), pointer :: decomp + ! local + integer :: global_size + integer :: local_size + integer :: ioff, joff + integer :: n + integer, dimension(:), allocatable :: dof + integer :: isc, iec, jsc, jec + integer :: ni, nj, nig, njg + integer :: i,j,k,l + integer :: unlim_dim_index + + unlim_dim_index = get_variable_unlimited_dimension_index(fileobj, variable_name) + + if (unlim_dim_index /= -1) then + if (unlim_dim_index /= ndims) call error("unlimited dimension must be the slowest varying dimension") + endif + + call mpp_get_compute_domain(fileobj%domain, xbegin=isc, xend=iec, xsize=ni) + call mpp_get_compute_domain(fileobj%domain, ybegin=jsc, yend=jec, ysize=nj) + call mpp_get_global_domain (fileobj%domain, xsize = nig, ysize = njg) + + ! TODO: check if this subroutine works as expected when buffer_includes_halos + + ! initialize decomp and its iodesc member + allocate(decomp) + decomp%basetype_nf = basetype_nf + decomp%unlim_dim_index = unlim_dim_index + decomp%ndims = ndims + allocate(decomp%dshape(ndims)); decomp%dshape(:) = vdata_shape(:) + + + ioff = 1 - isc + joff = 1 - jsc + call mpp_max(ioff) + call mpp_max(joff) + + if (xdim_index /= 1 .or. ydim_index /= 2) then + call error("Dimension index ordering assumptions violated. Need to update FMS PIO code.") + endif + + if (ndims==2 .or. (ndims==3 .and. unlim_dim_index==3)) then + local_size = vdata_shape(1) * vdata_shape(2) + global_size = nig * njg + allocate(dof(local_size)) + + n = 1 + do j=jsc+joff,jec+joff + do i=isc+ioff,iec+ioff + dof(n) = i + (j-1)*nig + n = n+1 + enddo + enddo + + ! sanity check: + if (any(dof<1) .or. any(dof>global_size)) then + print *, "global size: ", global_size, & + "minval(dof) ", minval(dof), & + "maxval(dof) ", maxval(dof) + call error('error in dof construction') + endif + + call PIO_initdecomp(pio_iosystem, basetype_nf, (/nig, njg/), dof, decomp%iodesc) + deallocate(dof) + + else if (ndims==3 .or. (ndims==4 .and. unlim_dim_index==4)) then + local_size = vdata_shape(1) * vdata_shape(2) * vdata_shape(3) + global_size = nig * njg * vdata_shape(3) + allocate(dof(local_size)) + + n = 1 + do k=1, vdata_shape(3) + do j=jsc+joff,jec+joff + do i=isc+ioff,iec+ioff + dof(n) = i + (j-1)*nig + (k-1)*nig*njg + n = n+1 + enddo + enddo + enddo + + ! sanity check: + if (any(dof<1) .or. any(dof>global_size)) then + print *, "global size: ", global_size, & + "minval(dof) ", minval(dof), & + "maxval(dof) ", maxval(dof) + call error('error in dof construction') + endif + + call PIO_initdecomp(pio_iosystem, basetype_nf, (/nig, njg, vdata_shape(3)/), dof, decomp%iodesc) + deallocate(dof) + + else if (ndims==4 .or. (ndims==5 .and. unlim_dim_index==5)) then + local_size = vdata_shape(1) * vdata_shape(2) * vdata_shape(3) * vdata_shape(4) + global_size = nig * njg * vdata_shape(3) * vdata_shape(4) + allocate(dof(local_size)) + + n = 1 + do l=1, vdata_shape(4) + do k=1, vdata_shape(3) + do j=jsc+joff,jec+joff + do i=isc+ioff,iec+ioff + dof(n) = i + (j-1)*nig + (k-1)*nig*njg + (l-1)*nig*njg*vdata_shape(3) + n = n+1 + enddo + enddo + enddo + enddo + + ! sanity check: + if (any(dof<1) .or. any(dof>global_size)) then + print *, "global size: ", global_size, & + "minval(dof) ", minval(dof), & + "maxval(dof) ", maxval(dof) + call error('error in dof construction') + endif + + call PIO_initdecomp(pio_iosystem, basetype_nf, (/nig, njg, vdata_shape(3), vdata_shape(4) /), dof, decomp%iodesc) + deallocate(dof) + else + call error("Unsupported number of dimensions encountered in get_decomp.") + endif + +end function create_decomp + include "register_domain_restart_variable.inc" include "domain_read.inc" include "domain_write.inc" diff --git a/fms2_pio_io/include/domain_write.inc b/fms2_pio_io/include/domain_write.inc index f27c46900..f54a01918 100644 --- a/fms2_pio_io/include/domain_write.inc +++ b/fms2_pio_io/include/domain_write.inc @@ -140,8 +140,11 @@ subroutine domain_write_2d(fileobj, variable_name, vdata, unlim_dim_level, & integer :: xgmin !< Starting x index of the global io domain integer :: ygmax !< Ending y index of the global io domain integer :: ygmin !< Ending y index of the global io domain - - print *, "ERRORline domain_write_2d"; call error("PIO version not implemented!!!") + integer :: basetype_nf + type(PIO_decomp), pointer :: decomp + integer :: varid + type(var_desc_t) :: vardesc + integer :: err if (.not. is_variable_domain_decomposed(fileobj, variable_name, .true., & xdim_index, ydim_index, xpos, ypos)) then @@ -150,231 +153,53 @@ subroutine domain_write_2d(fileobj, variable_name, vdata, unlim_dim_level, & edge_lengths=edge_lengths) return endif - io_domain => mpp_get_io_domain(fileobj%domain) - call domain_offsets(size(vdata, xdim_index), size(vdata, ydim_index), fileobj%domain, & - xpos, ypos, isd, isc, xc_size, jsd, jsc, yc_size, & - buffer_includes_halos, msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name)) - c(:) = 1 - e(:) = shape(vdata) - - !I/O root gathers the data and writes it. - if (fileobj%is_root) then - allocate(pe_isc(size(fileobj%pelist))) - allocate(pe_iec(size(fileobj%pelist))) - allocate(pe_icsize(size(fileobj%pelist))) - call mpp_get_compute_domains(io_domain, xbegin=pe_isc, xend=pe_iec, xsize=pe_icsize, & - position=xpos) - allocate(pe_jsc(size(fileobj%pelist))) - allocate(pe_jec(size(fileobj%pelist))) - allocate(pe_jcsize(size(fileobj%pelist))) - call mpp_get_compute_domains(io_domain, ybegin=pe_jsc, yend=pe_jec, ysize=pe_jcsize, & - position=ypos) - - !< Determine the size of the global io domain - call mpp_get_global_domain(io_domain, ybegin=ygmin, yend=ygmax, position=ypos) - call mpp_get_global_domain(io_domain, xbegin=xgmin, xend=xgmax, position=xpos) - !Write the out the data. - !< Set e to equal the size of the global io domain - e(xdim_index) = xgmax - xgmin + 1 - e(ydim_index) = ygmax - ygmin + 1 - - !< Allocate a global buffer, get the fill value if it exists in the file, and initialize - !! the buffer to the fill value - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(global_buf_i4_kind, e) - global_buf_i4_kind = 0 - if (get_fill_value(fileobj, variable_name, fill_i4_kind, broadcast=.false.)) then - global_buf_i4_kind = fill_i4_kind - endif - type is (integer(kind=i8_kind)) - call allocate_array(global_buf_i8_kind, e) - global_buf_i8_kind = 0 - if (get_fill_value(fileobj, variable_name, fill_i8_kind, broadcast=.false.)) then - global_buf_i8_kind = fill_i8_kind - endif - type is (real(kind=r4_kind)) - call allocate_array(global_buf_r4_kind, e) - global_buf_r4_kind = 0. - if (get_fill_value(fileobj, variable_name, fill_r4_kind, broadcast=.false.)) then - global_buf_r4_kind = fill_r4_kind - endif - type is (real(kind=r8_kind)) - call allocate_array(global_buf_r8_kind, e) - global_buf_r8_kind = 0. - if (get_fill_value(fileobj, variable_name, fill_r8_kind, broadcast=.false.)) then - global_buf_r8_kind = fill_r8_kind - endif - class default - call error("unsupported variable type: domain_write_2d: file: "//trim(fileobj%path)//" variable:"// & + select type(vdata) + type is (integer(kind=i4_kind)) + basetype_nf = PIO_INT + type is (real(kind=r4_kind)) + basetype_nf = PIO_REAL + type is (real(kind=r8_kind)) + basetype_nf = PIO_DOUBLE + class default + call error("unsupported variable type: domain_write_4d: file: "//trim(fileobj%path)//" variable:"// & & trim(variable_name)) - end select - - do i = 1, size(fileobj%pelist) - !< Set c relative to the starting global io domain index - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - e(xdim_index) = pe_icsize(i) - e(ydim_index) = pe_jcsize(i) - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(buf_i4_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_i4_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_i4_kind, size(buf_i4_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_i4_kind, global_buf_i4_kind, c, e) - deallocate(buf_i4_kind) - type is (integer(kind=i8_kind)) - call allocate_array(buf_i8_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_i8_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_i8_kind, size(buf_i8_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_i8_kind, global_buf_i8_kind, c, e) - deallocate(buf_i8_kind) - type is (real(kind=r4_kind)) - call allocate_array(buf_r4_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_r4_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_r4_kind, size(buf_r4_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_r4_kind, global_buf_r4_kind, c, e) - deallocate(buf_r4_kind) - type is (real(kind=r8_kind)) - call allocate_array(buf_r8_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_r8_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_r8_kind, size(buf_r8_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_r8_kind, global_buf_r8_kind, c, e) - deallocate(buf_r8_kind) - end select - enddo - deallocate(pe_isc) - deallocate(pe_iec) - deallocate(pe_icsize) - deallocate(pe_jsc) - deallocate(pe_jec) - deallocate(pe_jcsize) + end select + + ! get varid and construct a temporary vardesc to pass to PIO routines + varid = get_variable_id(fileobj%ncid, trim(variable_name), & + msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name)) + vardesc%varID = varid + vardesc%ncid = fileobj%ncid + + decomp => get_decomp(fileobj, variable_name, vardesc, shape(vdata), basetype_nf, 2, xdim_index, ydim_index, unlim_dim_level) + + err = 0 + select type(vdata) + type is (integer(kind=i4_kind)) + if (get_fill_value(fileobj, variable_name, fill_i4_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_i4_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + type is (real(kind=r4_kind)) + if (get_fill_value(fileobj, variable_name, fill_r4_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_r4_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + type is (real(kind=r8_kind)) + if (get_fill_value(fileobj, variable_name, fill_r8_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_r8_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + class default + call error("unsupported variable type: domain_write_4d: file: "//trim(fileobj%path)//" variable:"// & + & trim(variable_name)) + end select + call check_netcdf_code(err, "domain_write.inc") - !Write the out the data. - select type(vdata) - type is (integer(kind=i4_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_i4_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_i4_kind) - type is (integer(kind=i8_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_i8_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_i8_kind) - type is (real(kind=r4_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_r4_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_r4_kind) - type is (real(kind=r8_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_r8_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_r8_kind) - end select - else - if (buffer_includes_halos) then - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - endif - e(xdim_index) = xc_size - e(ydim_index) = yc_size - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(buf_i4_kind, e) - call get_array_section(buf_i4_kind, vdata, c, e) - call mpp_send(buf_i4_kind, size(buf_i4_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_i4_kind) - type is (integer(kind=i8_kind)) - call allocate_array(buf_i8_kind, e) - call get_array_section(buf_i8_kind, vdata, c, e) - call mpp_send(buf_i8_kind, size(buf_i8_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_i8_kind) - type is (real(kind=r4_kind)) - call allocate_array(buf_r4_kind, e) - call get_array_section(buf_r4_kind, vdata, c, e) - call mpp_send(buf_r4_kind, size(buf_r4_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_r4_kind) - type is (real(kind=r8_kind)) - call allocate_array(buf_r8_kind, e) - call get_array_section(buf_r8_kind, vdata, c, e) - call mpp_send(buf_r8_kind, size(buf_r8_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_r8_kind) - class default - call error("unsupported variable type: domain_write_2d: file: "//trim(fileobj%path)//" variable:"// & - & trim(variable_name)) - end select - endif end subroutine domain_write_2d @@ -438,8 +263,11 @@ subroutine domain_write_3d(fileobj, variable_name, vdata, unlim_dim_level, & integer :: xgmin !< Starting x index of the global io domain integer :: ygmax !< Ending y index of the global io domain integer :: ygmin !< Ending y index of the global io domain - - print *, "ERRORline domain_write_3d"; call error("PIO version not implemented!!!") + integer :: basetype_nf + integer :: varid + integer :: err + type(var_desc_t) :: vardesc + type(PIO_decomp), pointer :: decomp if (.not. is_variable_domain_decomposed(fileobj, variable_name, .true., & xdim_index, ydim_index, xpos, ypos)) then @@ -448,232 +276,54 @@ subroutine domain_write_3d(fileobj, variable_name, vdata, unlim_dim_level, & edge_lengths=edge_lengths) return endif - io_domain => mpp_get_io_domain(fileobj%domain) - call domain_offsets(size(vdata, xdim_index), size(vdata, ydim_index), fileobj%domain, & - xpos, ypos, isd, isc, xc_size, jsd, jsc, yc_size, & - buffer_includes_halos, msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name)) - c(:) = 1 - e(:) = shape(vdata) - !I/O root gathers the data and writes it. - if (fileobj%is_root) then - allocate(pe_isc(size(fileobj%pelist))) - allocate(pe_iec(size(fileobj%pelist))) - allocate(pe_icsize(size(fileobj%pelist))) - call mpp_get_compute_domains(io_domain, xbegin=pe_isc, xend=pe_iec, xsize=pe_icsize, & - position=xpos) - allocate(pe_jsc(size(fileobj%pelist))) - allocate(pe_jec(size(fileobj%pelist))) - allocate(pe_jcsize(size(fileobj%pelist))) - call mpp_get_compute_domains(io_domain, ybegin=pe_jsc, yend=pe_jec, ysize=pe_jcsize, & - position=ypos) + select type(vdata) + type is (integer(kind=i4_kind)) + basetype_nf = PIO_INT + type is (real(kind=r4_kind)) + basetype_nf = PIO_REAL + type is (real(kind=r8_kind)) + basetype_nf = PIO_DOUBLE + class default + call error("unsupported variable type: domain_write_4d: file: "//trim(fileobj%path)//" variable:"// & + & trim(variable_name)) + end select + + ! get varid and construct a temporary vardesc to pass to PIO routines + varid = get_variable_id(fileobj%ncid, trim(variable_name), & + msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name)) + vardesc%varID = varid + vardesc%ncid = fileobj%ncid + + decomp => get_decomp(fileobj, variable_name, vardesc, shape(vdata), basetype_nf, 3, xdim_index, ydim_index, unlim_dim_level) + + err = 0 + select type(vdata) + type is (integer(kind=i4_kind)) + if (get_fill_value(fileobj, variable_name, fill_i4_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_i4_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + type is (real(kind=r4_kind)) + if (get_fill_value(fileobj, variable_name, fill_r4_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_r4_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + type is (real(kind=r8_kind)) + if (get_fill_value(fileobj, variable_name, fill_r8_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_r8_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + class default + call error("unsupported variable type: domain_write_4d: file: "//trim(fileobj%path)//" variable:"// & + & trim(variable_name)) + end select + call check_netcdf_code(err, "domain_write.inc") - !< Determine the size of the global io domain - call mpp_get_global_domain(io_domain, ybegin=ygmin, yend=ygmax, position=ypos) - call mpp_get_global_domain(io_domain, xbegin=xgmin, xend=xgmax, position=xpos) - - !Write the out the data. - !< Set e to equal the size of the global io domain - e(xdim_index) = xgmax - xgmin + 1 - e(ydim_index) = ygmax - ygmin + 1 - - !< Allocate a global buffer, get the fill value if it exists in the file, and initialize - !! the buffer to the fill value - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(global_buf_i4_kind, e) - global_buf_i4_kind = 0 - if (get_fill_value(fileobj, variable_name, fill_i4_kind, broadcast=.false.)) then - global_buf_i4_kind = fill_i4_kind - endif - type is (integer(kind=i8_kind)) - call allocate_array(global_buf_i8_kind, e) - global_buf_i8_kind = 0 - if (get_fill_value(fileobj, variable_name, fill_i8_kind, broadcast=.false.)) then - global_buf_i8_kind = fill_i8_kind - endif - type is (real(kind=r4_kind)) - call allocate_array(global_buf_r4_kind, e) - global_buf_r4_kind = 0. - if (get_fill_value(fileobj, variable_name, fill_r4_kind, broadcast=.false.)) then - global_buf_r4_kind = fill_r4_kind - endif - type is (real(kind=r8_kind)) - call allocate_array(global_buf_r8_kind, e) - global_buf_r8_kind = 0. - if (get_fill_value(fileobj, variable_name, fill_r8_kind, broadcast=.false.)) then - global_buf_r8_kind = fill_r8_kind - endif - class default - call error("unsupported variable type: domain_write_3d: file: "//trim(fileobj%path)//" variable:"// & - & trim(variable_name)) - end select - - do i = 1, size(fileobj%pelist) - !< Set c relative to the starting global io domain index - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - e(xdim_index) = pe_icsize(i) - e(ydim_index) = pe_jcsize(i) - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(buf_i4_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_i4_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_i4_kind, size(buf_i4_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_i4_kind, global_buf_i4_kind, c, e) - deallocate(buf_i4_kind) - type is (integer(kind=i8_kind)) - call allocate_array(buf_i8_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_i8_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_i8_kind, size(buf_i8_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_i8_kind, global_buf_i8_kind, c, e) - deallocate(buf_i8_kind) - type is (real(kind=r4_kind)) - call allocate_array(buf_r4_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_r4_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_r4_kind, size(buf_r4_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_r4_kind, global_buf_r4_kind, c, e) - deallocate(buf_r4_kind) - type is (real(kind=r8_kind)) - call allocate_array(buf_r8_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_r8_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_r8_kind, size(buf_r8_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_r8_kind, global_buf_r8_kind, c, e) - deallocate(buf_r8_kind) - end select - enddo - deallocate(pe_isc) - deallocate(pe_iec) - deallocate(pe_icsize) - deallocate(pe_jsc) - deallocate(pe_jec) - deallocate(pe_jcsize) - - !Write the out the data. - select type(vdata) - type is (integer(kind=i4_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_i4_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_i4_kind) - type is (integer(kind=i8_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_i8_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_i8_kind) - type is (real(kind=r4_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_r4_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_r4_kind) - type is (real(kind=r8_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_r8_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_r8_kind) - end select - else - if (buffer_includes_halos) then - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - endif - e(xdim_index) = xc_size - e(ydim_index) = yc_size - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(buf_i4_kind, e) - call get_array_section(buf_i4_kind, vdata, c, e) - call mpp_send(buf_i4_kind, size(buf_i4_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_i4_kind) - type is (integer(kind=i8_kind)) - call allocate_array(buf_i8_kind, e) - call get_array_section(buf_i8_kind, vdata, c, e) - call mpp_send(buf_i8_kind, size(buf_i8_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_i8_kind) - type is (real(kind=r4_kind)) - call allocate_array(buf_r4_kind, e) - call get_array_section(buf_r4_kind, vdata, c, e) - call mpp_send(buf_r4_kind, size(buf_r4_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_r4_kind) - type is (real(kind=r8_kind)) - call allocate_array(buf_r8_kind, e) - call get_array_section(buf_r8_kind, vdata, c, e) - call mpp_send(buf_r8_kind, size(buf_r8_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_r8_kind) - class default - call error("unsupported variable type: domain_write_3d: file: "//trim(fileobj%path)//" variable:"// & - & trim(variable_name)) - end select - endif -end subroutine domain_write_3d +end subroutine domain_write_3d !> @brief Gather "compute" domain data on the I/O root rank and then have @@ -736,6 +386,13 @@ subroutine domain_write_4d(fileobj, variable_name, vdata, unlim_dim_level, & integer :: xgmin !< Starting x index of the global io domain integer :: ygmax !< Ending y index of the global io domain integer :: ygmin !< Ending y index of the global io domain + integer :: basetype_nf + integer :: varid + integer :: err + integer :: ndims + type(var_desc_t) :: vardesc + type(PIO_decomp), pointer :: decomp + if (.not. is_variable_domain_decomposed(fileobj, variable_name, .true., & xdim_index, ydim_index, xpos, ypos)) then @@ -745,233 +402,63 @@ subroutine domain_write_4d(fileobj, variable_name, vdata, unlim_dim_level, & return endif - print *, "ERRORline domain_write_4d"; call error("PIO version not implemented!!!") - - io_domain => mpp_get_io_domain(fileobj%domain) - call domain_offsets(size(vdata, xdim_index), size(vdata, ydim_index), fileobj%domain, & - xpos, ypos, isd, isc, xc_size, jsd, jsc, yc_size, & - buffer_includes_halos, msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name)) - c(:) = 1 - e(:) = shape(vdata) - - !I/O root gathers the data and writes it. - if (fileobj%is_root) then - allocate(pe_isc(size(fileobj%pelist))) - allocate(pe_iec(size(fileobj%pelist))) - allocate(pe_icsize(size(fileobj%pelist))) - call mpp_get_compute_domains(io_domain, xbegin=pe_isc, xend=pe_iec, xsize=pe_icsize, & - position=xpos) - allocate(pe_jsc(size(fileobj%pelist))) - allocate(pe_jec(size(fileobj%pelist))) - allocate(pe_jcsize(size(fileobj%pelist))) - call mpp_get_compute_domains(io_domain, ybegin=pe_jsc, yend=pe_jec, ysize=pe_jcsize, & - position=ypos) - - !< Determine the size of the global io domain - call mpp_get_global_domain(io_domain, ybegin=ygmin, yend=ygmax, position=ypos) - call mpp_get_global_domain(io_domain, xbegin=xgmin, xend=xgmax, position=xpos) + varid = get_variable_id(fileobj%ncid, trim(variable_name), & + msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name)) - !Write the out the data. - !< Set e to equal the size of the global io domain - e(xdim_index) = xgmax - xgmin + 1 - e(ydim_index) = ygmax - ygmin + 1 + ! The diagnostics manager always calls domain_write_4d regardless of the actual ndims of var + ! So, make sure to execute the correct routine: + err = pio_inq_varndims(fileobj%ncid, varid, ndims) + if (ndims == 2) then + call domain_write_2d(fileobj, variable_name, vdata(:,:,1,1), unlim_dim_level) ; return + elseif (ndims == 3) then + call domain_write_3d(fileobj, variable_name, vdata(:,:,1,:), unlim_dim_level) ; return + endif - !< Allocate a global buffer, get the fill value if it exists in the file, and initialize - !! the buffer to the fill value - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(global_buf_i4_kind, e) - global_buf_i4_kind = 0 - if (get_fill_value(fileobj, variable_name, fill_i4_kind, broadcast=.false.)) then - global_buf_i4_kind = fill_i4_kind - endif - type is (integer(kind=i8_kind)) - call allocate_array(global_buf_i8_kind, e) - global_buf_i8_kind = 0 - if (get_fill_value(fileobj, variable_name, fill_i8_kind, broadcast=.false.)) then - global_buf_i8_kind = fill_i8_kind - endif - type is (real(kind=r4_kind)) - call allocate_array(global_buf_r4_kind, e) - global_buf_r4_kind = 0. - if (get_fill_value(fileobj, variable_name, fill_r4_kind, broadcast=.false.)) then - global_buf_r4_kind = fill_r4_kind - endif - type is (real(kind=r8_kind)) - call allocate_array(global_buf_r8_kind, e) - global_buf_r8_kind = 0. - if (get_fill_value(fileobj, variable_name, fill_r8_kind, broadcast=.false.)) then - global_buf_r8_kind = fill_r8_kind - endif - class default - call error("unsupported variable type: domain_write_4d: file: "//trim(fileobj%path)//" variable:"// & - & trim(variable_name)) - end select - do i = 1, size(fileobj%pelist) - !< Set c relative to the starting global io domain index - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - e(xdim_index) = pe_icsize(i) - e(ydim_index) = pe_jcsize(i) - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(buf_i4_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_i4_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_i4_kind, size(buf_i4_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_i4_kind, global_buf_i4_kind, c, e) - deallocate(buf_i4_kind) - type is (integer(kind=i8_kind)) - call allocate_array(buf_i8_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_i8_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_i8_kind, size(buf_i8_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_i8_kind, global_buf_i8_kind, c, e) - deallocate(buf_i8_kind) - type is (real(kind=r4_kind)) - call allocate_array(buf_r4_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_r4_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_r4_kind, size(buf_r4_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_r4_kind, global_buf_r4_kind, c, e) - deallocate(buf_r4_kind) - type is (real(kind=r8_kind)) - call allocate_array(buf_r8_kind, e) - !Get the data for fileobj%pelist(i)'s portion of the compute domain. - if (i .eq. 1) then - !Root rank gets the data directly. - if (buffer_includes_halos) then - !Adjust if the input buffer has room for halos. - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - else - c(xdim_index) = 1 - c(ydim_index) = 1 - endif - call get_array_section(buf_r8_kind, vdata, c, e) - c(xdim_index) = pe_isc(i) - xgmin + 1 - c(ydim_index) = pe_jsc(i) - ygmin + 1 - else - !Receive data from non-root ranks. - call mpp_recv(buf_r8_kind, size(buf_r8_kind), fileobj%pelist(i), block=.true.) - endif - !Put local data into the global buffer. - call put_array_section(buf_r8_kind, global_buf_r8_kind, c, e) - deallocate(buf_r8_kind) - end select - enddo - deallocate(pe_isc) - deallocate(pe_iec) - deallocate(pe_icsize) - deallocate(pe_jsc) - deallocate(pe_jec) - deallocate(pe_jcsize) + select type(vdata) + type is (integer(kind=i4_kind)) + basetype_nf = PIO_INT + type is (real(kind=r4_kind)) + basetype_nf = PIO_REAL + type is (real(kind=r8_kind)) + basetype_nf = PIO_DOUBLE + class default + call error("unsupported variable type: domain_write_4d: file: "//trim(fileobj%path)//" variable:"// & + & trim(variable_name)) + end select + + ! set the members of temporary vardesc instance, which is passed when calling PIO_write_darray + vardesc%varID = varid + vardesc%ncid = fileobj%ncid + + decomp => get_decomp(fileobj, variable_name, vardesc, shape(vdata), basetype_nf, 4, xdim_index, ydim_index, unlim_dim_level) + + err = 0 + select type(vdata) + type is (integer(kind=i4_kind)) + if (get_fill_value(fileobj, variable_name, fill_i4_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_i4_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + type is (real(kind=r4_kind)) + if (get_fill_value(fileobj, variable_name, fill_r4_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_r4_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + type is (real(kind=r8_kind)) + if (get_fill_value(fileobj, variable_name, fill_r8_kind, broadcast=.false.)) then + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err, fillval=fill_r8_kind) + else + call PIO_write_darray(fileobj%file_desc, vardesc, decomp%iodesc, vdata, err) + endif + class default + call error("unsupported variable type: domain_write_4d: file: "//trim(fileobj%path)//" variable:"// & + & trim(variable_name)) + end select + call check_netcdf_code(err, "domain_write.inc") - !Write the out the data. - select type(vdata) - type is (integer(kind=i4_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_i4_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_i4_kind) - type is (integer(kind=i8_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_i8_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_i8_kind) - type is (real(kind=r4_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_r4_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_r4_kind) - type is (real(kind=r8_kind)) - call netcdf_write_data(fileobj, variable_name, global_buf_r8_kind, & - unlim_dim_level=unlim_dim_level) - deallocate(global_buf_r8_kind) - end select - else - if (buffer_includes_halos) then - c(xdim_index) = isc - isd + 1 - c(ydim_index) = jsc - jsd + 1 - endif - e(xdim_index) = xc_size - e(ydim_index) = yc_size - select type(vdata) - type is (integer(kind=i4_kind)) - call allocate_array(buf_i4_kind, e) - call get_array_section(buf_i4_kind, vdata, c, e) - call mpp_send(buf_i4_kind, size(buf_i4_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_i4_kind) - type is (integer(kind=i8_kind)) - call allocate_array(buf_i8_kind, e) - call get_array_section(buf_i8_kind, vdata, c, e) - call mpp_send(buf_i8_kind, size(buf_i8_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_i8_kind) - type is (real(kind=r4_kind)) - call allocate_array(buf_r4_kind, e) - call get_array_section(buf_r4_kind, vdata, c, e) - call mpp_send(buf_r4_kind, size(buf_r4_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_r4_kind) - type is (real(kind=r8_kind)) - call allocate_array(buf_r8_kind, e) - call get_array_section(buf_r8_kind, vdata, c, e) - call mpp_send(buf_r8_kind, size(buf_r8_kind), fileobj%io_root) - call mpp_sync_self(check=event_send) - deallocate(buf_r8_kind) - class default - call error("unsupported variable type: domain_write_4d: file: "//trim(fileobj%path)//" variable:"// & - & trim(variable_name)) - end select - endif end subroutine domain_write_4d diff --git a/fms2_pio_io/include/netcdf_write_data.inc b/fms2_pio_io/include/netcdf_write_data.inc index 446cd5b81..51b756500 100644 --- a/fms2_pio_io/include/netcdf_write_data.inc +++ b/fms2_pio_io/include/netcdf_write_data.inc @@ -66,13 +66,25 @@ subroutine netcdf_write_data_0d(fileobj, variable_name, variable_data, unlim_dim varid = get_variable_id(fileobj%ncid, trim(variable_name), msg=append_error_msg) select type(variable_data) type is (integer(kind=i4_kind)) - err = pio_put_var(fileobj%file_desc, varid, variable_data) + if (present(unlim_dim_level)) then + err = pio_put_var(fileobj%file_desc, varid, start=c(1:1), count=[1], ival=[variable_data]) + else + err = pio_put_var(fileobj%file_desc, varid, variable_data) + endif type is (integer(kind=i8_kind)) call error(trim(fileobj%path)//": 64 bit integers are not supported with PIO.") type is (real(kind=r4_kind)) - err = pio_put_var(fileobj%file_desc, varid, variable_data) + if (present(unlim_dim_level)) then + err = pio_put_var(fileobj%file_desc, varid, start=c(1:1), count=[1], ival=[variable_data]) + else + err = pio_put_var(fileobj%file_desc, varid, variable_data) + endif type is (real(kind=r8_kind)) - err = pio_put_var(fileobj%file_desc, varid, variable_data) + if (present(unlim_dim_level)) then + err = pio_put_var(fileobj%file_desc, varid, start=c(1:1), count=[1], ival=[variable_data]) + else + err = pio_put_var(fileobj%file_desc, varid, variable_data) + endif type is (character(len=*)) call error(trim(fileobj%path)//": character write not implemented yet (netcdf_write_data.inc)") ndims = get_variable_num_dimensions(fileobj, variable_name, broadcast=.false.) @@ -101,7 +113,6 @@ subroutine netcdf_write_data_0d(fileobj, variable_name, variable_data, unlim_dim end select call check_netcdf_code(err, append_error_msg) !endif - print*, "dbg-wrote0d ", trim(variable_name) end subroutine netcdf_write_data_0d @@ -136,6 +147,7 @@ subroutine netcdf_write_data_1d(fileobj, variable_name, variable_data, unlim_dim integer :: j integer :: tlen character(len=200) :: append_error_msg !< Msg to be appended to FATAL error message + integer :: true_ndim append_error_msg = "netcdf_write_data_1d: file:"//trim(fileobj%path)//" variable: "//trim(variable_name) @@ -150,9 +162,14 @@ subroutine netcdf_write_data_1d(fileobj, variable_name, variable_data, unlim_dim else e(1:1) = shape(variable_data) endif + true_ndim = 1 if (present(unlim_dim_level)) then unlim_dim_index = get_variable_unlimited_dimension_index(fileobj, variable_name, & broadcast=.false.) + ! FMS doesn't count unlim (time) dimension towards the total number of variables dimension. + ! Thus, when the variable has an unlimited dimension, we call PIO appropriately, i.e., by + ! incrementing the number ndims to account for the time dimension. + true_ndim = 2 if (unlim_dim_index .ne. 2) then call error("unlimited dimension must be the slowest varying dimension: "//trim(append_error_msg)) endif @@ -162,13 +179,13 @@ subroutine netcdf_write_data_1d(fileobj, variable_name, variable_data, unlim_dim varid = get_variable_id(fileobj%ncid, trim(variable_name), msg=append_error_msg) select type(variable_data) type is (integer(kind=i4_kind)) - err = pio_put_var(fileobj%file_desc, varid, variable_data(c(1):e(1))) + err = pio_put_var(fileobj%file_desc, varid, start=c(1:true_ndim), count=e(1:true_ndim), ival=variable_data) type is (integer(kind=i8_kind)) call error(trim(fileobj%path)//": 64 bit integers are not supported with PIO.") type is (real(kind=r4_kind)) - err = pio_put_var(fileobj%file_desc, varid, variable_data(c(1):e(1))) + err = pio_put_var(fileobj%file_desc, varid, start=c(1:true_ndim), count=e(1:true_ndim), ival=variable_data) type is (real(kind=r8_kind)) - err = pio_put_var(fileobj%file_desc, varid, variable_data(c(1):e(1))) + err = pio_put_var(fileobj%file_desc, varid, start=c(1:true_ndim), count=e(1:true_ndim), ival=variable_data) type is (character(len=*)) call error(trim(fileobj%path)//": character write not implemented yet (netcdf_write_data.inc)") ndims = get_variable_num_dimensions(fileobj, variable_name, broadcast=.false.) @@ -362,44 +379,71 @@ subroutine netcdf_write_data_4d(fileobj, variable_name, variable_data, unlim_dim integer,dimension(5) :: c integer, dimension(5) :: e character(len=200) :: append_error_msg !< Msg to be appended to FATAL error message + integer :: ndim - print *, "ERRORline", 368, "write_data.inc"; call error("PIO version not implemented!!!") + unlim_dim_index = -1 + c(:) = 1 + if (present(corner)) then + c(1:4) = corner(:) + endif + e(:) = 1 + if (present(edge_lengths)) then + e(1:4) = edge_lengths(:) + else + e(1:4) = shape(variable_data) + endif + if (present(unlim_dim_level)) then + unlim_dim_index = get_variable_unlimited_dimension_index(fileobj, variable_name, & + broadcast=.false.) + if (unlim_dim_index /= -1) then + c(unlim_dim_index) = unlim_dim_level + endif + endif - append_error_msg = "netcdf_write_data_4d: file:"//trim(fileobj%path)//" variable: "//trim(variable_name) + ndim = get_variable_num_dimensions(fileobj, variable_name) - if (fileobj%is_root) then - c(:) = 1 - if (present(corner)) then - c(1:4) = corner(:) - endif - e(:) = 1 - if (present(edge_lengths)) then - e(1:4) = edge_lengths(:) - else - e(1:4) = shape(variable_data) + if (unlim_dim_index == -1) then + ! The field has NO time dimension + if (ndim == 1) then + if (size(variable_data,2)/=1 .or. size(variable_data,3)/=1 .or. size(variable_data,4)/=1) then + call error("PIO data size assertion violation in netcdf_write_data_4d") + endif + call netcdf_write_data_1d(fileobj, variable_name, variable_data(:,1,1,1), unlim_dim_level, c(1), e(1)) + return + elseif (ndim == 2) then + if (size(variable_data,3)/=1 .or. size(variable_data,4)/=1) then + call error("PIO data size assertion violation in netcdf_write_data_4d") + endif + call netcdf_write_data_2d(fileobj, variable_name, variable_data(:,:,1,1), unlim_dim_level, c(1:2), e(1:2)) + return + elseif (ndim == 3) then + if (size(variable_data,4)/=1) then + call error("PIO data size assertion violation in netcdf_write_data_4d") + endif + call netcdf_write_data_3d(fileobj, variable_name, variable_data(:,:,:,1), unlim_dim_level, c(1:3), e(1:3)) + return endif - if (present(unlim_dim_level)) then - unlim_dim_index = get_variable_unlimited_dimension_index(fileobj, variable_name, & - broadcast=.false.) - if (unlim_dim_index /= -1) then - c(unlim_dim_index) = unlim_dim_level - endif + else + ! The field has time dimension + if (ndim == 1) then + if (size(variable_data,1)/=1 .or. size(variable_data,2)/=1 .or. size(variable_data,3)/=1 .or. size(variable_data,4)/=1) then + call error("PIO data size assertion violation in netcdf_write_data_4d") + endif + call netcdf_write_data_0d(fileobj, variable_name, variable_data(1,1,1,1), unlim_dim_level, c(1)) + return + elseif (ndim == 2) then + if (size(variable_data,2)/=1 .or. size(variable_data,3)/=1 .or. size(variable_data,4)/=1) then + call error("PIO data size assertion violation in netcdf_write_data_4d") + endif + call netcdf_write_data_1d(fileobj, variable_name, variable_data(:,1,1,1), unlim_dim_level, c(1:2), e(1:2)) + return + elseif (ndim == 3) then + if (size(variable_data,3)/=1 .or. size(variable_data,4)/=1) then + call error("PIO data size assertion violation in netcdf_write_data_4d") + endif + call netcdf_write_data_2d(fileobj, variable_name, variable_data(:,:,1,1), unlim_dim_level, c(1:3), e(1:3)) + return endif - call set_netcdf_mode(fileobj%ncid, data_mode) - varid = get_variable_id(fileobj%ncid, trim(variable_name), msg=append_error_msg) - select type(variable_data) - type is (integer(kind=i4_kind)) - err = nf90_put_var(fileobj%ncid, varid, variable_data, start=c, count=e) - type is (integer(kind=i8_kind)) - err = nf90_put_var(fileobj%ncid, varid, variable_data, start=c, count=e) - type is (real(kind=r4_kind)) - err = nf90_put_var(fileobj%ncid, varid, variable_data, start=c, count=e) - type is (real(kind=r8_kind)) - err = nf90_put_var(fileobj%ncid, varid, variable_data, start=c, count=e) - class default - call error("Unsupported variable type: "//trim(append_error_msg)) - end select - call check_netcdf_code(err, append_error_msg) endif end subroutine netcdf_write_data_4d diff --git a/fms2_pio_io/netcdf_io.F90 b/fms2_pio_io/netcdf_io.F90 index 698bf2406..0806e2a60 100644 --- a/fms2_pio_io/netcdf_io.F90 +++ b/fms2_pio_io/netcdf_io.F90 @@ -34,7 +34,7 @@ module netcdf_io_mod use platform_mod use mpp_mod, only: mpp_pe, mpp_npes, get_mpp_comm, mpp_root_pe, stdout -use pio, only : PIO_init, pio_createfile, PIO_initdecomp +use pio, only : PIO_init, pio_createfile, PIO_initdecomp, PIO_syncfile use pio, only : pio_createfile, pio_openfile, pio_file_is_open, PIO_closefile use pio, only : File_desc_t, var_desc_t, iosystem_desc_t, IO_desc_t use pio, only : pio_write, pio_clobber, pio_noclobber, pio_nowrite @@ -42,7 +42,7 @@ module netcdf_io_mod use pio, only : pio_iotype_netcdf, pio_iotype_pnetcdf use pio, only : PIO_set_log_level use pio, only : PIO_DOUBLE, PIO_REAL, PIO_INT, PIO_CHAR -use pio, only : PIO_64BIT_OFFSET, PIO_64BIT_DATA +use pio, only : PIO_64BIT_OFFSET, PIO_64BIT_DATA, PIO_NOERR use pio, only : pio_redef, pio_enddef, pio_global, pio_inq_dimid use pio, only : pio_seterrorhandling, pio_return_error, pio_def_dim use pio, only : pio_def_var, pio_inq_varid, pio_inquire_variable, pio_inquire_dimension @@ -246,6 +246,7 @@ module netcdf_io_mod public :: get_fill_value public :: get_variable_sense public :: get_variable_missing +public :: get_variable_id public :: get_variable_units public :: get_time_calendar public :: is_registered_to_restart @@ -357,10 +358,10 @@ module netcdf_io_mod integer :: pio_numiotasks, pio_rearranger, pio_root, pio_stride, pio_log_level integer :: pio_optbase ! Start index of I/O processors namelist / fms2_pio_nml / & - pio_netcdf_format, pio_numiotasks, pio_rearranger, pio_root, pio_stride, pio_typename,& - pio_log_level, pio_optbase + pio_netcdf_format, pio_numiotasks, pio_rearranger, pio_root, pio_stride, & + pio_typename, pio_optbase, pio_log_level -type(iosystem_desc_t) :: pio_iosystem ! The ParallelIO system set up by PIO_init +type(iosystem_desc_t), public :: pio_iosystem ! The ParallelIO system set up by PIO_init integer :: pio_iotype ! PIO_IOTYPE_NETCDF or PNETCDF integer, dimension(200) :: pio_ncids ! todo - remove @@ -415,17 +416,38 @@ subroutine fms2_pio_init () integer :: numAggregator = 0 !TODO integer :: ierr + pe = mpp_pe() + npes = mpp_npes() + localcomm = get_mpp_comm() + ! default values for pio_nml namelist vars: pio_netcdf_format = "64bit_offset" pio_typename = "netcdf" pio_numiotasks = 1 pio_rearranger = 1 pio_root = 1 - pio_stride = 2 + pio_stride = 128 pio_optbase = 1 - pio_log_level = 1 + pio_log_level = 0 READ (input_nml_file, NML=fms2_pio_nml, IOSTAT=mystat) + if (mystat /= 0) then + call error("Unable to read fms2_pio_nml from input.nml") + endif + + if (pe == mpp_root_pe()) then + write(stdout(),*) "--- fms2_pio_init namelist ---" + write(stdout(),*) " pio_netcdf_format: ", trim(pio_netcdf_format) + write(stdout(),*) " pio_typename: ", trim(pio_typename) + write(stdout(),*) " pio_numiotasks: ", pio_numiotasks + write(stdout(),*) " pio_rearranger: ", pio_rearranger + write(stdout(),*) " pio_root: ", pio_root + write(stdout(),*) " pio_stride: ", pio_stride + write(stdout(),*) " pio_optbase: ", pio_optbase + write(stdout(),*) " pio_log_level: ", pio_log_level + endif + + ! todo: here. add consistency checks for npes, numiotasks, stride, etc. select case(trim(pio_typename)) case ("netcdf") @@ -436,24 +458,6 @@ subroutine fms2_pio_init () call mpp_error(FATAL,'Unknown PIO filetype') end select - pe = mpp_pe() - npes = mpp_npes() - localcomm = get_mpp_comm() - - if (pe==mpp_root_pe()) then - write(stdout(),*) "fms2_pio_init namelist ---",& - ' pio_netcdf_format: ', trim(pio_netcdf_format), & - ' pio_typename: ', trim(pio_typename), & - ' pio_numiotasks: ', pio_numiotasks, & - ' pio_rearranger: ', pio_rearranger, & - ' pio_root: ', pio_root, & - ' pio_stride: ', pio_stride, & - ' pio_optbase: ', pio_optbase, & - ' pio_log_level: ', pio_log_level - endif - - ! todo: here. add consistency checks for npes, numiotasks, stride, etc. - !initialize PIO call PIO_init( & pe, & ! MPI rank @@ -465,7 +469,10 @@ subroutine fms2_pio_init () pio_iosystem, & ! The ParallelIO system set up by PIO_init pio_optbase ) ! Start index of I/O processors (optional) - !todo ierr = PIO_set_log_level(pio_log_level) + ierr = PIO_set_log_level(pio_log_level) + if(ierr /= PIO_NOERR) then + call error('Error setting PIO logging level') + end if if (pe==mpp_root_pe()) write(stdout(),*) "PIO initialized by FMS2 io module." pio_initialized = .true. @@ -830,7 +837,6 @@ function netcdf_file_open_pio(fileobj, path, mode, nc_format, pelist, is_restart ! Set the is_open flag to true for this file object. if (.not.allocated(fileobj%is_open)) allocate(fileobj%is_open) fileobj%is_open = .true. - print *, "pio-dbg -- opened: ", path fileobj%bc_dimensions%xlen = 0 fileobj%bc_dimensions%ylen = 0 @@ -878,7 +884,7 @@ function netcdf_file_open(fileobj, path, mode, nc_format, pelist, is_restart, do logical :: is_res logical :: dont_add_res !< flag indicated to not add ".res" to the filename - print *, "pio-dbg -- attempting to open", trim(path), mode + !print *, "pio-dbg -- attempting to open", trim(path), mode if (.not. string_compare(mode, "read", .true.)) then success = netcdf_file_open_pio(fileobj, path, mode, nc_format, pelist, is_restart, dont_add_res_to_filename) return @@ -1009,14 +1015,17 @@ subroutine netcdf_file_close_pio(fileobj) integer :: err integer :: i - print *, "pio-dbg - closing file", trim(fileobj%path), fileobj%is_open - call PIO_closefile(fileobj%file_desc) deallocate(fileobj%file_desc) if (allocated(fileobj%is_open)) fileobj%is_open = .false. fileobj%path = missing_path fileobj%ncid = missing_ncid + if (allocated(fileobj%pelist)) then + deallocate(fileobj%pelist) + endif + fileobj%io_root = missing_rank + fileobj%is_root = .false. if (allocated(fileobj%restart_vars)) then deallocate(fileobj%restart_vars) endif @@ -1033,7 +1042,6 @@ subroutine netcdf_file_close_pio(fileobj) if (allocated(fileobj%compressed_dims)) then deallocate(fileobj%compressed_dims) endif - print *, "pio-dbg -- closed: ", trim(fileobj%path), allocated(fileobj%is_open), fileobj%is_open end subroutine netcdf_file_close_pio @@ -1045,7 +1053,7 @@ subroutine netcdf_file_close(fileobj) integer :: err integer :: i - print *, "pio-dbg -- attempting to close", trim(fileobj%path) + !print *, "pio-dbg -- attempting to close", trim(fileobj%path) if (ncid_handled_by_pio(fileobj%ncid)) then call netcdf_file_close_pio(fileobj) return @@ -1299,8 +1307,7 @@ subroutine netcdf_add_variable(fileobj, variable_name, variable_type, dimensions if (string_compare(variable_type, "int", .true.)) then vtype = pio_int elseif (string_compare(variable_type, "int64", .true.)) then - if ( .not. fileobj%is_netcdf4) call error(trim(fileobj%path)// & - & ": 64 bit integers are not supported with PIO") + call error(trim(fileobj%path)//": 64 bit integers are not supported with PIO") !if ( .not. fileobj%is_netcdf4) call error(trim(fileobj%path)//& ! &": 64 bit integers are only supported with 'netcdf4' file format"//& ! &". Set netcdf_default_format='netcdf4' in the fms2_io namelist OR "//& @@ -2964,7 +2971,10 @@ subroutine flush_file(fileobj) integer :: err !< Netcdf error code - print *, "ERRORline", __LINE__, trim(__FILE__); call error("PIO version not implemented!!!") + if (ncid_handled_by_pio(fileobj%ncid)) then + call PIO_syncfile(fileobj%file_desc) + return + endif if (fileobj%is_root) then err = nf90_sync(fileobj%ncid)