From ee0152f2699a1f38b9ca95005c3b9287535e6c57 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 18 Mar 2021 17:16:41 -0400 Subject: [PATCH 01/10] +Separate MOM_interp_infra axistype from MOM_io Use axistype for MOM_interp_infra directly from mpp_io_mod and add a copy of get_axis_data to both copies of infra/FMS[12]/MOM_interp_infra.F90, and then use these in framework/MOM_horizontal_regridding.F90, to permit the MOM6 I/O calls to use the FMS2 interfaces without simultaneously requiring changes to the horizontal interpolation code. All answers are bitwise identical, but there are changes to the interfaces offered by a public module. --- config_src/infra/FMS1/MOM_interp_infra.F90 | 13 +++++++++++-- config_src/infra/FMS2/MOM_interp_infra.F90 | 13 +++++++++++-- src/framework/MOM_horizontal_regridding.F90 | 4 ++-- src/framework/MOM_interpolate.F90 | 1 - 4 files changed, 24 insertions(+), 7 deletions(-) diff --git a/config_src/infra/FMS1/MOM_interp_infra.F90 b/config_src/infra/FMS1/MOM_interp_infra.F90 index ca5b2b8516..170573f7ec 100644 --- a/config_src/infra/FMS1/MOM_interp_infra.F90 +++ b/config_src/infra/FMS1/MOM_interp_infra.F90 @@ -4,9 +4,9 @@ module MOM_interp_infra ! This file is part of MOM6. See LICENSE.md for the license. use MOM_domain_infra, only : MOM_domain_type, domain2d -use MOM_io_infra, only : axistype use MOM_time_manager, only : time_type use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use mpp_io_mod, only : axistype, mpp_get_axis_data use time_interp_external_mod, only : time_interp_external use time_interp_external_mod, only : init_external_field, time_interp_external_init use time_interp_external_mod, only : get_external_field_size @@ -16,7 +16,7 @@ module MOM_interp_infra public :: horiz_interp_type, horiz_interp_init public :: time_interp_extern, init_extern_field, time_interp_external_init -public :: get_external_field_info +public :: get_external_field_info, axistype, get_axis_data public :: run_horiz_interp, build_horiz_interp_weights !> Read a field based on model time, and rotate to the model domain. @@ -114,6 +114,15 @@ subroutine build_horiz_interp_weights_2d_to_2d(Interp, lon_in, lat_in, lon_out, end subroutine build_horiz_interp_weights_2d_to_2d +!> Extracts and returns the axis data stored in an axistype. +subroutine get_axis_data( axis, dat ) + type(axistype), intent(in) :: axis !< An axis type + real, dimension(:), intent(out) :: dat !< The data in the axis variable + + call mpp_get_axis_data( axis, dat ) +end subroutine get_axis_data + + !> get size of an external field from field index function get_extern_field_size(index) diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index ca5b2b8516..170573f7ec 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -4,9 +4,9 @@ module MOM_interp_infra ! This file is part of MOM6. See LICENSE.md for the license. use MOM_domain_infra, only : MOM_domain_type, domain2d -use MOM_io_infra, only : axistype use MOM_time_manager, only : time_type use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use mpp_io_mod, only : axistype, mpp_get_axis_data use time_interp_external_mod, only : time_interp_external use time_interp_external_mod, only : init_external_field, time_interp_external_init use time_interp_external_mod, only : get_external_field_size @@ -16,7 +16,7 @@ module MOM_interp_infra public :: horiz_interp_type, horiz_interp_init public :: time_interp_extern, init_extern_field, time_interp_external_init -public :: get_external_field_info +public :: get_external_field_info, axistype, get_axis_data public :: run_horiz_interp, build_horiz_interp_weights !> Read a field based on model time, and rotate to the model domain. @@ -114,6 +114,15 @@ subroutine build_horiz_interp_weights_2d_to_2d(Interp, lon_in, lat_in, lon_out, end subroutine build_horiz_interp_weights_2d_to_2d +!> Extracts and returns the axis data stored in an axistype. +subroutine get_axis_data( axis, dat ) + type(axistype), intent(in) :: axis !< An axis type + real, dimension(:), intent(out) :: dat !< The data in the axis variable + + call mpp_get_axis_data( axis, dat ) +end subroutine get_axis_data + + !> get size of an external field from field index function get_extern_field_size(index) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 8d02dfdf3f..73fb1f0a41 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -11,9 +11,9 @@ module MOM_horizontal_regridding use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_interpolate, only : time_interp_external, get_external_field_info, horiz_interp_init +use MOM_interpolate, only : time_interp_external, horiz_interp_init use MOM_interpolate, only : build_horiz_interp_weights, run_horiz_interp, horiz_interp_type -use MOM_io_infra, only : axistype, get_axis_data +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 diff --git a/src/framework/MOM_interpolate.F90 b/src/framework/MOM_interpolate.F90 index f282d03ff6..4a931d0bf3 100644 --- a/src/framework/MOM_interpolate.F90 +++ b/src/framework/MOM_interpolate.F90 @@ -9,7 +9,6 @@ module MOM_interpolate use MOM_interp_infra, only : time_interp_external_init, get_external_field_info use MOM_interp_infra, only : horiz_interp_type, horiz_interp_init use MOM_interp_infra, only : run_horiz_interp, build_horiz_interp_weights -use MOM_io_infra, only : axistype use MOM_time_manager, only : time_type implicit none ; private From f89ddeaa80113e0806cb8b3441e39bfde33e7a4a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 18 Mar 2021 18:21:19 -0400 Subject: [PATCH 02/10] +Add edge_axis argument to write_metadata_axis Added a new optional logical argument to write_metadata_axis to indicate when an axis is staggered at the edges of the tracer grid, and changed calls to get_file_info to stop requesting the number of global attributes. Also eliminated some unused optional arguments to the FMS1 version of write_metadata field. All answers are bitwise identical, but there are minor changes to some I/O related interfaces. --- config_src/infra/FMS1/MOM_io_infra.F90 | 33 +++++++++++--------------- config_src/infra/FMS2/MOM_io_infra.F90 | 14 ++++++----- src/framework/MOM_io.F90 | 18 +++++++------- 3 files changed, 31 insertions(+), 34 deletions(-) diff --git a/config_src/infra/FMS1/MOM_io_infra.F90 b/config_src/infra/FMS1/MOM_io_infra.F90 index 3ea201235a..0d4cc0deb5 100644 --- a/config_src/infra/FMS1/MOM_io_infra.F90 +++ b/config_src/infra/FMS1/MOM_io_infra.F90 @@ -312,13 +312,12 @@ subroutine get_filename_suffix(suffix) end subroutine get_filename_suffix -!> Get information about the number of dimensions, variables, global attributes and time levels +!> Get information about the number of dimensions, variables and time levels !! in the file associated with an open file unit -subroutine get_file_info(IO_handle, ndim, nvar, natt, ntime) +subroutine get_file_info(IO_handle, ndim, nvar, ntime) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O integer, optional, intent(out) :: ndim !< The number of dimensions in the file integer, optional, intent(out) :: nvar !< The number of variables in the file - integer, optional, intent(out) :: natt !< The number of global attributes in the file integer, optional, intent(out) :: ntime !< The number of time levels in the file ! Local variables @@ -328,7 +327,6 @@ subroutine get_file_info(IO_handle, ndim, nvar, natt, ntime) if (present(ndim)) ndim = ndims if (present(nvar)) nvar = nvars - if (present(natt)) natt = natts if (present(ntime)) ntime = ntimes end subroutine get_file_info @@ -683,7 +681,7 @@ subroutine write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:,:), intent(inout) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value @@ -697,7 +695,7 @@ subroutine write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:), intent(inout) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value @@ -711,7 +709,7 @@ subroutine write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:), intent(inout) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value @@ -724,7 +722,7 @@ subroutine write_field_1d(IO_handle, field_md, field, tstamp) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, dimension(:), intent(in) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field call mpp_write(IO_handle%unit, field_md, field, tstamp=tstamp) end subroutine write_field_1d @@ -734,7 +732,7 @@ subroutine write_field_0d(IO_handle, field_md, field, tstamp) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, intent(in) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field call mpp_write(IO_handle%unit, field_md, field, tstamp=tstamp) end subroutine write_field_0d @@ -750,7 +748,8 @@ end subroutine MOM_write_axis !> Store information about an axis in a previously defined axistype and write this !! information to the file indicated by unit. -subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian, sense, domain, data, calendar) +subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian, sense, domain, & + data, edge_axis, calendar) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(axistype), intent(inout) :: axis !< The axistype where this information is stored. character(len=*), intent(in) :: name !< The name in the file of this axis @@ -763,6 +762,7 @@ subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian !! -1 if they increase downward. type(domain1D), optional, intent(in) :: domain !< The domain decomposion for this axis real, dimension(:), optional, intent(in) :: data !< The coordinate values of the points on this axis + logical, optional, intent(in) :: edge_axis !< If true, this axis marks an edge of the tracer cells character(len=*), optional, intent(in) :: calendar !< The name of the calendar used with a time axis call mpp_write_meta(IO_handle%unit, axis, name, units, longname, cartesian=cartesian, sense=sense, & @@ -772,19 +772,13 @@ end subroutine write_metadata_axis !> Store information about an output variable in a previously defined fieldtype and write this !! information to the file indicated by unit. subroutine write_metadata_field(IO_handle, field, axes, name, units, longname, & - min, max, fill, scale, add, pack, standard_name, checksum) + pack, standard_name, checksum) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(inout) :: field !< The fieldtype where this information is stored type(axistype), dimension(:), intent(in) :: axes !< Handles for the axis used for this variable character(len=*), intent(in) :: name !< The name in the file of this variable character(len=*), intent(in) :: units !< The units of this variable character(len=*), intent(in) :: longname !< The long description of this variable - real, optional, intent(in) :: min !< The minimum valid value for this variable - real, optional, intent(in) :: max !< The maximum valid value for this variable - real, optional, intent(in) :: fill !< Missing data fill value - real, optional, intent(in) :: scale !< An multiplicative factor by which to scale - !! the variable before output - real, optional, intent(in) :: add !< An offset to add to the variable before output integer, optional, intent(in) :: pack !< A precision reduction factor with which the !! variable. The default, 1, has no reduction, !! but 2 is not uncommon. @@ -793,8 +787,9 @@ subroutine write_metadata_field(IO_handle, field, axes, name, units, longname, & optional, intent(in) :: checksum !< Checksum values that can be used to verify reads. - call mpp_write_meta(IO_handle%unit, field, axes, name, units, longname, min=min, max=max, & - fill=fill, scale=scale, add=add, pack=pack, standard_name=standard_name, checksum=checksum) + call mpp_write_meta(IO_handle%unit, field, axes, name, units, longname, & + pack=pack, standard_name=standard_name, checksum=checksum) + ! unused opt. args: min=min, max=max, fill=fill, scale=scale, add=add, & end subroutine write_metadata_field diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index 22548218d1..8115d4acfa 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -879,7 +879,7 @@ subroutine write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:,:), intent(inout) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value @@ -893,7 +893,7 @@ subroutine write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:), intent(inout) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value @@ -907,7 +907,7 @@ subroutine write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:), intent(inout) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value @@ -920,7 +920,7 @@ subroutine write_field_1d(IO_handle, field_md, field, tstamp) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, dimension(:), intent(in) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field call mpp_write(IO_handle%unit, field_md, field, tstamp=tstamp) end subroutine write_field_1d @@ -930,7 +930,7 @@ subroutine write_field_0d(IO_handle, field_md, field, tstamp) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, intent(in) :: field !< Field to write - real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: tstamp !< Model time of this field call mpp_write(IO_handle%unit, field_md, field, tstamp=tstamp) end subroutine write_field_0d @@ -946,7 +946,8 @@ end subroutine MOM_write_axis !> Store information about an axis in a previously defined axistype and write this !! information to the file indicated by unit. -subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian, sense, domain, data, calendar) +subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian, sense, domain, data, & + edge_axis, calendar) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(axistype), intent(inout) :: axis !< The axistype where this information is stored. character(len=*), intent(in) :: name !< The name in the file of this axis @@ -959,6 +960,7 @@ subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian !! -1 if they increase downward. type(domain1D), optional, intent(in) :: domain !< The domain decomposion for this axis real, dimension(:), optional, intent(in) :: data !< The coordinate values of the points on this axis + logical, optional, intent(in) :: edge_axis !< If true, this axis marks an edge of the tracer cells character(len=*), optional, intent(in) :: calendar !< The name of the calendar used with a time axis call mpp_write_meta(IO_handle%unit, axis, name, units, longname, cartesian=cartesian, sense=sense, & diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 1159ac87e1..247a0a9678 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -254,11 +254,11 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim if (use_latq) & call write_metadata(IO_handle, axis_latq, name="latq", units=y_axis_units, longname="Latitude", & - cartesian='Y', domain=y_domain, data=gridLatB(JsgB:JegB)) + cartesian='Y', domain=y_domain, data=gridLatB(JsgB:JegB), edge_axis=.true.) if (use_lonq) & call write_metadata(IO_handle, axis_lonq, name="lonq", units=x_axis_units, longname="Longitude", & - cartesian='X', domain=x_domain, data=gridLonB(IsgB:IegB)) + cartesian='X', domain=x_domain, data=gridLonB(IsgB:IegB), edge_axis=.true.) if (use_layer) & call write_metadata(IO_handle, axis_layer, name="Layer", units=trim(GV%zAxisUnits), & @@ -383,7 +383,7 @@ subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, tim type(MOM_domain_type), pointer :: Domain => NULL() character(len=200) :: check_name, mesg - integer :: length, ndim, nvar, natt, ntime, thread + integer :: length, nvar, thread logical :: exists, one_file, domain_set thread = SINGLE_FILE @@ -418,7 +418,7 @@ subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, tim endif if (.not.file_is_open(IO_handle)) return - call get_file_info(IO_handle, ndim, nvar, natt, ntime) + call get_file_info(IO_handle, nvar=nvar) if (nvar == -1) then write (mesg,*) "Reopening file ",trim(filename)," apparently had ",nvar,& @@ -1343,7 +1343,7 @@ end subroutine query_vardesc !> Write a 4d field to an output file, potentially with rotation subroutine MOM_write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & fill_value, turns, scale) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:,:), intent(inout) :: field !< Unrotated field to write @@ -1378,7 +1378,7 @@ end subroutine MOM_write_field_4d !> Write a 3d field to an output file, potentially with rotation subroutine MOM_write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & fill_value, turns, scale) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:), intent(inout) :: field !< Unrotated field to write @@ -1413,7 +1413,7 @@ end subroutine MOM_write_field_3d !> Write a 2d field to an output file, potentially with rotation subroutine MOM_write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & fill_value, turns, scale) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:), intent(inout) :: field !< Unrotated field to write @@ -1447,7 +1447,7 @@ end subroutine MOM_write_field_2d !> Write a 1d field to an output file subroutine MOM_write_field_1d(IO_handle, field_md, field, tstamp, fill_value, scale) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, dimension(:), intent(in) :: field !< Field to write real, optional, intent(in) :: tstamp !< Model timestamp @@ -1476,7 +1476,7 @@ end subroutine MOM_write_field_1d !> Write a 0d field to an output file subroutine MOM_write_field_0d(IO_handle, field_md, field, tstamp, fill_value, scale) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, intent(in) :: field !< Field to write real, optional, intent(in) :: tstamp !< Model timestamp From 0fdc5c43816c716509a97b2ecaa9506375fffcca Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 18 Mar 2021 18:27:08 -0400 Subject: [PATCH 03/10] Correct warnings from categorize_axis Corrected the logic of a warning message in categorize_axis by adding parentheses. All answers are bitwise identical, and spurious warnings are no longer being issued. --- config_src/infra/FMS2/MOM_read_data_fms2.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config_src/infra/FMS2/MOM_read_data_fms2.F90 b/config_src/infra/FMS2/MOM_read_data_fms2.F90 index 83a10e7e30..4732c019f4 100644 --- a/config_src/infra/FMS2/MOM_read_data_fms2.F90 +++ b/config_src/infra/FMS2/MOM_read_data_fms2.F90 @@ -233,7 +233,7 @@ subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t endif ; enddo endif - if (.not.(x_found .and. y_found) .and. (ndims>2) .or. ((ndims==2) .and. .not.is_t(ndims))) then + if (.not.(x_found .and. y_found) .and. ((ndims>2) .or. ((ndims==2) .and. .not.is_t(ndims)))) then ! This is a case where one would expect to find x-and y-dimensions, but none have been found. if (is_root_pe()) then dim_list = trim(dim_names(1))//", "//trim(dim_names(2)) From 255233ba3b1de88cc772887062bd6666fb1b5ddf Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 18 Mar 2021 18:44:25 -0400 Subject: [PATCH 04/10] Added code to write via FMS2 interfaces Added a large number of calls to handle all of the writes via the FMS2 interfaces to infra/FMS2/MOM_io_infra.F90. There are newly defined private types in MOM_io_infra to wrap the axistype and fieldtype that had previously been offered from mpp_io_mod. All answers are bitwise identical and it has been verified that output files do not change and the restarts are still working. --- config_src/infra/FMS2/MOM_io_infra.F90 | 520 +++++++++++++++++++++---- 1 file changed, 452 insertions(+), 68 deletions(-) diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index 8115d4acfa..d4fe0b5387 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -8,8 +8,13 @@ module MOM_io_infra use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, FATAL, WARNING use MOM_read_data_fms2, only : prepare_to_read_var -use fms2_io_mod, only : fms2_open_file => open_file, fms2_close_file => close_file -use fms2_io_mod, only : FmsNetcdfDomainFile_t, fms2_read_data => read_data, check_if_open +use fms2_io_mod, only : fms2_open_file => open_file, check_if_open, fms2_close_file => close_file +use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, fms2_read_data => read_data +use fms2_io_mod, only : get_unlimited_dimension_name, get_num_dimensions, get_num_variables +use fms2_io_mod, only : get_variable_names, variable_exists, get_variable_size +use fms2_io_mod, only : register_field, write_data, register_variable_attribute +use fms2_io_mod, only : variable_att_exists, get_variable_attribute, get_variable_num_dimensions +use fms2_io_mod, only : get_dimension_size, is_dimension_registered, register_axis, unlimited use fms_mod, only : write_version_number, open_namelist_file, check_nml_error use fms_io_mod, only : file_exist, field_exist, field_size, read_data @@ -17,8 +22,8 @@ module MOM_io_infra use mpp_io_mod, only : mpp_open, mpp_close, mpp_flush use mpp_io_mod, only : mpp_write_meta, mpp_write use mpp_io_mod, only : mpp_get_atts, mpp_attribute_exist -use mpp_io_mod, only : mpp_get_axes, axistype, mpp_get_axis_data -use mpp_io_mod, only : mpp_get_fields, fieldtype +use mpp_io_mod, only : mpp_get_axes, mpp_axistype=>axistype, mpp_get_axis_data +use mpp_io_mod, only : mpp_get_fields, mpp_fieldtype=>fieldtype use mpp_io_mod, only : mpp_get_info, mpp_get_times use mpp_io_mod, only : mpp_io_init ! These are encoding constants. @@ -38,7 +43,7 @@ module MOM_io_infra public :: io_infra_init, io_infra_end, MOM_namelist_file, check_namelist_error, write_version ! These types are inherited from underlying infrastructure code, to act as containers for ! information about fields and axes, respectively, and are opaque to this module. -public :: fieldtype, axistype +! public :: file_type, fieldtype, axistype ! These are encoding constant parmeters. public :: ASCII_FILE, NETCDF_FILE, SINGLE_FILE, MULTIPLE public :: APPEND_FILE, READONLY_FILE, OVERWRITE_FILE, WRITEONLY_FILE @@ -99,13 +104,37 @@ module MOM_io_infra !> Type for holding a handle to an open file and related information type, public :: file_type ; private integer :: unit = -1 !< The framework identfier or netCDF unit number of an output file + type(FmsNetcdfDomainFile_t), pointer :: fileobj => NULL() !< A domain-decomposed + !! file object that is open for writing character(len=:), allocatable :: filename !< The path to this file, if it is open logical :: open_to_read = .false. !< If true, this file or fileset can be read logical :: open_to_write = .false. !< If true, this file or fileset can be written to + integer :: num_times !< The number of time levels in this file + real :: file_time !< The time of the latest entry in the file. + logical :: FMS2_file !< If true, this file-type is to be used with FMS2 interfaces. end type file_type +!> This type is a container for information about a variable in a file. +type, public :: fieldtype ; private + character(len=256) :: name !< The name of this field in the files. + type(mpp_fieldtype) :: FT !< The FMS1 field-type that this type wraps + character(len=:), allocatable :: longname !< The long name for this field + character(len=:), allocatable :: units !< The units for this field + integer(kind=int64) :: chksum_read !< A checksum that has been read from a file + logical :: valid_chksum !< If true, this field has a valid checksum value. + logical :: FMS2_field !< If true, this field-type should be used with FMS2 interfaces. +end type fieldtype + +!> This type is a container for information about an axis in a file. +type, public :: axistype ; private + character(len=256) :: name !< The name of this axis in the files. + type(mpp_axistype) :: AT !< The FMS1 axis-type that this type wraps + real, allocatable, dimension(:) :: ax_data !< The values of the data on the axis. +end type axistype + !> For now, this is hard-coded to exercise the new FMS2 interfaces. logical :: FMS2_reads = .true. +logical :: FMS2_writes = .true. contains @@ -115,17 +144,11 @@ subroutine read_field_chksum(field, chksum, valid_chksum) type(fieldtype), intent(in) :: field !< The field whose checksum attribute is to be read. integer(kind=int64), intent(out) :: chksum !< The checksum for the field. logical, intent(out) :: valid_chksum !< If true, chksum has been successfully read. - ! Local variables - integer(kind=int64), dimension(3) :: checksum_file - checksum_file(:) = -1 - valid_chksum = mpp_attribute_exist(field, "checksum") - if (valid_chksum) then - call get_field_atts(field, checksum=checksum_file) - chksum = checksum_file(1) - else - chksum = -1 - endif + chksum = -1 + valid_chksum = field%valid_chksum + if (valid_chksum) chksum = field%chksum_read + end subroutine read_field_chksum !> Returns true if the named file or its domain-decomposed variant exists. @@ -156,7 +179,7 @@ end function FMS_file_exists logical function file_is_open(IO_handle) type(file_type), intent(in) :: IO_handle !< Handle to a file to inquire about - file_is_open = (IO_handle%unit >= 0) + file_is_open = ((IO_handle%unit >= 0) .or. associated(IO_handle%fileobj)) end function file_is_open !> closes a file (or fileset). If the file handle does not point to an open file, @@ -164,9 +187,16 @@ end function file_is_open subroutine close_file_type(IO_handle) type(file_type), intent(inout) :: IO_handle !< The I/O handle for the file to be closed - call mpp_close(IO_handle%unit) + if (associated(IO_handle%fileobj)) then + call fms2_close_file(IO_handle%fileobj) + deallocate(IO_handle%fileobj) + else + call mpp_close(IO_handle%unit) + endif if (allocated(IO_handle%filename)) deallocate(IO_handle%filename) IO_handle%open_to_read = .false. ; IO_handle%open_to_write = .false. + IO_handle%num_times = 0 ; IO_handle%file_time = 0.0 + IO_handle%FMS2_file = .false. end subroutine close_file_type !> closes a file. If the unit does not point to an open file, @@ -178,10 +208,14 @@ subroutine close_file_unit(unit) end subroutine close_file_unit !> Ensure that the output stream associated with a file handle is fully sent to disk. -subroutine flush_file_type(file) - type(file_type), intent(in) :: file !< The I/O handle for the file to flush +subroutine flush_file_type(IO_handle) + type(file_type), intent(in) :: IO_handle !< The I/O handle for the file to flush - call mpp_flush(file%unit) + if (associated(IO_handle%fileobj)) then + ! There does not appear to be an fms2 flush call. + else + call mpp_flush(IO_handle%unit) + endif end subroutine flush_file_type !> Ensure that the output stream associated with a unit is fully sent to disk. @@ -272,20 +306,70 @@ subroutine open_file_type(IO_handle, filename, action, MOM_domain, threading, fi !! to threading=MULTIPLE write to the same file (SINGLE_FILE) !! or to one file per PE (MULTIPLE, the default). - if (present(MOM_Domain)) then - call mpp_open(IO_handle%unit, filename, action=action, form=NETCDF_FILE, threading=threading, & + ! Local variables + type(FmsNetcdfDomainFile_t) :: fileobj_read ! A handle to a domain-decomposed file for obtaining information + ! about the exiting time axis entries in append mode. + logical :: success ! If true, the file was opened successfully + integer :: file_mode ! An integer that encodes whether the file is to be opened for + ! reading, writing or appending + character(len=40) :: mode ! A character string that encodes whether the file is to be opened for + ! reading, writing or appending + character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file + + if (IO_handle%open_to_write) then + call MOM_error(WARNING, "open_file_type called for file "//trim(filename)//& + " with an IO_handle that is already open to to write.") + return + endif + if (IO_handle%open_to_read) then + call MOM_error(FATAL, "open_file_type called for file "//trim(filename)//& + " with an IO_handle that is already open to to read.") + endif + + file_mode = WRITEONLY_FILE ; if (present(action)) file_mode = action + + if (FMS2_writes .and. present(MOM_Domain)) then + if (.not.associated(IO_handle%fileobj)) allocate (IO_handle%fileobj) + + if (file_mode == WRITEONLY_FILE) then ; mode = "write" + elseif (file_mode == APPEND_FILE) then ; mode = "append" + elseif (file_mode == OVERWRITE_FILE) then ; mode = "overwrite" + elseif (file_mode == READONLY_FILE) then ; mode = "read" + else + call MOM_error(FATAL, "open_file_type called with unrecognized action.") + endif + + IO_handle%num_times = 0 + IO_handle%file_time = 0.0 + if ((file_mode == APPEND_FILE) .and. file_exists(filename, MOM_Domain)) then + ! Determine the latest file time and number of records so far. + success = fms2_open_file(fileObj_read, trim(filename), "read", MOM_domain%mpp_domain) + call get_unlimited_dimension_name(fileObj_read, dim_unlim_name) + if (len_trim(dim_unlim_name) > 0) & + call get_dimension_size(fileObj_read, trim(dim_unlim_name), IO_handle%num_times) + if (IO_handle%num_times > 0) & + call fms2_read_data(fileObj_read, trim(dim_unlim_name), IO_handle%file_time, & + unlim_dim_level=IO_handle%num_times) + call fms2_close_file(fileObj_read) + endif + + success = fms2_open_file(IO_handle%fileobj, trim(filename), trim(mode), & + MOM_domain%mpp_domain, is_restart=.false.) + if (.not.success) call MOM_error(FATAL, "Unable to open file "//trim(filename)) + IO_handle%FMS2_file = .true. + elseif (present(MOM_Domain)) then + call mpp_open(IO_handle%unit, filename, action=file_mode, form=NETCDF_FILE, threading=threading, & fileset=fileset, domain=MOM_Domain%mpp_domain) + IO_handle%FMS2_file = .false. else - call mpp_open(IO_handle%unit, filename, action=action, form=NETCDF_FILE, threading=threading, & + call mpp_open(IO_handle%unit, filename, action=file_mode, form=NETCDF_FILE, threading=threading, & fileset=fileset) + IO_handle%FMS2_file = .false. endif IO_handle%filename = trim(filename) - if (present(action)) then - if (action == READONLY_FILE) then - IO_handle%open_to_read = .true. ; IO_handle%open_to_write = .false. - else - IO_handle%open_to_read = .false. ; IO_handle%open_to_write = .true. - endif + + if (file_mode == READONLY_FILE) then + IO_handle%open_to_read = .true. ; IO_handle%open_to_write = .false. else IO_handle%open_to_read = .false. ; IO_handle%open_to_write = .true. endif @@ -319,43 +403,59 @@ subroutine get_filename_suffix(suffix) end subroutine get_filename_suffix -!> Get information about the number of dimensions, variables, global attributes and time levels +!> Get information about the number of dimensions, variables and time levels !! in the file associated with an open file unit -subroutine get_file_info(IO_handle, ndim, nvar, natt, ntime) +subroutine get_file_info(IO_handle, ndim, nvar, ntime) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O integer, optional, intent(out) :: ndim !< The number of dimensions in the file integer, optional, intent(out) :: nvar !< The number of variables in the file - integer, optional, intent(out) :: natt !< The number of global attributes in the file integer, optional, intent(out) :: ntime !< The number of time levels in the file ! Local variables + character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file integer :: ndims, nvars, natts, ntimes - call mpp_get_info(IO_handle%unit, ndims, nvars, natts, ntimes ) + if (IO_handle%FMS2_file) then + if (present(ndim)) ndim = get_num_dimensions(IO_handle%fileobj) + if (present(nvar)) nvar = get_num_variables(IO_handle%fileobj) + if (present(ntime)) then + ntime = 0 + call get_unlimited_dimension_name(IO_handle%fileobj, dim_unlim_name) + if (len_trim(dim_unlim_name) > 0) & + call get_dimension_size(IO_handle%fileobj, trim(dim_unlim_name), ntime) + endif + else + call mpp_get_info(IO_handle%unit, ndims, nvars, natts, ntimes ) - if (present(ndim)) ndim = ndims - if (present(nvar)) nvar = nvars - if (present(natt)) natt = natts - if (present(ntime)) ntime = ntimes + if (present(ndim)) ndim = ndims + if (present(nvar)) nvar = nvars + if (present(ntime)) ntime = ntimes + endif end subroutine get_file_info !> Get the times of records from a file - !### Modify this to also convert to time_type, using information about the dimensions? subroutine get_file_times(IO_handle, time_values, ntime) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O real, allocatable, dimension(:), intent(inout) :: time_values !< The real times for the records in file. integer, optional, intent(out) :: ntime !< The number of time levels in the file + character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file integer :: ntimes + !### Modify this routine to optionally convert to time_type, using information about the dimensions? if (allocated(time_values)) deallocate(time_values) call get_file_info(IO_handle, ntime=ntimes) if (present(ntime)) ntime = ntimes if (ntimes > 0) then allocate(time_values(ntimes)) - call mpp_get_times(IO_handle%unit, time_values) + if (IO_handle%FMS2_file) then + call get_unlimited_dimension_name(IO_handle%fileobj, dim_unlim_name) + call fms2_read_data(IO_handle%fileobj, trim(dim_unlim_name), time_values) + else + call mpp_get_times(IO_handle%unit, time_values) + endif endif end subroutine get_file_times @@ -365,7 +465,45 @@ subroutine get_file_fields(IO_handle, fields) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O type(fieldtype), dimension(:), intent(inout) :: fields !< Field-type descriptions of all of !! the variables in a file. - call mpp_get_fields(IO_handle%unit, fields) + type(mpp_fieldtype), dimension(size(fields)) :: mpp_fields + character(len=256), dimension(size(fields)) :: var_names + character(len=256) :: units + character(len=2048) :: longname + integer(kind=int64), dimension(3) :: checksum_file + integer :: i, nvar + + nvar = size(fields) + ! Local variables + if (IO_handle%FMS2_file) then + call get_variable_names(IO_handle%fileobj, var_names) + do i=1,nvar + fields(i)%name = trim(var_names(i)) + longname = "" + call get_variable_attribute(IO_handle%fileobj, var_names(i), 'long_name', longname) + fields(i)%longname = trim(longname) + units = "" + call get_variable_attribute(IO_handle%fileobj, var_names(i), 'units', units) + fields(i)%units = trim(units) + + fields(i)%valid_chksum = variable_att_exists(IO_handle%fileobj, var_names(i), "checksum") + if (fields(i)%valid_chksum) then + call get_variable_attribute(IO_handle%fileobj, var_names(i), 'checksum', checksum_file) + fields(i)%chksum_read = checksum_file(1) + endif + enddo + else + call mpp_get_fields(IO_handle%unit, mpp_fields) + do i=1,nvar + fields(i)%FT = mpp_fields(i) + call mpp_get_atts(fields(i)%FT, name=fields(i)%name, units=units, longname=longname, & + checksum=checksum_file) + fields(i)%longname = trim(longname) + fields(i)%units = trim(units) + fields(i)%valid_chksum = mpp_attribute_exist(fields(i)%FT, "checksum") + if (fields(i)%valid_chksum) fields(i)%chksum_read = checksum_file(1) + enddo + endif + end subroutine get_file_fields !> Extract information from a field type, as stored or as found in a file @@ -376,7 +514,12 @@ subroutine get_field_atts(field, name, units, longname, checksum) character(len=*), optional, intent(out) :: longname !< The long name of the variable integer(kind=int64), dimension(:), & optional, intent(out) :: checksum !< The checksums of the variable in a file - call mpp_get_atts(field, name=name, units=units, longname=longname, checksum=checksum) + + if (present(name)) name = trim(field%name) + if (present(units)) units = trim(field%units) + if (present(longname)) longname = trim(field%longname) + if (present(checksum)) checksum = field%chksum_read + end subroutine get_field_atts !> Field_exists returns true if the field indicated by field_name is present in the @@ -389,7 +532,44 @@ function field_exists(filename, field_name, domain, no_domain, MOM_domain) type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition logical :: field_exists !< True if filename exists and field_name is in filename - if (present(MOM_domain)) then + ! Local variables + type(FmsNetcdfDomainFile_t) :: fileObj_dd ! A handle to a domain-decomposed file for obtaining information + ! about the exiting time axis entries in append mode. + type(FmsNetcdfFile_t) :: fileObj_simple ! A handle to a non-domain-decomposed file for obtaining information + ! about the exiting time axis entries in append mode. + logical :: success ! If true, the file was opened successfully + logical :: domainless ! If true, this file does not use a domain-decomposed file. + + domainless = .not.(present(MOM_domain) .or. present(domain)) + if (present(no_domain)) then + if (domainless .and. .not.no_domain) call MOM_error(FATAL, & + "field_exists: When no_domain is present and false, a domain must be supplied in query about "//& + trim(field_name)//" in file "//trim(filename)) + domainless = no_domain + endif + + if (FMS2_reads) then + field_exists = .false. + if (file_exists(filename)) then + if (domainless) then + success = fms2_open_file(fileObj_simple, trim(filename), "read") + if (success) then + field_exists = variable_exists(fileObj_simple, field_name) + call fms2_close_file(fileObj_simple) + endif + else + if (present(MOM_domain)) then + success = fms2_open_file(fileObj_dd, trim(filename), "read", MOM_domain%mpp_domain) + else + success = fms2_open_file(fileObj_dd, trim(filename), "read", domain) + endif + if (success) then + field_exists = variable_exists(fileobj_dd, field_name) + call fms2_close_file(fileObj_dd) + endif + endif + endif + elseif (present(MOM_domain)) then field_exists = field_exist(filename, field_name, domain=MOM_domain%mpp_domain, no_domain=no_domain) else field_exists = field_exist(filename, field_name, domain=domain, no_domain=no_domain) @@ -408,7 +588,32 @@ subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain) logical, optional, intent(in) :: no_domain !< If present and true, do not check for file !! names with an appended tile number - call field_size(filename, fieldname, sizes, field_found=field_found, no_domain=no_domain) + ! Local variables + type(FmsNetcdfFile_t) :: fileobj_read ! A handle to a non-domain-decomposed file for obtaining information + ! about the exiting time axis entries in append mode. + logical :: success ! If true, the file was opened successfully + logical :: field_exists ! True if filename exists and field_name is in filename + integer :: i, ndims + + if (FMS2_reads) then + field_exists = .false. + if (file_exists(filename)) then + success = fms2_open_file(fileObj_read, trim(filename), "read") + if (success) then + field_exists = variable_exists(fileobj_read, fieldname) + if (field_exists) then + ndims = get_variable_num_dimensions(fileobj_read, fieldname) + if (ndims > size(sizes)) call MOM_error(FATAL, & + "get_field_size called with too few sizes for "//trim(fieldname)//" in "//trim(filename)) + call get_variable_size(fileobj_read, fieldname, sizes(1:ndims)) + do i=ndims+1,size(sizes) ; sizes(i) = 0 ; enddo + endif + endif + endif + if (present(field_found)) field_found = field_exists + else + call field_size(filename, fieldname, sizes, field_found=field_found, no_domain=no_domain) + endif end subroutine get_field_size @@ -417,7 +622,16 @@ subroutine get_axis_data( axis, dat ) type(axistype), intent(in) :: axis !< An axis type real, dimension(:), intent(out) :: dat !< The data in the axis variable - call mpp_get_axis_data( axis, dat ) + integer :: i + + if (allocated(axis%ax_data)) then + if (size(axis%ax_data) > size(dat)) call MOM_error(FATAL, & + "get_axis_data called with too small of an output data array for "//trim(axis%name)) + do i=1,size(axis%ax_data) ; dat(i) = axis%ax_data(i) ; enddo + elseif (.not.FMS2_writes) then + call mpp_get_axis_data( axis%AT, dat ) + endif + end subroutine get_axis_data !> This routine uses the fms_io subroutine read_data to read a scalar named @@ -587,7 +801,7 @@ subroutine MOM_read_data_2d_region(filename, fieldname, data, start, nread, MOM_ real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. - ! This subroutine does not have an FMS-2 variant yet. + !### This subroutine does not have an FMS-2 variant yet. if (present(MOM_Domain)) then call read_data(filename, fieldname, data, start, nread, domain=MOM_Domain%mpp_domain, & @@ -716,6 +930,7 @@ subroutine MOM_read_data_0d_int(filename, fieldname, data, timelevel) integer, intent(inout) :: data !< The 1-dimensional array into which the data integer, optional, intent(in) :: timelevel !< The time level in the file to read + !### This needs an FMS2 variant, eventually. call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) end subroutine MOM_read_data_0d_int @@ -728,6 +943,7 @@ subroutine MOM_read_data_1d_int(filename, fieldname, data, timelevel) integer, dimension(:), intent(inout) :: data !< The 1-dimensional array into which the data integer, optional, intent(in) :: timelevel !< The time level in the file to read + !### This needs an FMS2 variant, eventually. call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) end subroutine MOM_read_data_1d_int @@ -875,7 +1091,7 @@ end subroutine MOM_read_vector_3d !> Write a 4d field to an output file. subroutine write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:,:), intent(inout) :: field !< Field to write @@ -883,13 +1099,24 @@ subroutine write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value - call mpp_write(IO_handle%unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & - tile_count=tile_count, default_data=fill_value) + + ! Local variables + integer :: time_index + + if (IO_handle%FMS2_file .and. present(tstamp)) then + time_index = write_time_if_later(IO_handle, tstamp) + call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index) + elseif (IO_handle%FMS2_file) then + call write_data(IO_handle%fileobj, trim(field_md%name), field) + else + call mpp_write(IO_handle%unit, field_md%FT, MOM_domain%mpp_domain, field, tstamp=tstamp, & + tile_count=tile_count, default_data=fill_value) + endif end subroutine write_field_4d !> Write a 3d field to an output file. subroutine write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:,:), intent(inout) :: field !< Field to write @@ -897,13 +1124,23 @@ subroutine write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value - call mpp_write(IO_handle%unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & + ! Local variables + integer :: time_index + + if (IO_handle%FMS2_file .and. present(tstamp)) then + time_index = write_time_if_later(IO_handle, tstamp) + call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index) + elseif (IO_handle%FMS2_file) then + call write_data(IO_handle%fileobj, trim(field_md%name), field) + else + call mpp_write(IO_handle%unit, field_md%FT, MOM_domain%mpp_domain, field, tstamp=tstamp, & tile_count=tile_count, default_data=fill_value) + endif end subroutine write_field_3d !> Write a 2d field to an output file. subroutine write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition real, dimension(:,:), intent(inout) :: field !< Field to write @@ -911,36 +1148,92 @@ subroutine write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_c integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) real, optional, intent(in) :: fill_value !< Missing data fill value - call mpp_write(IO_handle%unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & + ! Local variables + integer :: time_index + + if (IO_handle%FMS2_file .and. present(tstamp)) then + time_index = write_time_if_later(IO_handle, tstamp) + call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index) + elseif (IO_handle%FMS2_file) then + call write_data(IO_handle%fileobj, trim(field_md%name), field) + else + call mpp_write(IO_handle%unit, field_md%FT, MOM_domain%mpp_domain, field, tstamp=tstamp, & tile_count=tile_count, default_data=fill_value) + endif end subroutine write_field_2d !> Write a 1d field to an output file. subroutine write_field_1d(IO_handle, field_md, field, tstamp) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, dimension(:), intent(in) :: field !< Field to write real, optional, intent(in) :: tstamp !< Model time of this field - call mpp_write(IO_handle%unit, field_md, field, tstamp=tstamp) + ! Local variables + integer :: time_index + + if (IO_handle%FMS2_file .and. present(tstamp)) then + time_index = write_time_if_later(IO_handle, tstamp) + call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index) + elseif (IO_handle%FMS2_file) then + call write_data(IO_handle%fileobj, trim(field_md%name), field) + else + call mpp_write(IO_handle%unit, field_md%FT, field, tstamp=tstamp) + endif end subroutine write_field_1d !> Write a 0d field to an output file. subroutine write_field_0d(IO_handle, field_md, field, tstamp) - type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(in) :: field_md !< Field type with metadata real, intent(in) :: field !< Field to write real, optional, intent(in) :: tstamp !< Model time of this field - call mpp_write(IO_handle%unit, field_md, field, tstamp=tstamp) + ! Local variables + integer :: time_index + + if (IO_handle%FMS2_file .and. present(tstamp)) then + time_index = write_time_if_later(IO_handle, tstamp) + call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index) + elseif (IO_handle%FMS2_file) then + call write_data(IO_handle%fileobj, trim(field_md%name), field) + else + call mpp_write(IO_handle%unit, field_md%FT, field, tstamp=tstamp) + endif end subroutine write_field_0d +!> Returns the integer time index for a write in this file, also writing the time variable to +!! the file if this time is later than what is already in the file. +integer function write_time_if_later(IO_handle, field_time) + type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing + real, intent(in) :: field_time !< Model time of this field + + ! Local variables + character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file + + if ((field_time > IO_handle%file_time) .or. (IO_handle%num_times == 0)) then + IO_handle%file_time = field_time + IO_handle%num_times = IO_handle%num_times + 1 + if (IO_handle%FMS2_file) then + call get_unlimited_dimension_name(IO_handle%fileobj, dim_unlim_name) + call write_data(IO_handle%fileobj, trim(dim_unlim_name), (/field_time/), & + corner=(/IO_handle%num_times/), edge_lengths=(/1/)) + endif + endif + + write_time_if_later = IO_handle%num_times +end function write_time_if_later + !> Write the data for an axis subroutine MOM_write_axis(IO_handle, axis) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(axistype), intent(in) :: axis !< An axis type variable with information to write - call mpp_write(IO_handle%unit, axis) + if (IO_handle%FMS2_file) then + call write_data(IO_handle%fileobj, trim(axis%name), axis%ax_data) + else + call mpp_write(IO_handle%unit, axis%AT) + endif end subroutine MOM_write_axis @@ -963,26 +1256,80 @@ subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian logical, optional, intent(in) :: edge_axis !< If true, this axis marks an edge of the tracer cells character(len=*), optional, intent(in) :: calendar !< The name of the calendar used with a time axis - call mpp_write_meta(IO_handle%unit, axis, name, units, longname, cartesian=cartesian, sense=sense, & - domain=domain, data=data, calendar=calendar) + character(len=:), allocatable :: cart ! A left-adjusted and trimmed copy of cartesian + logical :: is_x, is_y, is_t ! If true, this is a domain-decomposed axis in one of the directions. + integer :: position ! A flag indicating the axis staggering position. + + if (IO_handle%FMS2_file) then + if (is_dimension_registered(IO_handle%fileobj, trim(name))) then + call MOM_error(FATAL, "write_metadata_axis was called more than once for axis "//trim(name)//& + " in file "//trim(IO_handle%filename)) + return + endif + endif + + axis%name = trim(name) + if (present(data)) then + if (allocated(axis%ax_data)) call MOM_error(FATAL, & + "Data is already allocated in a call to write_metadata_axis for axis "//& + trim(name)//" in file "//trim(IO_handle%filename)) + allocate(axis%ax_data(size(data))) ; axis%ax_data(:) = data(:) + endif + + if (IO_handle%FMS2_file) then + is_x = .false. ; is_y = .false. ; is_t = .false. + position = CENTER + if (present(cartesian)) then + cart = trim(adjustl(cartesian)) + if ((index(cart, "X") == 1) .or. (index(cart, "x") == 1)) is_x = .true. + if ((index(cart, "Y") == 1) .or. (index(cart, "y") == 1)) is_y = .true. + if ((index(cart, "T") == 1) .or. (index(cart, "t") == 1)) is_t = .true. + endif + + if (is_x) then + if (present(edge_axis)) then ; if (edge_axis) position = EAST_FACE ; endif + call register_axis(IO_handle%fileobj, trim(name), 'x', domain_position=position) + elseif (is_y) then + if (present(edge_axis)) then ; if (edge_axis) position = NORTH_FACE ; endif + call register_axis(IO_handle%fileobj, trim(name), 'y', domain_position=position) + elseif (is_t .and. .not.present(data)) then + ! This is the unlimited (time) dimension. + call register_axis(IO_handle%fileobj, trim(name), unlimited) + else + if (.not.present(data)) call MOM_error(FATAL,"MOM_io:register_diagnostic_axis: "//& + "An axis_length argument is required to register the axis "//trim(name)) + call register_axis(IO_handle%fileobj, trim(name), size(data)) + endif + + ! Now create the variable that describes this axis. + call register_field(IO_handle%fileobj, trim(name), "double", dimensions=(/name/)) + if (len_trim(units) > 0) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'units', & + trim(units), len_trim(units)) + if (len_trim(longname) > 0) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'long_name', & + trim(longname), len_trim(longname)) + if (present(cartesian)) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'cartesian_axis', & + trim(cartesian), len_trim(cartesian)) + if (present(sense)) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'sense', sense) + else + call mpp_write_meta(IO_handle%unit, axis%AT, name, units, longname, cartesian=cartesian, sense=sense, & + domain=domain, data=data, calendar=calendar) + endif end subroutine write_metadata_axis !> Store information about an output variable in a previously defined fieldtype and write this !! information to the file indicated by unit. subroutine write_metadata_field(IO_handle, field, axes, name, units, longname, & - min, max, fill, scale, add, pack, standard_name, checksum) + pack, standard_name, checksum) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(fieldtype), intent(inout) :: field !< The fieldtype where this information is stored type(axistype), dimension(:), intent(in) :: axes !< Handles for the axis used for this variable character(len=*), intent(in) :: name !< The name in the file of this variable character(len=*), intent(in) :: units !< The units of this variable character(len=*), intent(in) :: longname !< The long description of this variable - real, optional, intent(in) :: min !< The minimum valid value for this variable - real, optional, intent(in) :: max !< The maximum valid value for this variable - real, optional, intent(in) :: fill !< Missing data fill value - real, optional, intent(in) :: scale !< An multiplicative factor by which to scale - !! the variable before output - real, optional, intent(in) :: add !< An offset to add to the variable before output integer, optional, intent(in) :: pack !< A precision reduction factor with which the !! variable. The default, 1, has no reduction, !! but 2 is not uncommon. @@ -990,9 +1337,46 @@ subroutine write_metadata_field(IO_handle, field, axes, name, units, longname, & integer(kind=int64), dimension(:), & optional, intent(in) :: checksum !< Checksum values that can be used to verify reads. + ! Local variables + character(len=256), dimension(size(axes)) :: dim_names ! The names of the dimensions + type(mpp_axistype), dimension(size(axes)) :: mpp_axes ! The array of mpp_axistypes for this variable + character(len=16) :: prec_string ! A string specifying the precision with which to save this variable + character(len=64) :: checksum_string ! checksum character array created from checksum argument + integer :: i, ndims + + ndims = size(axes) + if (IO_handle%FMS2_file) then + do i=1,ndims ; dim_names(i) = trim(axes(i)%name) ; enddo + prec_string = "double" ; if (present(pack)) then ; if (pack > 1) prec_string = "float" ; endif + call register_field(IO_handle%fileobj, trim(name), trim(prec_string), dimensions=dim_names) + if (len_trim(units) > 0) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'units', & + trim(units), len_trim(units)) + if (len_trim(longname) > 0) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'long_name', & + trim(longname), len_trim(longname)) + if (present(standard_name)) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'standard_name', & + trim(standard_name), len_trim(standard_name)) + if (present(checksum)) then + write (checksum_string,'(Z16)') checksum(1) ! Z16 is the hexadecimal format code + call register_variable_attribute(IO_handle%fileobj, trim(name), "checksum", & + trim(checksum_string), len_trim(checksum_string)) + endif + !### Add more attributes if they are present; remove attributes that are never used. + else + do i=1,ndims ; mpp_axes(i) = axes(i)%AT ; enddo + call mpp_write_meta(IO_handle%unit, field%FT, mpp_axes, name, units, longname, & + pack=pack, standard_name=standard_name, checksum=checksum) + ! unused opt. args: min=min, max=max, fill=fill, scale=scale, add=add, & + endif - call mpp_write_meta(IO_handle%unit, field, axes, name, units, longname, min=min, max=max, & - fill=fill, scale=scale, add=add, pack=pack, standard_name=standard_name, checksum=checksum) + ! Store information in the field-type, regardless of which interfaces are used. + field%name = trim(name) + field%longname = trim(longname) + field%units = trim(units) + field%chksum_read = -1 + field%valid_chksum = .false. end subroutine write_metadata_field From 03e174e1434b139c72ce53b6b53585953cb8d3aa Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 18 Mar 2021 18:52:15 -0400 Subject: [PATCH 05/10] Removed FMS2 MOM_axis and MOM_write_field_fms2 Deleted the unused MOM_axis.F90 and MOM_write_field_fms2.F90 modules, now that MOM_io_infra.F90 has been updated to be able to use FMS2 interfaces for both reading and writing to files. As these modules had never been used, they do not change any answers. --- config_src/infra/FMS2/MOM_axis.F90 | 616 ------- .../infra/FMS2/MOM_write_field_fms2.F90 | 1503 ----------------- 2 files changed, 2119 deletions(-) delete mode 100644 config_src/infra/FMS2/MOM_axis.F90 delete mode 100644 config_src/infra/FMS2/MOM_write_field_fms2.F90 diff --git a/config_src/infra/FMS2/MOM_axis.F90 b/config_src/infra/FMS2/MOM_axis.F90 deleted file mode 100644 index b5d2b3ed88..0000000000 --- a/config_src/infra/FMS2/MOM_axis.F90 +++ /dev/null @@ -1,616 +0,0 @@ -!> This module contains routines that define and register axes to files -module MOM_axis - -! This file is part of MOM6. See LICENSE.md for the license. -use MOM_domains, only : MOM_domain_type -use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING -use MOM_grid, only : ocean_grid_type -use MOM_dyn_horgrid, only : dyn_horgrid_type -use MOM_string_functions, only : lowercase -use MOM_verticalGrid, only : verticalGrid_type -use fms2_io_mod, only : is_dimension_registered, register_axis, is_dimension_unlimited -use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, unlimited -use fms2_io_mod, only : get_variable_size, get_variable_num_dimensions, check_if_open -use fms2_io_mod, only : fms2_open_file=>open_file, fms2_close_file=>close_file -use fms2_io_mod, only : get_variable_dimension_names, read_data, get_unlimited_dimension_name -use fms2_io_mod, only : get_dimension_size -use mpp_domains_mod, only : domain2d, CENTER, CORNER, NORTH_FACE=>NORTH, EAST_FACE=>EAST -use mpp_domains_mod, only : mpp_get_compute_domain -use netcdf -implicit none ; private - -public MOM_register_diagnostic_axis, get_var_dimension_metadata, get_time_units -public MOM_get_diagnostic_axis_data, MOM_register_variable_axes, get_time_index -public convert_checksum_to_string -!> A type for making arrays of pointers to real 1-d arrays -type p1d - real, dimension(:), pointer :: p => NULL() !< A pointer to a 1d array -end type p1d - -!> A structure with information about a single axis variable -type axis_atts - character(len=64) :: name !< Names of the axis - character(len=48) :: units !< Physical dimensions of the axis - character(len=240) :: longname !< Long name of the axis - character(len=8) :: positive !< Positive-definite direction: up, down, east, west, north, south - integer :: horgrid_position !< Horizontal grid position - logical :: is_domain_decomposed !< if .true. the axis data are domain-decomposed - !! and need to be indexed by the compute domain - !! before passing to write_data -end type axis_atts - -!> Type for describing an axis variable (e.g., lath, lonh, Time) -type, public :: axis_data_type - !> An array of descriptions of the registered axes - type(axis_atts), pointer :: axis(:) => NULL() !< structure with axis attributes - type(p1d), pointer :: data(:) => NULL() !< pointer to the axis data -end type axis_data_type - -!> interface for registering axes associated with a variable to a netCDF file object -interface MOM_register_variable_axes - module procedure MOM_register_variable_axes_subdomain - module procedure MOM_register_variable_axes_full -end interface MOM_register_variable_axes - -contains - -!> register a MOM diagnostic axis to a domain-decomposed file -subroutine MOM_register_diagnostic_axis(fileObj, axisName, axisLength) - type(FmsNetcdfDomainFile_t), intent(inout) :: fileObj !< netCDF file object returned by call to open_file - character(len=*), intent(in) :: axisName !< name of the axis to register to file - integer, intent(in), optional :: axisLength !< length of axis/dimension ;only needed for Layer, Interface, Time, - !! Period - select case (trim(lowercase(axisName))) - case ('latq'); call register_axis(fileObj,'latq','y', domain_position=NORTH_FACE) - case ('lath'); call register_axis(fileObj,'lath','y', domain_position=CENTER) - case ('lonq'); call register_axis(fileObj,'lonq','x', domain_position=EAST_FACE) - case ('lonh'); call register_axis(fileObj,'lonh','x', domain_position=CENTER) - case default - if (.not. present(axisLength)) call MOM_error(FATAL,"MOM_io:register_diagnostic_axis: "//& - "An axis_length argument is required to register the axis "//trim(axisName)) - call register_axis(fileObj, trim(axisName), axisLength) - end select -end subroutine MOM_register_diagnostic_axis - - -!> Get the horizontal grid, vertical grid, and/or time dimension names and lengths -!! for a single variable from the hor_grid, t_grid, and z_grid values returned by a prior call to query_vardesc -subroutine get_var_dimension_metadata(hor_grid, z_grid, t_grid_in, & - dim_names, dim_lengths, num_dims, G, dG, GV) - - character(len=*), intent(in) :: hor_grid !< horizontal grid - character(len=*), intent(in) :: z_grid !< vertical grid - character(len=*), intent(in) :: t_grid_in !< time grid - character(len=*), dimension(:), intent(inout) :: dim_names !< array of dimension names - integer, dimension(:), intent(inout) :: dim_lengths !< array of dimension sizes - integer, intent(inout) :: num_dims !< number of axes to register in the restart file - type(ocean_grid_type), optional, intent(in) :: G !< The ocean's grid structure - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure - - ! local - logical :: use_lath - logical :: use_lonh - logical :: use_latq - logical :: use_lonq - character(len=8) :: t_grid - character(len=8) :: t_grid_read - integer :: isg, ieg, jsg, jeg, IsgB, IegB, JsgB, JegB - !integer :: npes - real, pointer, dimension(:) :: gridLatT => NULL(), & ! The latitude or longitude of T or B points for - gridLatB => NULL(), & ! the purpose of labeling the output axes. - gridLonT => NULL(), & - gridLonB => NULL() - type(MOM_domain_type), pointer :: domain => NULL() ! Domain used to get the pe count - - use_lath = .false. - use_lonh = .false. - use_latq = .false. - use_lonq = .false. - - ! set the ocean grid coordinates - - if (present(G)) then - gridLatT => G%gridLatT ; gridLatB => G%gridLatB - gridLonT => G%gridLonT ; gridLonB => G%gridLonB - isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg - IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB - - call get_horizontal_grid_logic(hor_grid, use_lath, use_lonh, use_latq, use_lonq) - elseif (present(dG)) then - gridLatT => dG%gridLatT ; gridLatB => dG%gridLatB - gridLonT => dG%gridLonT ; gridLonB => dG%gridLonB - isg = dG%isg ; ieg = dG%ieg ; jsg = dG%jsg ; jeg = dG%jeg - IsgB = dG%IsgB ; IegB = dG%IegB ; JsgB = dG%JsgB ; JegB = dG%JegB - - call get_horizontal_grid_logic(hor_grid, use_lath, use_lonh, use_latq, use_lonq) - endif - - ! add longitude name to dimension name array - if (use_lonh) then - num_dims = num_dims+1 - dim_names(num_dims) = "" - dim_names(num_dims)(1:len_trim("lonh")) = "lonh" - dim_lengths(num_dims) = size(gridLonT(isg:ieg)) - elseif (use_lonq) then - num_dims = num_dims+1 - dim_names(num_dims) = "" - dim_names(num_dims)(1:len_trim("lonq")) = "lonq" - dim_lengths(num_dims) = size(gridLonB(IsgB:IegB)) - endif - ! add latitude name to dimension name array - if (use_lath) then - num_dims = num_dims+1 - dim_names(num_dims) = "" - dim_names(num_dims)(1:len_trim("lath")) = "lath" - dim_lengths(num_dims) = size(gridLatT(jsg:jeg)) - elseif (use_latq) then - num_dims = num_dims+1 - dim_names(num_dims) = "" - dim_names(num_dims)(1:len_trim("latq")) = "latq" - dim_lengths(num_dims) = size(gridLatB(JsgB:JegB)) - endif - - if (present(GV)) then - ! vertical grid - select case (trim(z_grid)) - case ('L') - num_dims = num_dims+1 - dim_names(num_dims) = "" - dim_names(num_dims)(1:len_trim("Layer")) = "Layer" - dim_lengths(num_dims) = GV%ke - case ('i') - num_dims = num_dims+1 - dim_names(num_dims) = "" - dim_names(num_dims)(1:len_trim("Interface")) = "Interface" - dim_lengths(num_dims) = GV%ke+1 - case ('1') ! Do nothing. - case default - call MOM_error(FATAL, "MOM_io: get_var_dimension_features: "//& - " has an unrecognized z_grid argument"//trim(z_grid)) - end select - endif - ! time - t_grid = adjustl(t_grid_in) - select case (t_grid(1:1)) - case ('s', 'a', 'm') - num_dims = num_dims+1 - dim_names(num_dims) = "" - dim_names(num_dims)(1:len_trim("Time")) = "Time" - dim_lengths(num_dims) = unlimited - case ('p') - if (len_trim(t_grid(2:8)) <= 0) then - call MOM_error(FATAL,"MOM_io:get_var_dimension_features: "//& - "No periodic axis length was specified in "//trim(t_grid)) - endif - num_dims = num_dims+1 - dim_names(num_dims) = "" - dim_names(num_dims)(1:len_trim("Period")) = "Period" - dim_lengths(num_dims) = unlimited - case ('1') ! Do nothing. - case default - call MOM_error(WARNING, "MOM_io: get_var_dimension_metadata: "//& - "Unrecognized t_grid "//trim(t_grid)) - end select -end subroutine get_var_dimension_metadata - - -!> Populate the axis_data structure with axis data and attributes for diagnostic and restart files -subroutine MOM_get_diagnostic_axis_data(axis_data_CS, axis_name, axis_number, G, dG, GV, time_val, time_units) - - type(axis_data_type), intent(inout) :: axis_data_CS !< structure containing the axis data and metadata - character(len=*), intent(in) :: axis_name !< name of the axis - integer, intent(in) :: axis_number !< positional value (wrt to file) of the axis to register - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the file uses any - !! horizontal grid axes. - type(verticalGrid_type), target, optional, intent(in) :: GV !< ocean vertical grid structure - real,dimension(:), target, optional, intent(in) :: time_val !< time value - character(len=*), optional,intent(in) :: time_units!< units for non-periodic time axis - ! local - character(len=40) :: x_axis_units='', y_axis_units='' - integer :: isg, ieg, jsg, jeg, IsgB, IegB, JsgB, JegB - real, pointer, dimension(:) :: gridLatT => NULL(), & ! The latitude or longitude of T or B points for - gridLatB => NULL(), & ! the purpose of labeling the output axes. - gridLonT => NULL(), & - gridLonB => NULL() - - ! initialize axis_data_CS elements - axis_data_CS%axis(axis_number)%name = '' - axis_data_CS%axis(axis_number)%longname = '' - axis_data_CS%axis(axis_number)%units = '' - axis_data_CS%axis(axis_number)%horgrid_position = 0 - axis_data_CS%axis(axis_number)%is_domain_decomposed = .false. - axis_data_CS%axis(axis_number)%positive = '' - axis_data_CS%data(axis_number)%p => NULL() - - ! set the ocean grid coordinates and metadata - if (present(G)) then - gridLatT => G%gridLatT ; gridLatB => G%gridLatB - gridLonT => G%gridLonT ; gridLonB => G%gridLonB - x_axis_units = G%x_axis_units ; y_axis_units = G%y_axis_units - isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg - IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB - elseif (present(dG)) then - gridLatT => dG%gridLatT ; gridLatB => dG%gridLatB - gridLonT => dG%gridLonT ; gridLonB => dG%gridLonB - x_axis_units = dG%x_axis_units ; y_axis_units = dG%y_axis_units - isg = dG%isg ; ieg = dG%ieg ; jsg = dG%jsg ; jeg = dG%jeg - IsgB = dG%IsgB ; IegB = dG%IegB ; JsgB = dG%JsgB ; JegB = dG%JegB - endif - - select case(trim(lowercase(axis_name))) - case('lath') - if (associated(gridLatT)) & - axis_data_CS%data(axis_number)%p=>gridLatT(jsg:jeg) - - axis_data_CS%axis(axis_number)%name = trim(axis_name) - axis_data_CS%axis(axis_number)%longname = 'Latitude' - axis_data_CS%axis(axis_number)%units = y_axis_units - axis_data_CS%axis(axis_number)%horgrid_position = CENTER - axis_data_CS%axis(axis_number)%is_domain_decomposed = .true. - case('lonh') - if (associated(gridLonT)) & - axis_data_CS%data(axis_number)%p=>gridLonT(isg:ieg) - - axis_data_CS%axis(axis_number)%name = trim(axis_name) - axis_data_CS%axis(axis_number)%horgrid_position = CENTER - axis_data_CS%axis(axis_number)%longname = 'Longitude' - axis_data_CS%axis(axis_number)%units = x_axis_units - axis_data_CS%axis(axis_number)%is_domain_decomposed = .true. - case('latq') - if (associated(gridLatB)) & - axis_data_CS%data(axis_number)%p=>gridLatB(JsgB:JegB) - - axis_data_CS%axis(axis_number)%name = trim(axis_name) - axis_data_CS%axis(axis_number)%longname = 'Latitude' - axis_data_CS%axis(axis_number)%units = y_axis_units - axis_data_CS%axis(axis_number)%horgrid_position = NORTH_FACE - axis_data_CS%axis(axis_number)%is_domain_decomposed = .true. - case('lonq') - if (associated(gridLonB)) & - axis_data_CS%data(axis_number)%p=>gridLonB(IsgB:IegB) - - axis_data_CS%axis(axis_number)%name = trim(axis_name) - axis_data_CS%axis(axis_number)%longname = 'Longitude' - axis_data_CS%axis(axis_number)%units = x_axis_units - axis_data_CS%axis(axis_number)%horgrid_position = EAST_FACE - axis_data_CS%axis(axis_number)%is_domain_decomposed = .true. - case('layer') - if (present(GV)) then - axis_data_CS%data(axis_number)%p=>GV%sLayer(1:GV%ke) - axis_data_CS%axis(axis_number)%name = trim(axis_name) - axis_data_CS%axis(axis_number)%longname = 'Layer pseudo-depth, -z*' - axis_data_CS%axis(axis_number)%units = GV%zAxisUnits - axis_data_CS%axis(axis_number)%positive = 'up' - endif - case('interface') - if (present(GV)) then - axis_data_CS%data(axis_number)%p=>GV%sInterface(1:GV%ke+1) - axis_data_CS%axis(axis_number)%name = trim(axis_name) - axis_data_CS%axis(axis_number)%longname = 'Interface pseudo-depth, -z*' - axis_data_CS%axis(axis_number)%units = GV%zAxisUnits - axis_data_CS%axis(axis_number)%positive = 'up' - endif - case('time') - if (.not.(present(time_val))) & - call MOM_error(FATAL, "MOM_io::get_diagnostic_axis_data: requires time_val"//& - " and time_units arguments for "//trim(axis_name)) - - axis_data_CS%data(axis_number)%p=>time_val - axis_data_CS%axis(axis_number)%name = trim(axis_name) - axis_data_CS%axis(axis_number)%longname = 'Time' - - if (present(time_units)) then - axis_data_CS%axis(axis_number)%units = time_units - else - axis_data_CS%axis(axis_number)%units = 'days' - endif - case('period') - if (.not.(present(time_val))) & - call MOM_error(FATAL, "MOM_axis::get_diagnostic_axis_data: requires a time_val argument "// & - "for "//trim(axis_name)) - axis_data_CS%data(axis_number)%p=>time_val - axis_data_CS%axis(axis_number)%name = trim(axis_name) - axis_data_CS%axis(axis_number)%longname = 'Periods for cyclical variables' - case default - call MOM_error(WARNING, "MOM_axis::get_diagnostic_axis_data:"//trim(axis_name)//" is an unrecognized axis") - end select - -end subroutine MOM_get_diagnostic_axis_data - - -!> set the logical variables that determine which diagnositic axes to use -subroutine get_horizontal_grid_logic(grid_string_id, use_lath, use_lonh, use_latq, use_lonq) - character(len=*), intent(in) :: grid_string_id !< horizontal grid string - logical, intent(out) :: use_lath !< if .true., y-axis is oriented in CENTER position - logical, intent(out) :: use_lonh !< if .true., x-axis is oriented in CENTER position - logical, intent(out) :: use_latq !< if .true., y-axis is oriented in NORTH_FACE position - logical, intent(out) :: use_lonq !< if .true., x-axis is oriented in EAST_FACE position - - use_lath = .false. - use_lonh = .false. - use_latq = .false. - use_lonq = .false. - select case (trim(grid_string_id)) - case ('h') ; use_lath = .true. ; use_lonh = .true. ! x=CENTER, y=CENTER - case ('q') ; use_latq = .true. ; use_lonq = .true. ! x=EAST_FACE, y=NORTH_FACE - case ('u') ; use_lath = .true. ; use_lonq = .true. ! x=EAST_FACE, y=CENTER - case ('v') ; use_latq = .true. ; use_lonh = .true. ! x=CENTER, y=NORTH_FACE - case ('T') ; use_lath = .true. ; use_lonh = .true. ! x=CENTER, y=CENTER - case ('Bu') ; use_latq = .true. ; use_lonq = .true. ! x=EAST_FACE, y=NORTH_FACE - case ('Cu') ; use_lath = .true. ; use_lonq = .true. ! x=EAST_FACE, y=CENTER - case ('Cv') ; use_latq = .true. ; use_lonh = .true. ! x=CENTER, y=NORTH_FACE - case ('1') ; ! x=0, y=0 - case default - call MOM_error(FATAL, "MOM_axis:get_var_dimension_features "//& - "Unrecognized hor_grid argument "//trim(grid_string_id)) - end select -end subroutine get_horizontal_grid_logic - -!> Define the time units for the input time value -function get_time_units(time_value) result(time_units_out) - real, intent(in) :: time_value !< numerical time value in seconds - !! i.e., before dividing by 86400. - ! local - character(len=10) :: time_units ! time units - character(len=10) :: time_units_out ! time units trimmed - time_units = '' - time_units_out = '' - if (time_value < 0.0) then - time_units = "days" ! The default value. - elseif (mod(time_value,86400.0)==0.0) then - time_units = "days" - elseif ((time_value >= 0.99) .and. (time_value < 1.01)) then - time_units = "seconds" - elseif ((time_value >= 3599.0) .and. (time_value < 3601.0)) then - time_units = "hours" - elseif ((time_value >= 86399.0) .and. (time_value < 86401.0)) then - time_units = "days" - elseif ((time_value >= 3.0e7) .and. (time_value < 3.2e7)) then - time_units = "years" - else - write(time_units,'(es8.2," s")') time_value - endif - time_units_out = trim(time_units) -end function get_time_units - -!> function to get the index of a time_value from a netCDF file -function get_time_index(filename, time_to_find) result (time_index) - character(len=*) :: filename ! name of the file to read in - real, intent(in) :: time_to_find ! time value to search for in file - ! local - type(fmsNetcdfFile_t) :: fileobj ! netCDF file object returned by open_file - real, allocatable, dimension(:) :: file_times ! array of time values read from file - integer :: dim_unlim_size, i, time_index - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - logical :: file_open_success - - time_index = 1 - dim_unlim_size = 0 - dim_unlim_name = "" - file_open_success = .false. - - if (.not. check_if_open(fileobj)) & - !call MOM_error(FATAL, "get_time_index_nodd: netcdf file object must be open.") - file_open_success=fms2_open_file(fileobj, trim(filename), "read", is_restart=.false.) - - call get_unlimited_dimension_name(fileobj, dim_unlim_name) - call get_dimension_size(fileObj, trim(dim_unlim_name), dim_unlim_size) - ! time index will be one more than the unlimited dimension size if the time_to_find is not in the file - if (dim_unlim_size .gt. 0) then - time_index = dim_unlim_size+1 - allocate(file_times(dim_unlim_size)) - call read_data(fileobj,trim(dim_unlim_name), file_times) - - do i=1,dim_unlim_size - if (ABS(file_times(i)-time_to_find) .gt. TINY(time_to_find)) then - continue - else - time_index = i - exit - endif - enddo - deallocate(file_times) - endif - if (check_if_open(fileobj)) call fms2_close_file(fileobj) -end function get_time_index - -!> register axes associated with a variable from a domain-decomposed netCDF file that are mapped to -!! a sub-domain (e.g., a supergrid). -!> \note The user must specify units for variables with longitude/x-axis and/or latitude/y-axis axes to obtain -!! the correct domain decomposition for the data buffer. -subroutine MOM_register_variable_axes_subdomain(fileObj, variableName, io_domain, position) - type(FmsNetcdfFile_t), intent(inout) :: fileObj !< netCDF file object returned by call to open_file - character(len=*), intent(in) :: variableName !< name of the variable - type(domain2d), intent(in) :: io_domain !< type that contains the mpp io domain - integer, optional, intent(in) :: position !< A flag indicating where this data is discretized - - ! Local variables - character(len=40) :: units ! units corresponding to a specific variable dimension - character(len=40), allocatable, dimension(:) :: dim_names ! variable dimension names - integer :: i, isg, ieg, isc, iec, jsg, jeg, jsc, jec, xlen, ylen - integer :: ndims ! number of dimensions - integer :: pos ! Discrete variable position. Default is CENTER - integer, allocatable, dimension(:) :: dimSizes ! variable dimension sizes - - if (.not. check_if_open(fileObj)) call MOM_error(FATAL,"MOM_axis:register_variable_axes_subdomain: The fileObj "// & - " has not been opened. Call fms2_open_file(fileObj,...) "// & - "before passing the fileObj argument to this function.") - - ! get variable dimension names and lengths - ndims = get_variable_num_dimensions(fileObj, trim(variableName)) - allocate(dimSizes(ndims)) - allocate(dim_names(ndims)) - call get_variable_size(fileObj, trim(variableName), dimSizes, broadcast=.true.) - call get_variable_dimension_names(fileObj, trim(variableName), dim_names) - - ! Get the lengths of the global indicies, using the discrete position of this variable - pos = CORNER ; if (present(position)) pos = position - call mpp_get_compute_domain(io_domain, xsize=xlen, ysize=ylen, position=pos) - ! register the axes - !>\note: This is not a comprehensive check for all possible supported horizontal axes associated with variables - !! read from netCDF files. Developers should add/remove cases as needed. - do i=1,ndims - !if (.not.(is_dimension_registered(fileObj, trim(dim_names(i))))) then - select case(trim(lowercase(dim_names(i)))) - case ("grid_x_t") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case ("nx") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case("nxp") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case("longitude") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case("long") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case("lon") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case("lonh") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case("lonq") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case("xh") - call register_axis(fileObj, trim(dim_names(i)), xlen) - case ("grid_y_t") - call register_axis(fileObj, trim(dim_names(i)), ylen) - case ("ny") - call register_axis(fileObj, trim(dim_names(i)), ylen) - case("nyp") - call register_axis(fileObj, trim(dim_names(i)), ylen) - case("latitude") - call register_axis(fileObj, trim(dim_names(i)), ylen) - case("lat") - call register_axis(fileObj, trim(dim_names(i)), ylen) - case("lath") - call register_axis(fileObj, trim(dim_names(i)), ylen) - case("latq") - call register_axis(fileObj, trim(dim_names(i)), ylen) - case("yh") - call register_axis(fileObj, trim(dim_names(i)), ylen) - case default ! assumes that the axis is not domain-decomposed - if (.not. is_dimension_unlimited(fileObj, trim(dim_names(i)))) & - call MOM_error(WARNING,"MOM_register_variable_axes_subdomain: the axis "//trim(dim_names(i))//& - "is not included in the valid x and y dimension cases. If the code hangs, check the whether "//& - "an x or y axis is being registered as a non-domain-decomposed variable, "//& - "and add it to the accepted cases if necessary.") - call register_axis(fileObj, trim(dim_names(i)), dimSizes(i)) - end select - ! endif - enddo - - if (allocated(dimSizes)) deallocate(dimSizes) - if (allocated(dim_names)) deallocate(dim_names) -end subroutine MOM_register_variable_axes_subdomain - -!> register axes associated with a variable from a domain-decomposed netCDF file -!> @note The user must specify units for variables with longitude/x-axis and/or latitude/y-axis axes -!! to obtain the correct domain decomposition for the data buffer. -subroutine MOM_register_variable_axes_full(fileObj, variableName, position) - type(FmsNetcdfDomainFile_t), intent(inout) :: fileObj !< netCDF file object returned by call to open_file - character(len=*), intent(in) :: variableName !< name of the variable - integer, optional, intent(in) :: position !< A flag indicating where this data is discretized - - ! Local variables - character(len=40) :: units ! units corresponding to a specific variable dimension - character(len=40), allocatable, dimension(:) :: dim_names ! variable dimension names - integer :: i - integer :: ndims ! number of dimensions - integer :: xPos, yPos ! domain positions for x and y axes. Default is CENTER - integer, allocatable, dimension(:) :: dimSizes ! variable dimension sizes - - if (.not. check_if_open(fileObj)) call MOM_error(FATAL,"MOM_axis:register_variable_axes: The fileObj has "// & - "not been opened. Call fms2_open_file(fileObj,...) before "// & - "passing the fileObj argument to this function.") - xpos = CENTER ; ypos = CENTER - if (present(position)) then - if ((position == CORNER) .or. (position == EAST_FACE)) xpos = EAST_FACE - if ((position == CORNER) .or. (position == NORTH_FACE)) ypos = NORTH_FACE - endif - - ! get variable dimension names and lengths - ndims = get_variable_num_dimensions(fileObj, trim(variableName)) - allocate(dimSizes(ndims)) - allocate(dim_names(ndims)) - call get_variable_size(fileObj, trim(variableName), dimSizes) - call get_variable_dimension_names(fileObj, trim(variableName), dim_names) - ! register the axes - !>@note: This is not a comprehensive check for all possible supported horizontal axes associated with variables - !! read from netCDF files. Developers should add/remove cases as needed. - do i=1,ndims - if (.not.(is_dimension_registered(fileobj, trim(dim_names(i))))) then - select case(trim(lowercase(dim_names(i)))) - case ("grid_x_t") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case ("nx") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case("nxp") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case("longitude") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case("long") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case("lon") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case("lonh") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case("lonq") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case("xh") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case("i") - call register_axis(fileObj, trim(dim_names(i)),"x", domain_position=xPos) - case ("grid_y_t") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case ("ny") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case("nyp") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case("latitude") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case("lat") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case("lath") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case("latq") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case("yh") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case("j") - call register_axis(fileObj, trim(dim_names(i)),"y", domain_position=yPos) - case default ! assumes that the axis is not domain-decomposed - if (.not. is_dimension_unlimited(fileObj, trim(dim_names(i)))) & - call MOM_error(WARNING,"MOM_register_variable_axes_full: the axis "//trim(dim_names(i))//" is not "//& - "included in the valid x and y dimension cases. If the code hangs, check the whether "//& - "an x or y axis is being registered as a non-domain-decomposed variable, "//& - "and add it to the accepted cases if necessary.") - call register_axis(fileObj, trim(dim_names(i)), dimSizes(i)) - end select - endif - enddo - - deallocate(dimSizes) - deallocate(dim_names) -end subroutine MOM_register_variable_axes_full - - -!> convert the variable checksum integer(s) to a single string -!! If there is more than 1 checksum, commas are inserted between -!! each checksum value in the output string -function convert_checksum_to_string(checksum_int) result (checksum_string) - integer(kind=8), intent(in) :: checksum_int !< checksum integer values -! local - character(len=64) :: checksum_string - integer :: i - - checksum_string = '' - - write (checksum_string,'(Z16)') checksum_int ! Z16 is the hexadecimal format code - -end function convert_checksum_to_string - - -end module MOM_axis diff --git a/config_src/infra/FMS2/MOM_write_field_fms2.F90 b/config_src/infra/FMS2/MOM_write_field_fms2.F90 deleted file mode 100644 index 24ba5ebb50..0000000000 --- a/config_src/infra/FMS2/MOM_write_field_fms2.F90 +++ /dev/null @@ -1,1503 +0,0 @@ -!> This module contains wrapper functions to write data to netcdf files -module MOM_write_field_fms2 - -! This file is part of MOM6. See LICENSE.md for the license. - -use MOM_axis, only : MOM_get_diagnostic_axis_data, MOM_register_diagnostic_axis -use MOM_axis, only : axis_data_type, get_time_index, get_var_dimension_metadata -use MOM_axis, only : get_time_units, convert_checksum_to_string -use MOM_coms_infra, only : PE_here, root_PE, num_PEs -use MOM_domain_infra, only : MOM_domain_type -use MOM_domain_infra, only : domain2d, CENTER, CORNER, NORTH_FACE, EAST_FACE -use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING -use MOM_grid, only : ocean_grid_type -use MOM_dyn_horgrid, only : dyn_horgrid_type -use MOM_string_functions, only : lowercase, append_substring -use MOM_verticalGrid, only : verticalGrid_type - -use netcdf, only : nf90_max_name -! fms2_io -use fms2_io_mod, only : check_if_open, get_dimension_size -use fms2_io_mod, only : get_num_dimensions, get_num_variables, get_variable_names -use fms2_io_mod, only : get_unlimited_dimension_name, get_variable_dimension_names -use fms2_io_mod, only : get_variable_num_dimensions, get_variable_size, get_variable_units -use fms2_io_mod, only : get_variable_unlimited_dimension_index, is_dimension_unlimited -use fms2_io_mod, only : is_dimension_registered, register_axis -use fms2_io_mod, only : register_field, register_variable_attribute, fms2_open_file => open_file -use fms2_io_mod, only : fms2_close_file => close_file, write_data, variable_exists -use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, unlimited - -implicit none; private - -public write_field - -! CAUTION: The following variables are saved by default, and are only necessary for consecutive calls to -! write_field with the same file name. The user should ensure that fms2_close_file on -! the fileobj_write_field structures are called at every requisite time step at after the last -! variable is written to the file by omitting the optional leave_file_open argument, or setting it to .false. - -!> netCDF non-domain-decomposed file object returned by call to open_file in write_field calls -type(FmsNetcdfFile_t), private :: fileobj_write_field - -!> netCDF domain-decomposed file object returned by call to open_file in write_field calls -type(FmsNetcdfDomainFile_t), private :: fileobj_write_field_dd - -!> index of the time_level value that is written to netCDF file by the write_field routines -integer, private :: write_field_time_index - -!> interface to write data to a netcdf file generated by create_file -interface write_field - module procedure write_field_4d_DD - module procedure write_field_3d_DD - module procedure write_field_2d_DD - module procedure write_field_1d_DD - module procedure write_scalar - module procedure write_field_4d_noDD - module procedure write_field_3d_noDD - module procedure write_field_2d_noDD - module procedure write_field_1d_noDD -end interface - -contains -!> This function uses the fms_io function write_data to write a 1-D domain-decomposed data field named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_field_1d_DD(filename, fieldname, data, mode, domain, hor_grid, z_grid, & - start_index, edge_lengths, time_level, time_units, & - checksums, G, dG, GV, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, target, dimension(:), intent(in) :: data !< The 1-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - type(MOM_domain_type),intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition - character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor - character(len=*), intent(in) :: z_grid !< vertical grid descriptor - integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is - !! the variable size - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day. - integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any - !! vertical grid axes. - logical, optional, intent(in) :: leave_file_open !< if .true., leave the file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success !.true. if call to open_file is successful - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - integer :: num_dims, substring_index - integer :: dim_unlim_size! size of the unlimited dimension - integer, dimension(1) :: start, nwrite ! indices for first data value and number of values to write - character(len=20) :: t_units ! time units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=64) :: checksum_char ! checksum character array created from checksum argument - character(len=1024) :: filename_temp - character(len=48), dimension(2) :: dim_names !< variable dimension names (or name, in the 1-D case); 1 extra - !! dimension in case appending along the time axis - integer, dimension(2) :: dim_lengths !< variable dimension lengths (or length, in the 1-D case) - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - dim_unlim_size=0 - dim_unlim_name="" - dim_names(:) = "" - dim_lengths(:) = 0 - num_dims = 0 - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - ! get the dimension names and lengths - ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified - ! and not assumed from the t_grid value - if (present(G)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, G=G) - elseif(present(dG)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, dG=dG) - endif - - if (present(GV)) & - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, GV=GV) - ! define the start and edge_length arguments - start(:) = 1 - nwrite(:) = dim_lengths(1) - if (present(start_index)) then - start(1) = max(1, start_index(1)) - endif - - if (present(edge_lengths)) then - nwrite(1) = max(dim_lengths(1),edge_lengths(1)) - endif - - if (.not.(check_if_open(fileobj_write_field_dd))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_write_field_fms2:write_1d_DD:mode argument must be write, overwrite, or append") - ! get the time_level index - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field_dd, trim(filename_temp), lowercase(trim(mode)), & - domain%mpp_domain, is_restart=.false.) - ! register the diagnostic axis associated with the variable - call MOM_register_diagnostic_axis(fileobj_write_field_dd, trim(dim_names(1)), dim_lengths(1)) - endif - ! register and write the time_level - if (present(time_level)) then - if (.not. (variable_exists(fileobj_write_field_dd, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field_dd, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field_dd, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) - else - ! write the time_level if it is larger than the most recent file time - if (write_field_time_index .gt. dim_unlim_size) & - call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), & - corner=(/write_field_time_index/), edge_lengths=(/1/)) - endif - endif - ! register the field if it is not already in the file - if (.not.(variable_exists(fileobj_write_field_dd, trim(fieldname)))) then - call register_field(fileobj_write_field_dd, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) - if (present(units)) & - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'long_name', trim(longname)) - ! write the checksum attribute - if (present(checksums)) then - ! convert the checksum to a string - checksum_char = "" - checksum_char = convert_checksum_to_string(checksums(1,1)) - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), "checksum", checksum_char) - endif - endif - ! write the variable - if (present(time_level)) then - call write_data(fileobj_write_field_dd, trim(fieldname), data, corner=start, edge_lengths=nwrite, & - unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field_dd, trim(fieldname), data, corner=start, edge_lengths=nwrite) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field_dd)) call fms2_close_file(fileobj_write_field_dd) - write_field_time_index = 0 - endif - -end subroutine write_field_1d_DD - -!> This function uses the fms_io function write_data to write a 2-D domain-decomposed data field named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_field_2d_DD(filename, fieldname, data, mode, domain, hor_grid, z_grid, & - start_index, edge_lengths, time_level, time_units, & - checksums, G, dG, GV, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, target, dimension(:,:), intent(in) :: data !< The 1-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - type(MOM_domain_type),intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition - character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor - character(len=*), intent(in) :: z_grid !< vertical grid descriptor - integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is - !! the variable size - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day. - integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any - !! vertical grid axes. - logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success !.true. if call to open_file is successful - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - integer :: i, is, ie, js, je, j, ndims, num_dims, substring_index - integer, allocatable, dimension(:) :: x_inds, y_inds - integer :: dim_unlim_size ! size of the unlimited dimension - integer :: file_dim_length - integer, dimension(2) :: start, nwrite ! indices for starting points and number of values to write - character(len=20) :: t_units ! time units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=1024) :: filename_temp - character(len=64) :: checksum_char ! checksum character array created from checksum argument - character(len=48), dimension(3) :: dim_names ! variable dimension names; 1 extra dimension in case appending - ! along the time axis - character(len=48), allocatable, dimension(:) :: file_dim_names - integer, dimension(3) :: dim_lengths ! variable dimension lengths - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - dim_lengths(:) = 0 - dim_names(:) = "" - dim_unlim_size = 0 - dim_unlim_name = "" - ndims = 2 - num_dims = 0 - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - ! get the dimension names and lengths - ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension - ! is user-specified rather than derived from the t_grid value - if (present(G)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, G=G) - elseif(present(dG)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, dG=dG) - endif - if (present(GV)) & - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, GV=GV) - ! set the start (start_index) and nwrite (edge_lengths) values - start(:) = 1 - nwrite(:) = dim_lengths(1:ndims) - - if (present(start_index)) then - do i=1,ndims - start(i) = max(1,start_index(i)) - enddo - endif - - if (present(edge_lengths)) then - do i=1,ndims - nwrite(i) = max(dim_lengths(i),edge_lengths(i)) - enddo - endif - - if (.not.(check_if_open(fileobj_write_field_dd))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_write_field_fms2:write_2d_DD:mode argument must be write, overwrite, or append") - ! get the time_level index - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field_dd, trim(filename_temp), lowercase(trim(mode)), & - domain%mpp_domain, is_restart=.false.) - endif - ! register the horizontal diagnostic axes associated with the variable - do i=1,num_dims - if (.not.(is_dimension_registered(fileobj_write_field_dd, trim(dim_names(i))))) & - call MOM_register_diagnostic_axis(fileobj_write_field_dd, trim(dim_names(i)), dim_lengths(i)) - enddo - ! register and write the time_level - if (present(time_level)) then - call get_unlimited_dimension_name(fileobj_write_field_dd, dim_unlim_name) - call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) - num_dims=num_dims+1 - dim_names(num_dims) = trim(dim_unlim_name) - - if (.not. (variable_exists(fileobj_write_field_dd, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field_dd, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field_dd, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) - else - ! write the time_level if it is larger than the most recent file time - if (write_field_time_index .gt. dim_unlim_size) & - call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), & - corner=(/write_field_time_index/), edge_lengths=(/1/)) - endif - endif - ! register the variable if it is not already in the file - if (.not.(variable_exists(fileobj_write_field_dd, trim(fieldname)))) then - call register_field(fileobj_write_field_dd, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) - if (present(units)) & - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'long_name', trim(longname)) - ! write the checksum attribute - if (present(checksums)) then - ! convert the checksum to a string - checksum_char = "" - checksum_char = convert_checksum_to_string(checksums(1,1)) - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), "checksum", checksum_char) - endif - endif - ! write the variable - if (present(time_level)) then - call write_data(fileobj_write_field_dd, trim(fieldname), data, corner=start, edge_lengths=nwrite, & - unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field_dd, trim(fieldname), data, corner=start, edge_lengths=nwrite) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field_dd)) call fms2_close_file(fileobj_write_field_dd) - write_field_time_index=0 - if (allocated(file_dim_names)) deallocate(file_dim_names) - endif - -end subroutine write_field_2d_DD - -!> This function uses the fms_io function write_data to write a 3-D domain-decomposed data field named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_field_3d_DD(filename, fieldname, data, mode, domain, hor_grid, z_grid, & - start_index, edge_lengths, time_level, time_units, & - checksums, G, dG, GV, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, target, dimension(:,:,:), intent(in) :: data !< The 1-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - type(MOM_domain_type),intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition - character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor - character(len=*), intent(in) :: z_grid !< vertical grid descriptor - integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is - !! the variable size - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day. - integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any - !! vertical grid axes. - logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success !.true. if call to open_file is successful - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - integer :: i, is, ie, js, je, ndims, num_dims, substring_index - integer :: dim_unlim_size ! size of the unlimited dimension - integer, dimension(3) :: start, nwrite ! indices for first data value and number of values to write - character(len=20) :: t_units ! time units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=1024) :: filename_temp - character(len=64) :: checksum_char ! checksum character array created from checksum argument - character(len=48), dimension(4) :: dim_names !< variable dimension names; 1 extra dimension in case appending - !! along the time axis - integer, dimension(4) :: dim_lengths !< variable dimension lengths - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - dim_unlim_size = 0 - dim_unlim_name = "" - dim_names(:) = "" - dim_lengths(:) = 0 - num_dims = 0 - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - ! get the dimension names and lengths - ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified - ! and not assumed from the t_grid value - if (present(G)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, G=G) - elseif(present(dG)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, dG=dG) - endif - - if (present(GV)) & - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, GV=GV) - ! set the start (start_index) and nwrite (edge_lengths) values - ndims = 3 - start(:) = 1 - nwrite(:) = dim_lengths(1:3) - if (present(start_index)) then - do i=1,ndims - start(i) = max(1,start_index(i)) - enddo - endif - - if (present(edge_lengths)) then - do i=1,ndims - nwrite(i) = max(dim_lengths(i), edge_lengths(i)) - enddo - endif - - ! open the file - if (.not.(check_if_open(fileobj_write_field_dd))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_write_field_fms2:write_3d_DD:mode argument must be write, overwrite, or append") - ! get the time_level index - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field_dd, trim(filename_temp), lowercase(trim(mode)), & - domain%mpp_domain, is_restart=.false.) - ! register the horizontal and vertical diagnostic axes associated with the variable - do i=1,ndims - call MOM_register_diagnostic_axis(fileobj_write_field_dd, trim(dim_names(i)), dim_lengths(i)) - enddo - endif - ! register and write the time_level - if (present(time_level)) then - call get_unlimited_dimension_name(fileobj_write_field_dd ,dim_unlim_name) - call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) - num_dims=num_dims+1 - dim_names(num_dims) = trim(dim_unlim_name) - - if (.not. (variable_exists(fileobj_write_field_dd, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field_dd, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field_dd, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) - else - ! write the time_level if it is larger than the most recent file time - if (write_field_time_index .gt. dim_unlim_size ) & - call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), & - corner=(/write_field_time_index/), edge_lengths=(/1/)) - endif - endif - ! register the field if it is not already in the file - if (.not.(variable_exists(fileobj_write_field_dd, trim(fieldname)))) then - call register_field(fileobj_write_field_dd, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) - if (present(units)) & - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'long_name', trim(longname)) - ! write the checksum attribute - if (present(checksums)) then - ! convert the checksum to a string - checksum_char = "" - checksum_char = convert_checksum_to_string(checksums(1,1)) - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), "checksum", checksum_char) - endif - endif - ! write the data - if (present(time_level)) then - call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) - call write_data(fileobj_write_field_dd, trim(fieldname), data, corner=start, edge_lengths=nwrite, & - unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field_dd, trim(fieldname), data, corner=start, edge_lengths=nwrite) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field_dd)) call fms2_close_file(fileobj_write_field_dd) - write_field_time_index=0 - endif - - -end subroutine write_field_3d_DD - -!> This function uses the fms_io function write_data to write a 4-D domain-decomposed data field named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_field_4d_DD(filename, fieldname, data, mode, domain, hor_grid, z_grid, t_grid, & - start_index, edge_lengths, time_level, time_units, & - checksums, G, dG, GV, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, target, dimension(:,:,:,:), intent(in) :: data !< The 1-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - type(MOM_domain_type),intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition - character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor - character(len=*), intent(in) :: z_grid !< vertical grid descriptor - character(len=*), intent(in) :: t_grid !< time descriptor - integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is - !! the variable size - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day. - integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any - !! vertical grid axes. - logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success !.true. if call to open_file is successful - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - real :: file_time ! most recent time currently written to file - integer :: i, ndims, num_dims, substring_index - integer :: dim_unlim_size ! size of the unlimited dimension - integer, dimension(4) :: start, nwrite ! indices for first data value and number of values to write - character(len=20) :: t_units ! time units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=1024) :: filename_temp - character(len=64) :: checksum_char ! checksum character array created from checksum argument - character(len=48), dimension(4) :: dim_names ! variable dimension names - integer, dimension(4) :: dim_lengths ! variable dimension lengths - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - num_dims = 0 - dim_unlim_size = 0 - dim_unlim_name = "" - dim_names(:) = "" - dim_lengths(:) = 0 - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - ! get the dimension names and lengths - if (present(G)) then - call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & - dim_lengths, num_dims, G=G) - elseif(present(dG)) then - call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & - dim_lengths, num_dims, dG=dG) - endif - - if (present(GV)) & - call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & - dim_lengths, num_dims, GV=GV) - ! set the start (start_index) and nwrite (edge_lengths) values - ndims = 4 - start(:) = 1 - nwrite(:) = dim_lengths(:) - if (present(start_index)) then - do i=1,ndims - start(i) = max(1,start_index(i)) - enddo - endif - - if (present(edge_lengths)) then - do i=1,ndims - nwrite(i) = max(dim_lengths(i), edge_lengths(i)) - enddo - endif - - ! open the file - if (.not.(check_if_open(fileobj_write_field_dd))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_write_field_fms2:write_4d_DD:mode argument must be write, overwrite, or append") - ! get the index of the corresponding time_level the first time the file is opened - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field_dd, trim(filename_temp), lowercase(trim(mode)), & - domain%mpp_domain, is_restart=.false.) - ! register the horizontal and vertical diagnostic axes associated with the variable - do i=1,ndims - call MOM_register_diagnostic_axis(fileobj_write_field_dd, trim(dim_names(i)), dim_lengths(i)) - enddo - endif - ! register the time dimension and write the time_level - if (present(time_level)) then - call get_unlimited_dimension_name(fileobj_write_field_dd, dim_unlim_name) - call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) - num_dims=num_dims+1 - dim_names(num_dims) = trim(dim_unlim_name) - - if (.not. (variable_exists(fileobj_write_field_dd, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field_dd, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field_dd, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) - else - ! write the time_level if it is larger than the most recent file time - if (write_field_time_index .gt. dim_unlim_size) & - call write_data(fileobj_write_field_dd, trim(dim_unlim_name), (/time_level/), & - corner=(/write_field_time_index/), edge_lengths=(/1/)) - endif - endif - ! register the variable if it is not already in the file - if (.not.(variable_exists(fileobj_write_field_dd, trim(fieldname)))) then - call register_field(fileobj_write_field_dd, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) - if (present(units)) & - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), 'long_name', trim(longname)) - ! write the checksum attribute - if (present(checksums)) then - ! convert the checksum to a string - checksum_char = "" - checksum_char = convert_checksum_to_string(checksums(1,1)) - call register_variable_attribute(fileobj_write_field_dd, trim(fieldname), "checksum", checksum_char) - endif - endif - ! write the data - if (present(time_level)) then - call get_dimension_size(fileobj_write_field_dd, trim(dim_unlim_name), dim_unlim_size) - call write_data(fileobj_write_field_dd, trim(fieldname), data, corner=start, edge_lengths=nwrite, & - unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field_dd, trim(fieldname), data, corner=start, edge_lengths=nwrite) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field_dd)) call fms2_close_file(fileobj_write_field_dd) - write_field_time_index=0 - endif - -end subroutine write_field_4d_DD - -!> This routine uses the fms_io function write_data to write a scalar variable named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_scalar(filename, fieldname, data, mode, time_level, time_units, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, intent(in) :: data !< The 1-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day - logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success !.true. if call to open_file is successful - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - character(len=20) :: t_units ! time units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=1024) :: filename_temp - character(len=48), dimension(1) :: dim_names ! variable dimension names - integer :: i, num_dims, substring_index - integer :: dim_unlim_size ! size of the unlimited dimension - real, allocatable, dimension(:) :: file_times - integer, dimension(1) :: dim_lengths ! variable dimension lengths - integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file - - dim_unlim_size = 0 - dim_unlim_name= "" - dim_names(:) = "" - dim_lengths(:) = 0 - num_dims = 0 - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - - if (.not.(check_if_open(fileobj_write_field))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_write_field_fms2:write_scaler:mode argument must be write, overwrite, or append") - ! get the index of the corresponding time_level the first time the file is opened - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! get the pes associated with the file. - !>\note this is required so that only pe(1) is identified as the root pe to create the file - !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure - if (.not.(allocated(pelist))) then - allocate(pelist(num_PEs())) - pelist(:) = 0 - do i=1,size(pelist) - pelist(i) = i-1 - enddo - endif - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), trim(mode), is_restart=.false., & - pelist=pelist) - endif - if (present(time_level)) then - call get_unlimited_dimension_name(fileobj_write_field, dim_unlim_name) - call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) - ! write the time value if it is not already written to the file - if (.not.(variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/)) - else - ! write the next time value if it is larger than the most recent file time - if (write_field_time_index .gt. dim_unlim_size) & - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & - edge_lengths=(/1/)) - endif - endif - ! register the variable - if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then - if (present(time_level)) then - call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=(/trim(dim_unlim_name)/)) - else - call register_field(fileobj_write_field, trim(fieldname), "double") - endif - if (present(units)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) - endif - ! write the data - if (present(time_level)) then - call write_data(fileobj_write_field, trim(fieldname), data, unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field, trim(fieldname), data) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) - if (allocated(pelist)) deallocate(pelist) - write_field_time_index=0 - endif -end subroutine write_scalar - -!> This function uses the fms_io function write_data to write a 1-D non-domain-decomposed data field named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_field_1d_noDD(filename, fieldname, data, mode, hor_grid, z_grid, & - start_index, edge_lengths, time_level, time_units, & - checksums, G, dG, GV, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, target, dimension(:), intent(in) :: data !< The 1-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor - character(len=*), intent(in) :: z_grid !< vertical grid descriptor - integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is the - !! variable size - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day. - integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any - !! vertical grid axes. - logical, optional, intent(in) :: leave_file_open !< if .true., leave file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success !.true. if call to open_file is successful - integer :: i, ndims, num_dims, substring_index - integer :: dim_unlim_size ! size of the unlimited dimension - integer, dimension(1) :: start, nwrite ! indices for first data value and number of values to write - character(len=20) :: t_units ! time units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=1024) :: filename_temp - character(len=64) :: checksum_char ! checksum character array created from checksum argument - character(len=48), dimension(2) :: dim_names ! variable dimension names (up to 2 if appended at time level) - integer, dimension(2) :: dim_lengths ! variable dimension lengths - integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - dim_unlim_size = 0 - dim_unlim_name= "Time" - dim_names(:) = "" - dim_lengths(:) = 0 - num_dims = 0 - - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - ! get the dimension names and lengths - ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified - ! and not assumed from the t_grid value. - if (present(G)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, G=G) - elseif(present(dG)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, dG=dG) - endif - - if (present(GV)) & - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, GV=GV) - ! set the start (start_index) and nwrite (edge_lengths) values - start(:) = 1 - nwrite(:) = dim_lengths(1) - if (present(start_index)) then - start(1) = max(1,start_index(1)) - endif - - if (present(edge_lengths)) then - nwrite(1) = max(dim_lengths(1),edge_lengths(1)) - endif - - if (.not.(check_if_open(fileobj_write_field))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_write_field_fms2:write_1d_noDD:mode argument must be write, overwrite, or append") - ! get the index of the corresponding time_level the first time the file is opened - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! get the pes associated with the file. - !>\note this is required so that only pe(1) is identified as the root pe to create the file - !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure - if (.not.(allocated(pelist))) then - allocate(pelist(num_PEs())) - pelist(:) = 0 - do i=1,size(pelist) - pelist(i) = i-1 - enddo - endif - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), lowercase(trim(mode)), & - is_restart=.false., pelist=pelist) - endif - ! write the data, and the time value if it is not already written to the file - if (present(time_level)) then - call get_unlimited_dimension_name(fileobj_write_field,dim_unlim_name) - call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) - num_dims=num_dims+1 - dim_names(num_dims) = trim(dim_unlim_name) - - if (.not. (variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) - else - ! write the time value if it is larger than the most recent file time - if (write_field_time_index .gt. dim_unlim_size) & - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & - edge_lengths=(/1/)) - endif - endif - ! register the field if it is not already in the file - if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then - call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) - if (present(units)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) - ! write the checksum attribute - if (present(checksums)) then - ! convert the checksum to a string - checksum_char = '' - checksum_char = convert_checksum_to_string(checksums(1,1)) - call register_variable_attribute(fileobj_write_field, trim(fieldname), "checksum", checksum_char) - endif - endif - ! write the variable to the file - if (present(time_level)) then - call write_data(fileobj_write_field, trim(fieldname), data, corner=start, edge_lengths=nwrite, & - unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field, trim(fieldname), data, corner=start, edge_lengths=nwrite) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) - if (allocated(pelist)) deallocate(pelist) - write_field_time_index = 0 - endif - - -end subroutine write_field_1d_noDD - -!> This function uses the fms_io function write_data to write a scalar variable named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_field_2d_noDD(filename, fieldname, data, mode, hor_grid, z_grid, & - start_index, edge_lengths, time_level, time_units, & - checksums, G, dG, GV, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, target, dimension(:,:), intent(in) :: data !< The 2-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor - character(len=*), intent(in) :: z_grid !< vertical grid descriptor - integer, dimension(2), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(2), optional, intent(in) :: edge_lengths !< number of data values to read in; default is the - !! variable size - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day. - integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any - !! vertical grid axes. - logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success ! .true. if call to open_file is successful - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - integer :: i, ndims, num_dims, substring_index - integer :: dim_unlim_size ! size of the unlimited dimension - integer, dimension(2) :: start, nwrite ! indices for starting points and number of values to write - character(len=20) :: t_units ! time units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=1024) :: filename_temp - character(len=64) :: checksum_char ! checksum character array created from checksum argument - character(len=48), dimension(3) :: dim_names ! variable dimension names - integer, dimension(3) :: dim_lengths ! variable dimension lengths - integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - dim_unlim_size = 0 - dim_unlim_name = "" - dim_names(:) = "" - dim_lengths(:) = 0 - num_dims = 0 - - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - - ! get the dimension names and lengths - ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified - ! and not assumed from the t_grid value - if (present(G)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, G=G) - elseif(present(dG)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, dG=dG) - endif - - if (present(GV)) & - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, GV=GV) - - ! set the start (start_index) and nwrite (edge_lengths) values - ndims=2 - start(:) = 1 - nwrite(:) = dim_lengths(1:2) - if (present(start_index)) then - do i=1,ndims - start(i) = max(1,start_index(i)) - enddo - endif - - if (present(edge_lengths)) then - do i=1,ndims - nwrite(i) = max(dim_lengths(i),edge_lengths(i)) - enddo - endif - - if (.not.(check_if_open(fileobj_write_field))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_write_field_fms2:write_2d_noDD:mode argument must be write, overwrite, or append") - ! get the time_level index - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! get the pes associated with the file. - !>\note this is required so that only pe(1) is identified as the root pe to create the file - !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure - if(.not.(allocated(pelist))) then - allocate(pelist(num_PEs())) - pelist(:) = 0 - do i=1,size(pelist) - pelist(i) = i-1 - enddo - endif - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), lowercase(trim(mode)), & - is_restart=.false., pelist=pelist) - endif - - if (present(time_level)) then - call get_unlimited_dimension_name(fileobj_write_field,dim_unlim_name) - call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) - num_dims=num_dims+1 - dim_names(num_dims) = trim(dim_unlim_name) - - if (.not. (variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) - else - ! write the time value if it is larger than the most recent file time - if (write_field_time_index .gt. dim_unlim_size) & - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & - edge_lengths=(/1/)) - endif - endif - - ! register the variable to the file - if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then - call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) - if (present(units)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) - ! write the checksum attribute - if (present(checksums)) then - ! convert the checksum to a string - checksum_char = "" - checksum_char = convert_checksum_to_string(checksums(1,1)) - call register_variable_attribute(fileobj_write_field, trim(fieldname), "checksum", checksum_char) - endif - endif - ! write the variable to the file - if (present(time_level)) then - call write_data(fileobj_write_field, trim(fieldname), data, corner=start, edge_lengths=nwrite, & - unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field, trim(fieldname), data, corner=start, edge_lengths=nwrite) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) - if (allocated(pelist)) deallocate(pelist) - write_field_time_index=0 - endif - -end subroutine write_field_2d_noDD - -!> This function uses the fms_io function write_data to write a 3-D non-domain-decomposed data field named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_field_3d_noDD(filename, fieldname, data, mode, hor_grid, z_grid, & - start_index, edge_lengths, time_level, time_units, & - checksums, G, dG, GV, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, target, dimension(:,:,:), intent(in) :: data !< The 3-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor - character(len=*), intent(in) :: z_grid !< vertical grid descriptor - integer, dimension(4), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(4), optional, intent(in) :: edge_lengths !< number of data values to read in; default is the - !! variable size - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day. - integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any - !! vertical grid axes. - logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success !.true. if call to open_file is successful - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - integer :: i, ndims, num_dims, substring_index - integer :: dim_unlim_size ! size of the unlimited dimension - integer, dimension(3) :: start, nwrite ! indices for first data value and number of values to write - character(len=20) :: t_units ! time_units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=1024) :: filename_temp - character(len=64) :: checksum_char ! checksum character array created from checksum argument - character(len=48), dimension(4) :: dim_names ! variable dimension names - integer, dimension(4) :: dim_lengths ! variable dimension lengths - integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - dim_unlim_size = 0 - dim_unlim_name = "" - dim_names(:) = "" - dim_lengths(:) = 0 - num_dims = 0 - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - ! get the dimension names and lengths - ! NOTE: the t_grid argument is set to '1' (do nothing) because the presence of a time dimension is user-specified - ! and not assumed from the t_grid value - if (present(G)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, G=G) - elseif(present(dG)) then - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, dG=dG) - endif - - if (present(GV)) & - call get_var_dimension_metadata(hor_grid, z_grid, '1', dim_names, & - dim_lengths, num_dims, GV=GV) - ! set the start (start_index) and nwrite (edge_lengths) values - ndims = 3 - start(:) = 1 - nwrite(:) = dim_lengths(1:3) - if (present(start_index)) then - do i=1,ndims - start(i) = max(1,start_index(i)) - enddo - endif - - if (present(edge_lengths)) then - do i=1,ndims - nwrite(i) = max(dim_lengths(i), edge_lengths(i)) - enddo - endif - - ! open the file - if (.not.(check_if_open(fileobj_write_field))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_io:write_3d_noDD:mode argument must be write, overwrite, or append") - ! get the time_level index - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! get the pes associated with the file. - !>\note this is required so that only pe(1) is identified as the root pe to create the file - !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure - if (.not.(allocated(pelist))) then - allocate(pelist(num_PEs())) - pelist(:) = 0 - do i=1,size(pelist) - pelist(i) = i-1 - enddo - endif - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), lowercase(trim(mode)), & - is_restart=.false., pelist=pelist) - endif - ! register and write the time_level - if (present(time_level)) then - call get_unlimited_dimension_name(fileobj_write_field,dim_unlim_name) - call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) - num_dims=num_dims+1 - dim_names(num_dims) = trim(dim_unlim_name) - - if (.not. (variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) - else - ! write the time_level if it is larger than the most recent file time - if (write_field_time_index .gt. dim_unlim_size) & - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & - edge_lengths=(/1/)) - endif - endif - ! register the field if it is not already in the file - if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then - call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) - if (present(units)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) - ! write the checksum attribute - if (present(checksums)) then - ! convert the checksum to a string - checksum_char = "" - checksum_char = convert_checksum_to_string(checksums(1,1)) - call register_variable_attribute(fileobj_write_field, trim(fieldname), "checksum", checksum_char) - endif - endif - - if (present(time_level)) then - call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) - call write_data(fileobj_write_field, trim(fieldname), data, corner=start, edge_lengths=nwrite, & - unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field, trim(fieldname), data, corner=start, edge_lengths=nwrite) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) - if (allocated(pelist)) deallocate(pelist) - write_field_time_index=0 - endif - -end subroutine write_field_3d_noDD - -!> This function uses the fms_io function write_data to write a 4-D non-domain-decomposed data field named "fieldname" -!! to the file "filename" in "write", "overwrite", or "append" mode. It should be called after create_file in the MOM -!! file write procedure. -subroutine write_field_4d_noDD(filename, fieldname, data, mode, hor_grid, z_grid, t_grid, & - start_index, edge_lengths, time_level, time_units, & - checksums, G, dG, GV, leave_file_open, units, longname) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, target, dimension(:,:,:,:), intent(in) :: data !< The 1-dimensional data array to pass to read_data - character(len=*), intent(in) :: mode !< "write", "overwrite", or "append" - character(len=*), intent(in) :: hor_grid !< horizontal grid descriptor - character(len=*), intent(in) :: z_grid !< vertical grid descriptor - character(len=*), intent(in) :: t_grid !< time descriptor - integer, dimension(4), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(4), optional, intent(in) :: edge_lengths !< number of data values to read in; default is the - !! variable size - real, optional, intent(in) :: time_level !< time value to write - real, optional, intent(in) :: time_units !< length of the units for time [s]. The - !! default value is 86400.0, for 1 day. - integer(kind=8), dimension(:,:), optional, intent(in) :: checksums !< variable checksum - type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG - !! is required if the new file uses any - !! horizontal grid axes. - type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any - !! vertical grid axes. - logical, optional, intent(in) :: leave_file_open !< flag indicating whether to leave the file open - character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: longname !< long name variable attribute - ! local - logical :: file_open_success !.true. if call to open_file is successful - logical :: close_the_file ! indicates whether to close the file after write_data is called; default is .true. - integer :: i, ndims, num_dims, substring_index - integer :: dim_unlim_size ! size of the unlimited dimension - integer, dimension(4) :: start, nwrite ! indices for first data value and number of values to write - character(len=20) :: t_units ! time units - character(len=nf90_max_name) :: dim_unlim_name ! name of the unlimited dimension in the file - character(len=1024) :: filename_temp - character(len=64) :: checksum_char ! checksum character array created from checksum argument - character(len=48), dimension(4) :: dim_names ! variable dimension names - integer, dimension(4) :: dim_lengths ! variable dimension lengths - integer, allocatable, dimension(:) :: pelist ! list of pes associated with the netCDF file - - close_the_file = .true. - if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - - dim_unlim_size = 0 - dim_unlim_name = "" - dim_names(:) = "" - dim_lengths(:) = 0 - ndims = 4 - num_dims = 0 - ! append '.nc' to the file name if it is missing - filename_temp = "" - substring_index = 0 - substring_index = index(trim(filename), ".nc") - if (substring_index <= 0) then - filename_temp = append_substring(filename,".nc") - else - filename_temp = filename - endif - - ! get the dimension names and lengths - if (present(G)) then - call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & - dim_lengths, num_dims, G=G) - elseif(present(dG)) then - call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & - dim_lengths, num_dims, dG=dG) - endif - if (present(GV)) & - call get_var_dimension_metadata(hor_grid, z_grid, t_grid, dim_names, & - dim_lengths, num_dims, GV=GV) - ! set the start (start_index) and nwrite (edge_lengths) values - start(:) = 1 - nwrite(:) = dim_lengths(:) - if (present(start_index)) then - do i=1,ndims - start(i) = max(1, start_index(i)) - enddo - endif - - if (present(edge_lengths)) then - do i=1,ndims - nwrite(i) = max(dim_lengths(i), edge_lengths(i)) - enddo - endif - - ! open the file - if (.not.(check_if_open(fileobj_write_field))) then - if ((lowercase(trim(mode)) .ne. "write") .and. (lowercase(trim(mode)) .ne. "append") .and. & - (lowercase(trim(mode)) .ne. "overwrite")) & - call MOM_error(FATAL,"MOM_write_field_fms2:write_4d_noDD:mode argument must be write, overwrite, or append") - ! get the time_level index - if (present(time_level)) write_field_time_index = get_time_index(trim(filename_temp), time_level) - ! get the pes associated with the file. - !>\note this is required so that only pe(1) is identified as the root pe to create the file - !! Otherwise, multiple pes may try to open the file in write (NC_NOCLOBBER) mode, leading to failure - if (.not.(allocated(pelist))) then - allocate(pelist(num_PEs())) - pelist(:) = 0 - do i=1,size(pelist) - pelist(i) = i-1 - enddo - endif - ! open the file in write or append mode - file_open_success = fms2_open_file(fileobj_write_field, trim(filename_temp), lowercase(trim(mode)), & - is_restart=.false., pelist=pelist) - endif - ! register and write the time_level - if (present(time_level)) then - call get_unlimited_dimension_name(fileobj_write_field,dim_unlim_name) - call get_dimension_size(fileobj_write_field, trim(dim_unlim_name), dim_unlim_size) - num_dims=num_dims+1 - dim_names(num_dims) = trim(dim_unlim_name) - ! write the time value if it is not already written to the file - if (.not. (variable_exists(fileobj_write_field, trim(dim_unlim_name)))) then - ! set the time units - t_units = "" - if (present(time_units)) then - t_units = get_time_units(time_units) - else - t_units = "days" - endif - - call register_field(fileobj_write_field, trim(dim_unlim_name), "double", dimensions=(/trim(dim_unlim_name)/)) - call register_variable_attribute(fileobj_write_field, trim(dim_unlim_name), 'units', trim(t_units)) - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/)) - else - if (write_field_time_index .gt. dim_unlim_size) & - call write_data(fileobj_write_field, trim(dim_unlim_name), (/time_level/), corner=(/write_field_time_index/), & - edge_lengths=(/1/)) - endif - endif - ! register the variable - if (.not.(variable_exists(fileobj_write_field, trim(fieldname)))) then - call register_field(fileobj_write_field, trim(fieldname), "double", dimensions=dim_names(1:num_dims)) - if (present(units)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'units', trim(units)) - if (present(longname)) & - call register_variable_attribute(fileobj_write_field, trim(fieldname), 'long_name', trim(longname)) - ! write the checksum attribute - if (present(checksums)) then - ! convert the checksum to a string - checksum_char = "" - checksum_char = convert_checksum_to_string(checksums(1,1)) - call register_variable_attribute(fileobj_write_field, trim(fieldname), "checksum", checksum_char) - endif - endif - ! write the variable to the file - if (present(time_level)) then - call write_data(fileobj_write_field, trim(fieldname), data, corner=start, edge_lengths=nwrite, & - unlim_dim_level=write_field_time_index) - else - call write_data(fileobj_write_field, trim(fieldname), data, corner=start, edge_lengths=nwrite) - endif - ! close the file - if (close_the_file) then - if (check_if_open(fileobj_write_field)) call fms2_close_file(fileobj_write_field) - deallocate(pelist) - write_field_time_index=0 - endif -end subroutine write_field_4d_nodd - -end module MOM_write_field_fms2 From 84d5e217653047eab320fce54a02ad6afcfc260a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 19 Mar 2021 12:43:03 -0400 Subject: [PATCH 06/10] +FMS2 reads for ints & MOM_read_data_fms2 cleanup Added variants of MOM_read_data_0d, MOM_read_data_0d_int, MOM_read_data_1d, MOM_read_data_1d_int, and MOM_read_data_2d_region that use the FMS2 interfaces to read data. Also altered prepare_to_read_var so that it does not open a file (which had been an option before); one argument to prepare_to_read_var was no longer needed and so it was removed. Also added a public interface to find_varname_in_file that does the same as prepare_to_read_var but works on the other FMS2 file type without domain decomposition. Several unused and now redundant routines were removed from MOM_read_data_fms2.F90. Comments describing a number of variables were also added. All of these changes are confined to config_src/infra/FMS2, and all answers are bitwise identical. --- config_src/infra/FMS2/MOM_io_infra.F90 | 227 ++++++++++--- config_src/infra/FMS2/MOM_read_data_fms2.F90 | 335 +++++-------------- 2 files changed, 254 insertions(+), 308 deletions(-) diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index d4fe0b5387..c545462ad8 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -6,8 +6,8 @@ module MOM_io_infra use MOM_domain_infra, only : MOM_domain_type, rescale_comp_data, AGRID, BGRID_NE, CGRID_NE use MOM_domain_infra, only : domain2d, domain1d, CENTER, CORNER, NORTH_FACE, EAST_FACE use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, FATAL, WARNING +use MOM_read_data_fms2, only : prepare_to_read_var, find_varname_in_file -use MOM_read_data_fms2, only : prepare_to_read_var use fms2_io_mod, only : fms2_open_file => open_file, check_if_open, fms2_close_file => close_file use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, fms2_read_data => read_data use fms2_io_mod, only : get_unlimited_dimension_name, get_num_dimensions, get_num_variables @@ -41,9 +41,9 @@ module MOM_io_infra public :: MOM_read_data, MOM_read_vector, write_metadata, write_field public :: field_exists, get_field_atts, get_field_size, get_axis_data, read_field_chksum public :: io_infra_init, io_infra_end, MOM_namelist_file, check_namelist_error, write_version -! These types are inherited from underlying infrastructure code, to act as containers for -! information about fields and axes, respectively, and are opaque to this module. -! public :: file_type, fieldtype, axistype +! These types act as containers for information about files, fields and axes, respectively, +! and may also wrap opaque types from the underlying infrastructure. +public :: file_type, fieldtype, axistype ! These are encoding constant parmeters. public :: ASCII_FILE, NETCDF_FILE, SINGLE_FILE, MULTIPLE public :: APPEND_FILE, READONLY_FILE, OVERWRITE_FILE, WRITEONLY_FILE @@ -102,7 +102,7 @@ module MOM_io_infra end interface flush_file !> Type for holding a handle to an open file and related information -type, public :: file_type ; private +type :: file_type ; private integer :: unit = -1 !< The framework identfier or netCDF unit number of an output file type(FmsNetcdfDomainFile_t), pointer :: fileobj => NULL() !< A domain-decomposed !! file object that is open for writing @@ -115,7 +115,7 @@ module MOM_io_infra end type file_type !> This type is a container for information about a variable in a file. -type, public :: fieldtype ; private +type :: fieldtype ; private character(len=256) :: name !< The name of this field in the files. type(mpp_fieldtype) :: FT !< The FMS1 field-type that this type wraps character(len=:), allocatable :: longname !< The long name for this field @@ -126,14 +126,14 @@ module MOM_io_infra end type fieldtype !> This type is a container for information about an axis in a file. -type, public :: axistype ; private +type :: axistype ; private character(len=256) :: name !< The name of this axis in the files. type(mpp_axistype) :: AT !< The FMS1 axis-type that this type wraps real, allocatable, dimension(:) :: ax_data !< The values of the data on the axis. end type axistype -!> For now, this is hard-coded to exercise the new FMS2 interfaces. -logical :: FMS2_reads = .true. +!> For now, these module-variables are hard-coded to exercise the new FMS2 interfaces. +logical :: FMS2_reads = .true. logical :: FMS2_writes = .true. contains @@ -353,8 +353,7 @@ subroutine open_file_type(IO_handle, filename, action, MOM_domain, threading, fi call fms2_close_file(fileObj_read) endif - success = fms2_open_file(IO_handle%fileobj, trim(filename), trim(mode), & - MOM_domain%mpp_domain, is_restart=.false.) + success = fms2_open_file(IO_handle%fileobj, trim(filename), trim(mode), MOM_domain%mpp_domain) if (.not.success) call MOM_error(FATAL, "Unable to open file "//trim(filename)) IO_handle%FMS2_file = .true. elseif (present(MOM_Domain)) then @@ -442,7 +441,8 @@ subroutine get_file_times(IO_handle, time_values, ntime) integer, optional, intent(out) :: ntime !< The number of time levels in the file character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file - integer :: ntimes + integer :: ntimes ! The number of time levels in the file + !### Modify this routine to optionally convert to time_type, using information about the dimensions? if (allocated(time_values)) deallocate(time_values) @@ -465,12 +465,13 @@ subroutine get_file_fields(IO_handle, fields) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O type(fieldtype), dimension(:), intent(inout) :: fields !< Field-type descriptions of all of !! the variables in a file. - type(mpp_fieldtype), dimension(size(fields)) :: mpp_fields - character(len=256), dimension(size(fields)) :: var_names - character(len=256) :: units - character(len=2048) :: longname - integer(kind=int64), dimension(3) :: checksum_file - integer :: i, nvar + type(mpp_fieldtype), dimension(size(fields)) :: mpp_fields ! Fieldtype structures for the variables + character(len=256), dimension(size(fields)) :: var_names ! The names of all variables + character(len=256) :: units ! The units of a variable as recorded in the file + character(len=2048) :: longname ! The long-name of a variable as recorded in the file + integer(kind=int64), dimension(3) :: checksum_file ! The checksums for a variable in the file + integer :: nvar ! The number of variables in the file + integer :: i nvar = size(fields) ! Local variables @@ -647,20 +648,40 @@ subroutine MOM_read_data_0d(filename, fieldname, data, timelevel, scale, MOM_Dom optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition ! Local variables - type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object - logical :: has_time_dim ! True if the variable has an unlimited time axis. + type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file + type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object character(len=96) :: var_to_read ! Name of variable to read from the netcdf file - logical :: success ! True if the file was successfully opened + logical :: has_time_dim ! True if the variable has an unlimited time axis. + logical :: success ! True if the file was successfully opened if (present(MOM_Domain) .and. FMS2_reads) then ! Open the FMS2 file-set. - success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain) if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) ! Find the matching case-insensitive variable name in the file and prepare to read it. - call prepare_to_read_var(fileobj, fieldname, MOM_domain, "MOM_read_data_1d: ", filename, & + call prepare_to_read_var(fileobj_DD, fieldname, "MOM_read_data_0d: ", filename, & var_to_read, has_time_dim, timelevel) + ! Read the data. + if (present(timelevel) .and. has_time_dim) then + call fms2_read_data(fileobj_DD, var_to_read, data, unlim_dim_level=timelevel) + else + call fms2_read_data(fileobj_DD, var_to_read, data) + endif + + ! Close the file-set. + if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD) + elseif (FMS2_reads) then + ! Open the FMS2 file-set. + success = fms2_open_file(fileObj, trim(filename), "read") + if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) + + ! Find the matching case-insensitive variable name in the file, and determine whether it + ! has a time dimension. + call find_varname_in_file(fileObj, fieldname, "MOM_read_data_0d: ", filename, & + var_to_read, has_time_dim, timelevel) + ! Read the data. if (present(timelevel) .and. has_time_dim) then call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel) @@ -695,20 +716,40 @@ subroutine MOM_read_data_1d(filename, fieldname, data, timelevel, scale, MOM_Dom optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition ! Local variables - type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object - logical :: has_time_dim ! True if the variable has an unlimited time axis. + type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file + type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object character(len=96) :: var_to_read ! Name of variable to read from the netcdf file - logical :: success ! True if the file was successfully opened + logical :: has_time_dim ! True if the variable has an unlimited time axis. + logical :: success ! True if the file was successfully opened if (present(MOM_Domain) .and. FMS2_reads) then ! Open the FMS2 file-set. - success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain) if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) ! Find the matching case-insensitive variable name in the file and prepare to read it. - call prepare_to_read_var(fileobj, fieldname, MOM_domain, "MOM_read_data_1d: ", filename, & + call prepare_to_read_var(fileobj_DD, fieldname, "MOM_read_data_1d: ", filename, & var_to_read, has_time_dim, timelevel) + ! Read the data. + if (present(timelevel) .and. has_time_dim) then + call fms2_read_data(fileobj_DD, var_to_read, data, unlim_dim_level=timelevel) + else + call fms2_read_data(fileobj_DD, var_to_read, data) + endif + + ! Close the file-set. + if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD) + elseif (FMS2_reads) then + ! Open the FMS2 file-set. + success = fms2_open_file(fileObj, trim(filename), "read") + if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) + + ! Find the matching case-insensitive variable name in the file, and determine whether it + ! has a time dimension. + call find_varname_in_file(fileObj, fieldname, "MOM_read_data_1d: ", filename, & + var_to_read, has_time_dim, timelevel) + ! Read the data. if (present(timelevel) .and. has_time_dim) then call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel) @@ -747,17 +788,17 @@ subroutine MOM_read_data_2d(filename, fieldname, data, MOM_Domain, & ! Local variables type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object - logical :: has_time_dim ! True if the variable has an unlimited time axis. character(len=96) :: var_to_read ! Name of variable to read from the netcdf file - logical :: success ! True if the file was successfully opened + logical :: has_time_dim ! True if the variable has an unlimited time axis. + logical :: success ! True if the file was successfully opened if (FMS2_reads) then ! Open the FMS2 file-set. - success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain) if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) ! Find the matching case-insensitive variable name in the file and prepare to read it. - call prepare_to_read_var(fileobj, fieldname, MOM_domain, "MOM_read_data_2d: ", filename, & + call prepare_to_read_var(fileobj, fieldname, "MOM_read_data_2d: ", filename, & var_to_read, has_time_dim, timelevel, position) ! Read the data. @@ -801,9 +842,41 @@ subroutine MOM_read_data_2d_region(filename, fieldname, data, start, nread, MOM_ real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. - !### This subroutine does not have an FMS-2 variant yet. + ! Local variables + type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file + type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object + character(len=96) :: var_to_read ! Name of variable to read from the netcdf file + logical :: success ! True if the file was successfully opened - if (present(MOM_Domain)) then + if (present(MOM_Domain) .and. FMS2_reads) then + ! Open the FMS2 file-set. + success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain) + if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) + + ! Find the matching case-insensitive variable name in the file and prepare to read it. + call prepare_to_read_var(fileobj_DD, fieldname, "MOM_read_data_2d_region: ", & + filename, var_to_read) + + ! Read the data. + call fms2_read_data(fileobj_DD, var_to_read, data, corner=start(1:2), edge_lengths=nread(1:2)) + + ! Close the file-set. + if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD) + elseif (FMS2_reads) then + ! Open the FMS2 file-set. + success = fms2_open_file(fileObj, trim(filename), "read") + if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) + + ! Find the matching case-insensitive variable name in the file, and determine whether it + ! has a time dimension. + call find_varname_in_file(fileObj, fieldname, "MOM_read_data_2d_region: ", filename, var_to_read) + + ! Read the data. + call fms2_read_data(fileobj, var_to_read, data, corner=start(1:2), edge_lengths=nread(1:2)) + + ! Close the file-set. + if (check_if_open(fileobj)) call fms2_close_file(fileobj) + elseif (present(MOM_Domain)) then ! Read the variable using the FMS-1 interface. call read_data(filename, fieldname, data, start, nread, domain=MOM_Domain%mpp_domain, & no_domain=no_domain) else @@ -838,17 +911,17 @@ subroutine MOM_read_data_3d(filename, fieldname, data, MOM_Domain, & ! Local variables type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object - logical :: has_time_dim ! True if the variable has an unlimited time axis. character(len=96) :: var_to_read ! Name of variable to read from the netcdf file - logical :: success ! True if the file was successfully opened + logical :: has_time_dim ! True if the variable has an unlimited time axis. + logical :: success ! True if the file was successfully opened if (FMS2_reads) then ! Open the FMS2 file-set. - success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain) if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) ! Find the matching case-insensitive variable name in the file and prepare to read it. - call prepare_to_read_var(fileobj, fieldname, MOM_domain, "MOM_read_data_3d: ", filename, & + call prepare_to_read_var(fileobj, fieldname, "MOM_read_data_3d: ", filename, & var_to_read, has_time_dim, timelevel, position) ! Read the data. @@ -895,11 +968,11 @@ subroutine MOM_read_data_4d(filename, fieldname, data, MOM_Domain, & if (FMS2_reads) then ! Open the FMS2 file-set. - success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain) if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) ! Find the matching case-insensitive variable name in the file and prepare to read it. - call prepare_to_read_var(fileobj, fieldname, MOM_domain, "MOM_read_data_4d: ", filename, & + call prepare_to_read_var(fileobj, fieldname, "MOM_read_data_4d: ", filename, & var_to_read, has_time_dim, timelevel, position) ! Read the data. @@ -930,8 +1003,34 @@ subroutine MOM_read_data_0d_int(filename, fieldname, data, timelevel) integer, intent(inout) :: data !< The 1-dimensional array into which the data integer, optional, intent(in) :: timelevel !< The time level in the file to read - !### This needs an FMS2 variant, eventually. - call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) + ! Local variables + type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file + logical :: has_time_dim ! True if the variable has an unlimited time axis. + character(len=96) :: var_to_read ! Name of variable to read from the netcdf file + logical :: success ! If true, the file was opened successfully + + if (FMS2_reads) then + ! Open the FMS2 file-set. + success = fms2_open_file(fileObj, trim(filename), "read") + if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) + + ! Find the matching case-insensitive variable name in the file, and determine whether it + ! has a time dimension. + call find_varname_in_file(fileObj, fieldname, "MOM_read_data_0d_int: ", filename, & + var_to_read, has_time_dim, timelevel) + + ! Read the data. + if (present(timelevel) .and. has_time_dim) then + call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel) + else + call fms2_read_data(fileobj, var_to_read, data) + endif + + ! Close the file-set. + if (check_if_open(fileobj)) call fms2_close_file(fileobj) + else + call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) + endif end subroutine MOM_read_data_0d_int @@ -943,8 +1042,35 @@ subroutine MOM_read_data_1d_int(filename, fieldname, data, timelevel) integer, dimension(:), intent(inout) :: data !< The 1-dimensional array into which the data integer, optional, intent(in) :: timelevel !< The time level in the file to read - !### This needs an FMS2 variant, eventually. - call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) + ! Local variables + type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file for obtaining information + ! about the exiting time axis entries in append mode. + logical :: has_time_dim ! True if the variable has an unlimited time axis. + character(len=96) :: var_to_read ! Name of variable to read from the netcdf file + logical :: success ! If true, the file was opened successfully + + if (FMS2_reads) then + ! Open the FMS2 file-set. + success = fms2_open_file(fileObj, trim(filename), "read") + if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) + + ! Find the matching case-insensitive variable name in the file, and determine whether it + ! has a time dimension. + call find_varname_in_file(fileObj, fieldname, "MOM_read_data_1d_int: ", filename, & + var_to_read, has_time_dim, timelevel) + + ! Read the data. + if (present(timelevel) .and. has_time_dim) then + call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel) + else + call fms2_read_data(fileobj, var_to_read, data) + endif + + ! Close the file-set. + if (check_if_open(fileobj)) call fms2_close_file(fileobj) + else + call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) + endif end subroutine MOM_read_data_1d_int @@ -983,13 +1109,13 @@ subroutine MOM_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data if (FMS2_reads) then ! Open the FMS2 file-set. - success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain) if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) ! Find the matching case-insensitive u- and v-variable names in the file and prepare to read them. - call prepare_to_read_var(fileobj, u_fieldname, MOM_domain, "MOM_read_vector_2d: ", filename, & + call prepare_to_read_var(fileobj, u_fieldname, "MOM_read_vector_2d: ", filename, & u_var, has_time_dim, timelevel, position=u_pos) - call prepare_to_read_var(fileobj, v_fieldname, MOM_domain, "MOM_read_vector_2d: ", filename, & + call prepare_to_read_var(fileobj, v_fieldname, "MOM_read_vector_2d: ", filename, & v_var, has_time_dim, timelevel, position=v_pos) ! Read the u-data and v-data. There would already been an error message for one @@ -1053,13 +1179,13 @@ subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data if (FMS2_reads) then ! Open the FMS2 file-set. - success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain, is_restart=.false.) + success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain) if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename)) ! Find the matching case-insensitive u- and v-variable names in the file and prepare to read them. - call prepare_to_read_var(fileobj, u_fieldname, MOM_domain, "MOM_read_vector_3d: ", filename, & + call prepare_to_read_var(fileobj, u_fieldname, "MOM_read_vector_3d: ", filename, & u_var, has_time_dim, timelevel, position=u_pos) - call prepare_to_read_var(fileobj, v_fieldname, MOM_domain, "MOM_read_vector_3d: ", filename, & + call prepare_to_read_var(fileobj, v_fieldname, "MOM_read_vector_3d: ", filename, & v_var, has_time_dim, timelevel, position=v_pos) ! Read the u-data and v-data, dangerously assuming either both or neither have time dimensions. @@ -1363,7 +1489,6 @@ subroutine write_metadata_field(IO_handle, field, axes, name, units, longname, & call register_variable_attribute(IO_handle%fileobj, trim(name), "checksum", & trim(checksum_string), len_trim(checksum_string)) endif - !### Add more attributes if they are present; remove attributes that are never used. else do i=1,ndims ; mpp_axes(i) = axes(i)%AT ; enddo call mpp_write_meta(IO_handle%unit, field%FT, mpp_axes, name, units, longname, & diff --git a/config_src/infra/FMS2/MOM_read_data_fms2.F90 b/config_src/infra/FMS2/MOM_read_data_fms2.F90 index 4732c019f4..c632e95177 100644 --- a/config_src/infra/FMS2/MOM_read_data_fms2.F90 +++ b/config_src/infra/FMS2/MOM_read_data_fms2.F90 @@ -1,36 +1,31 @@ -!> This module contains routines that wrap the fms2 read_data calls +!> This module contains routines that encapsulate common preparatory work for FMS2 read_data calls module MOM_read_data_fms2 ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, FATAL, WARNING, is_root_PE -use MOM_domain_infra, only : MOM_domain_type, AGRID, BGRID_NE, CGRID_NE -use MOM_domain_infra, only : domain2d, CENTER, CORNER, NORTH_FACE, EAST_FACE +use MOM_domain_infra, only : CENTER, CORNER, NORTH_FACE, EAST_FACE use MOM_string_functions, only : lowercase -use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t -use fms2_io_mod, only : fms2_open_file => open_file, fms2_close_file => close_file -use fms2_io_mod, only : get_num_variables, get_variable_names, check_if_open -use fms2_io_mod, only : read_data, variable_exists, get_variable_size, get_variable_units -use fms2_io_mod, only : get_variable_attribute, attribute_exists => variable_att_exists +use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, check_if_open +use fms2_io_mod, only : variable_exists, get_num_variables, get_variable_names +use fms2_io_mod, only : get_variable_size, get_variable_units +use fms2_io_mod, only : get_variable_attribute, variable_att_exists use fms2_io_mod, only : get_variable_num_dimensions, get_variable_dimension_names use fms2_io_mod, only : is_dimension_unlimited, get_dimension_size use fms2_io_mod, only : is_dimension_registered, register_axis implicit none ; private -public prepare_to_read_var -! public MOM_read_data_scalar, MOM_read_data_2d_noDD, MOM_read_data_1d_noDD +public prepare_to_read_var, find_varname_in_file contains -!> Find the case-insensitive name match with a variable in a domain-decomposed file-set -!! opening the file(s) as necessary, prepare FMS2 to read this variable, and return some -!! information needed to call read_data correctly for this variable and file. -subroutine prepare_to_read_var(fileobj, fieldname, domain, err_header, filename, var_to_read, & +!> Find the case-insensitive name match with a variable in an open domain-decomposed file-set, +!! prepare FMS2 to read this variable, and return some information needed to call fms2_read_data +!! correctly for this variable and file. +subroutine prepare_to_read_var(fileobj, fieldname, err_header, filename, var_to_read, & has_time_dim, timelevel, position) - type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< A handle to an FMS2 file object, that - !! will be opened if necessary + type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< An FMS2 handle to an open domain-decomposed file character(len=*), intent(in) :: fieldname !< The variable name to seek in the file - type(MOM_domain_type), intent(in) :: domain !< MOM domain attribute with the mpp_domain decomposition character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(out) :: var_to_read !< The variable name to read from the file @@ -39,33 +34,30 @@ subroutine prepare_to_read_var(fileobj, fieldname, domain, err_header, filename, integer, optional, intent(in) :: position !< A flag indicating where this variable is discretized ! Local variables - logical :: file_open_success !.true. if call to open_file is successful logical :: variable_found ! Is a case-insensitive version of the variable found in the netCDF file? - character(len=96), allocatable, dimension(:) :: var_names !< array for names of variables in a netCDF - !! file opened to read - character(len=96), allocatable :: dim_names(:) ! variable dimension names - integer :: nvars ! The number of variables in the file. - integer :: i, dim_unlim_size, num_var_dims, time_dim + character(len=256), allocatable, dimension(:) :: var_names ! The names of all the variables in the netCDF file + character(len=256), allocatable :: dim_names(:) ! The names of a variable's dimensions + integer :: nvars ! The number of variables in the file. + integer :: dim_unlim_size ! The current size of the unlimited (time) dimension in the file. + integer :: num_var_dims ! The number of dimensions a variable has in the file. + integer :: time_dim ! The position of the unlimited (time) dimension for a variable, or -1 + ! if it has no unlimited dimension. + integer :: i ! Open the file if necessary - if (.not.(check_if_open(fileobj))) then - file_open_success = fms2_open_file(fileobj, filename, "read", domain%mpp_domain, is_restart=.false.) - if (.not.file_open_success) call MOM_error(FATAL, trim(err_header)//" failed to open "//trim(filename)) - endif + if (.not.check_if_open(fileobj)) & + call MOM_error(FATAL, trim(err_header)//trim(filename)//" was not open in call to prepare_to_read_var.") ! Search for the variable in the file, looking for the case-sensitive name first. if (variable_exists(fileobj, trim(fieldname))) then var_to_read = trim(fieldname) - variable_found = .true. else ! Look for case-insensitive variable name matches. - var_to_read = "" - variable_found = .false. - nvars = get_num_variables(fileobj) if (nvars < 1) call MOM_error(FATAL, "nvars is less than 1 for file "//trim(filename)) allocate(var_names(nvars)) call get_variable_names(fileobj, var_names) + variable_found = .false. do i=1,nvars if (lowercase(trim(var_names(i))) == lowercase(trim(fieldname))) then variable_found = .true. @@ -116,8 +108,6 @@ subroutine prepare_to_read_var(fileobj, fieldname, domain, err_header, filename, end subroutine prepare_to_read_var !> register axes associated with a variable from a domain-decomposed netCDF file -!> @note The user must specify units for variables with longitude/x-axis and/or latitude/y-axis axes -!! to obtain the correct domain decomposition for the data buffer. subroutine MOM_register_variable_axes(fileObj, variableName, filename, position) type(FmsNetcdfDomainFile_t), intent(inout) :: fileObj !< Handle to an open FMS2 netCDF file object character(len=*), intent(in) :: variableName !< name of the variable @@ -125,19 +115,15 @@ subroutine MOM_register_variable_axes(fileObj, variableName, filename, position) integer, optional, intent(in) :: position !< A flag indicating where this data is discretized ! Local variables - character(len=40) :: units ! units corresponding to a specific variable dimension - character(len=40), allocatable, dimension(:) :: dim_names ! variable dimension names + character(len=256), allocatable, dimension(:) :: dim_names ! variable dimension names integer, allocatable, dimension(:) :: dimSizes ! variable dimension sizes logical, allocatable, dimension(:) :: is_x ! Is this a (likely domain-decomposed) x-axis logical, allocatable, dimension(:) :: is_y ! Is this a (likely domain-decomposed) y-axis logical, allocatable, dimension(:) :: is_t ! Is this a time axis or another unlimited axis integer :: ndims ! number of dimensions + integer :: xPos, yPos ! Discrete positions for x and y axes. Default is CENTER integer :: i - integer :: xPos, yPos ! domain positions for x and y axes. Default is CENTER - if (.not. check_if_open(fileObj)) call MOM_error(FATAL,"MOM_axis:register_variable_axes: The fileObj has "// & - "not been opened. Call fms2_open_file(fileObj,...) before "// & - "passing the fileObj argument to this function.") xPos = CENTER ; yPos = CENTER if (present(position)) then if ((position == CORNER) .or. (position == EAST_FACE)) xPos = EAST_FACE @@ -168,9 +154,7 @@ subroutine MOM_register_variable_axes(fileObj, variableName, filename, position) endif enddo - deallocate(dimSizes) - deallocate(dim_names) - deallocate(is_x, is_y, is_t) + deallocate(dimSizes, dim_names, is_x, is_y, is_t) end subroutine MOM_register_variable_axes !> Determine whether a variable's axes are associated with x-, y- or time-dimensions. Other @@ -184,11 +168,12 @@ subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t logical, dimension(ndims), intent(out) :: is_y !< Indicates if each dimension a (likely decomposed) y-axis logical, dimension(ndims), intent(out) :: is_t !< Indicates if each dimension unlimited (usually time) axis - integer :: i - character(len=256) :: cartesian ! A flag indicating a Cartesian direction - usually a single character. + ! Local variables + character(len=128) :: cartesian ! A flag indicating a Cartesian direction - usually a single character. character(len=512) :: dim_list ! A concatenated list of dimension names. - character(len=40) :: units ! units corresponding to a specific variable dimension + character(len=128) :: units ! units corresponding to a specific variable dimension logical :: x_found, y_found ! Indicate whether an x- or y- dimension have been found. + integer :: i x_found = .false. ; y_found = .false. is_x(:) = .false. ; is_y(:) = .false. @@ -197,16 +182,12 @@ subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t ! First look for indicative variable attributes if (.not.is_t(i)) then if (variable_exists(fileobj, trim(dim_names(i)))) then - if (attribute_exists(fileobj, trim(dim_names(i)), "cartesian_axis")) then + if (variable_att_exists(fileobj, trim(dim_names(i)), "cartesian_axis")) then call get_variable_attribute(fileobj, trim(dim_names(i)), "cartesian_axis", cartesian) cartesian = adjustl(cartesian) if ((index(cartesian, "X") == 1) .or. (index(cartesian, "x") == 1)) is_x(i) = .true. if ((index(cartesian, "Y") == 1) .or. (index(cartesian, "y") == 1)) is_y(i) = .true. if ((index(cartesian, "T") == 1) .or. (index(cartesian, "t") == 1)) is_t(i) = .true. - ! if (is_root_pe() .and. is_x(i)) & - ! call MOM_error(NOTE, "X-dimension determined from cartesian_axis for "//trim(dim_names(i))) - ! if (is_root_pe() .and. is_y(i)) & - ! call MOM_error(NOTE, "Y-dimension determined from cartesian_axis for "//trim(dim_names(i))) endif endif endif @@ -308,199 +289,44 @@ subroutine categorize_axis_from_name(dimname, is_x, is_y) end subroutine categorize_axis_from_name -!===== Everything below this pertains to reading non-decomposed variables ===! -!===== using FMS2 interfaces will probably be discarded eventually. =========! - -!!> This routine calls the fms_io read_data subroutine to read a scalar (0-D) field named "fieldname" -!! from file "filename". -subroutine MOM_read_data_scalar(filename, fieldname, data, timelevel, scale, leave_file_open) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, intent(inout) :: data !< The variable to read from read_data - integer, optional, intent(in) :: timelevel !< time level to read - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by - logical, optional, intent(in) :: leave_file_open !< if .true., leave file open - - ! Local variables - type(FmsNetcdfFile_t) :: fileobj ! A handle to a simple netCDF file - logical :: close_the_file ! indicates whether to close the file after read_data is called. - character(len=96) :: var_to_read ! variable to read from the netcdf file - character(len=48) :: err_header ! A preamble for error messages - - err_header = "MOM_read_data_fms2:MOM_read_data_scalar: " - - ! Find the matching variable name in the file, opening it and reading metadata if necessary. - call find_varname_in_file(fileobj, fieldname, err_header, filename, var_to_read) - - ! read the data - if (present(timelevel)) then - call read_data(fileobj, trim(var_to_read), data, unlim_dim_level=timelevel) - else - call read_data(fileobj, trim(var_to_read), data) - endif - - ! Close the file, if necessary - close_the_file = .true. ; if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - if (close_the_file .and. check_if_open(fileobj)) call fms2_close_file(fileobj) - - ! Rescale the data that was read if necessary. - if (present(scale)) then ; if (scale /= 1.0) then - data = scale*data - endif ; endif - -end subroutine MOM_read_data_scalar - -!> This routine calls the fms_io read_data subroutine to read 1-D non-domain-decomposed data field named "fieldname" -!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. -subroutine MOM_read_data_1d_noDD(filename, fieldname, data, start_index, & - edge_lengths, timelevel, scale, leave_file_open) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, dimension(:), intent(inout) :: data !< The 1-dimensional data array to pass to read_data - integer, dimension(1), optional, intent(in) :: start_index !< starting index of data buffer. Default is 1 - integer, dimension(1), optional, intent(in) :: edge_lengths !< number of data values to read in; default is - !! the variable size - integer, optional, intent(in) :: timelevel !< time level to read - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by - logical, optional, intent(in) :: leave_file_open !< if .true., leave file open - - ! Local variables - type(FmsNetcdfFile_t) :: fileobj ! A handle to a simple netCDF file - logical :: close_the_file ! indicates whether to close the file after read_data is called. - integer :: time_dim ! The dimension position of a variables unlimited time axis, or -1 if it has none. - integer, parameter :: ndim = 1 ! The dimensionality of the array being read - integer, dimension(ndim) :: start, nread ! indices for first data value and number of values to read - character(len=96) :: var_to_read ! variable to read from the netcdf file - character(len=48) :: err_header ! A preamble for error messages - - err_header = "MOM_read_data_fms2:MOM_read_data_1d_noDD: " - - ! Find the matching case-insensitive variable name in the file, opening the file if necessary. - call find_varname_in_file(fileobj, fieldname, err_header, filename, var_to_read) - - ! set the start and nread values that will be passed as the read_data corner and edge_lengths arguments - start(:) = 1 ; if (present(start_index)) start(:) = start_index(:) - nread(:) = shape(data) ; if (present(edge_lengths)) nread(:) = edge_lengths(:) - - time_dim = -1 - if (present(timelevel)) then - time_dim = get_time_dim(fileobj, var_to_read, err_header, filename, timelevel) - if (time_dim == ndim) then ; nread(ndim) = 1 ; start(ndim) = timelevel ; endif - endif - - ! read the data - if (time_dim > 0) then - call read_data(fileobj, trim(var_to_read), data, corner=start, edge_lengths=nread, & - unlim_dim_level=timelevel) - else - call read_data(fileobj, trim(var_to_read), data, corner=start, edge_lengths=nread) - endif - - ! Close the file, if necessary - close_the_file = .true. ; if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - if (close_the_file .and. check_if_open(fileobj)) call fms2_close_file(fileobj) - - ! Rescale the data that was read if necessary. - if (present(scale)) then ; if (scale /= 1.0) then - data(:) = scale*data(:) - endif ; endif - -end subroutine MOM_read_data_1d_noDD - -!> This routine calls the fms_io read_data subroutine to read a 2-D non-domain-decomposed data field named "fieldname" -!! from file "filename". The routine multiplies the data by "scale" if the optional argument is included in the call. -subroutine MOM_read_data_2d_noDD(filename, fieldname, data, start_index, & - edge_lengths, timelevel, position, scale, leave_file_open) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, dimension(:,:), intent(inout) :: data !< The 2-dimensional data array to pass to read_data - integer, dimension(2), optional, intent(in) :: start_index !< starting indices of data buffer. Default is 1 - integer, dimension(2), optional, intent(in) :: edge_lengths !< number of data values to read in. - !! Default values are the variable dimension sizes - integer, optional, intent(in) :: timelevel !< time level to read - integer, optional, intent(in) :: position !< A flag indicating where this data is located - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by - logical, optional, intent(in) :: leave_file_open !< if .true., leave file open - - ! Local variables - type(FmsNetcdfFile_t) :: fileobj ! A handle to a simple netCDF file - logical :: close_the_file ! indicates whether to close the file after read_data is called. - integer :: time_dim ! The dimension position of a variables unlimited time axis, or -1 if it has none. - integer, parameter :: ndim = 2 ! The dimensionality of the array being read - integer, dimension(ndim) :: start, nread ! indices for first data value and number of values to read - character(len=96) :: var_to_read ! variable to read from the netcdf file - character(len=48) :: err_header ! A preamble for error messages - - err_header = "MOM_read_data_fms2:MOM_read_data_2d_DD: " - - ! Find the matching case-insensitive variable name in the file, opening the file if necessary. - call find_varname_in_file(fileobj, fieldname, err_header, filename, var_to_read) - - ! set the start and nread values that will be passed as the read_data corner and edge_lengths arguments - start(:) = 1 ; if (present(start_index)) start(:) = start_index(:) - nread(:) = shape(data) ; if (present(edge_lengths)) nread(:) = edge_lengths(:) - - time_dim = -1 - if (present(timelevel)) then - time_dim = get_time_dim(fileobj, var_to_read, err_header, filename, timelevel) - if (time_dim == ndim) then ; nread(ndim) = 1 ; start(ndim) = timelevel ; endif - endif - - ! read the data - if (time_dim > 0) then - call read_data(fileobj, trim(var_to_read), data, corner=start, edge_lengths=nread, & - unlim_dim_level=timelevel) - else - call read_data(fileobj, trim(var_to_read), data, corner=start, edge_lengths=nread) - endif - - ! Close the file, if necessary - close_the_file = .true. ; if (present(leave_file_open)) close_the_file = .not.(leave_file_open) - if (close_the_file .and. check_if_open(fileobj)) call fms2_close_file(fileobj) - - ! Rescale the data that was read if necessary. - if (present(scale)) then ; if (scale /= 1.0) then - data(:,:) = scale*data(:,:) - endif ; endif - -end subroutine MOM_read_data_2d_noDD - !> Find the case-sensitive name of the variable in a netCDF file with a case-insensitive name match. -subroutine find_varname_in_file(fileobj, fieldname, err_header, filename, var_to_read) - type(FmsNetcdfFile_t), intent(inout) :: fileobj !< A handle to a file object, that - !! will be opened if necessary +!! Optionally also determine whether this variable has an unlimited time dimension. +subroutine find_varname_in_file(fileobj, fieldname, err_header, filename, var_to_read, has_time_dim, timelevel) + type(FmsNetcdfFile_t), intent(inout) :: fileobj !< An FMS2 handle to an open NetCDF file character(len=*), intent(in) :: fieldname !< The variable name to seek in the file character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(out) :: var_to_read !< The variable name to read from the file + logical, optional, intent(out) :: has_time_dim !< Indicates whether fieldname has a time dimension + integer, optional, intent(in) :: timelevel !< A time level to read ! Local variables - logical :: file_open_success !.true. if call to open_file is successful logical :: variable_found ! Is a case-insensitive version of the variable found in the netCDF file? - character(len=96), allocatable, dimension(:) :: var_names !< array for names of variables in a netCDF - !! file opened to read - integer :: nvars ! The number of variables in the file. + character(len=256), allocatable, dimension(:) :: var_names ! The names of all the variables in the netCDF file + character(len=256), allocatable :: dim_names(:) ! The names of a variable's dimensions + integer :: nvars ! The number of variables in the file + integer :: dim_unlim_size ! The current size of the unlimited (time) dimension in the file. + integer :: num_var_dims ! The number of dimensions a variable has in the file. + integer :: time_dim ! The position of the unlimited (time) dimension for a variable, or -1 + ! if it has no unlimited dimension. integer :: i - var_to_read = "" - ! Open the file if necessary - if (.not.(check_if_open(fileobj))) then - file_open_success = fms2_open_file(fileobj, filename, "read", is_restart=.false.) - if (.not.file_open_success) call MOM_error(FATAL, trim(err_header)//" failed to open "//trim(filename)) - endif + if (.not.check_if_open(fileobj)) & + call MOM_error(FATAL, trim(err_header)//trim(filename)//" was not open in call to find_varname_in_file.") - if (variable_exists(fileobj, fieldname)) then - var_to_read = fieldname - else - variable_found = .false. + ! Search for the variable in the file, looking for the case-sensitive name first. + if (variable_exists(fileobj, trim(fieldname))) then + var_to_read = trim(fieldname) + else ! Look for case-insensitive variable name matches. nvars = get_num_variables(fileobj) if (nvars < 1) call MOM_error(FATAL, "nvars is less than 1 for file "//trim(filename)) allocate(var_names(nvars)) call get_variable_names(fileobj, var_names) ! search for the variable in the file + variable_found = .false. do i=1,nvars if (lowercase(trim(var_names(i))) == lowercase(trim(fieldname))) then variable_found = .true. @@ -513,43 +339,38 @@ subroutine find_varname_in_file(fileobj, fieldname, err_header, filename, var_to deallocate(var_names) endif -end subroutine find_varname_in_file + ! FMS2 can not handle a timelevel argument if the variable does not have one in the file, + ! so some error checking and logic are required. + if (present(has_time_dim) .or. present(timelevel)) then + time_dim = -1 -!> Return the number of the time dimension for a variable in an open non-domain-decomposed file, -!! or -1 if it has no time (or other unlimited) dimension. -integer function get_time_dim(fileobj, var_to_read, err_header, filename, timelevel) - type(FmsNetcdfFile_t), intent(in) :: fileobj !< A handle to an open file object - character(len=*), intent(in) :: var_to_read !< The variable name to read from the file - character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages - character(len=*), intent(in) :: filename !< The name of the file to read - integer, optional, intent(in) :: timelevel !< A time level to read + num_var_dims = get_variable_num_dimensions(fileobj, trim(var_to_read)) + allocate(dim_names(num_var_dims)) ; dim_names(:) = "" + call get_variable_dimension_names(fileobj, trim(var_to_read), dim_names) - ! Local variables - integer :: i, dim_unlim_size, num_var_dims - character(len=96), allocatable :: dim_names(:) ! variable dimension names - - num_var_dims = get_variable_num_dimensions(fileobj, trim(var_to_read)) - allocate(dim_names(num_var_dims)) ; dim_names(:) = "" - call get_variable_dimension_names(fileobj, trim(var_to_read), dim_names) - - get_time_dim = -1 - do i=1,num_var_dims - if (is_dimension_unlimited(fileobj, dim_names(i))) then - get_time_dim = i - if (present(timelevel)) then - call get_dimension_size(fileobj, dim_names(i), dim_unlim_size) - if (timelevel > dim_unlim_size) call MOM_error(FATAL, trim(err_header)//& - "Attempting to read a time level of "//trim(var_to_read)//& - " that exceeds the size of "//trim(filename)) + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj, dim_names(i))) then + time_dim = i + if (present(timelevel)) then + call get_dimension_size(fileobj, dim_names(i), dim_unlim_size) + if ((timelevel > dim_unlim_size) .and. is_root_PE()) call MOM_error(FATAL, & + trim(err_header)//"Attempting to read a time level of "//trim(var_to_read)//& + " that exceeds the size of the time dimension in "//trim(filename)) + endif + exit endif - exit - endif - enddo - if (get_time_dim < 0) & - call MOM_error(WARNING, trim(err_header)//"time level specified, but the variable "//& + enddo + deallocate(dim_names) + + if (present(timelevel) .and. (time_dim < 0) .and. is_root_PE()) & + call MOM_error(WARNING, trim(err_header)//"time level specified, but the variable "//& trim(var_to_read)//" does not have an unlimited dimension in "//trim(filename)) - deallocate(dim_names) + if ((.not.present(timelevel)) .and. (time_dim > 0) .and. is_root_PE()) & + call MOM_error(WARNING, trim(err_header)//"The variable "//trim(var_to_read)//& + " has an unlimited dimension in "//trim(filename)//" but no time level is specified.") + if (present(has_time_dim)) has_time_dim = (time_dim > 0) + endif -end function get_time_dim +end subroutine find_varname_in_file end module MOM_read_data_fms2 From 3fe07d4c66918249e91fa192fae76e4572655d00 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 23 Mar 2021 05:55:12 -0400 Subject: [PATCH 07/10] +Move FMS2 read helper routines into MOM_io_infra Moved the routines prepare_to_read_var and find_varname_in_file, along with four other routines that prepare_to_read_var calls from MOM_read_data_fms2 to MOM_io_infra, and eliminated the file MOM_read_data_fms2.F90. All answers are bitwise identical, but this rearrangement of identical code does eliminate one public module, so that the file structure of infra/FMS2 is identical to that of infra/FMS1. --- config_src/infra/FMS2/MOM_io_infra.F90 | 364 +++++++++++++++++- config_src/infra/FMS2/MOM_read_data_fms2.F90 | 376 ------------------- 2 files changed, 360 insertions(+), 380 deletions(-) delete mode 100644 config_src/infra/FMS2/MOM_read_data_fms2.F90 diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index c545462ad8..07a7c798d6 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -5,16 +5,17 @@ module MOM_io_infra use MOM_domain_infra, only : MOM_domain_type, rescale_comp_data, AGRID, BGRID_NE, CGRID_NE use MOM_domain_infra, only : domain2d, domain1d, CENTER, CORNER, NORTH_FACE, EAST_FACE -use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, FATAL, WARNING -use MOM_read_data_fms2, only : prepare_to_read_var, find_varname_in_file +use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, FATAL, WARNING, is_root_PE +use MOM_string_functions, only : lowercase use fms2_io_mod, only : fms2_open_file => open_file, check_if_open, fms2_close_file => close_file use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, fms2_read_data => read_data use fms2_io_mod, only : get_unlimited_dimension_name, get_num_dimensions, get_num_variables -use fms2_io_mod, only : get_variable_names, variable_exists, get_variable_size +use fms2_io_mod, only : get_variable_names, variable_exists, get_variable_size, get_variable_units use fms2_io_mod, only : register_field, write_data, register_variable_attribute use fms2_io_mod, only : variable_att_exists, get_variable_attribute, get_variable_num_dimensions -use fms2_io_mod, only : get_dimension_size, is_dimension_registered, register_axis, unlimited +use fms2_io_mod, only : get_variable_dimension_names, is_dimension_registered, get_dimension_size +use fms2_io_mod, only : is_dimension_unlimited, register_axis, unlimited use fms_mod, only : write_version_number, open_namelist_file, check_nml_error use fms_io_mod, only : file_exist, field_exist, field_size, read_data @@ -1215,6 +1216,361 @@ subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data end subroutine MOM_read_vector_3d +!> Find the case-sensitive name of the variable in a netCDF file with a case-insensitive name match. +!! Optionally also determine whether this variable has an unlimited time dimension. +subroutine find_varname_in_file(fileobj, fieldname, err_header, filename, var_to_read, has_time_dim, timelevel) + type(FmsNetcdfFile_t), intent(inout) :: fileobj !< An FMS2 handle to an open NetCDF file + character(len=*), intent(in) :: fieldname !< The variable name to seek in the file + character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(out) :: var_to_read !< The variable name to read from the file + logical, optional, intent(out) :: has_time_dim !< Indicates whether fieldname has a time dimension + integer, optional, intent(in) :: timelevel !< A time level to read + + ! Local variables + logical :: variable_found ! Is a case-insensitive version of the variable found in the netCDF file? + character(len=256), allocatable, dimension(:) :: var_names ! The names of all the variables in the netCDF file + character(len=256), allocatable :: dim_names(:) ! The names of a variable's dimensions + integer :: nvars ! The number of variables in the file + integer :: dim_unlim_size ! The current size of the unlimited (time) dimension in the file. + integer :: num_var_dims ! The number of dimensions a variable has in the file. + integer :: time_dim ! The position of the unlimited (time) dimension for a variable, or -1 + ! if it has no unlimited dimension. + integer :: i + + ! Open the file if necessary + if (.not.check_if_open(fileobj)) & + call MOM_error(FATAL, trim(err_header)//trim(filename)//" was not open in call to find_varname_in_file.") + + ! Search for the variable in the file, looking for the case-sensitive name first. + if (variable_exists(fileobj, trim(fieldname))) then + var_to_read = trim(fieldname) + else ! Look for case-insensitive variable name matches. + nvars = get_num_variables(fileobj) + if (nvars < 1) call MOM_error(FATAL, "nvars is less than 1 for file "//trim(filename)) + allocate(var_names(nvars)) + call get_variable_names(fileobj, var_names) + + ! search for the variable in the file + variable_found = .false. + do i=1,nvars + if (lowercase(trim(var_names(i))) == lowercase(trim(fieldname))) then + variable_found = .true. + var_to_read = trim(var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) & + call MOM_error(FATAL, trim(err_header)//trim(fieldname)//" not found in "//trim(filename)) + deallocate(var_names) + endif + + ! FMS2 can not handle a timelevel argument if the variable does not have one in the file, + ! so some error checking and logic are required. + if (present(has_time_dim) .or. present(timelevel)) then + time_dim = -1 + + num_var_dims = get_variable_num_dimensions(fileobj, trim(var_to_read)) + allocate(dim_names(num_var_dims)) ; dim_names(:) = "" + call get_variable_dimension_names(fileobj, trim(var_to_read), dim_names) + + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj, dim_names(i))) then + time_dim = i + if (present(timelevel)) then + call get_dimension_size(fileobj, dim_names(i), dim_unlim_size) + if ((timelevel > dim_unlim_size) .and. is_root_PE()) call MOM_error(FATAL, & + trim(err_header)//"Attempting to read a time level of "//trim(var_to_read)//& + " that exceeds the size of the time dimension in "//trim(filename)) + endif + exit + endif + enddo + deallocate(dim_names) + + if (present(timelevel) .and. (time_dim < 0) .and. is_root_PE()) & + call MOM_error(WARNING, trim(err_header)//"time level specified, but the variable "//& + trim(var_to_read)//" does not have an unlimited dimension in "//trim(filename)) + if ((.not.present(timelevel)) .and. (time_dim > 0) .and. is_root_PE()) & + call MOM_error(WARNING, trim(err_header)//"The variable "//trim(var_to_read)//& + " has an unlimited dimension in "//trim(filename)//" but no time level is specified.") + if (present(has_time_dim)) has_time_dim = (time_dim > 0) + endif + +end subroutine find_varname_in_file + + +!> Find the case-insensitive name match with a variable in an open domain-decomposed file-set, +!! prepare FMS2 to read this variable, and return some information needed to call fms2_read_data +!! correctly for this variable and file. +subroutine prepare_to_read_var(fileobj, fieldname, err_header, filename, var_to_read, & + has_time_dim, timelevel, position) + type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< An FMS2 handle to an open domain-decomposed file + character(len=*), intent(in) :: fieldname !< The variable name to seek in the file + character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(out) :: var_to_read !< The variable name to read from the file + logical, optional, intent(out) :: has_time_dim !< Indicates whether fieldname has a time dimension + integer, optional, intent(in) :: timelevel !< A time level to read + integer, optional, intent(in) :: position !< A flag indicating where this variable is discretized + + ! Local variables + logical :: variable_found ! Is a case-insensitive version of the variable found in the netCDF file? + character(len=256), allocatable, dimension(:) :: var_names ! The names of all the variables in the netCDF file + character(len=256), allocatable :: dim_names(:) ! The names of a variable's dimensions + integer :: nvars ! The number of variables in the file. + integer :: dim_unlim_size ! The current size of the unlimited (time) dimension in the file. + integer :: num_var_dims ! The number of dimensions a variable has in the file. + integer :: time_dim ! The position of the unlimited (time) dimension for a variable, or -1 + ! if it has no unlimited dimension. + integer :: i + + ! Open the file if necessary + if (.not.check_if_open(fileobj)) & + call MOM_error(FATAL, trim(err_header)//trim(filename)//" was not open in call to prepare_to_read_var.") + + ! Search for the variable in the file, looking for the case-sensitive name first. + if (variable_exists(fileobj, trim(fieldname))) then + var_to_read = trim(fieldname) + else ! Look for case-insensitive variable name matches. + nvars = get_num_variables(fileobj) + if (nvars < 1) call MOM_error(FATAL, "nvars is less than 1 for file "//trim(filename)) + allocate(var_names(nvars)) + call get_variable_names(fileobj, var_names) + + variable_found = .false. + do i=1,nvars + if (lowercase(trim(var_names(i))) == lowercase(trim(fieldname))) then + variable_found = .true. + var_to_read = trim(var_names(i)) + exit + endif + enddo + if (.not.(variable_found)) & + call MOM_error(FATAL, trim(err_header)//trim(fieldname)//" not found in "//trim(filename)) + deallocate(var_names) + endif + + ! FMS2 can not handle a timelevel argument if the variable does not have one in the file, + ! so some error checking and logic are required. + if (present(has_time_dim) .or. present(timelevel)) then + time_dim = -1 + + num_var_dims = get_variable_num_dimensions(fileobj, trim(var_to_read)) + allocate(dim_names(num_var_dims)) ; dim_names(:) = "" + call get_variable_dimension_names(fileobj, trim(var_to_read), dim_names) + + do i=1,num_var_dims + if (is_dimension_unlimited(fileobj, dim_names(i))) then + time_dim = i + if (present(timelevel)) then + call get_dimension_size(fileobj, dim_names(i), dim_unlim_size) + if ((timelevel > dim_unlim_size) .and. is_root_PE()) call MOM_error(FATAL, & + trim(err_header)//"Attempting to read a time level of "//trim(var_to_read)//& + " that exceeds the size of the time dimension in "//trim(filename)) + endif + exit + endif + enddo + deallocate(dim_names) + + if (present(timelevel) .and. (time_dim < 0) .and. is_root_PE()) & + call MOM_error(WARNING, trim(err_header)//"time level specified, but the variable "//& + trim(var_to_read)//" does not have an unlimited dimension in "//trim(filename)) + if ((.not.present(timelevel)) .and. (time_dim > 0) .and. is_root_PE()) & + call MOM_error(WARNING, trim(err_header)//"The variable "//trim(var_to_read)//& + " has an unlimited dimension in "//trim(filename)//" but no time level is specified.") + if (present(has_time_dim)) has_time_dim = (time_dim > 0) + endif + + ! Registering the variable axes essentially just specifies the discrete position of this variable. + call MOM_register_variable_axes(fileobj, var_to_read, filename, position) + +end subroutine prepare_to_read_var + +!> register axes associated with a variable from a domain-decomposed netCDF file +subroutine MOM_register_variable_axes(fileObj, variableName, filename, position) + type(FmsNetcdfDomainFile_t), intent(inout) :: fileObj !< Handle to an open FMS2 netCDF file object + character(len=*), intent(in) :: variableName !< name of the variable + character(len=*), intent(in) :: filename !< The name of the file to read + integer, optional, intent(in) :: position !< A flag indicating where this data is discretized + + ! Local variables + character(len=256), allocatable, dimension(:) :: dim_names ! variable dimension names + integer, allocatable, dimension(:) :: dimSizes ! variable dimension sizes + logical, allocatable, dimension(:) :: is_x ! Is this a (likely domain-decomposed) x-axis + logical, allocatable, dimension(:) :: is_y ! Is this a (likely domain-decomposed) y-axis + logical, allocatable, dimension(:) :: is_t ! Is this a time axis or another unlimited axis + integer :: ndims ! number of dimensions + integer :: xPos, yPos ! Discrete positions for x and y axes. Default is CENTER + integer :: i + + xPos = CENTER ; yPos = CENTER + if (present(position)) then + if ((position == CORNER) .or. (position == EAST_FACE)) xPos = EAST_FACE + if ((position == CORNER) .or. (position == NORTH_FACE)) yPos = NORTH_FACE + endif + + ! get variable dimension names and lengths + ndims = get_variable_num_dimensions(fileObj, trim(variableName)) + allocate(dimSizes(ndims)) + allocate(dim_names(ndims)) + allocate(is_x(ndims)) ; is_x(:) = .false. + allocate(is_y(ndims)) ; is_y(:) = .false. + allocate(is_t(ndims)) ; is_t(:) = .false. + call get_variable_size(fileObj, trim(variableName), dimSizes) + call get_variable_dimension_names(fileObj, trim(variableName), dim_names) + call categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t) + + ! register the axes + do i=1,ndims + if ( .not.is_dimension_registered(fileobj, trim(dim_names(i))) ) then + if (is_x(i)) then + call register_axis(fileObj, trim(dim_names(i)), "x", domain_position=xPos) + elseif (is_y(i)) then + call register_axis(fileObj, trim(dim_names(i)), "y", domain_position=yPos) + else + call register_axis(fileObj, trim(dim_names(i)), dimSizes(i)) + endif + endif + enddo + + deallocate(dimSizes, dim_names, is_x, is_y, is_t) +end subroutine MOM_register_variable_axes + +!> Determine whether a variable's axes are associated with x-, y- or time-dimensions. Other +!! unlimited dimensions are also labeled as time axes for these purposes. +subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t) + type(FmsNetcdfDomainFile_t), intent(in) :: fileObj !< Handle to an open FMS2 netCDF file object + character(len=*), intent(in) :: filename !< The name of the file to read + integer, intent(in) :: ndims !< The number of dimensions associated with a variable + character(len=*), dimension(ndims), intent(in) :: dim_names !< Names of the dimensions associated with a variable + logical, dimension(ndims), intent(out) :: is_x !< Indicates if each dimension a (likely decomposed) x-axis + logical, dimension(ndims), intent(out) :: is_y !< Indicates if each dimension a (likely decomposed) y-axis + logical, dimension(ndims), intent(out) :: is_t !< Indicates if each dimension unlimited (usually time) axis + + ! Local variables + character(len=128) :: cartesian ! A flag indicating a Cartesian direction - usually a single character. + character(len=512) :: dim_list ! A concatenated list of dimension names. + character(len=128) :: units ! units corresponding to a specific variable dimension + logical :: x_found, y_found ! Indicate whether an x- or y- dimension have been found. + integer :: i + + x_found = .false. ; y_found = .false. + is_x(:) = .false. ; is_y(:) = .false. + do i=1,ndims + is_t(i) = is_dimension_unlimited(fileObj, trim(dim_names(i))) + ! First look for indicative variable attributes + if (.not.is_t(i)) then + if (variable_exists(fileobj, trim(dim_names(i)))) then + if (variable_att_exists(fileobj, trim(dim_names(i)), "cartesian_axis")) then + call get_variable_attribute(fileobj, trim(dim_names(i)), "cartesian_axis", cartesian) + cartesian = adjustl(cartesian) + if ((index(cartesian, "X") == 1) .or. (index(cartesian, "x") == 1)) is_x(i) = .true. + if ((index(cartesian, "Y") == 1) .or. (index(cartesian, "y") == 1)) is_y(i) = .true. + if ((index(cartesian, "T") == 1) .or. (index(cartesian, "t") == 1)) is_t(i) = .true. + endif + endif + endif + if (is_x(i)) x_found = .true. + if (is_y(i)) y_found = .true. + enddo + + if (.not.(x_found .and. y_found)) then + ! Next look for hints from axis names for uncharacterized axes + do i=1,ndims ; if (.not.(is_x(i) .or. is_y(i) .or. is_t(i))) then + call categorize_axis_from_name(dim_names(i), is_x(i), is_y(i)) + if (is_x(i)) x_found = .true. + if (is_y(i)) y_found = .true. + endif ; enddo + endif + + if (.not.(x_found .and. y_found)) then + ! Look for hints from CF-compliant axis units for uncharacterized axes + do i=1,ndims ; if (.not.(is_x(i) .or. is_y(i) .or. is_t(i))) then + call get_variable_units(fileobj, trim(dim_names(i)), units) + call categorize_axis_from_units(units, is_x(i), is_y(i)) + if (is_x(i)) x_found = .true. + if (is_y(i)) y_found = .true. + endif ; enddo + endif + + if (.not.(x_found .and. y_found) .and. ((ndims>2) .or. ((ndims==2) .and. .not.is_t(ndims)))) then + ! This is a case where one would expect to find x-and y-dimensions, but none have been found. + if (is_root_pe()) then + dim_list = trim(dim_names(1))//", "//trim(dim_names(2)) + do i=3,ndims ; dim_list = trim(dim_list)//", "//trim(dim_names(i)) ; enddo + call MOM_error(WARNING, "categorize_axes: Failed to identify x- and y- axes in the axis list ("//& + trim(dim_list)//") of a variable being read from "//trim(filename)) + endif + endif + +end subroutine categorize_axes + +!> Determine whether an axis is associated with the x- or y-directions based on a comparison of +!! its units with CF-compliant variants of latitude or longitude units. +subroutine categorize_axis_from_units(unit_string, is_x, is_y) + character(len=*), intent(in) :: unit_string !< string of units + logical, intent(out) :: is_x !< Indicates if the axis units are associated with an x-direction axis + logical, intent(out) :: is_y !< Indicates if the axis units are associated with an y-direction axis + + is_x = .false. ; is_y = .false. + select case (lowercase(trim(unit_string))) + case ("degrees_north"); is_y = .true. + case ("degree_north") ; is_y = .true. + case ("degrees_n") ; is_y = .true. + case ("degree_n") ; is_y = .true. + case ("degreen") ; is_y = .true. + case ("degreesn") ; is_y = .true. + case ("degrees_east") ; is_x = .true. + case ("degree_east") ; is_x = .true. + case ("degreese") ; is_x = .true. + case ("degreee") ; is_x = .true. + case ("degree_e") ; is_x = .true. + case ("degrees_e") ; is_x = .true. + case default ; is_x = .false. ; is_y = .false. + end select + +end subroutine categorize_axis_from_units + +!> Tries to determine whether the axis name is commonly associated with an x- or y- axis. This +!! approach is fragile and unreliable, but it a backup to reading a CARTESIAN file attribute. +subroutine categorize_axis_from_name(dimname, is_x, is_y) + character(len=*), intent(in) :: dimname !< A dimension name + logical, intent(out) :: is_x !< Indicates if the axis name is associated with an x-direction axis + logical, intent(out) :: is_y !< Indicates if the axis name is associated with an y-direction axis + + is_x = .false. ; is_y = .false. + select case(trim(lowercase(dimname))) + case ("grid_x_t") ; is_x = .true. + case ("nx") ; is_x = .true. + case ("nxp") ; is_x = .true. + case ("longitude") ; is_x = .true. + case ("long") ; is_x = .true. + case ("lon") ; is_x = .true. + case ("lonh") ; is_x = .true. + case ("lonq") ; is_x = .true. + case ("xh") ; is_x = .true. + case ("xq") ; is_x = .true. + case ("i") ; is_x = .true. + + case ("grid_y_t") ; is_y = .true. + case ("ny") ; is_y = .true. + case ("nyp") ; is_y = .true. + case ("latitude") ; is_y = .true. + case ("lat") ; is_y = .true. + case ("lath") ; is_y = .true. + case ("latq") ; is_y = .true. + case ("yh") ; is_y = .true. + case ("yq") ; is_y = .true. + case ("j") ; is_y = .true. + + case default ; is_x = .false. ; is_y = .false. + end select + +end subroutine categorize_axis_from_name + + !> Write a 4d field to an output file. subroutine write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value) type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing diff --git a/config_src/infra/FMS2/MOM_read_data_fms2.F90 b/config_src/infra/FMS2/MOM_read_data_fms2.F90 deleted file mode 100644 index c632e95177..0000000000 --- a/config_src/infra/FMS2/MOM_read_data_fms2.F90 +++ /dev/null @@ -1,376 +0,0 @@ -!> This module contains routines that encapsulate common preparatory work for FMS2 read_data calls -module MOM_read_data_fms2 - -! This file is part of MOM6. See LICENSE.md for the license. -use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, FATAL, WARNING, is_root_PE -use MOM_domain_infra, only : CENTER, CORNER, NORTH_FACE, EAST_FACE -use MOM_string_functions, only : lowercase -use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, check_if_open -use fms2_io_mod, only : variable_exists, get_num_variables, get_variable_names -use fms2_io_mod, only : get_variable_size, get_variable_units -use fms2_io_mod, only : get_variable_attribute, variable_att_exists -use fms2_io_mod, only : get_variable_num_dimensions, get_variable_dimension_names -use fms2_io_mod, only : is_dimension_unlimited, get_dimension_size -use fms2_io_mod, only : is_dimension_registered, register_axis - -implicit none ; private - -public prepare_to_read_var, find_varname_in_file - -contains - -!> Find the case-insensitive name match with a variable in an open domain-decomposed file-set, -!! prepare FMS2 to read this variable, and return some information needed to call fms2_read_data -!! correctly for this variable and file. -subroutine prepare_to_read_var(fileobj, fieldname, err_header, filename, var_to_read, & - has_time_dim, timelevel, position) - type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< An FMS2 handle to an open domain-decomposed file - character(len=*), intent(in) :: fieldname !< The variable name to seek in the file - character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(out) :: var_to_read !< The variable name to read from the file - logical, optional, intent(out) :: has_time_dim !< Indicates whether fieldname has a time dimension - integer, optional, intent(in) :: timelevel !< A time level to read - integer, optional, intent(in) :: position !< A flag indicating where this variable is discretized - - ! Local variables - logical :: variable_found ! Is a case-insensitive version of the variable found in the netCDF file? - character(len=256), allocatable, dimension(:) :: var_names ! The names of all the variables in the netCDF file - character(len=256), allocatable :: dim_names(:) ! The names of a variable's dimensions - integer :: nvars ! The number of variables in the file. - integer :: dim_unlim_size ! The current size of the unlimited (time) dimension in the file. - integer :: num_var_dims ! The number of dimensions a variable has in the file. - integer :: time_dim ! The position of the unlimited (time) dimension for a variable, or -1 - ! if it has no unlimited dimension. - integer :: i - - ! Open the file if necessary - if (.not.check_if_open(fileobj)) & - call MOM_error(FATAL, trim(err_header)//trim(filename)//" was not open in call to prepare_to_read_var.") - - ! Search for the variable in the file, looking for the case-sensitive name first. - if (variable_exists(fileobj, trim(fieldname))) then - var_to_read = trim(fieldname) - else ! Look for case-insensitive variable name matches. - nvars = get_num_variables(fileobj) - if (nvars < 1) call MOM_error(FATAL, "nvars is less than 1 for file "//trim(filename)) - allocate(var_names(nvars)) - call get_variable_names(fileobj, var_names) - - variable_found = .false. - do i=1,nvars - if (lowercase(trim(var_names(i))) == lowercase(trim(fieldname))) then - variable_found = .true. - var_to_read = trim(var_names(i)) - exit - endif - enddo - if (.not.(variable_found)) & - call MOM_error(FATAL, trim(err_header)//trim(fieldname)//" not found in "//trim(filename)) - deallocate(var_names) - endif - - ! FMS2 can not handle a timelevel argument if the variable does not have one in the file, - ! so some error checking and logic are required. - if (present(has_time_dim) .or. present(timelevel)) then - time_dim = -1 - - num_var_dims = get_variable_num_dimensions(fileobj, trim(var_to_read)) - allocate(dim_names(num_var_dims)) ; dim_names(:) = "" - call get_variable_dimension_names(fileobj, trim(var_to_read), dim_names) - - do i=1,num_var_dims - if (is_dimension_unlimited(fileobj, dim_names(i))) then - time_dim = i - if (present(timelevel)) then - call get_dimension_size(fileobj, dim_names(i), dim_unlim_size) - if ((timelevel > dim_unlim_size) .and. is_root_PE()) call MOM_error(FATAL, & - trim(err_header)//"Attempting to read a time level of "//trim(var_to_read)//& - " that exceeds the size of the time dimension in "//trim(filename)) - endif - exit - endif - enddo - deallocate(dim_names) - - if (present(timelevel) .and. (time_dim < 0) .and. is_root_PE()) & - call MOM_error(WARNING, trim(err_header)//"time level specified, but the variable "//& - trim(var_to_read)//" does not have an unlimited dimension in "//trim(filename)) - if ((.not.present(timelevel)) .and. (time_dim > 0) .and. is_root_PE()) & - call MOM_error(WARNING, trim(err_header)//"The variable "//trim(var_to_read)//& - " has an unlimited dimension in "//trim(filename)//" but no time level is specified.") - if (present(has_time_dim)) has_time_dim = (time_dim > 0) - endif - - ! Registering the variable axes essentially just specifies the discrete position of this variable. - call MOM_register_variable_axes(fileobj, var_to_read, filename, position) - -end subroutine prepare_to_read_var - -!> register axes associated with a variable from a domain-decomposed netCDF file -subroutine MOM_register_variable_axes(fileObj, variableName, filename, position) - type(FmsNetcdfDomainFile_t), intent(inout) :: fileObj !< Handle to an open FMS2 netCDF file object - character(len=*), intent(in) :: variableName !< name of the variable - character(len=*), intent(in) :: filename !< The name of the file to read - integer, optional, intent(in) :: position !< A flag indicating where this data is discretized - - ! Local variables - character(len=256), allocatable, dimension(:) :: dim_names ! variable dimension names - integer, allocatable, dimension(:) :: dimSizes ! variable dimension sizes - logical, allocatable, dimension(:) :: is_x ! Is this a (likely domain-decomposed) x-axis - logical, allocatable, dimension(:) :: is_y ! Is this a (likely domain-decomposed) y-axis - logical, allocatable, dimension(:) :: is_t ! Is this a time axis or another unlimited axis - integer :: ndims ! number of dimensions - integer :: xPos, yPos ! Discrete positions for x and y axes. Default is CENTER - integer :: i - - xPos = CENTER ; yPos = CENTER - if (present(position)) then - if ((position == CORNER) .or. (position == EAST_FACE)) xPos = EAST_FACE - if ((position == CORNER) .or. (position == NORTH_FACE)) yPos = NORTH_FACE - endif - - ! get variable dimension names and lengths - ndims = get_variable_num_dimensions(fileObj, trim(variableName)) - allocate(dimSizes(ndims)) - allocate(dim_names(ndims)) - allocate(is_x(ndims)) ; is_x(:) = .false. - allocate(is_y(ndims)) ; is_y(:) = .false. - allocate(is_t(ndims)) ; is_t(:) = .false. - call get_variable_size(fileObj, trim(variableName), dimSizes) - call get_variable_dimension_names(fileObj, trim(variableName), dim_names) - call categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t) - - ! register the axes - do i=1,ndims - if ( .not.is_dimension_registered(fileobj, trim(dim_names(i))) ) then - if (is_x(i)) then - call register_axis(fileObj, trim(dim_names(i)), "x", domain_position=xPos) - elseif (is_y(i)) then - call register_axis(fileObj, trim(dim_names(i)), "y", domain_position=yPos) - else - call register_axis(fileObj, trim(dim_names(i)), dimSizes(i)) - endif - endif - enddo - - deallocate(dimSizes, dim_names, is_x, is_y, is_t) -end subroutine MOM_register_variable_axes - -!> Determine whether a variable's axes are associated with x-, y- or time-dimensions. Other -!! unlimited dimensions are also labeled as time axes for these purposes. -subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t) - type(FmsNetcdfDomainFile_t), intent(in) :: fileObj !< Handle to an open FMS2 netCDF file object - character(len=*), intent(in) :: filename !< The name of the file to read - integer, intent(in) :: ndims !< The number of dimensions associated with a variable - character(len=*), dimension(ndims), intent(in) :: dim_names !< Names of the dimensions associated with a variable - logical, dimension(ndims), intent(out) :: is_x !< Indicates if each dimension a (likely decomposed) x-axis - logical, dimension(ndims), intent(out) :: is_y !< Indicates if each dimension a (likely decomposed) y-axis - logical, dimension(ndims), intent(out) :: is_t !< Indicates if each dimension unlimited (usually time) axis - - ! Local variables - character(len=128) :: cartesian ! A flag indicating a Cartesian direction - usually a single character. - character(len=512) :: dim_list ! A concatenated list of dimension names. - character(len=128) :: units ! units corresponding to a specific variable dimension - logical :: x_found, y_found ! Indicate whether an x- or y- dimension have been found. - integer :: i - - x_found = .false. ; y_found = .false. - is_x(:) = .false. ; is_y(:) = .false. - do i=1,ndims - is_t(i) = is_dimension_unlimited(fileObj, trim(dim_names(i))) - ! First look for indicative variable attributes - if (.not.is_t(i)) then - if (variable_exists(fileobj, trim(dim_names(i)))) then - if (variable_att_exists(fileobj, trim(dim_names(i)), "cartesian_axis")) then - call get_variable_attribute(fileobj, trim(dim_names(i)), "cartesian_axis", cartesian) - cartesian = adjustl(cartesian) - if ((index(cartesian, "X") == 1) .or. (index(cartesian, "x") == 1)) is_x(i) = .true. - if ((index(cartesian, "Y") == 1) .or. (index(cartesian, "y") == 1)) is_y(i) = .true. - if ((index(cartesian, "T") == 1) .or. (index(cartesian, "t") == 1)) is_t(i) = .true. - endif - endif - endif - if (is_x(i)) x_found = .true. - if (is_y(i)) y_found = .true. - enddo - - if (.not.(x_found .and. y_found)) then - ! Next look for hints from axis names for uncharacterized axes - do i=1,ndims ; if (.not.(is_x(i) .or. is_y(i) .or. is_t(i))) then - call categorize_axis_from_name(dim_names(i), is_x(i), is_y(i)) - if (is_x(i)) x_found = .true. - if (is_y(i)) y_found = .true. - endif ; enddo - endif - - if (.not.(x_found .and. y_found)) then - ! Look for hints from CF-compliant axis units for uncharacterized axes - do i=1,ndims ; if (.not.(is_x(i) .or. is_y(i) .or. is_t(i))) then - call get_variable_units(fileobj, trim(dim_names(i)), units) - call categorize_axis_from_units(units, is_x(i), is_y(i)) - if (is_x(i)) x_found = .true. - if (is_y(i)) y_found = .true. - endif ; enddo - endif - - if (.not.(x_found .and. y_found) .and. ((ndims>2) .or. ((ndims==2) .and. .not.is_t(ndims)))) then - ! This is a case where one would expect to find x-and y-dimensions, but none have been found. - if (is_root_pe()) then - dim_list = trim(dim_names(1))//", "//trim(dim_names(2)) - do i=3,ndims ; dim_list = trim(dim_list)//", "//trim(dim_names(i)) ; enddo - call MOM_error(WARNING, "categorize_axes: Failed to identify x- and y- axes in the axis list ("//& - trim(dim_list)//") of a variable being read from "//trim(filename)) - endif - endif - -end subroutine categorize_axes - -!> Determine whether an axis is associated with the x- or y-directions based on a comparison of -!! its units with CF-compliant variants of latitude or longitude units. -subroutine categorize_axis_from_units(unit_string, is_x, is_y) - character(len=*), intent(in) :: unit_string !< string of units - logical, intent(out) :: is_x !< Indicates if the axis units are associated with an x-direction axis - logical, intent(out) :: is_y !< Indicates if the axis units are associated with an y-direction axis - - is_x = .false. ; is_y = .false. - select case (lowercase(trim(unit_string))) - case ("degrees_north"); is_y = .true. - case ("degree_north") ; is_y = .true. - case ("degrees_n") ; is_y = .true. - case ("degree_n") ; is_y = .true. - case ("degreen") ; is_y = .true. - case ("degreesn") ; is_y = .true. - case ("degrees_east") ; is_x = .true. - case ("degree_east") ; is_x = .true. - case ("degreese") ; is_x = .true. - case ("degreee") ; is_x = .true. - case ("degree_e") ; is_x = .true. - case ("degrees_e") ; is_x = .true. - case default ; is_x = .false. ; is_y = .false. - end select - -end subroutine categorize_axis_from_units - -!> Tries to determine whether the axis name is commonly associated with an x- or y- axis. This -!! approach is fragile and unreliable, but it a backup to reading a CARTESIAN file attribute. -subroutine categorize_axis_from_name(dimname, is_x, is_y) - character(len=*), intent(in) :: dimname !< A dimension name - logical, intent(out) :: is_x !< Indicates if the axis name is associated with an x-direction axis - logical, intent(out) :: is_y !< Indicates if the axis name is associated with an y-direction axis - - is_x = .false. ; is_y = .false. - select case(trim(lowercase(dimname))) - case ("grid_x_t") ; is_x = .true. - case ("nx") ; is_x = .true. - case ("nxp") ; is_x = .true. - case ("longitude") ; is_x = .true. - case ("long") ; is_x = .true. - case ("lon") ; is_x = .true. - case ("lonh") ; is_x = .true. - case ("lonq") ; is_x = .true. - case ("xh") ; is_x = .true. - case ("xq") ; is_x = .true. - case ("i") ; is_x = .true. - - case ("grid_y_t") ; is_y = .true. - case ("ny") ; is_y = .true. - case ("nyp") ; is_y = .true. - case ("latitude") ; is_y = .true. - case ("lat") ; is_y = .true. - case ("lath") ; is_y = .true. - case ("latq") ; is_y = .true. - case ("yh") ; is_y = .true. - case ("yq") ; is_y = .true. - case ("j") ; is_y = .true. - - case default ; is_x = .false. ; is_y = .false. - end select - -end subroutine categorize_axis_from_name - - -!> Find the case-sensitive name of the variable in a netCDF file with a case-insensitive name match. -!! Optionally also determine whether this variable has an unlimited time dimension. -subroutine find_varname_in_file(fileobj, fieldname, err_header, filename, var_to_read, has_time_dim, timelevel) - type(FmsNetcdfFile_t), intent(inout) :: fileobj !< An FMS2 handle to an open NetCDF file - character(len=*), intent(in) :: fieldname !< The variable name to seek in the file - character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(out) :: var_to_read !< The variable name to read from the file - logical, optional, intent(out) :: has_time_dim !< Indicates whether fieldname has a time dimension - integer, optional, intent(in) :: timelevel !< A time level to read - - ! Local variables - logical :: variable_found ! Is a case-insensitive version of the variable found in the netCDF file? - character(len=256), allocatable, dimension(:) :: var_names ! The names of all the variables in the netCDF file - character(len=256), allocatable :: dim_names(:) ! The names of a variable's dimensions - integer :: nvars ! The number of variables in the file - integer :: dim_unlim_size ! The current size of the unlimited (time) dimension in the file. - integer :: num_var_dims ! The number of dimensions a variable has in the file. - integer :: time_dim ! The position of the unlimited (time) dimension for a variable, or -1 - ! if it has no unlimited dimension. - integer :: i - - ! Open the file if necessary - if (.not.check_if_open(fileobj)) & - call MOM_error(FATAL, trim(err_header)//trim(filename)//" was not open in call to find_varname_in_file.") - - ! Search for the variable in the file, looking for the case-sensitive name first. - if (variable_exists(fileobj, trim(fieldname))) then - var_to_read = trim(fieldname) - else ! Look for case-insensitive variable name matches. - nvars = get_num_variables(fileobj) - if (nvars < 1) call MOM_error(FATAL, "nvars is less than 1 for file "//trim(filename)) - allocate(var_names(nvars)) - call get_variable_names(fileobj, var_names) - - ! search for the variable in the file - variable_found = .false. - do i=1,nvars - if (lowercase(trim(var_names(i))) == lowercase(trim(fieldname))) then - variable_found = .true. - var_to_read = trim(var_names(i)) - exit - endif - enddo - if (.not.(variable_found)) & - call MOM_error(FATAL, trim(err_header)//trim(fieldname)//" not found in "//trim(filename)) - deallocate(var_names) - endif - - ! FMS2 can not handle a timelevel argument if the variable does not have one in the file, - ! so some error checking and logic are required. - if (present(has_time_dim) .or. present(timelevel)) then - time_dim = -1 - - num_var_dims = get_variable_num_dimensions(fileobj, trim(var_to_read)) - allocate(dim_names(num_var_dims)) ; dim_names(:) = "" - call get_variable_dimension_names(fileobj, trim(var_to_read), dim_names) - - do i=1,num_var_dims - if (is_dimension_unlimited(fileobj, dim_names(i))) then - time_dim = i - if (present(timelevel)) then - call get_dimension_size(fileobj, dim_names(i), dim_unlim_size) - if ((timelevel > dim_unlim_size) .and. is_root_PE()) call MOM_error(FATAL, & - trim(err_header)//"Attempting to read a time level of "//trim(var_to_read)//& - " that exceeds the size of the time dimension in "//trim(filename)) - endif - exit - endif - enddo - deallocate(dim_names) - - if (present(timelevel) .and. (time_dim < 0) .and. is_root_PE()) & - call MOM_error(WARNING, trim(err_header)//"time level specified, but the variable "//& - trim(var_to_read)//" does not have an unlimited dimension in "//trim(filename)) - if ((.not.present(timelevel)) .and. (time_dim > 0) .and. is_root_PE()) & - call MOM_error(WARNING, trim(err_header)//"The variable "//trim(var_to_read)//& - " has an unlimited dimension in "//trim(filename)//" but no time level is specified.") - if (present(has_time_dim)) has_time_dim = (time_dim > 0) - endif - -end subroutine find_varname_in_file - -end module MOM_read_data_fms2 From 7bdecbc62be1c636fd49071f4d7830a76ffa79af Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 23 Mar 2021 20:19:56 -0400 Subject: [PATCH 08/10] Add missing ".nc" to FMS2 output filenames Add a missing ".nc" suffix to the output filename with FMS2_io, while also issuing a warning, following the practice of FMS1. Also reordered the calls to add the longname and axis attributes to FMS2 files, to follow the order used in MOM6 calls to FMS1. All answers are bitwise identical, but there are some changes to output filenames and orders of attributes in files (to revert to traditional behavior). --- config_src/infra/FMS2/MOM_io_infra.F90 | 71 +++++++++++++++++++++----- 1 file changed, 57 insertions(+), 14 deletions(-) diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index 07a7c798d6..009050985d 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -20,6 +20,7 @@ module MOM_io_infra use fms_mod, only : write_version_number, open_namelist_file, check_nml_error use fms_io_mod, only : file_exist, field_exist, field_size, read_data use fms_io_mod, only : fms_io_exit, get_filename_appendix +use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_global_domain use mpp_io_mod, only : mpp_open, mpp_close, mpp_flush use mpp_io_mod, only : mpp_write_meta, mpp_write use mpp_io_mod, only : mpp_get_atts, mpp_attribute_exist @@ -315,7 +316,9 @@ subroutine open_file_type(IO_handle, filename, action, MOM_domain, threading, fi ! reading, writing or appending character(len=40) :: mode ! A character string that encodes whether the file is to be opened for ! reading, writing or appending + character(len=:), allocatable :: filename_tmp ! A copy of filename with .nc appended if necessary. character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file + integer :: index_nc if (IO_handle%open_to_write) then call MOM_error(WARNING, "open_file_type called for file "//trim(filename)//& @@ -332,6 +335,15 @@ subroutine open_file_type(IO_handle, filename, action, MOM_domain, threading, fi if (FMS2_writes .and. present(MOM_Domain)) then if (.not.associated(IO_handle%fileobj)) allocate (IO_handle%fileobj) + ! The FMS1 interface automatically appends .nc if necessary, but FMS2 interface does not. + index_nc = index(trim(filename), ".nc") + if (index_nc > 0) then + filename_tmp = trim(filename) + else + filename_tmp = trim(filename)//".nc" + if (is_root_PE()) call MOM_error(WARNING, "Open_file is appending .nc to the filename "//trim(filename)) + endif + if (file_mode == WRITEONLY_FILE) then ; mode = "write" elseif (file_mode == APPEND_FILE) then ; mode = "append" elseif (file_mode == OVERWRITE_FILE) then ; mode = "overwrite" @@ -342,9 +354,9 @@ subroutine open_file_type(IO_handle, filename, action, MOM_domain, threading, fi IO_handle%num_times = 0 IO_handle%file_time = 0.0 - if ((file_mode == APPEND_FILE) .and. file_exists(filename, MOM_Domain)) then + if ((file_mode == APPEND_FILE) .and. file_exists(filename_tmp, MOM_Domain)) then ! Determine the latest file time and number of records so far. - success = fms2_open_file(fileObj_read, trim(filename), "read", MOM_domain%mpp_domain) + success = fms2_open_file(fileObj_read, trim(filename_tmp), "read", MOM_domain%mpp_domain) call get_unlimited_dimension_name(fileObj_read, dim_unlim_name) if (len_trim(dim_unlim_name) > 0) & call get_dimension_size(fileObj_read, trim(dim_unlim_name), IO_handle%num_times) @@ -354,8 +366,8 @@ subroutine open_file_type(IO_handle, filename, action, MOM_domain, threading, fi call fms2_close_file(fileObj_read) endif - success = fms2_open_file(IO_handle%fileobj, trim(filename), trim(mode), MOM_domain%mpp_domain) - if (.not.success) call MOM_error(FATAL, "Unable to open file "//trim(filename)) + success = fms2_open_file(IO_handle%fileobj, trim(filename_tmp), trim(mode), MOM_domain%mpp_domain) + if (.not.success) call MOM_error(FATAL, "Unable to open file "//trim(filename_tmp)) IO_handle%FMS2_file = .true. elseif (present(MOM_Domain)) then call mpp_open(IO_handle%unit, filename, action=file_mode, form=NETCDF_FILE, threading=threading, & @@ -626,6 +638,7 @@ subroutine get_axis_data( axis, dat ) integer :: i + ! This routine might not be needed for MOM6. if (allocated(axis%ax_data)) then if (size(axis%ax_data) > size(dat)) call MOM_error(FATAL, & "get_axis_data called with too small of an output data array for "//trim(axis%name)) @@ -1010,6 +1023,7 @@ subroutine MOM_read_data_0d_int(filename, fieldname, data, timelevel) character(len=96) :: var_to_read ! Name of variable to read from the netcdf file logical :: success ! If true, the file was opened successfully + ! This routine might not be needed for MOM6. if (FMS2_reads) then ! Open the FMS2 file-set. success = fms2_open_file(fileObj, trim(filename), "read") @@ -1050,6 +1064,7 @@ subroutine MOM_read_data_1d_int(filename, fieldname, data, timelevel) character(len=96) :: var_to_read ! Name of variable to read from the netcdf file logical :: success ! If true, the file was opened successfully + ! This routine might not be needed for MOM6. if (FMS2_reads) then ! Open the FMS2 file-set. success = fms2_open_file(fileObj, trim(filename), "read") @@ -1741,6 +1756,7 @@ subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian character(len=:), allocatable :: cart ! A left-adjusted and trimmed copy of cartesian logical :: is_x, is_y, is_t ! If true, this is a domain-decomposed axis in one of the directions. integer :: position ! A flag indicating the axis staggering position. + integer :: i, isc, iec, global_size if (IO_handle%FMS2_file) then if (is_dimension_registered(IO_handle%fileobj, trim(name))) then @@ -1751,12 +1767,9 @@ subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian endif axis%name = trim(name) - if (present(data)) then - if (allocated(axis%ax_data)) call MOM_error(FATAL, & + if (present(data) .and. allocated(axis%ax_data)) call MOM_error(FATAL, & "Data is already allocated in a call to write_metadata_axis for axis "//& trim(name)//" in file "//trim(IO_handle%filename)) - allocate(axis%ax_data(size(data))) ; axis%ax_data(:) = data(:) - endif if (IO_handle%FMS2_file) then is_x = .false. ; is_y = .false. ; is_t = .false. @@ -1783,20 +1796,50 @@ subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian call register_axis(IO_handle%fileobj, trim(name), size(data)) endif + if (present(data)) then + ! With FMS2, the data for the axis labels has to match the computational domain on this PE. + if (present(domain)) then + ! The commented-out code on the next ~11 lines runs but there is missing data in the output file + ! call mpp_get_compute_domain(domain, isc, iec) + ! call mpp_get_global_domain(domain, size=global_size) + ! if (size(data) == global_size) then + ! allocate(axis%ax_data(iec+1-isc)) ; axis%ax_data(:) = data(isc:iec) + ! ! A simpler set of labels: do i=1,iec-isc ; axis%ax_data(i) = real(isc + i) - 1.0 ; enddo + ! elseif (size(data) == global_size+1) then + ! ! This is an edge axis. Note the effective SW indexing convention here. + ! allocate(axis%ax_data(iec+2-isc)) ; axis%ax_data(:) = data(isc:iec+1) + ! ! A simpler set of labels: do i=1,iec+1-isc ; axis%ax_data(i) = real(isc + i) - 1.5 ; enddo + ! else + ! call MOM_error(FATAL, "Unexpected size of data for "//trim(name)//" in write_metadata_axis.") + ! endif + + ! This works for a simple 1x1 IO layout, but gives errors for nontrivial IO layouts + allocate(axis%ax_data(size(data))) ; axis%ax_data(:) = data(:) + + else ! Store the entire array of axis labels. + allocate(axis%ax_data(size(data))) ; axis%ax_data(:) = data(:) + endif + endif + + ! Now create the variable that describes this axis. call register_field(IO_handle%fileobj, trim(name), "double", dimensions=(/name/)) - if (len_trim(units) > 0) & - call register_variable_attribute(IO_handle%fileobj, trim(name), 'units', & - trim(units), len_trim(units)) if (len_trim(longname) > 0) & call register_variable_attribute(IO_handle%fileobj, trim(name), 'long_name', & trim(longname), len_trim(longname)) + if (len_trim(units) > 0) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'units', & + trim(units), len_trim(units)) if (present(cartesian)) & call register_variable_attribute(IO_handle%fileobj, trim(name), 'cartesian_axis', & trim(cartesian), len_trim(cartesian)) if (present(sense)) & call register_variable_attribute(IO_handle%fileobj, trim(name), 'sense', sense) else + if (present(data)) then + allocate(axis%ax_data(size(data))) ; axis%ax_data(:) = data(:) + endif + call mpp_write_meta(IO_handle%unit, axis%AT, name, units, longname, cartesian=cartesian, sense=sense, & domain=domain, data=data, calendar=calendar) endif @@ -1831,12 +1874,12 @@ subroutine write_metadata_field(IO_handle, field, axes, name, units, longname, & do i=1,ndims ; dim_names(i) = trim(axes(i)%name) ; enddo prec_string = "double" ; if (present(pack)) then ; if (pack > 1) prec_string = "float" ; endif call register_field(IO_handle%fileobj, trim(name), trim(prec_string), dimensions=dim_names) - if (len_trim(units) > 0) & - call register_variable_attribute(IO_handle%fileobj, trim(name), 'units', & - trim(units), len_trim(units)) if (len_trim(longname) > 0) & call register_variable_attribute(IO_handle%fileobj, trim(name), 'long_name', & trim(longname), len_trim(longname)) + if (len_trim(units) > 0) & + call register_variable_attribute(IO_handle%fileobj, trim(name), 'units', & + trim(units), len_trim(units)) if (present(standard_name)) & call register_variable_attribute(IO_handle%fileobj, trim(name), 'standard_name', & trim(standard_name), len_trim(standard_name)) From ba643bde228221d07623c9a3dd3f4b71b7758854 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 25 Mar 2021 16:19:33 -0400 Subject: [PATCH 09/10] Explicit domain decomposition of horizontal axes FMS2 restart routines expect axes to be domain-decomposed. However, the domain_write_1d function does not apply this decomposition and instead routes this operation to compressed_write_1d. In order to accommodate this, we explicitly slice the 1d arrays of any axes into its domain-decomposed segment before passing to write_data. We have also introduced a control flag to MOM's FMS2 axistype to direct MOM_write_axis when this needs to be applied. We currently apply the domain decomposition flag to all horizontal axes regardless of circumstances. For now this is probably sufficient, but may need further testing (e.g. cube sphere). --- config_src/infra/FMS2/MOM_io_infra.F90 | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index 009050985d..ddc85da570 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -16,6 +16,7 @@ module MOM_io_infra use fms2_io_mod, only : variable_att_exists, get_variable_attribute, get_variable_num_dimensions use fms2_io_mod, only : get_variable_dimension_names, is_dimension_registered, get_dimension_size use fms2_io_mod, only : is_dimension_unlimited, register_axis, unlimited +use fms2_io_mod, only : get_global_io_domain_indices use fms_mod, only : write_version_number, open_namelist_file, check_nml_error use fms_io_mod, only : file_exist, field_exist, field_size, read_data @@ -132,6 +133,7 @@ module MOM_io_infra character(len=256) :: name !< The name of this axis in the files. type(mpp_axistype) :: AT !< The FMS1 axis-type that this type wraps real, allocatable, dimension(:) :: ax_data !< The values of the data on the axis. + logical :: domain_decomposed = .false. !< True if axis is domain-decomposed end type axistype !> For now, these module-variables are hard-coded to exercise the new FMS2 interfaces. @@ -1726,8 +1728,16 @@ subroutine MOM_write_axis(IO_handle, axis) type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing type(axistype), intent(in) :: axis !< An axis type variable with information to write + integer :: is, ie + if (IO_handle%FMS2_file) then - call write_data(IO_handle%fileobj, trim(axis%name), axis%ax_data) + if (axis%domain_decomposed) then + ! FMS2 does not domain-decompose 1d arrays, so we explicitly slice it + call get_global_io_domain_indices(IO_handle%fileobj, trim(axis%name), is, ie) + call write_data(IO_handle%fileobj, trim(axis%name), axis%ax_data(is:ie)) + else + call write_data(IO_handle%fileobj, trim(axis%name), axis%ax_data) + endif else call mpp_write(IO_handle%unit, axis%AT) endif @@ -1781,6 +1791,10 @@ subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian if ((index(cart, "T") == 1) .or. (index(cart, "t") == 1)) is_t = .true. endif + ! For now, we assume that all horizontal axes are domain-decomposed. + if (is_x .or. is_y) & + axis%domain_decomposed = .true. + if (is_x) then if (present(edge_axis)) then ; if (edge_axis) position = EAST_FACE ; endif call register_axis(IO_handle%fileobj, trim(name), 'x', domain_position=position) From 95ad9372b312476d2f741b5de38e8bf4f0cbb43e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 26 Mar 2021 07:34:11 -0400 Subject: [PATCH 10/10] Fix distributed reads of checksums using FMS2_io Corrects reads of hexadecimal checksum attributes from distributed files when using the FMS2 IO interfaces. All answers are bitwise identical, and reads and writes of distributed sets of restart files are now working with the FMS2 IO interfaces. --- config_src/infra/FMS2/MOM_io_infra.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index ddc85da570..df9d6dc7ca 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -484,6 +484,7 @@ subroutine get_file_fields(IO_handle, fields) character(len=256), dimension(size(fields)) :: var_names ! The names of all variables character(len=256) :: units ! The units of a variable as recorded in the file character(len=2048) :: longname ! The long-name of a variable as recorded in the file + character(len=64) :: checksum_char ! The hexadecimal checksum read from the file integer(kind=int64), dimension(3) :: checksum_file ! The checksums for a variable in the file integer :: nvar ! The number of variables in the file integer :: i @@ -503,8 +504,9 @@ subroutine get_file_fields(IO_handle, fields) fields(i)%valid_chksum = variable_att_exists(IO_handle%fileobj, var_names(i), "checksum") if (fields(i)%valid_chksum) then - call get_variable_attribute(IO_handle%fileobj, var_names(i), 'checksum', checksum_file) - fields(i)%chksum_read = checksum_file(1) + call get_variable_attribute(IO_handle%fileobj, var_names(i), 'checksum', checksum_char) + ! If there are problems, there might need to be code added to handle commas. + read (checksum_char(1:16), '(Z16)') fields(i)%chksum_read endif enddo else