From f7a2254f3652727217d578f3aac4966a0d4da864 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Sat, 8 Jan 2022 06:12:14 -0500 Subject: [PATCH] Rewrite horizontal regridding to use netCDF wrapper functions (#48) * Refresh attempt to get rid of NetCDF calls * Fix comments * Set netCDF attrs in MOM_horizontal_regridding This patch sets the following netCDF attributes in the function `horiz_interp_and_extrap_tracer_record` via `read_attribute`. * `missing_value` (as `_FillValue`) * `scale_factor` * `add_offset` This resolves some issues related to uninitialized values. * read_variable_2d in horizontal remapping This patch extends the `read_variable` interface to include 2d array support, in order to facilitate domainless I/O via netCDF calls. This is far from the best implementation (e.g. read_variable_2d introduces another `broadcast` alongside the original one in the horizontal regridding) but it addresses the immediate issues with `MOM_read_data()`. * set default scale factor to 1 * add missing start/count arguments * Update MOM_io.F90 * Manage optional args in read_variable_2d This patch modifies read_variable_2d so that the size() tests of the optional arguments are applied before the call to nf90_get_var. The tests are also wrapped inside present() flags to avoid checking unassigned variables. Thanks to Robert Hallberg for the suggestions. Co-authored-by: Matthew Harrison --- src/framework/MOM_horizontal_regridding.F90 | 112 +++++--------- src/framework/MOM_io.F90 | 163 +++++++++++++++++++- 2 files changed, 204 insertions(+), 71 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 0f16a5b301..de511688a9 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -16,9 +16,8 @@ module MOM_horizontal_regridding use MOM_interpolate, only : build_horiz_interp_weights, run_horiz_interp, horiz_interp_type use MOM_interp_infra, only : axistype, get_external_field_info, get_axis_data use MOM_time_manager, only : time_type - -use netcdf, only : NF90_OPEN, NF90_NOWRITE, NF90_GET_ATT, NF90_GET_VAR -use netcdf, only : NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, NF90_INQUIRE_DIMENSION +use MOM_io, only : axis_info, get_axis_info, get_var_axes_info, MOM_read_data +use MOM_io, only : read_attribute, read_variable implicit none ; private @@ -304,10 +303,12 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real :: max_lat, min_lat, pole, max_depth, npole real :: roundoff ! The magnitude of roundoff, usually ~2e-16. real :: add_offset, scale_factor + logical :: found_attr logical :: add_np logical :: is_ongrid character(len=8) :: laynum type(horiz_interp_type) :: Interp + type(axis_info), dimension(4) :: axes_info ! Axis information used for regridding integer :: is, ie, js, je ! compute domain indices integer :: isc, iec, jsc, jec ! global compute domain indices integer :: isg, ieg, jsg, jeg ! global extent @@ -334,6 +335,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, is_ongrid = .false. if (present(ongrid)) is_ongrid = ongrid + if (allocated(tr_z)) deallocate(tr_z) + if (allocated(mask_z)) deallocate(mask_z) + PI_180 = atan(1.0)/45. ! Open NetCDF file and if present, extract data and spatial coordinate information @@ -341,64 +345,23 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, call cpu_clock_begin(id_clock_read) - rcode = NF90_OPEN(filename, NF90_NOWRITE, ncid) - if (rcode /= 0) call MOM_error(FATAL,"error opening file "//trim(filename)//& - " in hinterp_extrap") - rcode = NF90_INQ_VARID(ncid, varnam, varid) - if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(varnam)//& - " in file "//trim(filename)//" in hinterp_extrap") - - rcode = NF90_INQUIRE_VARIABLE(ncid, varid, ndims=ndims, dimids=dims) - if (rcode /= 0) call MOM_error(FATAL, "Error inquiring about the dimensions of "//trim(varnam)//& - " in file "//trim(filename)//" in hinterp_extrap") - if (ndims < 3) call MOM_error(FATAL,"Variable "//trim(varnam)//" in file "//trim(filename)// & - " has too few dimensions to be read as a 3-d array.") - - rcode = NF90_INQUIRE_DIMENSION(ncid, dims(1), dim_name(1), len=id) - if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 data for "// & - trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap") - rcode = NF90_INQ_VARID(ncid, dim_name(1), dim_id(1)) - if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(1))//& - " in file "//trim(filename)//" in hinterp_extrap") - rcode = NF90_INQUIRE_DIMENSION(ncid, dims(2), dim_name(2), len=jd) - if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 2 data for "// & - trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap") - rcode = NF90_INQ_VARID(ncid, dim_name(2), dim_id(2)) - if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(2))//& - " in file "//trim(filename)//" in hinterp_extrap") - rcode = NF90_INQUIRE_DIMENSION(ncid, dims(3), dim_name(3), len=kd) - if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 3 data for "// & - trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap") - rcode = NF90_INQ_VARID(ncid, dim_name(3), dim_id(3)) - if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(3))//& - " in file "//trim(filename)//" in hinterp_extrap") - - missing_value=0.0 - rcode = NF90_GET_ATT(ncid, varid, "_FillValue", missing_value) - if (rcode /= 0) call MOM_error(FATAL,"error finding missing value for "//trim(varnam)//& - " in file "// trim(filename)//" in hinterp_extrap") - - rcode = NF90_GET_ATT(ncid, varid, "add_offset", add_offset) - if (rcode /= 0) add_offset = 0.0 - - rcode = NF90_GET_ATT(ncid, varid, "scale_factor", scale_factor) - if (rcode /= 0) scale_factor = 1.0 + call get_var_axes_info(trim(filename), trim(varnam), axes_info) + + if (allocated(z_in)) deallocate(z_in) + if (allocated(z_edges_in)) deallocate(z_edges_in) + if (allocated(tr_z)) deallocate(tr_z) + if (allocated(mask_z)) deallocate(mask_z) + + call get_axis_info(axes_info(1),ax_size=id) + call get_axis_info(axes_info(2),ax_size=jd) + call get_axis_info(axes_info(3),ax_size=kd) allocate(lon_in(id), lat_in(jd), z_in(kd), z_edges_in(kd+1)) allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd)) - start = 1 ; count = 1 ; count(1) = id - rcode = NF90_GET_VAR(ncid, dim_id(1), lon_in, start, count) - if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 values for var_name "// & - trim(varnam)//",dim_name "//trim(dim_name(1))//" in file "// trim(filename)//" in hinterp_extrap") - start = 1 ; count = 1 ; count(1) = jd - rcode = NF90_GET_VAR(ncid, dim_id(2), lat_in, start, count) - if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 2 values for var_name "// & - trim(varnam)//",dim_name "//trim(dim_name(2))//" in file "// trim(filename)//" in hinterp_extrap") - start = 1 ; count = 1 ; count(1) = kd - rcode = NF90_GET_VAR(ncid, dim_id(3), z_in, start, count) - if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 3 values for var_name "// & - trim(varnam//",dim_name "//trim(dim_name(3)))//" in file "// trim(filename)//" in hinterp_extrap") + call get_axis_info(axes_info(1),ax_data=lon_in) + call get_axis_info(axes_info(2),ax_data=lat_in) + call get_axis_info(axes_info(3),ax_data=z_in) call cpu_clock_end(id_clock_read) @@ -422,6 +385,21 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, endif ! construct level cell boundaries as the mid-point between adjacent centers + ! Set the I/O attributes + call read_attribute(trim(filename), "_FillValue", missing_value, & + varname=trim(varnam), found=found_attr) + if (.not. found_attr) call MOM_error(FATAL, & + "error finding missing value for " // trim(varnam) // & + " in file " // trim(filename) // " in hinterp_extrap") + + call read_attribute(trim(filename), "scale_factor", scale_factor, & + varname=trim(varnam), found=found_attr) + if (.not. found_attr) scale_factor = 1. + + call read_attribute(trim(filename), "add_offset", add_offset, & + varname=trim(varnam), found=found_attr) + if (.not. found_attr) add_offset = 0. + z_edges_in(1) = 0.0 do K=2,kd z_edges_in(K) = 0.5*(z_in(k-1)+z_in(k)) @@ -458,12 +436,8 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, mask_in = 0.0 if (is_ongrid) then start(1) = is+G%HI%idg_offset ; start(2) = js+G%HI%jdg_offset ; start(3) = k - count(1) = ie-is+1 ; count(2) = je-js+1; count(3) = 1 - rcode = NF90_GET_VAR(ncid,varid, tr_in, start, count) - if (rcode /= 0) call MOM_error(FATAL,"horiz_interp_and_extrap_tracer_record: "//& - "error reading level "//trim(laynum)//" of variable "//& - trim(varnam)//" in file "// trim(filename)) - + count(1) = ie-is+1 ; count(2) = je-js+1; count(3) = 1; start(4) = 1; count(4) = 1 + call MOM_read_data(trim(filename), trim(varnam), tr_in, G%Domain, timelevel=1) do j=js,je do i=is,ie if (abs(tr_in(i,j)-missing_value) > abs(roundoff*missing_value)) then @@ -474,15 +448,11 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, endif enddo enddo - else + start(:) = 1 ; start(3) = k + count(:) = 1 ; count(1) = id ; count(2) = jd + call read_variable(trim(filename), trim(varnam), tr_in, start=start, nread=count) if (is_root_pe()) then - start = 1 ; start(3) = k ; count(:) = 1 ; count(1) = id ; count(2) = jd - rcode = NF90_GET_VAR(ncid,varid, tr_in, start, count) - if (rcode /= 0) call MOM_error(FATAL,"horiz_interp_and_extrap_tracer_record: "//& - "error reading level "//trim(laynum)//" of variable "//& - trim(varnam)//" in file "// trim(filename)) - if (add_np) then pole = 0.0 ; npole = 0.0 do i=1,id @@ -603,6 +573,8 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, enddo ! kd + deallocate(lon_in, lat_in) + end subroutine horiz_interp_and_extrap_tracer_record !> Extrapolate and interpolate using a FMS time interpolation handle diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 563f9f9f8a..2b8fb210d5 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -51,6 +51,8 @@ module MOM_io public :: slasher, write_field, write_version_number public :: io_infra_init, io_infra_end public :: stdout_if_root +public :: get_var_axes_info +public :: get_axis_info ! This is used to set up information descibing non-domain-decomposed axes. public :: axis_info, set_axis_info, delete_axis_info ! This is used to set up global file attributes @@ -98,6 +100,7 @@ module MOM_io interface read_variable module procedure read_variable_0d, read_variable_0d_int module procedure read_variable_1d, read_variable_1d_int + module procedure read_variable_2d end interface read_variable !> Read a global or variable attribute from a named netCDF file using netCDF calls @@ -887,6 +890,65 @@ subroutine read_variable_1d_int(filename, varname, var, ncid_in) call broadcast(var, size(var), blocking=.true.) end subroutine read_variable_1d_int +!> Read a 2d array from a netCDF input file and save to a variable. +!! +!! Start and nread ranks may exceed var, but must match the rank of the +!! variable in the netCDF file. This allows for reading slices of larger +!! arrays. +!! +!! I/O occurs only on the root PE, and data is broadcast to other ranks. +!! Due to potentially large memory communication and storage, this subroutine +!! should only be used when domain-decomposition is unavaialable. +subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in) + character(len=*), intent(in) :: filename !< Name of file to be read + character(len=*), intent(in) :: varname !< Name of variable to be read + real, intent(out) :: var(:,:) !< Output array of variable + integer, optional, intent(in) :: start(:) !< Starting index on each axis. + integer, optional, intent(in) :: nread(:) !< Number of values to be read along each axis + integer, optional, intent(in) :: ncid_in !< netCDF ID of an opened file. + !! If absent, the file is opened and closed within this routine. + + integer :: ncid, varid, ndims, rc + character(len=*), parameter :: hdr = "read_variable_2d" + character(len=128) :: msg + logical :: size_mismatch + + if (is_root_pe()) then + if (present(ncid_in)) then + ncid = ncid_in + else + call open_file_to_Read(filename, ncid) + endif + + call get_varid(varname, ncid, filename, varid, match_case=.false.) + if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//& + " in "//trim(filename)) + ! Verify that start(:) and nread(:) ranks match variable's dimension count + rc = nf90_inquire_variable(ncid, varid, ndims=ndims) + if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& + " Difficulties reading "//trim(varname)//" from "//trim(filename)) + + size_mismatch = .false. + if (present(start)) size_mismatch = size_mismatch .or. size(start) /= ndims + if (present(nread)) size_mismatch = size_mismatch .or. size(nread) /= ndims + + if (size_mismatch) then + write (msg, '("'// hdr //': size(start) ", i0, " and/or size(nread) ", & + i0, " do not match ndims ", i0)') size(start), size(nread), ndims + call MOM_error(FATAL, trim(msg)) + endif + ! NOTE: We could check additional information here (type, size, ...) + + rc = nf90_get_var(ncid, varid, var, start, nread) + if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& + " Difficulties reading "//trim(varname)//" from "//trim(filename)) + + if (.not.present(ncid_in)) call close_file_to_read(ncid, filename) + endif + + call broadcast(var, size(var), blocking=.true.) +end subroutine read_variable_2d + !> Read a character-string global or variable attribute subroutine read_attribute_str(filename, attname, att_val, varname, found, all_read, ncid_in) character(len=*), intent(in) :: filename !< Name of the file to read @@ -1542,6 +1604,32 @@ subroutine delete_axis_info(axes) enddo end subroutine delete_axis_info + +!> Retrieve the information from an axis_info type. +subroutine get_axis_info(axis,name,longname,units,cartesian,ax_size,ax_data) + type(axis_info), intent(in) :: axis !< An axis type + character(len=*), intent(out), optional :: name !< The axis name. + character(len=*), intent(out), optional :: longname !< The axis longname. + character(len=*), intent(out), optional :: units !< The axis units. + character(len=*), intent(out), optional :: cartesian !< The cartesian attribute + !! of the axis [X,Y,Z,T]. + integer, intent(out), optional :: ax_size !< The size of the axis. + real, optional, allocatable, dimension(:), intent(out) :: ax_data !< The axis label data. + + if (present(ax_data)) then + if (allocated(ax_data)) deallocate(ax_data) + allocate(ax_data(axis%ax_size)) + ax_data(:)=axis%ax_data + endif + + if (present(name)) name=axis%name + if (present(longname)) longname=axis%longname + if (present(units)) units=axis%units + if (present(cartesian)) cartesian=axis%cartesian + if (present(ax_size)) ax_size=axis%ax_size + +end subroutine get_axis_info + !> Store information that can be used to create an attribute in a subsequent call to create_file. subroutine set_attribute_info(attribute, name, str_value) type(attribute_info), intent(inout) :: attribute !< A type with information about a named attribute @@ -2233,7 +2321,80 @@ subroutine MOM_io_init(param_file) call log_version(param_file, mdl, version) end subroutine MOM_io_init - +!> Returns the dimension variable information for a netCDF variable +subroutine get_var_axes_info(filename, fieldname, axes_info) + character(len=*), intent(in) :: filename !< A filename from which to read + character(len=*), intent(in) :: fieldname !< The name of the field to read + type(axis_info), dimension(4), intent(inout) :: axes_info !< A returned array of field axis information + + !! local variables + integer :: rcode + logical :: success + integer :: ncid, varid, ndims + integer :: id, jd, kd + integer, dimension(4) :: dims, dim_id + real :: missing_value + character(len=128) :: dim_name(4) + integer, dimension(1) :: start, count + !! cartesian axis data + real, allocatable, dimension(:) :: x + real, allocatable, dimension(:) :: y + real, allocatable, dimension(:) :: z + + + call open_file_to_read(filename, ncid, success=success) + + rcode = NF90_INQ_VARID(ncid, trim(fieldname), varid) + if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(fieldname)//& + " in file "//trim(filename)//" in hinterp_extrap") + + rcode = NF90_INQUIRE_VARIABLE(ncid, varid, ndims=ndims, dimids=dims) + if (rcode /= 0) call MOM_error(FATAL, "Error inquiring about the dimensions of "//trim(fieldname)//& + " in file "//trim(filename)//" in hinterp_extrap") + if (ndims < 3) call MOM_error(FATAL,"Variable "//trim(fieldname)//" in file "//trim(filename)// & + " has too few dimensions to be read as a 3-d array.") + rcode = NF90_INQUIRE_DIMENSION(ncid, dims(1), dim_name(1), len=id) + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 data for "// & + trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap") + rcode = NF90_INQ_VARID(ncid, dim_name(1), dim_id(1)) + if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(1))//& + " in file "//trim(filename)//" in hinterp_extrap") + rcode = NF90_INQUIRE_DIMENSION(ncid, dims(2), dim_name(2), len=jd) + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 2 data for "// & + trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap") + rcode = NF90_INQ_VARID(ncid, dim_name(2), dim_id(2)) + if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(2))//& + " in file "//trim(filename)//" in hinterp_extrap") + rcode = NF90_INQUIRE_DIMENSION(ncid, dims(3), dim_name(3), len=kd) + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 3 data for "// & + trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap") + rcode = NF90_INQ_VARID(ncid, dim_name(3), dim_id(3)) + if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(3))//& + " in file "//trim(filename)//" in hinterp_extrap") + allocate(x(id), y(jd), z(kd)) + + start = 1 ; count = 1 ; count(1) = id + rcode = NF90_GET_VAR(ncid, dim_id(1), x, start, count) + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 values for var_name "// & + trim(fieldname)//",dim_name "//trim(dim_name(1))//" in file "// trim(filename)//" in hinterp_extrap") + start = 1 ; count = 1 ; count(1) = jd + rcode = NF90_GET_VAR(ncid, dim_id(2), y, start, count) + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 2 values for var_name "// & + trim(fieldname)//",dim_name "//trim(dim_name(2))//" in file "// trim(filename)//" in hinterp_extrap") + start = 1 ; count = 1 ; count(1) = kd + rcode = NF90_GET_VAR(ncid, dim_id(3), z, start, count) + if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 3 values for var_name "// & + trim(fieldname//",dim_name "//trim(dim_name(3)))//" in file "// trim(filename)//" in hinterp_extrap") + + call set_axis_info(axes_info(1), name=trim(dim_name(1)), ax_size=id, ax_data=x,cartesian='X') + call set_axis_info(axes_info(2), name=trim(dim_name(2)), ax_size=jd, ax_data=y,cartesian='Y') + call set_axis_info(axes_info(3), name=trim(dim_name(3)), ax_size=kd, ax_data=z,cartesian='Z') + + call close_file_to_read(ncid, filename) + + deallocate(x,y,z) + +end subroutine get_var_axes_info !> \namespace mom_io !! !! This file contains a number of subroutines that manipulate