-
Notifications
You must be signed in to change notification settings - Fork 82
Rewrite horizontal regridding to use netCDF wrapper functions #48
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Changes from 9 commits
2400ea0
95f509e
2909e52
221df78
051e8a4
577c0c1
93b74ad
282395b
c664aa5
8fae1ee
9865334
ecfa52c
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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,60 @@ 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, start, nread) | ||
| 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 | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think that there a several problems with this test.
Could this work properly if the tests were changed to
Member
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Optional arguments are pass-through, so the call to |
||
| 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 | ||
|
|
@@ -1542,6 +1599,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 +2316,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 | ||
|
marshallward marked this conversation as resolved.
Outdated
|
||
| 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 | ||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Consider changing this to something like:
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If we're going to support reading 2d slices in a 3d+ field, I would prefer stronger dimensional testing.
In truth, I don't think we should be encouraging use of this function at all. It does not support parallel I/O and problems at high resolution are inevitable.