From 2400ea0a0b6ce261df718e919c6da2641bf38621 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 25 Oct 2021 10:22:23 -0400 Subject: [PATCH 1/8] Refresh attempt to get rid of NetCDF calls --- src/framework/MOM_horizontal_regridding.F90 | 90 +++++------------ src/framework/MOM_io.F90 | 102 +++++++++++++++++++- 2 files changed, 125 insertions(+), 67 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 0f16a5b301..046de24548 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -16,7 +16,7 @@ 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 MOM_io, only : axis_info, get_axis_info, get_var_axes_info, MOM_read_data 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 @@ -308,6 +308,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, 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) @@ -458,12 +421,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 @@ -478,11 +437,8 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, else 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)) - + start(4) = 1; count(4) = 1 + call MOM_read_data(trim(filename), trim(varnam), tr_in, start, count, no_domain=.true.) if (add_np) then pole = 0.0 ; npole = 0.0 do i=1,id @@ -603,6 +559,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 00eeb4cf89..18351de97c 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 @@ -1540,6 +1542,31 @@ 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 array with information about named axes + character(len=*), intent(out), optional :: name !< An associated name. + character(len=*), intent(out), optional :: longname !< An associated longname. + character(len=*), intent(out), optional :: units !< Units + character(len=*), intent(out), optional :: cartesian !< cartesian direction [X,Y,Z,T] + integer, intent(out), optional :: ax_size !< The size of the axis. + real, optional, allocatable, dimension(:), intent(out) :: ax_data !< 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 @@ -2231,7 +2258,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 + real :: 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 From 95f509e338c904c00a7d67c4d0e30626a8febd2d Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 1 Nov 2021 10:20:53 -0400 Subject: [PATCH 2/8] Fix comments --- src/framework/MOM_io.F90 | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 18351de97c..ba6c0394ea 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -1545,13 +1545,14 @@ 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 array with information about named axes - character(len=*), intent(out), optional :: name !< An associated name. - character(len=*), intent(out), optional :: longname !< An associated longname. - character(len=*), intent(out), optional :: units !< Units - character(len=*), intent(out), optional :: cartesian !< cartesian direction [X,Y,Z,T] - integer, intent(out), optional :: ax_size !< The size of the axis. - real, optional, allocatable, dimension(:), intent(out) :: ax_data !< axis label 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) From 2909e52a89d029caec1402ef56c62b08bb13a164 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Mon, 8 Nov 2021 13:21:23 -0500 Subject: [PATCH 3/8] 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. --- src/framework/MOM_horizontal_regridding.F90 | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 046de24548..f026c8aa88 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -17,8 +17,7 @@ module MOM_horizontal_regridding use MOM_interp_infra, only : axistype, get_external_field_info, get_axis_data use MOM_time_manager, only : time_type use MOM_io, only : axis_info, get_axis_info, get_var_axes_info, MOM_read_data -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 : read_attribute implicit none ; private @@ -304,6 +303,7 @@ 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 @@ -385,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 = 0. + + 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)) From 051e8a4ec8de783941ea85694c0f4d428996d767 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 18 Nov 2021 13:35:58 -0500 Subject: [PATCH 4/8] 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()`. --- src/framework/MOM_horizontal_regridding.F90 | 9 ++-- src/framework/MOM_io.F90 | 56 +++++++++++++++++++++ 2 files changed, 60 insertions(+), 5 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index f026c8aa88..ecc86ea397 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -17,7 +17,7 @@ module MOM_horizontal_regridding use MOM_interp_infra, only : axistype, get_external_field_info, get_axis_data use MOM_time_manager, only : time_type use MOM_io, only : axis_info, get_axis_info, get_var_axes_info, MOM_read_data -use MOM_io, only : read_attribute +use MOM_io, only : read_attribute, read_variable implicit none ; private @@ -448,12 +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 - start(4) = 1; count(4) = 1 - call MOM_read_data(trim(filename), trim(varnam), tr_in, start, count, no_domain=.true.) if (add_np) then pole = 0.0 ; npole = 0.0 do i=1,id diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index ba6c0394ea..ee074ce01f 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -100,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 +888,61 @@ 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 + + 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)) + + ! NOTE: We could check additional information here (type, size, ...) + + rc = nf90_get_var(ncid, varid, var) + if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& + " Difficulties reading "//trim(varname)//" from "//trim(filename)) + + if (size(start) /= ndims .or. size(nread) /= ndims) 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 + + 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 From 93b74ad25a6991bc6f090dc8e5ee9263cb52c115 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Wed, 29 Dec 2021 14:58:05 -0500 Subject: [PATCH 5/8] set default scale factor to 1 --- src/framework/MOM_horizontal_regridding.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index ecc86ea397..de511688a9 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -394,7 +394,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, call read_attribute(trim(filename), "scale_factor", scale_factor, & varname=trim(varnam), found=found_attr) - if (.not. found_attr) scale_factor = 0. + if (.not. found_attr) scale_factor = 1. call read_attribute(trim(filename), "add_offset", add_offset, & varname=trim(varnam), found=found_attr) From c664aa51a8ba0cfd753c3fd440ab79d97f656122 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 3 Jan 2022 21:50:28 -0500 Subject: [PATCH 6/8] add missing start/count arguments --- src/framework/MOM_io.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 1fdb535d4d..f4c30e0237 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -928,8 +928,7 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in) " Difficulties reading "//trim(varname)//" from "//trim(filename)) ! NOTE: We could check additional information here (type, size, ...) - - rc = nf90_get_var(ncid, varid, var) + 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)) From 8fae1eef56086317ea0f0ee8f41099b93a1890ec Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 4 Jan 2022 10:05:33 -0500 Subject: [PATCH 7/8] Update MOM_io.F90 --- src/framework/MOM_io.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index f4c30e0237..85e66d5312 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -2323,7 +2323,7 @@ subroutine get_var_axes_info(filename, fieldname, axes_info) type(axis_info), dimension(4), intent(inout) :: axes_info !< A returned array of field axis information !! local variables - real :: rcode + integer :: rcode logical :: success integer :: ncid, varid, ndims integer :: id, jd, kd From 9865334d20b9caaf98327be5bd9ac758f4ace00b Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 4 Jan 2022 10:29:49 -0500 Subject: [PATCH 8/8] 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. --- src/framework/MOM_io.F90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 85e66d5312..2b8fb210d5 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -911,6 +911,7 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in) 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 @@ -927,16 +928,20 @@ subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in) if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //& " Difficulties reading "//trim(varname)//" from "//trim(filename)) - ! 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)) + 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(start) /= ndims .or. size(nread) /= ndims) then + 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