Skip to content

Commit

Permalink
Rewrite horizontal regridding to use netCDF wrapper functions (#48)
Browse files Browse the repository at this point in the history 
* 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 <Matthew.Harrison@noaa.gov>
  • Loading branch information
@marshallward @MJHarrison-GFDL
marshallward and MJHarrison-GFDL authored Jan 8, 2022
1 parent df46be4 commit f7a2254
Show file tree
Hide file tree
Showing 2 changed files with 204 additions and 71 deletions.
112 changes: 42 additions & 70 deletions src/framework/MOM_horizontal_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -334,71 +335,33 @@ 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
! The convention adopted here requires that the data be written in (i,j,k) ordering.

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)

Expand All @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
163 changes: 162 additions & 1 deletion src/framework/MOM_io.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit f7a2254

Please sign in to comment.